library(data.table)
library(ggplot2)
gg_color_hue <- function(n) {
  hues = seq(15, 375, length = n + 1)
  hcl(h = hues, l = 65, c = 100)[1:n]
}
cols2 = gg_color_hue(2)

Introduction

By conditional, we could be referring to any multivariable linear predictor. However, here I am specifically looking at a random intercept model and integrating out the random effect to yield marginal estimates of a treatment effect.

Consider a two arm trial where we collect several measurements of a binary variable of interest. For a contrived and silly example, perhaps we are investigating the psychic effects of a new drug and the participants are asked to guess the pictures that we are looking at but that they can only see through their minds eye. For each picture, the assessor records the participants response as correct (1) or incorrect (0). We incorporate all the data into the model via a person level random intercept, i.e.

\[ \begin{equation} \begin{split} \mathsf{Pr}(Y_{ij} = 1 | u_i) = \beta_0^\prime + u_i + \beta_1^\prime x_{ij} \end{split} \tag{1} \end{equation} \]

We simulate data from this experiment with parameters \(\beta_0^\prime = -2\), \(\beta_1^\prime = 0.4\) and \(\mathsf{Var}(u) = 2\) as follows

set.seed(1)
N <- 600
m <- 10
sdu <- sqrt(2)
id <- 1:N
id <- rep(id, each = m)
trt <- rbinom(N, 1, prob = 0.5)
trt <- rep(trt, each = m)
u <- rnorm(N, 0, sdu)
u <- rep(u, each = m)
b <- c(-2, 0.4)
d1 <- data.table(
  id = id, 
  trt = trt,
  u = u
)
d1[, eta := b[1] + u + b[2]* trt]
d1[, eta0 := b[1] + u]
d1[, eta1 := b[1] + u + b[2]]
d1[, y := rbinom(.N, 1, plogis(eta))]
d1[trt == 1, treatment := "Active"]
d1[trt == 0, treatment := "Control"]

This set of parameters imply that a participant with \(u_i = 0\) has a 0.12 probability of success when in the control arm and a 0.17 probability of success when in the active arm.

The following plot shows the probability of success obtained from 600 simulated participants assuming that we could see their probability of success under both the control and active arms simultaneously. The dashed lines showing the risks for a participant with a zero value from the random effect for each treatment group.

Note that while the odds ratio is, by definition, constant across the whole sample, the participant level change in risk varies depending on their baseline risk, which is a sum of the intercept term and the cohort heterogeneity represented by the random effect.

d2 <- unique(d1[, .(id, u, p0 = plogis(eta0), p1 = plogis(eta1))][order(u)])
d2 <- melt(d2, id.vars = c("id", "u"))
d2[variable == "p0", treatment := "Control"]
d2[variable == "p1", treatment := "Active"]
d2[, variable := NULL]
ggplot(d2, aes(x = u, y = value, group = treatment, col = treatment)) +
  geom_line() +
  geom_hline(yintercept = plogis(b[1]), col = cols2[2], lty = 2) +
  geom_hline(yintercept = plogis(sum(b)), col = cols2[1], lty = 2) +
  geom_vline(xintercept = 0, lwd = 0.25) +
  scale_y_continuous("Probability of success")+
  scale_x_continuous("u", limits = c(-4, 4)) +
  theme_bw() + 
  theme(legend.pos = "bottom")
## Warning: Removed 8 row(s) containing missing values (geom_path).

The marginal risk can be computed by averaging out the random effect per

\[ \begin{equation} \begin{split} \mathsf{Pr}(Y_{ij} = 1) = \int \frac{\exp(\beta_0^\prime + u_i + \beta_1^\prime)}{1 + \exp(\beta_0^\prime + u_i + \beta_1^\prime)} f(u_i; \nu^2) \mathsf{du} \end{split} \tag{2} \end{equation} \]

As we have assumed a normal random effect with zero mean and variance 2 we can approximate this quantity to an arbitrary level of precision via random sampling.

u <- rnorm(100000, 0, sqrt(2))
b0_marg <- exp(b[1] + u)/ ( 1 + exp(b[1] + u) )
b1_marg <- exp(b[1] + u + b[2])/ ( 1 + exp(b[1] + u + b[2]) )
c(mean(b0_marg), mean(b1_marg))
## [1] 0.1836537 0.2344197

These quantities imply a marginal log-odds of response in the control group equal to -1.49 and a log-odds-ratio for the effect of treatment equal to 0.31.

Thus in the marginal and random effects models, the parameters differ when using logistic regression. The former describes the ratio of population odds; the latter describes the ratio of an individual’s odds. The conditional parameter values will be larger in absolute terms than their random effect analogues, but (typically) the conditional parameter will also be less precise.

I reiterate - the marginal treatment effect relates to the effect on the entire population whereas the conditional treatment effect relates to the estimated effect on an individual.

What this means is that if the conditional OR associated with treatment is 1.4 then you can say that giving this intervention will increase the odds of the event in the person sitting in front of you by a factor of 1.4 relative to the control. However, for a marginal OR 1.4, you could say that if you gave the entire population the treatment then you will increase the odds of the event in the population by a factor of 1.4 relative to the control. Arguably, the former is of more interest to the individual receiving treatment whereas the latter is of more interest to public health professionals, planners and policy makers.

Note that this characteristic property arises both with odds-ratios and hazard ratios, but it does not apply for linear regression or risk ratios (log link).

For further information, see section 7.4 of Diggle’s text on the Analysis of Longitudinal Data.