首页 文章

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

提问于
浏览 次
3

我有一个混合效果模型，我想在我的随机效应协方差矩阵中删除一些相关性以减少我的自由度 . 要做到这一点，我想我应该使用 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

任何帮助赞赏 . 谢谢

r covariance mixed-models nlme

1 回答

  • 3

    将 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

    我没有't think there'以任何方式构建你的第二个例子（ age 和 age^3 不相关，所有其他相关性非零）与 pdBlocked - 没有办法重新排列术语的顺序（矩阵的行/列），以便这个矩阵是块对角线 . 原则上你可以编写自己的 pdMatrix 类，但这不会太容易......

    我开始在 lme4 中找到如何做到这一点，它采用模块化设计，可以让你轻松地做到这一点，但发现了你的模型的另一个问题;它's overdetermined for this data set (I don'知道它是否适用于您的真实数据集） . 因为 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)

    但是，我建议你仔细考虑这样做是否有意义 . 我不认为你会以这种方式丢弃术语的精确度/功率方面实际获得很大的收益（一般来说，从事后模型简化得出的假设检验能力的明显增益是虚幻的 - 参见Harrell回归建模策略），除非你有机械或基于主题的理由来期待这种特殊的协方差结构，我认为我不会打扰 .

    回复于

相关问题