---
title: "Measures of Assocation using epiR"
author: "Mark Stevenson and Evan Sergeant"
date: "`r Sys.Date()`"
bibliography: epiR_measures_of_assoc.bib
link-citations: yes
output:
  html_document:
    fig_caption: true
    css: styles.css
    toc: false
vignette: >
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteIndexEntry{Measures of association}
  %\VignetteEncoding{UTF-8}
  
  
---
```{r, echo = FALSE, message = FALSE}

# Use tinytex (instead of MiKTex) to generate PDFs. See -Tools -Global options -Sweave. Use tinytex.
# tinytex::install_tinytex()

library(dplyr); library(flextable); library(knitr); library(officer);  library(tidyr)
knitr::opts_chunk$set(collapse = T, comment = "#>")
options(tibble.print_min = 4L, tibble.print_max = 4L)
```
A common task in epidemiology is to quantify the strength of association between exposures ('risk factors') and disease outcomes. In this context the term 'exposure' is taken to mean a variable whose association with the outcome is to be estimated. 

Exposures can be harmful, beneficial or both harmful and beneficial (e.g., if an immunisable disease is circulating, exposure to immunising agents helps most recipients but may harm those who experience adverse reactions). The term 'outcome' is used to describe all the possible results that may arise from exposure to a causal factor or from preventive or therapeutic interventions [@porta:2014]. In human and animal health an 'outcome-positive' individual is an individual who has experienced a given disease of interest. 

In this vignette we describe how `epiR` can be used to compute the various measures of association used in epidemiology notably the risk ratio, odds ratio, attributable risk in the exposed, attributable fraction in the exposed, attributable risk in the population and attributable fraction in the population. Examples are provided to demonstrate how the package can be used to deal with exposure-outcome data presented in various formats.  

This vignette has been written assuming the reader routinely formats their 2 $\times$ 2 table data with the outcome status as columns and exposure status as rows. If this is not the case the argument `outcome = "as.columns"` (the default) can be changed to `outcome = "as.rows"`.

## Measures of association strength

### The incidence risk ratio

Consider a study where study subjects are monitored for the presence of disease over a defined follow-up period. At the start of the follow-up period study subjects are classified according to exposure to a hypothesised risk factor. If the exposure and outcome are binary variables (yes or no) we can present the counts of subjects in each of the four exposure-outcome categories in a 2 $\times$ 2 table.

```{r echo = FALSE, results = 'asis'}
twobytwo.df <- data.frame("exp" = c("Exposure+","Exposure-","Total"), "dpos" = c("a","c","a + c"), "dneg" = c("b","d","b + d"), "total" = c("a + b","c + d","a + b + c + d"))

# Create a header key data frame:
hkey.df <- data.frame(col_keys = c("exp","dpos","dneg","total"),
  h1 = c("", "Outcome+", "Outcome-", "Total"), stringsAsFactors = FALSE)

# Create table:
caption.t <- flextable::as_paragraph(as_chunk("Table 1: A 2 × 2 contingency table.",
    props = fp_text(font.size = 9, font.family = "Arial", bold = FALSE)))

border_h = fp_border(color = "black", width = 1)

ft <- flextable(twobytwo.df) %>%
  width(j = 1, width = 1.00) %>%
  width(j = 2, width = 1.00) %>%
  width(j = 3, width = 1.00) %>%
  width(j = 4, width = 1.00) %>%
  
  set_header_df(mapping = hkey.df, key = "col_keys") %>%
  fontsize(size = 9, part = "all") %>%
  
  bg(bg = "grey80", part = "header") %>%
  hline_top(border = border_h, part = "all" ) %>%
  align(align = "left", part = "all") %>%
  set_caption(caption = caption.t)
ft
```
When our data are in this format we can calculate the incidence risk of the outcome in those that were exposed $R_{[D+|E+]}$, the incidence risk in those that were not exposed $R_{[D+|E-]}$ and finally, the incidence risk in both the exposed and the unexposed (i.e., the entire study population) $R_{[D+|E\pm]}$:

```{r echo = FALSE, results = 'asis'}
irr.df <- data.frame(
  "exp" = c("Exposure+","Exposure-","Total"), 
  "dpos" = c("a","c","a + c"), 
  "dneg" = c("b","d","b + d"), 
  "total" = c("a + b","c + d","a + b + c + d"), 
  risk = c("R[D+|E+] = a \U00F7 (a + b)", "R[D+|E-] = c \U00F7 (c + d)", "R[D+|E±] = (a + c) \U00F7 (a + b + c + d)"))

# Create a header key data frame:
hkey.df <- data.frame(col_keys = c("exp","dpos","dneg","total","risk"),
  h1 = c("", "Outcome+", "Outcome-", "Total", "Risk"), stringsAsFactors = FALSE)

# Create table:
caption.t <- flextable::as_paragraph(as_chunk("Table 2: A 2 × 2 contingency table with outcome risks calculated for the exposure positive, the exposure negative and the entire study population.",
    props = fp_text(font.size = 9, font.family = "Arial", bold = FALSE)))

border_h = fp_border(color = "black", width = 1)

ft <- flextable(irr.df) %>%
  width(j = 1, width = 1.00) %>%
  width(j = 2, width = 1.00) %>%
  width(j = 3, width = 1.00) %>%
  width(j = 4, width = 1.00) %>%
  width(j = 5, width = 2.00) %>%
  
  set_header_df(mapping = hkey.df, key = "col_keys") %>%
  fontsize(size = 9, part = "all") %>%
  
  bg(bg = "grey80", part = "header") %>%
  hline_top(border = border_h, part = "all" ) %>%
  align(align = "left", part = "all") %>%
  set_caption(caption = caption.t)

vfixa <- c("[D+|E+]","[D+|E-]","[D+|E±]")
vfixb <- c(" = a \U00F7 (a + b)"," = c \U00F7 (c + d)"," = (a + c) \U00F7 (a + b + c + d)")

for (i in 1:3) {
  ft <- compose(
    ft, 
    j = "risk", 
    i = i,
    value = as_paragraph(
      as_chunk("R", props = fp_text(font.size = 9)),
      as_chunk(vfixa[i], props = fp_text(vertical.align = "subscript", font.size = 9)),
      as_chunk(vfixb[i], props = fp_text(font.size = 9))
    )
  )
}

ft
```

In Figure 1 we present the outcome incidence risk in the exposed, unexposed and exposed plus unexposed groups as a bar chart. The outcome incidence risk ratio equals the outcome incidence risk in the exposed divided by the outcome incidence risk in the unexposed:

```{r massoc01, echo = FALSE, out.width = "70%", fig.cap = "\\label{fig:massoc01}Figure 1: The incidence risk ratio.", out.width = "80%", fig.align = "center"}
knitr::include_graphics("figures/risk_ratio.png")
```

The outcome incidence risk ratio provides an estimate of how many times more likely exposed individuals are to experience the outcome of interest, compared with non-exposed individuals. 

If the outcome incidence risk ratio equals 1, then the risk of the outcome in both the exposed and non-exposed groups are the same. If the outcome incidence risk ratio is greater than 1, then exposure increases the outcome incidence risk with greater departures from 1 indicative of a stronger association If the outcome incidence risk ratio is less than 1, exposure reduces the outcome incidence risk and exposure is said to be protective.

### The odds ratio --- cohort studies

In a cohort study definition of exposure status (exposure-positive, exposure-negative) comes first. Subjects are then followed over time to determine their outcome status (outcome-positive, outcome-negative). The outcome odds in the exposed and unexposed populations are calculated as follows:

```{r echo = FALSE, results = 'asis'}
orcohort.df <- data.frame(
  "exp" = c("Exposure+","Exposure-","Total"), 
  "dpos" = c("a","c","a + c"), "dneg" = c("b","d","b + d"), 
  "total" = c("a + b","c + d","a + b + c + d"), 
  odds = c("OE+ = a \U00F7 b","OE- = c \U00F7 d", "OT = (a + c) \U00F7 (b + d)"))

# Create a header key data frame:
hkey.df <- data.frame(col_keys = c("exp","dpos","dneg","total","odds"),
  h1 = c("", "Outcome+", "Outcome-", "Total", "Odds"), stringsAsFactors = FALSE)

# Create table:
caption.t <- flextable::as_paragraph(as_chunk("Table 3: A 2 × 2 contigency table with outcome odds calculated for the exposure positive, the exposure negative and the entire study population.",
    props = fp_text(font.size = 9, font.family = "Arial", bold = FALSE)))

border_h = fp_border(color = "black", width = 1)

ft <- flextable(orcohort.df) %>%
  width(j = 1, width = 1.00) %>%
  width(j = 2, width = 1.00) %>%
  width(j = 3, width = 1.00) %>%
  width(j = 4, width = 1.00) %>%
  width(j = 5, width = 2.00) %>%
  
  set_header_df(mapping = hkey.df, key = "col_keys") %>%
  fontsize(size = 9, part = "all") %>%
  
  bg(bg = "grey80", part = "header") %>%
  hline_top(border = border_h, part = "all" ) %>%
  align(align = "left", part = "all") %>%
  set_caption(caption = caption.t)

vfixa <- c("[D+|E+]","[D+|E-]","[D+|E±]")
vfixb <- c(" = a \U00F7 b"," = c \U00F7 d"," = (a + c) \U00F7 (b + d)")

for (i in 1:3) {
  ft <- compose(
    ft, 
    j = "odds", 
    i = i,
    value = as_paragraph(
      as_chunk("O", props = fp_text(font.size = 9)),
      as_chunk(vfixa[i], props = fp_text(vertical.align = "subscript", font.size = 9)),
      as_chunk(vfixb[i], props = fp_text(font.size = 9))
    )
  )
}
ft
```
The outcome odds ratio for a cohort study is the outcome odds in the exposed divided by the outcome odds in the unexposed.

### The odds ratio --- case-control studies

In a case-control study outcome status (outcome-positive, outcome-negative) is defined first. The history provided by each study subject then provides information about exposure status (exposure-positive, exposure-negative). For case-control studies, instead of talking about the outcome odds in the exposed and unexposed groups (as we did when we were working with data from a cohort study) we talk about the odds of *exposure odds* in the case (outcome-positive) and control (outcome-negative) groups. To emphasise this difference we list outcome status (case, control) as rows and exposure status (exposure-positive, exposure-negative) in our 2 $\times$ 2 table. A way to remember this is to think about the temporal order of events in a study. In a case-control study outcome status is defined first: exposure status is defined second. Because outcome status comes first it is listed as the first column of the 2 $\times$ 2 table.

```{r echo = FALSE, results = 'asis'}
orcc.df <- data.frame(
  "out" = c("Outcome+","Outcome-","Total"), 
  "epos" = c("a","c","a + c"), 
  "eneg" = c("b","d","b + d"), 
  "total" = c("a + b","c + d","a + b + c + d"), 
  odds = c("OE+ = a \U00F7 b","OE- = c \U00F7 d", "OT = (a + c) \U00F7 (b + d)"))

# Create a header key data frame:
hkey.df <- data.frame(col_keys = c("out","epos","eneg","total","odds"),
  h1 = c("", "Exposure+", "Exposure-", "Total", "Odds"), stringsAsFactors = FALSE)

# Create table:
caption.t <- flextable::as_paragraph(as_chunk("Table 4: A 2 × 2 contingency table with exposure odds calculated for the outcome positive, the outcome negative and the entire study population.",
    props = fp_text(font.size = 9, font.family = "Arial", bold = FALSE)))

border_h = fp_border(color = "black", width = 1)

ft <- flextable(orcc.df) %>%
  width(j = 1, width = 1.00) %>%
  width(j = 2, width = 1.00) %>%
  width(j = 3, width = 1.00) %>%
  width(j = 4, width = 1.00) %>%
  width(j = 5, width = 2.00) %>%
  
  set_header_df(mapping = hkey.df, key = "col_keys") %>%
  fontsize(size = 9, part = "all") %>%
  
  bg(bg = "grey80", part = "header") %>%
  hline_top(border = border_h, part = "all" ) %>%
  align(align = "left", part = "all") %>%
  set_caption(caption = caption.t)

vfixa <- c("[E+|D+]","[E+|D-]","[E+|D±]")
vfixb <- c(" = a \U00F7 b"," = c \U00F7 d"," = (a + c) \U00F7 (b + d)")

for (i in 1:3) {
  ft <- compose(
    ft, 
    j = "odds", 
    i = i,
    value = as_paragraph(
      as_chunk("O", props = fp_text(font.size = 9)),
      as_chunk(vfixa[i], props = fp_text(vertical.align = "subscript", font.size = 9)),
      as_chunk(vfixb[i], props = fp_text(font.size = 9))
    )
  )
}
ft
```
The odds ratio for a case-control study is defined as exposure odds in the cases $O_{[E+|D+]}$ divided by the exposure odds in the controls $O_{[E+|D-]}$. 

Note that the numeric estimate of the odds ratio for a case-contrl study is exactly the same as that calculated for a cohort study. The expression of the result is the only thing that differs. Key points:

* In a *cohort study* we talk about the outcome odds being $x$ times greater (or less) in *the exposed*, compared with the unexposed. 
* In a *case-control study* we talk about the exposure odds being $x$ times greater (or less) in *cases*, compared with controls.

## Measures of effect in the exposed

### The attributable risk in the exposed

The attributable risk is defined as the increase or decrease in the incidence risk of the outcome in the exposed that is attributable to exposure (Figure 2). The attributable risk (unlike the incidence risk ratio) provides a measure of the absolute frequency of the outcome that is associated with exposure.

```{r massoc02, echo = FALSE, out.width = "70%", fig.cap = "\\label{fig:massoc02}Figure 2: The attributable risk in the exposed.", out.width = "80%", fig.align = "center"}
knitr::include_graphics("figures/attributable_risk.png")
```

A useful way of expressing attributable risk in a clinical setting is in terms of the number needed to treat, NNT. NNT equals the inverse of the attributable risk. Depending on the outcome of interest we use different labels for NNT. 

When dealing with an outcome that is 'desirable' (e.g., treatment success) we call NNT the number needed to treat for benefit, NNTB. NNTB equals the number of subjects who would have to be exposed to result in a single (desirable) outcome. When dealing with an outcome that is 'undesirable' (e.g., death) we call NNT the number needed to treat for harm, NNTH. NNTH equals the number of subjects who would have to be exposed to result in a single (undesirable) outcome.

### The attributable fraction in the exposed

The attributable fraction in the exposed is the proportion of outcome-positive subjects in the exposed group that is due to exposure (Figure 3).

```{r massoc03, echo = FALSE, out.width = "70%", fig.cap = "\\label{fig:massoc03}Figure 3: The attributable fraction in the exposed.", out.width = "80%", fig.align = "center"}
knitr::include_graphics("figures/attributable_fraction.png")
```

## Measures of effect in the population

### The attributable risk in the population

The population attributable risk is the increase or decrease in incidence risk of the outcome in the study population that is attributable to exposure (Figure 4).

```{r massoc04, echo = FALSE, out.width = "70%", fig.cap = "\\label{fig:massoc04}Figure 4: The attributable risk in the population.", out.width = "80%", fig.align = "center"}
knitr::include_graphics("figures/population_attributable_risk.png")
```

### The attributable fraction in the population

The population attributable fraction (also known as the aetiologic fraction) is the proportion of outcome-positive subjects in the study population that is due to the exposure (Figure 5).

```{r massoc05, echo = FALSE, out.width = "70%", fig.cap = "\\label{fig:massoc05}Figure 5: The attributable fraction in the population.", out.width = "80%", fig.align = "center"}
knitr::include_graphics("figures/population_attributable_fraction.png")
```

On the condition that the exposure of interest is a cause of the disease outcome, the population attributable fraction represents the proportional reduction in average disease risk over a specified period of time that would be achieved by eliminating the exposure while the distribution of other risk factors in the population remained unchanged. For this reason, PAFs are particularly useful to guide policy makers when planning public health interventions. 

In summary, if you're going to use PAFs as a means for informing policy, make sure that: (1) the exposure of interest is causally related to the outcome; and (2) the exposure of interest is something amenable to intervention.

## Theory to practice: Calculating measures of association using R

### Direct entry of 2 by 2 table contingency table cell frequencies

Firstly, a 2 $\times$ 2 table can be created by listing the contingency table cell frequencies in vector format. Take the following example. 

A cohort study investigating the relationship between dry cat food (DCF) and presence or absence of feline lower urinary tract disease (FLUTD) was conducted [@willeberg:1977]. Counts of individuals in each group were as follows. DCF-exposed cats (outcome-positive, outcome-negative) 13, 2163. Non DCF-exposed cats (outcome-positive, outcome-negative) 5, 3349. We can enter these data directly into R as a vector of length four. Check that your counts have been entered in the correct order by viewing the data as a matrix.

```{r}
dat.v01 <- c(13,2163,5,3349); dat.v01

# View the data in the usual 2 by 2 table format:
matrix(dat.v01, nrow = 2, byrow = TRUE)
```

Calculate the outcome incidence risk ratio, incidence odds ratio, attributable prevalence in the exposed, attributable fraction in the exposed, attributable prevalence in the population and the attributable fraction in the population using function `epi.2by2`. Note that we use the terms incidence risk ratio and incidence odds ratio here because we're dealing with data from a cohort study --- following determination of dry food exposure status cats have been followed over time to determine the incidence risk of FLUTD. 

```{r, message = FALSE}
library(epiR)

epi.2by2(dat = dat.v01, method = "cohort.count", elab = "Dry food", olab = "FLUTD", digits = 2, conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
```

The incidence risk of FLUTD among those that were dry food positive was 4.0 (95% CI 1.4 to 11) times the incidence risk of FLUTD among those that were dry food negative. 

In DCF exposed cats, 75% of FLUTD was attributable to DCF (95% CI 33% to 91%). Fifty-four percent of FLUTD cases in this cat population were attributable to DCF (95% CI 32% to 75%).

Need a hand to get the correct wording to explain each of the listed measures of association and measures of effect? Set `interpret = TRUE` in `epi.2by2`:

```{r, message = FALSE}
epi.2by2(dat = dat.v01, method = "cohort.count", elab = "Dry food", olab = "FLUTD", digits = 2, conf.level = 0.95, units = 100, interpret = TRUE, outcome = "as.columns")
```

### Data frame with one row per observation

Here we provide examples where you have exposure status and outcome status listed for each member of your study population. There are two options for contingency table preparation in this situation: (1) using base R's table function; or (2) using functions available in the `tidyverse` group of packages.

For this example we use the low infant birth weight data presented by @hosmer_lemeshow:2000 and available in the `MASS` package in R. The `birthwt` data frame has 189 rows and 10 columns. The data were collected at Baystate Medical Center, Springfield, Massachusetts USA in 1986. 

## Two by two table preparation using the `table` function in base R

```{r, message = FALSE}
library(MASS)

# Load and view the data:
dat.df02 <- birthwt; head(dat.df02)
```
Each row of this data set provides details of a child born to a single mother. We're interested in the association between `smoke` (the mother's smoking status during pregnancy) and `low` (delivery of a baby less than 2.5 kg bodyweight). 

Its important that the table you present to `epi.2by2` is in the correct format: Outcome positives in the first column, outcome negatives in the second column, exposure positives in the first row and exposure negatives in the second row. If we run the `table` function on the `bwt` data the output table is in the incorrect format. Outcome negatives are listed first in the columns of the table and exposure negatives are listed first in the rows:

```{r, message = FALSE}
dat.tab02 <- table(dat.df02$smoke, dat.df02$low, dnn = c("Smoke", "Low BW")); dat.tab02
```
There are two ways to fix this problem. The quick fix is to simply ask R to switch the order of the rows and columns in the output table:

```{r, message = FALSE}
dat.tab02 <- table(dat.df02$smoke, dat.df02$low, dnn = c("Smoke", "Low BW")); dat.tab02
dat.tab02 <- dat.tab02[2:1,2:1]; dat.tab02
```
The second approach is to set the exposure variable and the outcome variable as a factor and to define the levels of each factor using `levels = c(1,0)`:

```{r, message = FALSE}
# Variables low, smoke and race as factors. Put an 'f' in front of the variable names to remind you that they're factors:

dat.df02$flow <- factor(dat.df02$low, levels = c(1,0))
dat.df02$fsmoke <- factor(dat.df02$smoke, levels = c(1,0))
dat.df02$frace <- factor(dat.df02$race, levels = c(1,2,3))

dat.tab02 <- table(dat.df02$fsmoke, dat.df02$flow, dnn = c("Smoke", "Low BW")); dat.tab02
```
With the table in the correct format we can compute the odds ratio for smoking and delivery of a low birth weight baby:

```{r, message = FALSE}
dat.epi02 <- epi.2by2(dat = dat.tab02, method = "cohort.count", elab = "Smoke", olab = "Low BW", digits = 2, conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
dat.epi02
```
The incidence odds of low BW among those mothers that were smoke positive was 2.0 (95% CI 1.1 to 3.8) times the incidence odds of low BW among those mothers that were smoke negative. All of the calculated measures of association computed by `epi.2by2` can be listed using:

```{r, message = FALSE}
names(dat.epi02$massoc.detail)
```
Compare the Wald confidence intervals for the odds ratio with the score confidence intervals:

```{r, message = FALSE}
dat.epi02$massoc.detail$OR.strata.wald
# Wald confidence intervals: 2.0 (95% CI 1.1 to 3.8)

dat.epi02$massoc.detail$OR.strata.score
# Score confidence intervals: 2.0 (95% CI 1.1 to 3.8)
```

## Two by two table preparation using `tidyverse`

The `dplyr` and `tidyr` packages (available within the `tidyverse` group of packages) can also be used to prepare your data in the required format:

```{r, message = FALSE}
library(dplyr); library(tidyr)

dat.df03 <- birthwt; head(dat.df03)

# Here we set the factor levels and tabulate the data in a single call using pipe operators:
dat.tab03 <- dat.df03 %>%
  mutate(flow = factor(low, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(fsmoke = factor(smoke, levels = c(1,0), labels = c("yes","no"))) %>%
  group_by(fsmoke, flow) %>%
  summarise(n = n()) 

# View the data:
dat.tab03

# View the data in conventional 2 by 2 table format:
pivot_wider(dat.tab03, id_cols = c(fsmoke), 
   names_from = flow, values_from = n)
```
As before, compute the odds ratio for smoking and delivery of a low birth weight baby:

```{r, message = FALSE}
dat.epi03 <- epi.2by2(dat = dat.tab03, method = "cohort.count", elab = "Smoke", olab = "Low BW", digits = 2, conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
dat.epi03
```

The incidence odds of low BW among those mothers that were smoke positive was 2.0 (95% CI 1.1 to 3.8) times the incidence odds of low BW among those mothers that were smoke negative.

## Measures of association for a series of candidate risk factors

The code that follows provides an approach for calculating measures of association for a series of candidate risk factors, writing the results of each analysis to a data frame. Here we set up a loop to calculate the odds ratio for age, smoke and race:

```{r, message = FALSE}
dat.df04 <- birthwt; head(dat.df04)

dat.df04$flow <- factor(dat.df04$low, levels = c(1,0))
dat.df04$fage <- ifelse(dat.df04$age > 23, 0,1)
dat.df04$fage <- factor(dat.df04$fage, levels = c(1,0))
dat.df04$fsmoke <- factor(dat.df04$smoke, levels = c(1,0))

# Race is coded 1 = white, 2 = black and 3 = other. Set white as the reference level:
dat.df04$frace <- ifelse(dat.df04$race == 1, 0, 1)
dat.df04$frace <- factor(dat.df04$frace, levels = c(1,0))

# Empty vectors to collect results:
rfactor <- ref <- or.est <- or.low <- or.upp <- c() 

# The candidate risk factors are in columns 12 to 14 of data frame dat.df04:
for(i in 12:14){
  tdat.tab04 <- table(dat.df04[,i], dat.df04$flow)
  tdat.epi04 <- epi.2by2(dat = tdat.tab04, method = "cohort.count", 
   digits = 2, conf.level = 0.95, units = 100, interpret = FALSE, outcome = "as.columns")
  
  trfactor <- as.character(names(dat.df04)[i]) 
  rfactor <- c(rfactor, trfactor) 
  
  tref <- as.character(paste("Reference: ", trfactor, " - ", levels(dat.df04[,i])[2], sep = ""))
  ref <- c(ref, tref)
  
  tor.est <- as.numeric(tdat.epi04$massoc.detail$OR.strata.wald[1])
  or.est <- c(or.est, tor.est)
  
  tor.low <- as.numeric(tdat.epi04$massoc.detail$OR.strata.wald[2])
  or.low <- c(or.low, tor.low)
  
  tor.upp <- as.numeric(tdat.epi04$massoc.detail$OR.strata.wald[3])
  or.upp <- c(or.upp, tor.upp)
}

gdat.df04 <- data.frame(ybrk = 1:3, ylab = rfactor, ref, or.est, or.low, or.upp)
gdat.df04
```

Plot the point estimate of the odds ratio for each candidate risk factor and its 95% confidence interval. Annotate the plot with the reference category for each risk factor. Note that a log 2 scale is used for the horizontal axis.

```{r massoc06, warning = FALSE, echo = TRUE, fig.cap = "\\label{fig:massoc06}Figure 6: Risk factors for low birth weight babies. Error bar plot showing the point estimate of the odds ratio and its 95% confidence interval for maternal age, smoking and race.", out.width = "80%", fig.align = "center"}
library(ggplot2); library(scales)

xbrk <- seq(from = -2, to = 2, by = 1)
xlab <- 2^xbrk

ggplot(data = gdat.df04, aes(x = log2(or.est), y = ybrk)) +
  theme_bw() +
  geom_point() + 
  geom_errorbar(aes(xmin = log2(or.low), xmax = log2(or.upp), width = 0.2)) + 
  scale_x_continuous(breaks = xbrk, labels = xlab, limits = range(xbrk), 
   name = "Odds ratio") + 
  scale_y_continuous(breaks = gdat.df04$ybrk, labels = gdat.df04$ylab, name = "Risk factor") + 
  geom_vline(xintercept = log2(1), linetype = "dashed") + 
  annotate("text", x = log2(0.25), y = gdat.df04$ybrk, label = gdat.df04$ref, hjust = 0, size = 3) +
  coord_fixed(ratio = 0.75 / 1) + 
  theme(plot.margin = unit(c(0,0,0,0), "mm")) +
  theme(axis.title.y = element_text(vjust = 0))
```

## Confounding

We're concerned that the mother's race may confound the association between low birth weight and delivery of a low birth weight baby so we'll stratify the data by race and compute the Mantel-Haenszel adjusted odds ratio. As before, our tables can be prepared using either base R or `tidyverse`.

### Stratified two by two table preparation using the table function in base R

```{r, message = FALSE}
dat.df05 <- birthwt; head(dat.df05)

dat.df05$flow <- factor(dat.df05$low, levels = c(1,0))
dat.df05$fsmoke <- factor(dat.df05$smoke, levels = c(1,0))
dat.df05$frace <- factor(dat.df05$race, levels = c(1,2,3))

dat.tab05 <- table(dat.df05$fsmoke, dat.df05$flow, dat.df05$frace, 
   dnn = c("Smoke", "Low BW", "Race")); dat.tab05
```

Compute the Mantel-Haenszel adjusted odds ratio for smoking and delivery of a low birth weight baby, adjusting for the effect of race. Function `epi.2by2` automatically returns the Mantel-Haenszel odds ratio and risk ratio when presented with stratified contingency tables.

```{r, message = FALSE}
dat.epi05 <- epi.2by2(dat = dat.tab05, method = "cohort.count", elab = "smoke",
   olab = "Low BW", digits = 2, conf.level = 0.95, units = 100, 
   interpret = FALSE, outcome = "as.columns")
dat.epi05
```

The Mantel-Haenszel test of homogeneity of the strata odds ratios is not significant ($\chi^{2}$ test statistic 2.800; df 2; p-value = 0.25). We accept the null hypothesis and conclude that the odds ratios for each strata of race are the same. Because the stratum specific odds ratios are not statistically significantly different the Mantel-Haenszel adjusted odds ratio provides an appropriate summary of the association between smoking and low birth weight.

After accounting for the confounding effect of race, the odds of having a low birth weight baby for smokers was 3.1 (95% CI 1.5 to 6.4) times that of non-smokers.

### Stratified two by two table preparation using tidyverse

```{r, message = FALSE}
dat.df06 <- birthwt

dat.tab06 <- dat.df06 %>%
  mutate(flow = factor(low, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(fsmoke = factor(smoke, levels = c(1,0), labels = c("yes","no"))) %>%
  mutate(frace = factor(race)) %>%
  group_by(frace, fsmoke, flow) %>%
  summarise(n = n()) 
dat.tab06

# View the data in conventional 2 by 2 table format:
pivot_wider(dat.tab06, id_cols = c(frace, fsmoke), 
   names_from = flow, values_from = n)
```

Compute the Mantel-Haenszel adjusted odds ratio for smoking and delivery of a low birth weight baby, adjusting for the effect of race:

```{r, message = FALSE}
dat.epi06 <- epi.2by2(dat = dat.tab06, method = "cohort.count", 
   elab = "Smoke", olab = "Low BW", digits = 2, conf.level = 0.95, 
   units = 100, interpret = FALSE, outcome = "as.columns")

dat.epi06
```

Plot the individual strata odds ratios and the Mantel-Haenszel summary odds ratio as an error bar plot to better understand how the Mantel-Haenszel adjusted odds ratio relates to the individual strata odds ratios:

```{r massoc07, warning = FALSE, echo = TRUE, fig.cap = "\\label{fig:massoc07}Figure 7: Risk factors for low birth weight babies. Error bar plot showing the odds of having a low birth weight baby for smokers of maternal race categories 1, 2 and 3 and the Mantel-Haenszel odds of having a low birth weight baby for smokers, adjusted for maternal race.", out.width = "80%", fig.align = "center"}
xbrk <- seq(from = -5, to = 5, by = 1)
xlab <- 2^xbrk

nstrata <- dat.epi06$n.strata
ybrk <- c(1:nstrata, max(nstrata) + 1)
ylab <- c("M-H", paste("Strata ", 1:nstrata, sep = ""))

or.est <- c(dat.epi06$massoc.detail$OR.mh$est, 
   dat.epi06$massoc.detail$OR.strata.cfield$est)
or.low <- c(dat.epi06$massoc.detail$OR.mh$lower, 
   dat.epi06$massoc.detail$OR.strata.cfield$lower)
or.upp <- c(dat.epi06$massoc.detail$OR.mh$upper, 
   dat.epi06$massoc.detail$OR.strata.cfield$upper)
gdat.df06 <- data.frame(ybrk, ylab, or.est, or.low, or.upp)

ggplot(data = gdat.df06, aes(x = log2(or.est), y = ybrk)) +
  theme_bw() +
  geom_point() + 
  geom_errorbar(aes(xmin = log2(or.upp), xmax = log2(or.low), width = 0.2)) + 
  scale_x_continuous(breaks = xbrk, labels = xlab, limits = range(xbrk), 
   name = "Odds ratio") + 
  scale_y_continuous(breaks = gdat.df06$ybrk, labels = gdat.df06$ylab, name = "Risk factor") + 
  geom_vline(xintercept = log2(1), linetype = "dashed") + 
  coord_fixed(ratio = 0.75 / 1) +
  theme(plot.margin = unit(c(0,0,0,0), "mm")) +
  theme(axis.title.y = element_text(vjust = 0))
```

## References

---
nocite: '@*'
---
