R4SME 11 logistic interaction.Rmd

---
title: "R for SME 11: Logistic regression (interaction)"
author: Andrea Mazzella [link](https://github.com/andreamazzella)
output: html_notebook
---

## Basics
Load packages
```{r}
library(haven)
library(epiDisplay)
library(magrittr)
library(tidyverse)
```

# Part 1: mortality

## Data import, exploration, management

Make sure you have the mortality.dta dataset in the same folder as this .rmd

```{r Import data}
mortality<-read_dta("./mortality.dta")
mortality %<>% mutate_if(is.labelled,as_factor)
```

```{r Explore data}
glimpse(mortality)
summary(mortality)
View(mortality)
```

Let's change the outcome into a labelled factor.
```{r Factors and labelling}
mortality %<>%
  mutate(died = factor(
    died,
    levels = c(0, 1),
    labels = c("alive", "dead")
  ))
```

## Data analysis

1.

Analyse the association between visual impairment and death, stratifying on sex, using a Mantel-Haenszel approach.
```{r}
# 2x2 table
mortality %$% tabpct(died, vimp, percent = "col", graph = F)

# Crude OR
mortality %$% cc(died, vimp, graph = F)

# Stratified tables
mortality %$% table(died, vimp, sex)

# Stratum-specific OR and MH-OR
mortality %$% mhor(died, vimp, sex, graph = F)
```
Crude OR: 5.57 (3.78-8.2)

Stratum-specific OR (CI)
- male   3.94       2.15       6.98
- female 8.05       4.41      14.29

MH OR  5.43       3.69       8.00

Homogeneity test p = 0.07: there is weak evidence of interaction between sex and visual impairment; the OR in females is about twice the OR in males.

2.

Fit logistic regression models to estimate the same association, without interaction and with an interaction.
NB: in this output for the model with interaction the interpretation of adj OR changes - it's a *stratum-specific* OR (in the baseline stratum of the other covariate)
```{r}
# Model without interaction
glm(died ~ vimp + sex,
    data = mortality,
    family = binomial()) %>%
  logistic.display()

# Model with interaction
glm(died ~ vimp * sex, # the asterisk marks the interaction
    data = mortality,
    family = binomial()) %>%
  logistic.display()
```
- What is the OR for visual impairment adjusted for sex?
Model without interaction: Adj OR (died/vimp//sex) = 5.53 (3.75,8.16)  (Wald's p < 0.001)

Now use the above output to calculate, by hand, the odds and OR in the four groups.

- What is the OR for visual impairment *in males*?
3.95

- What is the OR for females (vs males) *in the visually unimpaired*?
0.77

- What is the interaction term?
2.04

- What's the OR for visual impairment *in females*?
_Issue:_ the logistic.display() output doesn't show the equivalent of STATA's "_cons", so you can't calculate the table for OR in the four groups. _Workaround_: you let R do the work, see below (Q3).

- Compare answers with MH analysis above.

3.

Calculate the other stratum-specific OR (for females) by making females the baseline group.
```{r}
mortality$sex2 <- factor(mortality$sex,
                         levels = c("Female", "Male")) # girl power
glm(died ~ vimp * sex2,
    data = mortality,
    family = binomial()) %>%
  logistic.display()
mortality$sex2 <- factor(mortality$sex,
                         levels = c("Male", "Female")) # the patriarchy strikes back
```


4.

Fit a logistic regression model with interaction between visual impairment and sex, and control for age
```{r}
glm(died ~ vimp * sex + agegrp,
    data = mortality,
    family = binomial()) %>%
  logistic.display()
```
- OR for visual impairment in males: ???
- OR for females in visually unimpaired: 0.85
- Interaction parameter between vimp and sex: 1.93
- OR for visual impairment in females: ??
- Impact on controlling for age: 

5.

I'm not sure what the equivalent to Stats's lincom is. I'm also not sure why we need it, as one can calculate all stratum-specific ORs by changing the baseline group.

# Part 2: Mwanza
## Data import, exploration & management

6.

Make sure you have the mwanza.dta dataset in the same folder as this .rmd, and load it. It contains data on HIV infection among women in Mwanza, Tanzania.
```{r}
# Import the dataset
mwanza <- read_dta("./mwanza.dta")
```

```{r}
#Familiarise yourself with the data
View(mwanza)
glimpse(mwanza)
summary(mwanza)
```


Use logistic regression to assess whether the association between education (ed) and HIV (case) is modified by age (age1). First, relevel these variables:
- ed into two groups: "none" and "any formal education";
- age1 into three groups: 15-24, 25-34, 34+.
```{r}
# Relevel "ed" and "age1"
mwanza %<>%
  mutate(ed2 = as.factor(
    case_when(
      ed == 1 ~ "none",
      TRUE ~ "any formal"
    )
  ))
mwanza %$% table(ed, ed2)

mwanza %<>%
  mutate(age3 = as.factor(
    case_when(
      age1 == 1 ~ "15-24",
      age1 == 2 ~ "15-24",
      age1 == 3 ~ "25-34",
      age1 == 4 ~ "25-34",
      TRUE ~ "34+"
    )
  ))
mwanza %$% table(age1, age3)
```
Now use MH and logistic regression to assess for interaction.
Is there evidence of interaction between age and education?
```{r}
# Crude OR
mwanza %$% cc(case, ed2, graph = F)

# Stratum-specific OR and MH-OR
mwanza %$% mhor(case, ed2, age3, graph = F)

# Logistic regression without interaction
logit_without <- glm(case ~ ed2 + age3,
                     data = mwanza,
                     family = binomial())

# Logistic regression with interaction
logit_interact <- glm(case ~ ed2 * age3,
                      data = mwanza,
                      family = binomial())
lrtest(logit_without, logit_interact)
```


Crude OR: 0.41 (0.29-0.60) p < 0.001
MH OR:    0.43 (0.29-0.63) p < 0.001
Homogeneity test, p-value = 0.006. There is strong evidence for interaction between age and education.
Interaction parameter: 0.2 in age group 25-34, 0.41 in age group 34+.


Calculate age-specific ORs for the association between education and HIV.

```{r}
logistic.display(logit_interact)
```
- In the 15-24 age group, odds of HIV given no education = 1.03 odds of HIV given any formal education.

```{r}
# 25-34
mwanza$age3 <- factor(mwanza$age3,
                      levels = c("25-34", "34+", "15-24"))
glm(case ~ ed2 * age3,
    data = mwanza,
    family = binomial()) %>% logistic.display()

# 34+
mwanza$age3 <- factor(mwanza$age3,
                      levels = c("34+", "15-24", "25-34"))
glm(case ~ ed2 * age3,
    data = mwanza,
    family = binomial()) %>% logistic.display()

# Change back
mwanza$age3 <- factor(mwanza$age3,
                      levels = c("15-24", "25-34", "34+"))
```
- In the 25-34 age group, odds of HIV given no education = 0.21 odds of HIV given any formal education.
- In the 35+ age group, odds of HIV given no education = 0.42 odds of HIV given any formal education.