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 迴歸建模策略），除非你有機械或基於主題的理由期望這種特殊的協變結構，我不認爲我會打擾。