混合効果モデルで共分散行列を指定するためのpdBlockedの構文nlme

Question

私は混合効果モデルを持っています。私は自分のdofを減らすために、ランダム効果共分散行列に相関関係の一部を落としたいと思います。これを行うには、私は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)との相関関係があります。このような何か：このscenrioにも

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の用語をまとめて1つの言葉で表しているようです。

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では、構築する方法はないと思う

- 用語の順序を再配置する方法はありません（行/行列の列）を使用して、この行列がブロック対角になるようにします。あなたがこれをやらせるだろうモジュラー設計を有しており、私はlme4でこれを行う方法を見つけ出すために始めた

...あなたはあなた自身のpdMatrixクラスを原則で書くことができますが、それは超簡単ではありません少しが簡単に見つかりましたが、モデルに問題がありました。このデータセットでは過度に判断されます（実際のデータセットの場合はわかりません）。 Orthodontデータセットには1件につき4つの観測値しかないので、個体あたり4つのランダム効果値（インターセプト+ 3多項式値）を当てはめると、ランダム効果の分散が残差分散と混同されるモデルが得られますこれらのモデルから）。試してみると、lme4はエラーとなります。

しかし、これをやりたければエラーを上書きすることができます（危険なロビンソン！）まず、下三角のコレスキー係数を掛ける線形代数をいくつか実行しなければなりません[lme4は、 がtheta[2]*theta[4]+theta[5]*theta[7]と等価であることを確信してください。ここで、thetaは、コレスキー係数の要素のベクトルです（下三角行列、列一次）。したがって、完全な10パラメータモデルではなく、theta[7]が-theta[2]*theta[4]/theta[5]に設定された9パラメータモデルをフィッティングすることで、これを行うことができます。

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 回帰モデリング戦略）この特定の共分散構造を期待する機械的理由や主題に基づく理由がない限り、私は気にしません。