反復措置のためにvarIdentを持つlme

Question

反復尺度とvarIdentでlmeを使用すると、私は奇妙な結果に遭遇しました。これに関する助けは非常に高く評価されます！

時系列に沿った葉の13Cシグナルが2種（AとB）で異なるかどうかをテストしています。私は基本的に特定の時点ではなく、種間の全体的な違いに興味があります。

Block Species time X13C 
1 B 2 0.775040865 
2 B 2 0.343913792 
3 B 2 0.381053614 
1 A 2 0.427101597 
2 A 2 0.097743662 
3 A 2 0.748345826 
1 B 24 0.416700446 
2 B 24 0.230773558 
3 B 24 0.681386484 
1 A 24 0.334026511 
2 A 24 0.866426406 
3 A 24 0.606346215 
1 B 48 0.263085491 
2 B 48 0.083323709 
3 B 48 0.534697801 
1 A 48 0.30594443 
2 A 48 0.024555489 
3 A 48 0.790670392 
1 B 96 0.158090804 
2 B 96 0.254880689 
3 B 96 0.082666799 
1 A 96 0.139189281 
2 A 96 0.300340119 
3 A 96 0.233149535 
1 B 192 0.055421148 
2 B 192 0.082582155 
3 B 192 0.136636735 
1 A 192 0.03641637 
2 A 192 0.06082544 
3 A 192 0.126029308 

私は次のモデル適用しています：ここで

は私のデータセットです時間の残差の不均一があるので

bulk<-lme(X13Catex ~ Species*time, random = ~1|Block/Species, method='REML', na.action=na.exclude, data=VacL, corAR1())

を、私はモデルを改良しvarIdentは、適用されますフィット（AIC）。正規化された残差プロットも良く見えました。

bulk.var<-lme(X13Catex ~ Species*time, random = ~1|Block/Species, method='REML', na.action=na.exclude, data=VacL, corAR1(), weights=varIdent(form=~1|time))

事は、このコードで、私は種のための重要なp値を得ることであるが、それは種が全く異なることは思えない私のデータを見て...私はそれが非常に奇妙だと思いますエラーバーが各時点で重なり、いくつかの時点AでBがBより大きく、他のいくつかの点では逆になるので、このような低いp値を得る。問題は、各サンプリング時点でのそれぞれの種のための低複製することができれば、私は他の類似の変数を分析した場合、それが再び起こった

​​

...

が、私は疑問に思う（N = 3）。このように複製量が少ないvarIdentと「比較的複雑な」モデルを適用すると、p値の有意な値がわかりますか？これに対処する方法に関する提案はありますか？

ありがとうございました！

Answer 1

よろしくお願いします。

まず、あなたの相関構造が正しいとは思われません。あなたはそこで共存する時間を使う必要があります。

fit0 <- lme(X13C ~ Species*time, random = ~1|Block/Species, method='REML', 
      na.action=na.exclude, data=VacL, 
      corAR1(0.9, form = ~ time | Block/Species)) 
summary(fit0) 

次に、ネストされたランダム効果の分散はやや小さいようです。このパラメータを削除しようとしましょう。

fit1 <- lme(X13C ~ Species*time, random = ~1|Block, method='REML', 
      na.action=na.exclude, data=VacL, 
      corAR1(0.9, form = ~ time | Block/Species)) 
summary(fit1) 

anova(fit0, fit1) 
#  Model df  AIC  BIC logLik Test  L.Ratio p-value 
#fit0  1 8 35.47003 45.5348 -9.735014       
#fit1  2 7 33.47003 42.2767 -9.735014 1 vs 2 8.37192e-09 0.9999 

plot(fit1) 

プロットは確かに強い異質性を示しています。この時点で、私は真剣にGLMMを検討するだろう。忘れないでください、デルタ13 Cは（変換された）フラクション13 C/12 Cです。通常の仮定はこれに少し疑わしいようです（私はたまにデルタ値のためにそれを使用します）。しかし、当てはめられた値に依存して分散をモデル化することができるように私には見えます。

fit2 <- lme(X13C ~ Species*time, random = ~1|Block, method='REML', 
      na.action=na.exclude, data=VacL, 
      corAR1(0.9, form = ~ time | Block/Species), 
      weights = varPower()) 

plot(fit2, resid(., type = "normalized") ~ fitted(.)) 

anova(fit1, fit2) 
#  Model df  AIC  BIC logLik Test L.Ratio p-value 
#fit1  1 7 33.47003 42.27670 -9.735014       
#fit2  2 8 11.34319 21.40796 2.328405 1 vs 2 24.12684 <.0001 

ない、あまりにも悪いです。 p値を調べましょう。

coef(summary(fit2)) 
#      Value Std.Error DF t-value  p-value 
#(Intercept) 0.3906798322 0.0640391495 24 6.1006405 2.661703e-06 
#SpeciesB  -0.0303078937 0.0777180616 24 -0.3899723 6.999965e-01 
#time   -0.0016541839 0.0003059863 24 -5.4060727 1.491893e-05 
#SpeciesB:time 0.0002578782 0.0004048368 24 0.6369930 5.301592e-01 

インターセプトもスロープも大きく異なりません。今、ANOVAを試してみましょう。分散構造のない比較のために

anova(fit2) 
      numDF denDF F-value p-value 
(Intercept)  1 24 9.0061 0.0062 
Species   1 24 525.7457 <.0001 
time    1 24 56.5135 <.0001 
Species:time  1 24 0.4058 0.5302 

：

anova(fit1) 
      numDF denDF F-value p-value 
(Intercept)  1 24 29.536428 <.0001 
Species   1 24 0.319802 0.5770 
time    1 24 17.602173 0.0003 
Species:time  1 24 0.041482 0.8403 

だから、あなたはより少ない4つのパラメータが使用するモデルと同じ問題を取得します。そして今、私は、対応するモデルパラメータが重要ではないが、分散構造が含まれている場合、なぜ種の影響が連続ANOVAにおいて重要であるのか分からない。あなたはピンヘイロ&ベイツ2000を勉強して、自分自身を見つけようとしたり、Cross Validatedで尋ねることができます。