## 7.12 Interactions

As with linear and logistic regression, include interactions in a Cox model to assess effect modification. Including an interaction allows you to assess if the association between a risk factor and time to an event depends on another variable, and to estimate the HR for one variable at different levels of another.

**Example 7.7:** Using the Natality teaching dataset, see if the association between previous preterm birth and time to preterm birth (adjusted for age, race/ethnicity, and marital status) depends on whether or not the mother had pre-pregnancy hypertension (`RF_PHYPE`

). As with any model with an interaction, include both the main effects and the interaction term.

```
cox.ex7.7.int <- coxph(Surv(gestage37, preterm01) ~ MAGER + MRACEHISP + DMAR +
RF_PPTERM + RF_PHYPE + # Main effects
RF_PPTERM:RF_PHYPE, # Interaction
data = natality)
round(
cbind("AHR" = exp(summary(cox.ex7.7.int)$coef[, "coef"]),
exp(confint(cox.ex7.7.int)),
"p-value" = summary(cox.ex7.7.int)$coef[, "Pr(>|z|)"])
, 3)
```

```
## AHR 2.5 % 97.5 % p-value
## MAGER 1.032 1.008 1.057 0.009
## MRACEHISPNH Black 1.699 1.194 2.416 0.003
## MRACEHISPNH Other 0.929 0.520 1.657 0.802
## MRACEHISPHispanic 1.304 0.918 1.852 0.139
## DMARUnmarried 1.785 1.320 2.413 0.000
## RF_PPTERMYes 2.660 1.634 4.328 0.000
## RF_PHYPEYes 1.287 0.526 3.148 0.580
## RF_PPTERMYes:RF_PHYPEYes 6.774 1.195 38.392 0.031
```

The interaction is statistically significant (the p-value for the `RF_PPTERMYes:RF_PHYPEYes`

row is .031). Had either of the terms in the interaction been categorical with more than two levels, we would have used `car::Anova(cox.ex7.7.int, type = 3, test = "Wald")`

to get the interaction p-value as it would have been a multiple degree of freedom test.

### 7.12.1 Overall test of a predictor involved in an interaction

As shown in Section 5.9.11, we can also carry out an overall test of previous preterm birth by comparing the full model to the model with both the main effect (`RF_PPTERM`

) and interaction (`RF_PPTERM:RF_PHYPE`

) removed. We did not remove missing data prior to fitting the full model, so we need to create a dataset with no missing data first so both models can be fit to the same observations, or use the `subset`

option. As with MLR, use `anova()`

to compare the models but for a Cox regression we must specify `test = "Chisq"`

to get a p-value. There is no option to specify a Wald test in this case; `anova()`

uses a likelihood ratio test when comparing Cox models.

```
# Fit the reduced model using the subset option
# to remove missing values for the variable
# being removed to create the reduced model
fit0 <- coxph(Surv(gestage37, preterm01) ~
MAGER + MRACEHISP + DMAR +
RF_PHYPE,
data = natality,
subset = complete.cases(RF_PPTERM))
# Compare to full model
anova(fit0, cox.ex7.7.int, test = "Chisq")
```

```
## Analysis of Deviance Table
## Cox model: response is Surv(gestage37, preterm01)
## Model 1: ~ MAGER + MRACEHISP + DMAR + RF_PHYPE
## Model 2: ~ MAGER + MRACEHISP + DMAR + RF_PPTERM + RF_PHYPE + RF_PPTERM:RF_PHYPE
## loglik Chisq Df Pr(>|Chi|)
## 1 -1587
## 2 -1578 19.1 2 0.000071 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
```

The 2 df test (2 df because the models differ in two terms – the main effect and the interaction) results in a p-value that is very small. Thus, after adjusting for age, race/ethnicity, marital status, and pre-pregnancy hypertension, previous preterm birth is significantly associated with time to preterm birth.

### 7.12.2 Estimating the HR at each level of the other variable

In addition to the test of significance, we also would like to know what is the magnitude of the AHR for previous preterm birth among those with and without pre-pregnancy hypertension. The main effect for `RF_PPTERM`

provides the AHR for previous preterm birth among those at the reference level of `RF_PHYPE`

(“No”). For linear and logistic regression models, we used `gmodels::estimable()`

to estimate the AHR at non-reference levels of the other predictor in the interaction. Unfortunately, `gmodels::estimable()`

does not work for a `coxph`

object. Instead, refit the model after changing the reference level of `RF_PHYPE`

to “Yes” which will result in a main effect for `RF_PPTERM`

at the new reference level. In this way, we can get the AHR, 95% CI, and p-value for `RF_PPTERM`

at each level of `RF_PHYPE`

.

First, using the original model, extract the results for the main effect `RF_PPTERM`

, which is the effect at `RF_PHYPE`

= “No” (the reference level for `RF_PHYPE`

).

```
round(
cbind("AHR" = exp(summary(cox.ex7.7.int)$coef[, "coef"]),
exp(confint(cox.ex7.7.int)),
"p-value" = summary(cox.ex7.7.int)$coef[, "Pr(>|z|)"])["RF_PPTERMYes", ]
, 3)
```

```
## AHR 2.5 % 97.5 % p-value
## 2.660 1.634 4.328 0.000
```

Next, re-level `RF_PHYPE`

(not `RF_PPTERM`

), re-fit the model, and extract the main effect results again to get the `RF_PPTERM`

effect at `RF_PHYPE`

= “Yes”.

```
natality_b <- natality %>%
mutate(RF_PHYPE = relevel(RF_PHYPE, ref = "Yes"))
# Refit with data = natality_b
cox.ex7.7.int_b <- coxph(Surv(gestage37, preterm01) ~
MAGER + MRACEHISP + DMAR +
RF_PPTERM + RF_PHYPE +
RF_PPTERM:RF_PHYPE,
data = natality_b)
round(
cbind("AHR" = exp(summary(cox.ex7.7.int_b)$coef[, "coef"]),
exp(confint(cox.ex7.7.int_b)),
"p-value" = summary(cox.ex7.7.int_b)$coef[, "Pr(>|z|)"])["RF_PPTERMYes",]
, 3)
```

```
## AHR 2.5 % 97.5 % p-value
## 18.016 3.406 95.289 0.001
```

The AHR for previous preterm birth is greater among those with pre-pregnancy hypertension (18.016) than among those without (2.660) (interaction p-value = .031, a test of the null hypothesis that these two AHRs are the same).

The AHR among those with pre-pregnancy hypertension is quite large and has a very wide CI. Why? It turns out that there were not many cases with pre-pregnancy hypertension, resulting in high variability when estimating the AHR in this subgroup and a wide CI.

```
##
## No Yes
## 1965 34
```

Also, among those with pre-pregnancy hypertension, only two mothers had a previous preterm birth and both experienced a preterm birth, which explains the large `RF_PPTERM`

AHR in this group.

```
# Previous preterm vs. preterm
# among those with hypertension
table(natality$RF_PPTERM[natality$RF_PHYPE == "Yes"],
natality$preterm01[natality$RF_PHYPE == "Yes"])
```

```
##
## 0 1
## No 27 5
## Yes 0 2
```

You may have noticed that there is a 0 in this table. With logistic regression, a 0 in the two-way table comparing a categorical predictor to the outcome was an indication that logistic regression would have convergence problems due to separation (Section 6.10). In Section 7.13 we will discuss this same issue in the context of a Cox regression model. It turns out this particular 0 does not cause a problem with separation, but others will.

### 7.12.3 Visualizing an interaction

Plot the survival functions illustrating the previous preterm birth effect among those without and with pre-pregnancy hypertension. The code is similar to that in Section 7.11, but now we need a `data.frame`

for each combination of levels of `RF_PPTERM`

and `RF_PHYPE`

.

```
# RF_PPTERM = "Yes", RF_PHYPE = "Yes"
DAT_PPTERM_Yes_PHYPE_Yes <- data.frame(
preterm01 = 1,
gestage37 = 17:36,
DMAR = "Married",
MRACEHISP = "NH White",
MAGER = mean(natality.complete$MAGER),
RF_PPTERM = "Yes",
RF_PHYPE = "Yes")
# Other 3 combinations
# Copy first data.frame
DAT_PPTERM_Yes_PHYPE_No <- DAT_PPTERM_Yes_PHYPE_Yes
DAT_PPTERM_No_PHYPE_Yes <- DAT_PPTERM_Yes_PHYPE_Yes
DAT_PPTERM_No_PHYPE_No <- DAT_PPTERM_Yes_PHYPE_Yes
# Update values of the 2 terms in the interaction
DAT_PPTERM_Yes_PHYPE_No$RF_PPTERM <- "Yes"
DAT_PPTERM_Yes_PHYPE_No$RF_PHYPE <- "No"
DAT_PPTERM_No_PHYPE_Yes$RF_PPTERM <- "No"
DAT_PPTERM_No_PHYPE_Yes$RF_PHYPE <- "Yes"
DAT_PPTERM_No_PHYPE_No$RF_PPTERM <- "No"
DAT_PPTERM_No_PHYPE_No$RF_PHYPE <- "No"
```

When computing the predictions, be careful to use the original fit for the combinations with `RF_PHYPE`

= “No” (the original reference level of the other term in the interaction) and the re-leveled fit for `RF_PHYPE`

= “Yes” (the reference level after re-leveling).

```
# Predictions
# Use re-leveled fit when RF_PHYPE = "Yes" (ref level after re-leveling)
PRED_PPTERM_Yes_PHYPE_Yes <-
predict(cox.ex7.7.int_b, DAT_PPTERM_Yes_PHYPE_Yes, type = "survival")
# Use original fit when RF_PHYPE = "No" (original ref level)
PRED_PPTERM_Yes_PHYPE_No <-
predict(cox.ex7.7.int, DAT_PPTERM_Yes_PHYPE_No, type = "survival")
# Use re-leveled fit when RF_PHYPE = "Yes" (ref level after re-leveling)
PRED_PPTERM_No_PHYPE_Yes <-
predict(cox.ex7.7.int_b, DAT_PPTERM_No_PHYPE_Yes, type = "survival")
# Use original fit when RF_PHYPE = "No" (original ref level)
PRED_PPTERM_No_PHYPE_No <-
predict(cox.ex7.7.int, DAT_PPTERM_No_PHYPE_No, type = "survival")
```

Finally, plot the predictions.

```
# RF_PHYPE = "No"
par(mfrow=c(1,2))
plot(0, 0, col = "white",
xlab = "Gestational Age (weeks)", ylab = "P(Not yet preterm)",
xlim = c(17, 36),
ylim = c(0, 1),
font.axis = 2, font.lab = 2,
main = "Without hypertension")
# NOTE: For lines() the arguments are x, y not y ~ x
lines( DAT_PPTERM_No_PHYPE_No$gestage37,
PRED_PPTERM_No_PHYPE_No, lwd = 2, lty = 1)
lines( DAT_PPTERM_Yes_PHYPE_No$gestage37,
PRED_PPTERM_Yes_PHYPE_No, lwd = 2, lty = 2)
legend(17, 0.55, c("No", "Yes"), title = "Previous preterm birth",
lty=c(1,2), lwd=c(2,2), bty = "n")
# RF_PHYPE = "Yes"
plot(0, 0, col = "white",
xlab = "Gestational Age (weeks)", ylab = "P(Not yet preterm)",
xlim = c(17, 36),
ylim = c(0, 1),
font.axis = 2, font.lab = 2,
main = "With hypertension")
# NOTE: For lines() the arguments are x, y not y ~ x
lines(DAT_PPTERM_No_PHYPE_Yes$gestage37,
PRED_PPTERM_No_PHYPE_Yes, lwd = 2, lty = 1)
lines(DAT_PPTERM_Yes_PHYPE_Yes$gestage37,
PRED_PPTERM_Yes_PHYPE_Yes, lwd = 2, lty = 2)
legend(17, 0.55, c("No", "Yes"), title = "Previous preterm birth",
lty=c(1,2), lwd=c(2,2), bty = "n")
```

It is a good idea to compare the visualization of an interaction to the estimated AHR at each level of the other variable in the interaction. If they are not consistent, then the estimation and/or plotting code have an error. In Section 7.12.2, we estimated that the AHR for `RF_PPTERM`

at `RF_PHYPE`

= “No” is 2.660 and at `RF_PHYPE`

= “Yes” is 18.016. These are consistent with Figure 7.15. Both AHRs are greater than 1 and in both plots the survival curve for `RF_PPTERM`

= “Yes” is lower than that for “No”. Also, the AHR for `RF_PPTERM`

at `RF_PHYPE`

= “Yes” is much larger than at `RF_PHYPE`

= “No” and the lines in the right panel are much further apart than those in the left panel.

The right panel in Figure 7.15 shows that, based on this dataset, we estimate that those with both a previous preterm birth and pre-pregnancy hypertension (and other predictors as specified above) have over 85% probability of preterm birth (the dashed survival function drops below 0.15). Be wary of this estimate, however, since, as seen above, it is based on a very small sample size. Using the methods from Section 7.10, compute a 95% CI for this quantity to get an idea of how imprecise it is.

```
unlist(summary(survfit(cox.ex7.7.int_b,
DAT_PPTERM_Yes_PHYPE_Yes[DAT_PPTERM_Yes_PHYPE_Yes$gestage37 == 36,],
se.fit =T, conf.int = 0.95),
times=36)[c("surv", "std.err", "lower", "upper")])
```

```
## surv std.err lower upper
## 0.141915 0.204637 0.008407 1.000000
```

**Conclusion:** After adjusting for mother’s age, race/ethnicity, marital status, and pre-pregnancy hypertension, previous preterm birth is significantly associated with time to preterm birth (p < .001 – from the test comparing the model with an interaction to the model with no `RF_PPTERM`

main effect or interaction). The association between previous preterm birth and time to preterm birth differs significantly between those without and with pre-pregnancy hypertension (interaction p = .031). Among mothers who did not have pre-pregnancy hypertension, those with a previous preterm birth have 2.7 times the hazard of preterm birth of those who did not (AHR = 2.66; 95% CI = 1.63, 4.33; p < .001). Among those who did have pre-pregnancy hypertension, the hazard is even larger (AHR = 18.0; 95% CI = 3.4, 95.3; p < .001).