Multilevel Quantile Regression

Generate data

Show the code

n <- 12
id <- 1:n
ln <- 1:2
roi <- 1:3
rep <- 1:50
d <- expand.grid(rep = rep, roi = roi, ln = ln, id = id)
d$trt <- ifelse(d$id < 7, 1, 0)

d$id_ln <- as.numeric(interaction(d$id, d$ln))
d$id_ln_roi <- as.numeric(interaction(d$id, d$ln, d$roi))
set.seed(1)
e <- rnorm(n = nrow(d), sd = ifelse(d$trt == 0,0.5,2*0.5))
e_id <- rnorm(n = n, sd = 0.8)
e_id_ln <- rnorm(n = n*2, sd = 0.5)
e_id_ln_roi <- rnorm(n = n*2*3, sd = 0.7)
d$y <- 4 + 0.4*(d$trt == 1) + 
  e_id[d$id] + 
  e_id_ln[d$id_ln] +
  e_id_ln_roi[d$id_ln_roi] +
  e
d$id_ln <- NULL
d$rep <- NULL
d$trt <- factor(d$trt)
d$roi <- factor(d$roi)
d <- d[,c("id", "ln", "roi", "trt", "y", "id_ln_roi")]

Inspect and visualize data

• id = subject id
• ln = lymph node
• roi = region of interest
• trt = treatment indicator (0 = control/1 = treated)
• y = dependent variable
• id_ln_roi = grouping indicator (id/ln/roi groups)

head(d)

  id ln roi trt        y id_ln_roi
1  1  1   1   1 4.964657         1
2  1  1   1   1 5.774754         1
3  1  1   1   1 4.755482         1
4  1  1   1   1 7.186392         1
5  1  1   1   1 5.920619         1
6  1  1   1   1 4.770643         1

Plot of y versus treatment group, colored by roi.

library(ggplot2)
ggplot(d) +
  aes(x = trt, y = y, color = roi) +
  geom_jitter(width = 0.2, height = 0)

Summary of y grouped by trt.

tapply(d$y, d$trt, summary)

$`0`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
0.02821 2.73509 3.82543 3.67880 4.73462 6.79213 

$`1`
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-0.4181  3.2678  4.2715  4.2853  5.2978  8.8127 

Multilevel Quantile Model

Model the 0.75 quantile (ie, 75th percentile) as a function of treatment, with a random intercept for the id/ln/roi group. To model a different quantile, change the tau argument.

library(lqmm)
fit.lqmm <- lqmm(y ~ trt, random = ~ 1 , 
                 group = id_ln_roi, 
                 data = d, tau = 0.75)
summary(fit.lqmm)

Call: lqmm(fixed = y ~ trt, random = ~1, group = id_ln_roi, tau = 0.75, 
    data = d)

Quantile 0.75 

Fixed effects:
              Value Std. Error lower bound upper bound Pr(>|t|)    
(Intercept) 3.94630    0.33451     3.27408      4.6185 6.31e-16 ***
trt1        1.06129    0.43329     0.19056      1.9320  0.01793 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
AIC:
[1] 9782 (df = 4)

We estimate the 75th quantile of treatment is about 1.06 higher than the 75th quantile of the control group.