我有一个混合效应模型，我想在我的随机效应协方差矩阵中删除一些相关性以减少我的自由度。为此，我认为我应该使用 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
任何帮助表示赞赏。谢谢
最佳答案
将 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 会给你一个错误。但是，如果您仍然真的想这样做，您可以覆盖错误(危险将罗宾逊!)您首先必须做一些线性代数，乘以下三角 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)
但是，我建议您仔细考虑这样做是否有意义。我认为您实际上不会通过以这种方式删除术语而在精度/功率方面获得很大 yield (一般来说，从事后模型简化中获得的假设检验能力的明显增益是虚幻的 - 请参阅 Harrell Regression Modeling Strategies)，除非您有机械或基于主题的理由来期待这种特殊的协方差结构，我认为我不会打扰。
关于r - 在混合效应模型 nlme 中指定协方差矩阵的 pdBlocked 语法，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/38976189/