マルチレベルモデリングのための共分散構造

Question

私はトロポニン上昇を予測する反復測定を最大10回繰り返す約300人の患者のマルチレベル反復測定データセットを持っています。データセットには他の変数がありますが、ここにはそれらを含めていません。 nlmeを使用して、ランダムな傾き、患者間で効果が異なるランダム切片モデルを作成しようとしています。異なる患者では時間の効果が異なります。時間のために測定値の相関関係を可能にする一次共分散構造を導入しようとすると、次のエラーメッセージが表示されます。

Error in `coef<-.corARMA`(`*tmp*`, value = value[parMap[, i]]) : Coefficient matrix not invertible 

私のコードとデータセットのサンプルが含まれており、私は知恵の言葉に非常に感謝しています。

#baseline model includes only the intercept. Random slopes - intercept varies across patients 
randomintercept <- lme(troponin ~ 1, 
         data = df, random = ~1|record_id, method = "ML", 
         na.action = na.exclude, 
         control = list(opt="optim")) 

#random intercept and time as fixed effect 
timeri <- update(randomintercept,.~. + day) 
#random slopes and intercept: effect of time is different in different people 
timers <- update(timeri, random = ~ day|record_id) 
#model covariance structure. corAR1() first order autoregressive covariance structure, timepoints equally spaced 
armodel <- update(timers, correlation = corAR1(0, form = ~day|record_id)) 
Error in `coef<-.corARMA`(`*tmp*`, value = value[parMap[, i]]) : Coefficient matrix not invertible 

データ：

record_id day troponin 
1 1 32 
2 0  NA 
2 1 NA 
2 2 NA 
2 3 8 
2 4 6 
2 5 7 
2 6 7 
2 7 7 
2 8 NA 
2 9 9 
3 0 14 
3 1 1167 
3 2 1935 
4 0 19 
4 1 16 
4 2 29 
5 0 NA 
5 1 17 
5 2 47 
5 3 684 
6 0 46 
6 1 45440 
6 2 47085 
7 0 48 
7 1 87 
7 2 44 
7 3 20 
7 4 15 
7 5 11 
7 6 10 
7 7 11 
7 8 197 
8 0 28 
8 1 31 
9 0 NA 
9 1 204 
10 0 NA 
10 1 19 

Answer 1

あなたが「nlminb」あなたのオプティマイザを変更する場合はこれに合うことができる（あるいは少なくともそれが低減されたデータは、あなたが投稿設定で動作します）。あなたがフィットモデルを見れば

armodel <- update(timers, 
       correlation = corAR1(0, form = ~day|record_id), 
       control=list(opt="nlminb")) 

しかし、あなたが問題を持って表示されます - 推定AR1パラメータである-1ランダム切片と傾き用語は、R = 0.998と相関しています。

問題はデータの性質によると思います。データの大部分は10-50の範囲内にあるように見えるが、1桁または2桁の変動（例えば、個人6、最大約45000）がある。このスパイクのデータにモデルを適合させるのは難しいかもしれません。私は強くあなたのデータを変換することを提案します。標準的な診断プロット（plot(randomintercept)）は次のようになります。フィッティングのに対し

を対数スケールで

rlog <- update(randomintercept,log10(troponin) ~ .) 
plot(rlog) 

はややより合理的で、まだいくつかの証拠があるが異分散性

AR +ランダム斜面モデルはOKフィット：intervals(ar.rlog)で

ar.rlog <- update(rlog, 
       random = ~day|record_id, 
       correlation = corAR1(0, form = ~day|record_id)) 
## Linear mixed-effects model fit by maximum likelihood 
## ... 
## Random effects: 
## Formula: ~day | record_id 
## Structure: General positive-definite, Log-Cholesky parametrization 
##    StdDev Corr 
## (Intercept) 0.1772409 (Intr) 
## day   0.6045765 0.992 
## Residual 0.4771523  
## 
## Correlation Structure: ARMA(1,0) 
## Formula: ~day | record_id 
## Parameter estimate(s): 
##  Phi1 
## 0.09181557 
## ... 

一目は自己回帰パラメータの信頼区間は（-0.52,0.65）であることを示しているので、それは維持する価値はないかもしれません...不均一は、もはや問題のようだモデルにランダムな斜面で

...

plot(rlog,sqrt(abs(resid(.)))~fitted(.),type=c("p","smooth"))