How to compute custom contrasts for linear quantile mixed model from lqmm

Question

I would like to get custom pairwise contrasts and Holm adjustment for a linear quantile mixed model generated using the lqmm and glht functions, but my attempt generates an error. For comparison, I show that this is possible with a linear mixed model using lmer.

## Load necessary packages
library(lqmm); library(lme4); library(multcomp)

## Orthodont data
data(Orthodont) # head(Orthodont) str(Orthodont)
Orthodont$age = as.factor(Orthodont$age)

# Random intercept model
fit_lmm <- lmer(distance ~ (1 | Subject ) + age, data = Orthodont) # summary(fitOi.lmm)
fit_lqmm <- lqmm(distance ~ age, random = ~ 1, group = Subject,
                   tau = c(0.1,0.5,0.9), data = Orthodont) # summary(fitOi.lqmm)

###  define contrast matrix 
contr <- rbind("10 - 8" = c(-1, 1, 0, 0),
               "12 - 8" = c(-1, 0, 1, 0), 
               "12 - 10" = c(0, -1, 1, 0))
contrast_lmm = summary(glht(fit_lmm, linfct = mcp(age = contr), test="holm")) # this works
contrast_lmm

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: User-defined Contrasts


Fit: lmer(formula = distance ~ (1 | Subject) + age, data = Orthodont)

Linear Hypotheses:
             Estimate Std. Error z value Pr(>|z|)    
10 - 8 == 0    0.9815     0.3924   2.501 0.033067 *  
12 - 8 == 0    2.4630     0.3924   6.277  < 1e-04 ***
12 - 10 == 0   1.4815     0.3924   3.776 0.000472 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
(Adjusted p values reported -- single-step method)

The following line of code generates an error.

contrast_lmm = summary(glht(fit_lqmm, linfct = mcp(age = contr), test="holm")) 

Error in terms.default(object) : no terms component nor attribute
Error in factor_contrasts(model) : 
  no ‘model.matrix’ method for ‘model’ found!

EDIT 1

The emmeans package also does not work with lqmm models, and the qdrg function does not work to create a reference grid. 

EDIT 2

Here I try out and expand on the suggestion on 6 Apr 2020 by Russ Lenth using the emmeans package.

dummy = 1
fit_lqmm_slopes = as.data.frame(t(c(dummy=dummy, fit_lqmm$theta_x[2:nrow(fit_lqmm$theta_x),"0.5"])))

fit_lqmm_intcpts_df = data.frame(dummy=dummy, intercept = fit_lqmm$y) # head(fit_lqmm_intcpts_df)

fit_lqmm_coef_mat = data.matrix(merge(fit_lqmm_intcpts_df,fit_lqmm_slopes)[,2:(ncol(fit_lqmm_slopes)+1)]) # head(fit_lqmm_coef_mat) # coef(fit_lmm)

fit_lqmm_qdrg = qdrg(formula = ~age, vcov = VarCorr(fit_lqmm)$'0.5', coef = fit_lqmm_coef_mat, df = fit_lqmm$rdf, data = Orthodont) 

This gives the error message:

Error in ref_grid(result, ...) : 
  Non-conformable elements in reference grid.
 Probably due to rank deficiency not handled as expected.

I am concerned that the variance-covariance entry (vcov argument) is not specified correctly. So far I can only extract the variance-covariance "matrix" for the random effect using lqmm:VarCorr, which yields a single value for each quantile. I cannot find a way to extract the full variance-covariance matrix from the lqmm object to match the structure of vcov(fit_lmm). I constructed the coefficient data frame to match the structure given by coef(fit_lmm). Should the formula include the random intercept (1 | Subject)?

Answer 1

Regarding emmeans::qdrg, it is noninformative to merely say "it doesn't work" because we have no way of knowing what was tried.

My suggestion would be to try harder to get qdrg() to work. It would not be at all surprising if its simpler use, with just object and data arguments, would fail, given that apparently such basic things as a terms method don't even appear to be implemented. But it seems likely that you can get it to work using something like

rq <- qdrg(formula = ~ age, data = Orthodont, coef = ???, vcov = ???, 
           df = fit_lqmm$rdf)
contrast(rg, as.data.frame(t(contr)))

where you somehow get the needed coef vector and vcov matrix from the fitted model. You may well need to extract just the parts of the coefficients and covariance matrix that pertain to the fixed effects (in this case the intercept and age coefficients). It's hard to believe it wouldn't be possible to get this information from the fitted model. It appears from the documentation for lqmm that this information can be extracted from fit_lqmm$theta (or from theta_x and theta_z?)

[Note that as.data.frame(t(contr)) turns your custom contrast coefficients into a data frame with the contrast coefficients as columns, as needed by emmeans::contrast().]

Addendum

It turns out to be tricky to get the vcov part. It requires a call to summary.lqmm. Here is code that works for your example:

require(emmeans)
require(lqmm)

bhat = coef(fit_lqmm)
vcv = summary(fit_lqmm, covariance = TRUE) $ Cov
rg0.5 = qdrg(formula = ~ age, data = Orthodont, df = fit_lqmm$rdf,
    coef = bhat[, "0.5"], vcov = vcv[1:4, 1:4, "0.5"])
cont = list("10 - 8"   = c(-1, 1, 0, 0),
             "12 - 8"  = c(-1, 0, 1, 0), 
             "12 - 10" = c(0, -1, 1, 0))

Note: The covariance matrices include some other parameters besides the coefficients, so I had to extract the ones that correspond only to the coefficienrts (elements 1--4 in this example).

Then we get:

> confint(rg0.5)
 age prediction    SE  df lower.CL upper.CL
 8         22.0 0.487 103     21.0     23.0
 10        22.9 0.479 103     22.0     23.9
 12        24.0 0.543 103     22.9     25.1
 14        25.8 0.444 103     24.9     26.7

Confidence level used: 0.95 

> contrast(rg0.5, cont, adjust = "mvt")
 contrast estimate    SE  df t.ratio p.value
 10 - 8      0.902 0.330 103 2.733   0.0198 
 12 - 8      2.000 0.259 103 7.719   <.0001 
 12 - 10     1.098 0.318 103 3.453   0.0023 

P value adjustment: mvt method for 3 tests 

The "mvt" adjustment is the same multiplicity adjustment used by default in summary.glht.

I did not need to call emmeans() here because the reference grid involves only the one factor, age. For more complicated models, you'd probably have to add something like emm = emmeans(rg, "treatment")