syntax of pdBlocked to specify covariance matrix in mixed-effects model nlme

Question

I have a mixed-effects model and I want to drop some of my correlations in my random-effects covariance matrix to reduce my dof. To do this I think I should use pdBlocked but can't get the correct syntax to get specifically what I want.

Example code:

library(nlme)
m3 <- lme(distance ~ age +I(age^2) + I(age^3), data = Orthodont, 
           random = list(Subject = pdBlocked(list(~ age,~0 + I(age^2),~0+I(age^3)))))

which gives the following covariance matrix:

getVarCov(m3)
Random effects variance covariance matrix
            (Intercept)      age         I(age^2)                     I(age^3)
(Intercept)      5.2217 -0.30418 0.00000000000000 0.00000000000000000000000000
age             -0.3042  0.04974 0.00000000000000 0.00000000000000000000000000
I(age^2)         0.0000  0.00000 0.00000000003593 0.00000000000000000000000000
I(age^3)         0.0000  0.00000 0.00000000000000 0.00000000000000000000002277
  Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

This is close to what I want but not quite. I would like to keep the correlation between I(age^3) and intercept, age zero but allow a correlation with I(age^2). Something like this:

getVarCov(m3)
Random effects variance covariance matrix
            (Intercept)      age         I(age^2)                     I(age^3)
(Intercept)      5.2217 -0.30418 0.00000000000000 0.00000000000000000000000000
age             -0.3042  0.04974 0.00000000000000 0.00000000000000000000000000
I(age^2)         0.0000  0.00000 0.00000000003593 a_value
I(age^3)         0.0000  0.00000 a_value          0.00000000000000000000002277
  Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

also for this scenrio

getVarCov(m3)
    Random effects variance covariance matrix
                (Intercept)      age         I(age^2)                     I(age^3)
    (Intercept)      5.2217 -0.30418 c_value          b_value          
    age             -0.3042  0.04974 d_value          0.00000000000000000000000000
    I(age^2)         c_value d_value 0.00000000003593 a_value
    I(age^3)         b_value 0.00000 a_value          0.00000000000000000000002277
      Standard Deviations: 2.285 0.223 0.000005994 0.000000000004772 

I'm just not sure how to make a flexible covariance matrix to able to pick which ones are zero. These links were very helpful but still cant figure it out exactly http://rpsychologist.com/r-guide-longitudinal-lme-lmer

https://stats.stackexchange.com/questions/87050/dropping-term-for-correlation-between-random-effects-in-lme-and-interpretting-su?rq=1

Any help appreciated. Thanks

Answer 1

Putting the age^2 andage^3 terms together in a single term seems to do it.

m4 <- lme(distance ~ age +I(age^2) + I(age^3), data = Orthodont, 
          random = list(Subject = pdBlocked(list(~ age,
                                                 ~0 + I(age^2)+I(age^3)))),
          control=lmeControl(opt="optim"))
getVarCov(m4)
## Random effects variance covariance matrix
##             (Intercept)       age    I(age^2)    I(age^3)
## (Intercept)     5.00960 -0.225450  0.0000e+00  0.0000e+00
## age            -0.22545  0.019481  0.0000e+00  0.0000e+00
## I(age^2)        0.00000  0.000000  4.1676e-04 -1.5164e-05
## I(age^3)        0.00000  0.000000 -1.5164e-05  5.5376e-07
##   Standard Deviations: 2.2382 0.13957 0.020415 0.00074415 

I don't think there's any way to construct your second example (age and age^3 uncorrelated, all other correlations non-zero) with pdBlocked - there's no way to re-arrange the order of terms (rows/columns of the matrix) so that this matrix is block diagonal. You could in principle write your own pdMatrix class, but that would not be super-easy ...

I started to figure out how to do this in lme4, which has a modular design that would let you do this slightly more easily, but discovered an additional problem with your model; it's overdetermined for this data set (I don't know if it is for your real data set). Because the Orthodont data set only has 4 observations per subject, fitting 4 random-effect values per individual (intercept plus 3 polynomial values) gives you a model where the random-effects variances are confounded with the residual variance (which can't be removed from these models). If you try, lme4 gives you an error.

However, if you still really want to do this, you can override the error (DANGER WILL ROBINSON!) You first have to do some linear algebra, multiplying the lower-triangular Cholesky factor [which is the way that lme4 parameterizes the variance-covariance matrix] to convince yourself that the Cov(age,age^3) is equivalent to theta[2]*theta[4]+theta[5]*theta[7], where theta is the vector of elements of the Cholesky factor (lower-triangular, column-first order). Thus we can do this by fitting a 9-parameter model instead of the full 10-parameter model, with theta[7] set equal to -theta[2]*theta[4]/theta[5] ...

lf <- lFormula(distance ~ age +I(age^2) + I(age^3) +
                 (age+ I(age^2) + I(age^3)|Subject), data = Orthodont,
               control=lmerControl(check.nobs.vs.nRE="ignore"))
devfun <- do.call(mkLmerDevfun,lf)
trans_theta <- function(theta)
  c(theta[1:6],-theta[2]*theta[4]/theta[5],theta[7:9])
devfun2 <- function(theta) {
  return(devfun(trans_theta(theta)))
}
diagval <- (lf$reTrms$lower==0)
opt <- minqa::bobyqa(fn=devfun2,par=ifelse(diagval,1,0)[-7],
             lower=lf$reTrms$lower[-7])
opt$par <- trans_theta(opt$par)
opt$conv <- 0
m1 <- mkMerMod(environment(devfun), opt, lf$reTrms, fr = lf$fr)
VarCorr(m1)

However, I would suggest you think carefully about whether it makes sense to do this. I don't think you will actually gain very much in terms of precision/power from dropping terms in this way (in general, apparent gains in hypothesis-testing power derived from post hoc model simplification are illusory - see Harrell Regression Modeling Strategies) unless you have mechanistic or subject-based reasons to expect this particular covariance structure, I don't think I would bother.