---
title: "Factors - handling categorical data"
output: 
  html_document:
    toc: TRUE
    toc_float: TRUE 
---

# Data

We are introducing a new dataset: `melanoma`. Remember to Restart your R session to clear up the Environment from a previous session.

```{r, message = FALSE}
library(tidyverse)
library(finalfit)
library(broom)
library(here)
library(knitr)

mydata <- boot::melanoma # cleverly loads it from the package - as it's an example dataset included with the boot package

# your data usually arrives in a csv file:
mydata2 <- read_csv(here("melanoma_data.csv"))

```

Click on `mydata` in your Environment and have a look at the values - you'll see that categorical variables are coded as numbers, rather than text. Because of this, R thinks that these are numeric, or continuous variables, rather than factories with distinct levels or categories. You will need to recode these numbers into proper factors. 

You can also check the data using the `summary`, `glimpse` or `skim` from the `skimr` package. Each provide a different level of detail. 

```{r}
mydata %>% summary() 

mydata %>% glimpse()

#install.packages(skimr)
skimr::skim(mydata)
```


# Recoding factors - `fct_recode()`

Press F1 on `boot::melanoma` or type `?boot::melanoma` into the console to see its description. Use the information from help to change the numbers (e.g. 0 - female, 1 - male) into proper factors.

```{r}
mydata <- mydata %>% 
  mutate(status.factor = status %>% 
           # changes the variable to a factor 
           factor() %>% 
           # recodes the numbers to be meaningful levels, notice the order is New Name = Old Name 
           fct_recode("Died"  = "1",
                      "Alive" = "2",
                      "Died - other causes" = "3"))
  

mydata <- mydata %>%
  mutate(sex.factor = sex %>% 
           # changes the variable to a factor 
           factor() %>% 
           # recodes the numbers to be meaningful levels
           fct_recode("Female" = "0",
                      "Male"   = "1"))
  
mydata <- mydata %>%
  mutate(ulcer.factor = ulcer %>% 
           #changes the variable to a factor
           factor() %>% 
           # recodes the numbers to be meaningful levels
           fct_recode("Present" = "1",
                      "Absent"  = "0"))
```

You will notice that if you run the summary code from before using `glimpse()` that these new variables say `<fct>` 

```{r}
glimpse(mydata)
```


# Change the order of values within a factor - `fct_relevel()`

While character columns are plotted or run through statistical tests in alphabetical order, with factors we can use any order we define.

For example, look at the order of levels of `mydata$status.factor` (levels are all the unique and ordered values a factor can have):

```{r}
mydata %>% 
  pull(status.factor) %>% 
  levels()
```

The reason "Alive" is currently second, and not alphabetical is, that in the original factor it was number "2". And when we recoded "2" to "Alive", it retained its order. If, however, we want to make statistical comparisons in the future, relative to "Alive", we need to move it to the first position by using `fct_relevel()`:


```{r}
mydata <- mydata %>% 
  mutate(status.factor = status.factor %>%
           fct_relevel("Alive"))
```

```{r}
mydata %>% 
  pull(status.factor) %>%
  levels()
```

# Group factor levels together - `fct_collapse()`

`status.factor` currently has 3 levels, but we might want to make it binary (two levels) so let's group "Died - other causes" with those "Alive", to give a disease-specific mortality variable:


```{r}
mydata <- mydata %>%
  mutate(status.factor = status.factor %>% 
	fct_collapse("Alive" = c("Alive", "Died - other causes")))
```

```{r}
mydata %>% 
  pull(status.factor) %>% 
  levels()
```

# Converting a continuous (numeric) variable into a categorical variable - `cut()`

Sometimes we want to cut a numeric variable into categories. This may be when dealing with known thresholds (e.g. it is known that 5 mg of something is completely harmless but that there is a step change after that). Another reason for cutting a continuous variable into categories is if it doesn't have a linear response to the outcome but we want to use a linear model.

```{r}
mydata <- mydata %>%
  mutate(age.factor = age %>% 
           cut(breaks = c(4,20,40,60,95), include.lowest = TRUE))
```

* `include.lowest = TRUE` is a necessary historical quirk of the `cut()` function. Without it, the lowest value (in this case, 4 years old) would get set to NA (missing). Without `include.lowest = TRUE` the lowest level would appear as "(4,20]" which doesn't include 4.  Instead it looks like `[4` and does (check this by looking at `mydata` again). 

### Example

Investigate `mydata$age.factor` and use `fct_recode()` to give it more human readable levels.

```{r}
mydata %>% 
  pull(age.factor) %>% 
  levels()
```

For example:

```{r, eval = FALSE}
mydata <- mydata %>%
  mutate(age.factor = age.factor %>% 
           fct_recode("4 to 20" = "[4,20]"))
```

```{r, echo = FALSE}
mydata <- mydata %>%
  mutate(age.factor = age.factor %>% 
           fct_recode(" 4 to 20" =  "[4,20]",
                      "21 to 40" = "(20,40]",
                      "41 to 60" = "(40,60]",
                      "61 to 95" = "(60,95]"))
```

*Note:* if you are running this notebook interactively (by clicking `ctrl+enter` in each chunk of code) you might see a warning message here telling you that there is an unknown level in `f`. All this means is that because you renamed the level the first time you used `fct_recode` it now can't find the original name (`[4,20]`). In this case we can safely ignore the message because we have already changed the value. The message is usually helpful if you accidentally mistype a factor level so that you can then change it.  We know it isn't an error message because the code ran successfully and the message didn't start with "Error:", but unhelpfully it is in red font which can seem alarming!

```{r}
mydata %>% 
  pull(age.factor) %>% 
  levels()
```


# Export object or table

We could save our data as a csv file using `write_csv()`:

```{r warning=FALSE}
write_csv(mydata, "melanoma_edited.csv")
```

You can open a csv file with Excel and copy the table into Word or PowerPoint for presenting.

But when you read in a spreadsheet, if it is in this format, it will not know which variables are factors, and which are just plain characters. And therefore, it won't know which order to plot the factors in - you will have to reapply the `fct_relevel()` line.

To overcome this issue, we can save our data as an R data file, an `.rda` file.  For saving tibbles as `.rda` files, we do:

```{r}
save(mydata, file = here::here("melanoma_factored.rda"))
```

A new data file should have appeared in your files tab. 

R data, or `.rda` files are read in using load: 

```{r}
load(here::here("melanoma_factored.rda"))
```

Whereas spreadsheets (or formats from other stats packages) are read:

```{r, eval = FALSE}
# read_csv(here::here("melanoma_data.csv"))
```


## Plotting

Always plot new data first! Data visualization is one of the best ways to get familiar with your dataset. 

```{r}
mydata %>% 
  ggplot(aes(x = ulcer.factor, fill=status.factor)) + 
  geom_bar(position = "stack") +
  theme_bw() +
  scale_fill_brewer(palette = "Paired")
```

Try changing the `position = "stack"` argument within `geom_bar()` to be `position = "dodge"` and see how that changes the figure. Another option is `position = fill`. Take note of the y-axis - how does this differ from `position = "stack"`?

`position = fill` shows relative proportions at each x by stacking the bars and then standardising each bar to have the same height, which is why the y-axis is now ranges from 0-1.

`position = "dodge2` adds some empty space between the side-by-side two bars.  

```{r}
mydata %>% 
  ggplot(aes(x = age.factor, fill = status.factor)) +
  geom_bar(position = "dodge2") +
  theme_bw() +
  scale_fill_brewer(palette = "Paired")
```

```{r}
mydata %>% 
  ggplot(aes(x = ulcer.factor, fill=status.factor)) + 
  geom_bar(position = "dodge") +
  theme_bw() +
  scale_fill_brewer(palette = "Paired") +
  facet_grid(sex.factor~age.factor)
```

# Bonus Content: Analysis

The following sections are optional and more advanced. You will not be assessed on this content.

## Using base R

Let's test mortality against sex. As both of these variables are categorical, a $\chi^{2}$ test is appropriate. 

```{r}
#first produce a contigency table 
table(mydata$status.factor, mydata$sex.factor)

# then run the test - the function automatically makes a contigency table 
chisq.test(mydata$status.factor, mydata$sex.factor)

```

Note that `chisq.test()` defaults to the Yates continuity correction when using a 2x2 contingency table. It is fine to use this, but if you have a particular need not to, turn if off with `chisq.test(mydata$status.factor, mydata$sex.factor, correct=FALSE)`.

There are 3 types of $\chi^{2}$ tests: 

1. $\chi^{2}$ test goodness of fit 

2. $\chi^{2}$ test of independence 

3. $\chi^{2}$ test for homogeneity

If you are interested in looking into $\chi^{2}$ tests further, this article nicely explains the differences between them (https://www.simplypsychology.org/chi-square.html)

## Summarising factors with `library(finalfit)`

* HealthyR package
* Contains functions to summarise data for publication tables and figures
* Easily run regression analyses


```{r}
library(finalfit)

# F1 summary_factorlist to see options
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = "sex.factor",
                     p = TRUE,
                     column = TRUE)
```



```{r}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = "age.factor",
                     p = TRUE,
                     column = TRUE)
```


`finalfit` also works with continuous variables


```{r}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = "thickness",
                     p = TRUE,
                     column = TRUE)
```



## Summarising multiple factors

`finalfit` is useful for summarising multiple variables. We often want to summarise more than one factor or continuous variable against our dependent variable of interest. Think of Table 1 in a journal article. 

Here is a quick way of doing so using the `finalfit` package.


```{r, message = FALSE}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = c("sex.factor", "thickness"),
                     p = TRUE,
                     column = TRUE) %>%
  kable()

```


## Summarising factors with `library(tidyverse)`

### Example

`Tidyverse` gives the flexibility and power to examine millions of rows of your data any way you wish. The following are intended as an extension to what you have already done. These demonstrate some more advanced approaches to combining `tidy` functions. 

```{r}

# Calculate number of patients in each group
counted_data <- mydata %>%
  count(ulcer.factor, status.factor)

counted_data %>% 
  kable(col.names = c("Ucler", "Status", "n")) 

# Add the total number of people in each status group
counted_data2 <- counted_data %>%
  group_by(status.factor) %>%
  mutate(total = sum(n))   

counted_data2 %>% 
  kable(col.names = c("Ucler", "Status", "n", "Total")) %>% 
  kable_styling(bootstrap_options = "striped", full_width = F)
  
```

```{r}

# Calculate the percentage of n to total
counted_data3 <- counted_data2 %>%
  mutate(percentage = round(100*n/total, 1))

counted_data3 %>% 
  gt(rowname_col = "ulcer.factor") %>%
  tab_stubhead(label = "Ulcer") %>%
  cols_label(total = "Total", percentage = "Percentage") %>%
  opt_row_striping()
```

```{r}

# Create a combined columns of both n and percentage,
# using paste to add brackets around the percentage
counted_data4 <- counted_data3 %>% 
  mutate(count_perc = paste0(n, " (", percentage, ")"))

counted_data4
```


Or combine everything together without the intermediate `counted_data` breaks.

```{r}

mydata %>%
  count(ulcer.factor, status.factor) %>%
  group_by(status.factor) %>%
  mutate(total = sum(n)) %>%
  mutate(percentage = round(100*n/total, 1)) %>% 
  mutate(count_perc = paste0(n, " (", percentage, ")")) %>% 
  select(-total, -n, -percentage) %>% 
  pivot_wider(names_from = status.factor, values_from = count_perc) %>%
  kable(col.names = c("Ulcer", "Alive", "Died"))

```



### Exercise 1

Using `finalfit` create a summary table with "status.factor" as the dependent variable and the following as explanatory variables:  

* `sex.factor`  
* `ulcer.factor`  
* `age.factor`   
* `thickness`   



Try changing the table to show `median` and `interquartile range` instead of mean and SD.  

```{r}
# Exercise 1 - your R code:

```


### Exercise 2
By changing one and only one line in the following block create a new table showing the breakdown of `status.factor` by age and secondly the breakdown of `status.factor` by sex:

```{r, eval=FALSE}

mydata %>%
  count(ulcer.factor, status.factor) %>%
  group_by(status.factor) %>%
  mutate(total = sum(n)) %>%
  mutate(percentage = round(100*n/total, 1)) %>% 
  mutate(count_perc = paste0(n, " (", percentage, ")")) %>% 
  select(-total, -n, -percentage) %>% 
  pivot_wider(names_from = status.factor, values_from = count_perc)

```

### Exercise 3
Reproduce the table from Exercise 1 but this time with row-wise percentages instead of col(column)-wise percentages.

```{r}
# Exercise 3 - your R code:

```



## Bonus Content: Solutions 


### Excercise 1 

```{r}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = c("sex.factor", "thickness",
                                     "age.factor", "ulcer.factor"),
                     p = TRUE,
                     column = TRUE)
```

To change the data summaries to show the `median` and `interquartile range` instead of mean and SD, you need to change the `cont` argument default 

```{r}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = c("sex.factor", "thickness",
                                     "age.factor", "ulcer.factor"),
                     cont = "median",
                     column = TRUE)
```

### Excercise 2

```{r, eval=FALSE}

mydata %>%
  count(ulcer.factor, status.factor, age.factor, sex.factor) %>%
  group_by(status.factor) %>%
  mutate(total = sum(n)) %>%
  mutate(percentage = round(100*n/total, 1)) %>% 
  mutate(count_perc = paste0(n, " (", percentage, ")")) %>% 
  select(-total, -n, -percentage) %>% 
  pivot_wider(names_from = status.factor, values_from = count_perc)

```

### Excercise 3 

```{r}
mydata %>% 
  summary_factorlist(dependent = "status.factor", 
                     explanatory = c("sex.factor", "thickness",
                                     "age.factor", "ulcer.factor"),
                     cont = "median",
                     column = FALSE)
```

