在混合效应模型 nlme 中指定协方差矩阵的 pdBlocked 语法

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

也适用于此场景

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 - 无法重新排列项的顺序（矩阵的 rows/columns）以便该矩阵是块对角线。您可以 原则上 编写您自己的 pdMatrix class，但这并不是超级容易...

我开始在 lme4 中弄清楚如何做到这一点，它有一个模块化设计，可以让你稍微更容易地做到这一点，但发现了一个额外的你的模型有问题；对于这个数据集是多定的（不知道是不是你的真实数据集）。因为 Orthodont 数据集每个受试者只有 4 个观察值，所以每个个体拟合 4 个随机效应值（截距加上 3 个多项式值）给你一个模型，其中随机效应方差与残差方差（可以' 从这些模型中删除）。如果您尝试，lme4 会给您一个错误。

但是，如果您仍然真的想这样做，您可以覆盖错误（DANGER WILL 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)

但是，我建议您仔细考虑这样做是否有意义。我认为您实际上不会从 precision/power 中以这种方式删除术语而获得太多收益（通常，从 post hoc[ 得出的假设检验能力明显增加=41=] 模型简化是虚幻的 - 请参阅 Harrell 回归建模策略 ) 除非你有机械或基于主题的理由期望这种特定的协方差结构，否则我认为我不会打扰。

在混合效应模型 nlme 中指定协方差矩阵的 pdBlocked 语法

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

r

covariance

mixed-models

nlme