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(age^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 

也为这个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 

我只是不知道如何做一个灵活的协方差矩阵能够挑选哪些是零。这些链接是非常有益的，但仍不能搞清楚究竟 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 - 有没有办法重新安排条款的顺序（行/矩阵的列）以便该矩阵是块对角线。你可以原则编写自己的pdMatrix类，但是这不会是超级容易...

我开始弄清楚如何在lme4，具有模块化设计做到这一点，它可以让你做到这一点稍微容易，但发现您的模型的另一个问题;它对这个数据集是否超定（我不知道它是否适用于你的真实数据集）。由于Orthodont数据集只有4个观测值，每个个体拟合4个随机效应值（截距加上3个多项式值）给出了一个模型，其中随机效应方差与剩余方差混淆（无法去除来自这些模型）。如果你尝试，lme4会给你一个错误。

然而，如果你仍然真的想这样做，你可以重写错误（危险将ROBINSON！）你首先必须做一些线性代数，乘以下三角Cholesky因子[这是lme4参数化的方式方差 - 协方差矩阵]来说服你自己Cov(age,age^3)相当于theta[2]*theta[4]+theta[5]*theta[7]，其中theta是Cholesky因子（下三角，列一阶）元素的向量。因此，我们可以通过拟合9参数模型而不是完整的10参数模型来完成此工作，其中theta[7]设置为-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) 

但是，我建议你仔细考虑一下这样做是否合理。我认为在这种方式下，通过下降条件，我认为你实际上不会获得非常多的精度/功率（一般来说，假设检验能力的明显提高源于事件模型简化是虚幻的 - 参见Harrell 回归建模策略），除非你有机械或基于主题的理由期望这种特殊的协变结构，我不认为我会打扰。