r - pdBlocked 的语法用于指定混合效应模型 nlme 中的协方差矩阵

Question

我有一个混合效应模型，我想删除随机效应协方差矩阵中的一些相关性以减少自由度。为此，我认为我应该使用 pdBlocked 但无法获得正确的语法来获得我想要的具体内容。

示例代码:

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

给出以下协方差矩阵:

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 

这接近我想要的，但不完全是。我想保持 I(age^3) 和 intercept、age 之间的相关性为零，但允许与 I(年龄^2)。像这样的事情:

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 

也适用于这种情况

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 

我只是不确定如何创建一个灵活的协方差矩阵来选择哪些为零。这些链接非常有帮助，但仍然无法准确弄清楚 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

任何帮助表示赞赏。谢谢

Answer 1

将 age^2 和 age^3 术语放在一个术语中似乎可以做到这一点。

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 

我认为没有任何方法可以构建您的第二个示例(age 和age^3 不相关，所有其他相关性非零)与 pdBlocked - 无法重新排列项的顺序(矩阵的行/列)，以便该矩阵是 block 对角矩阵。 原则上你可以编写自己的pdMatrix类，但这并不是 super 简单......

我开始弄清楚如何在 lme4 中执行此操作，它具有模块化设计，可以让您更轻松地执行此操作，但发现了另一个问题你的模型；对于这个数据集来说它是过度确定的(我不知道它是否适合你的真实数据集)。由于 Orthodont 数据集每个受试者只有 4 个观测值，因此对每个个体拟合 4 个随机效应值(截距加 3 个多项式值)会得到一个模型，其中随机效应方差与残差方差混杂在一起(无法从这些模型中删除)。如果您尝试，lme4 会给您一个错误。

但是，如果您仍然确实想这样做，您可以覆盖该错误(危险将罗宾逊!)您首先必须做一些线性代数，乘以下三角 Cholesky 因子 [这就是 lme4 参数化方差-协方差矩阵]，让自己相信 Cov(age,age^3) 等价于 theta[2]*theta[4]+theta[ 5]*theta[7]，其中 theta 是 Cholesky 因子(下三角、列优先)的元素向量。因此，我们可以通过拟合 9 参数模型而不是完整的 10 参数模型来实现此目的，并将 theta[7] 设置为等于 -theta[2]*theta[4]/θ[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)

但是，我建议您仔细考虑一下这样做是否有意义。我认为，通过这种方式删除项，您实际上不会在精度/功效方面获得太多 yield (一般来说，从事后模型简化中获得的假设检验能力的明显 yield 是虚幻的 -参见 Harrell 回归建模策略)除非您有机械或基于主题的原因来期望这种特定的协方差结构，否则我认为我不会打扰。