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

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
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)? • Re Edit 2: VarCorr is definitely the wrong thing to use. That would be estimates of the mixed-effects part of the model. You need the covariance matrix of the fixed effects. As I suggested, that's probably in theta_z. It may take some doing, as it appeared from the documentation that it saves only the lower (or upper?) triangle of that matrix (e.g., 1+2+3+4=10 elements of a 4x4 matrix) Apr 8, 2020 at 15:40 • I was wrong. the theta_z slot is also for the random effects, not the fixed effects. I'll add to my answer. Apr 8, 2020 at 18:52 ## 1 Answer 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().]

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")