steph_prac02_categorical.Rmd

# Comparing Two Categorical Variables: Practical 10

In this practical we focus on investigating the association between two categorical variables. We commence in the same way as the previous practical, loading libraries, setting options and themes, and reading in the data.


```{r, results = FALSE, warning = FALSE, message = FALSE}
#--- Load libraries - reminder to always load tidyverse last
library(epiDisplay)
library(magrittr)
library(foreign)
library(psych)
library(tidyverse)

#--- Change to show no scientific notation & round to 3 decimal places
options(scipen = 10, digits=3) 

#--- Set the plot theme to come out in black and white
theme_set(theme_bw())

#--- Read in the file
bab9 <- read.dta("./BAB9.dta", convert.factors = T) 
```


## Two-way frequency tables

We start by asking if birthweight (outcome) is associated with hypertension in the mother (exposure). To do this, we will want to group birthweight into two categories. Categorisation is a choice made by the researcher - you should be aware that you will lose information about this variable by making the simplifying assumption to categorise. However, in this case, there is a clinical definition of low birthweight, being <2500g, that may have real world relevance. We use **dplyr** commands from the **tidyverse** to generate this new variable, specifically mutate(). 

```{r, warnings = F, message = F}
#--- Generate low birth weight variable
bab9 <- bab9 %>% mutate(lbw = ifelse(bweight < 2500, "low", "normal"))
```

Let's briefly unpack this code. We pipe bab9 into the mutate() command to indicate we are working with the bab9 dataset. The mutate() command says we will create a new variable. We say that this new variable will be called lbw (low birthweight) and use an ifelse() statement to generate it. The ifelse statement says that _if_ birthweight is less than 2500g, then lbw should take the value "low", _else_ it should take the value "normal". Note that we do not need to generate variables and then label them as we do in Stata - R automatically detects that you are creating a factor (categorical) variable with two levels. Finally, we overwrite the original bab9 datafile with our new bab9 datafile that includes the created lbw variable. Note what has changed in the Environment pane.

The 2x2 frequency table and odds ratio calculations (including a chi-squared test and a Fisher's exact test for the null hypothesis that there is no association between the two variables) can be generated by the command cc() in the **epiDisplay** package. We add the option graph = F to suppress the default behaviour of the command to display one.

The **epiDisplay** command tabpct() also provides both row and column percentages for 2x2 tables. Note that the outcome variable has been put first, with the exposure variable second.

```{r}
#--- 2x2 table with OR
bab9 %$% cc(ht, lbw, graph = F)

#--- Reminder of same code, no pipe
# cc(bab9$ht, bab9$lbw, graph = F)

#--- 2x2 table with OR and percentages
bab9 %$% tabpct(ht, lbw, graph = F)
```


> Exercise 2.1: From the cc() table, how many low birthweight babies were born to women with hypertension? How many low birthweight babies were born to women with no hypertension?

> Exercise 2.2: From the tabpct() table, what proportion of babies born to women who were hypertensive were of low birthweight? How does this compare with the proportion of low birthweight babies born to women who were not hypertensive?

> Exercise 2.3: From the cc() table, what is your interpretation of the chi-squared test?


## Z-test for proportions

To carry out a z test to compare two proportions, we can use the prop.test() command, equivalent to _prtest_ in Stata. It is however, slightly more finicky to use. prop.test() requires two inputs: a vector of 'successes' (numerator) and a vector of 'counts' (denominator). A 'vector' in R is a sequence of objects that are of the same type (i.e. numerical or character...) - it is similar to a 'list', but a list may contain objects of any data type. 

You will thus need to manually extract these vectors from the output of the 2x2 table. For example, the vector of 'successes' (hypertension = yes) for these data is c(27, 62). The c() function creates a vector, the c standing for concatenate. The test is below, invoking the correct = F option to suppress the default behaviour of R of adding a continuity correction. 

Note that R does not provide the z score, but it can be derived by noting that the z score is the square root of the chi square -- the cc() command is much less work than running the actual test - thus demonstrating the utility of user-built packages! In the Environment pane, we can click on the ztest object and look at the other objects contained within it - one of which is 'statistic'. We use the $ symbol to show that we want the 'statistic' object from the 'ztest' object, called indexing. You have seen this before to extract variables from data frames. If you click on the 'ztest' object in the environment pane, a list of (nine) objects contained within it can be seen. Try and extract the p value from the test using this indexing strategy.

```{r}
#--- Run a z test
ztest <- prop.test(c(27, 53), n = c(89, 552), correct = F)
ztest

#--- Extract the z statistic
sqrt(ztest$statistic)
```


> Exercise 2.4: What is your interpretation of this result?

> Exercise 2.5: Check the relationship between the z and the chi-squared values obtained.

## Trends in proportions

Since birthweight is known to increase with length of gestation, we will use this fact to test for a linear trend in proportions. Gestational age was measured to the nearest week. 

Let's examine the histogram for gestational age.

```{r, warning = F, message = F}
#--- Plot gestational age
bab9 %>% ggplot(aes(x = gestwks)) + geom_histogram()
```

We notice that here we have left skew. This means that if we categorised the data by dividing gestwks into equally spaced groups of time (i.e. groups of 3.5 weeks) we would end up with unequal groups. So instead, we put an equal number of individuals into each group using the ntile() command and convert it to a factor. ntile() takes two arguments: the variable to divide into ntiles, and the number of ntiles desired.

The describeBy() command from the psych package gives a number of summary statistics for gestwks by values of gest5.

We can also get a 2x2 table from cc() in EpiDisplay. 

```{r}
#--- Get the quintiles
bab9$gest5 <- as.factor(ntile(bab9$gestwks, 5))

summary(bab9$gest5)

#--- Extract the top of each quintile by finding the maximum (and thus the cut points)
bab9 %>% group_by(gest5) %>% dplyr::summarise(max(gestwks))

#--- Examine gestational age by quintile
bab9 %$% describeBy(gestwks, gest5)

#--- Table for low birthweight and gestational age
bab9 %$% cc(lbw, gest5, graph= F)
```


## Further exercises

We read in the bab (original) dataset, replicating the command given in the first chunk and create a new variable for birthweight divided into those less than 2500g, those between 2500g and 3000g, and those over 3000g. We do this with the **tidyverse** function _case_when()_. This function uses the formula syntax we have seen before: we explicitly set out the condition on the left hand side of the tilde and the value we want on the right (e.g. the first line of case_when says that when birthweight is less than or equal to 2500, set the value of lbw3 to 1). To use multiple conditions, we use the & operator to stand for 'AND'. We can also use | to stand for 'OR' and ! to stand for 'NOT'. Note that after applying _case_when()_ we then convert this variable to a factor. 

We can also use a nested ifelse approach - if the bweight variable is less than 2500, give it the value 1, else continue to the next statement - we repeat this process. In general, _case_when()_ is more clear than using _ifelse()_ as it explicitly sets out each new level and generates an NA when a variable does not meet any conditions, while _ifelse()_ can lead to every value not specified being lumped into the final _else_ condition. In this case, we have only three possible categories and thus the ifelse approach is relatively straightforward, but nesting many statements can get rapidly confusing.


> Exercise 2.7: What percentage of birthweights were 3000g or more?


```{r}
#--- Read in the bab file
bab <- read.dta("./BAB.dta", convert.factors = T) 

#--- Create the low birthweight categorical variable
bab <- bab %>% mutate(lbw3 = case_when(bweight <= 2500 ~ 1,
                                       bweight < 3000 & bweight > 2500 ~ 2,
                                       bweight >= 3000 ~ 3)) %>% 
               mutate(lbw3 = as.factor(lbw3))

summary(bab$lbw3)

#--- Same as above using a nested ifelse approach                        
# bab <- bab %>% mutate(lbw3 = ifelse(bweight <= 2500, 1, ifelse(bweight < 3000, 2, 3))) %>% 
# mutate(lbw3 = as.factor(lbw3))

```

Now we conduct some further analyses.


> Exercise 2.8: Analyse whether there is any association between sex of infant and birthweight (using lbw3).
How many degrees of freedom is the chi-square test based on?


```{r}
#--- Get a 2x2 table with a chi-square test
bab %$% cc(sex, lbw3, graph = F)
```

Now we categorise maternal age into four groups (<30, 30-34, 35-39, 40+) to investigate the relationship with hypertension. We could do this by using three nested ifelse statements, or by using _case_when()_ but we demonstrate a straightforward way to categorise the continuous variable using the _cut()_ function.


> Exercise 2.9: Is there any evidence of a linear association between hypertension and maternal age?

```{r}
#--- Create the new variable and extract a summary
bab <- bab %>% mutate(matagegp = cut(matage,
    breaks = c(0, 29, 34, 39,+Inf),
    labels = c("<30", "30-34", "35-39", "40+")
  ))

summary(bab$matagegp)

#--- Investigate the relationship between maternal age group and hypertension
bab %$% cc(ht, matagegp, graph = F)
```

We create a new variable that codes which quarter of gestation age each baby falls into. This uses exactly the same code as presented above and the same ntile() function. 


> Exercise 2.10: Find the values of gestational age which define the quartiles.

> Exercise 2.11: What proportion of infants in the lowest quarter of gestational age were males?

```{r}
#--- Get the quintiles
bab$gest4 <- as.factor(ntile(bab$gestwks, 4))

summary(bab$gest4)

#--- Extract the top of each quintile (and thus the cut points)
bab %>% group_by(gest4) %>% dplyr::summarise(max(gestwks))

#--- Get the 2x2 table of gender and gestational quartile
bab %$% tabpct(sex, gest4, graph = F)
```