lme4: how to specify 2 correlations with random intercept, without adding correlations between random slopes

Question

(re-posted from stats.stackexchange.com)

I am trying to specify a model in R's lme4 package in which I have 2 correlations between random intercept and random slopes, but the random slopes are not allowed to correlate.

lmer (Y ~ A + B + (1+A+B|Subject), data=mydata)

is bad because it models correlation between the random slopes for A and B.

Whereas:

lmer (Y ~ A + B + (1+A|Subject) + (1+B|Subject), data=mydata)

is bad because the random intercept for Subject gets introduced into the model twice. Is there a third way, perhaps more hack-ish?

Answer 1

This turned out to be harder than I thought!

The variance-covariance matrices within lme4 are parameterized according to their Cholesky factors (essentially a matrix square root); therefore, if we want to set up a model with a particular correlation fixed to zero, we want

t1  0  0     t1 t2 t3      t1^2   t1*t2          t1*t3
t2 t4  0     0  t4 t5   =  t1*t2  t2^2 + t4^2    t2*t3 + t4*t5
t3 t5 t6     0   0 t6   =  t1*t3  t2*t3 + t4*t5  t3^2 + t5^2 + t6^2

and solve for the [3,2] element (the correlation between A and B) to be equal to zero; in other words, we will need t2 t3 + t4 t5 == 0, or a 6-element vector where t5 == -t2*t3/t4;

tfun <- function(theta) {
  theta5 <- -theta[2]*theta[3]/theta[4] 
  c(theta[1:4],theta5,theta[5])
}

Simulate some data:

set.seed(101)
dd <- data.frame(A=rnorm(1000),B=rnorm(1000),
                 Subject=factor(rep(1:20,50)))

library("lme4")
dd$Y <- simulate(~A+B+(1+A+B|Subject),
         newdata=dd,
         family=gaussian(),
         newparams=list(beta=c(1,2,3),
                        theta=tfun(c(1,0.2,0.3,2,3)),
                        sigma=1))[[1]]

Now follow the steps in ?modular:

lmod <- lFormula(Y ~ A + B + (1+A+B|Subject), data=dd)
devfun <- do.call(mkLmerDevfun, lmod)

A wrapper function for devfun() that will take a 5-element vector, compute the corresponding constrained theta vector, and pass it to devfun():

devfun2 <- function(theta) {
    devfun(tfun(theta))
}

Delete one term from the lower-bound vector:

lwr <- lmod$reTrms$lower
## [1]    0 -Inf -Inf    0 -Inf    0
lwr <- lwr[c(1:4,6)]
library("minqa")
## n.b. optwrap fails with minqa::bobyqa
opt <- lme4:::optwrap(optimizer=bobyqa,
                      par=ifelse(lwr==0,1,0),
                      fn=devfun2,
                      lower=lwr)

Now adjust the result according to the parameter transformation:

opt$par <- tfun(opt$par)
m1 <- mkMerMod(environment(devfun), opt, lmod$reTrms, fr = lmod$fr)

VarCorr(m1)    
##  Groups   Name        Std.Dev. Corr       
##  Subject  (Intercept) 1.41450             
##           A           1.49374  0.019      
##           B           2.47895  0.316 0.000
##  Residual             0.96617             

The desired correlation is now fixed to zero.