lme with varIdent for repeated measures

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值，但看着我的数据似乎并没有该物种的不同在所有...我认为这是很奇怪的得到如此低的p值，因为误差线在每个时间点重叠，并且在某些时间点A大于B，而在另一些时间点则相反。

​​

它再次发生了，当我分析了其他类似的变量...

我不知道一个问题可能是每个品种在每个取样时间的低复制（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是（转化的）13C/12C分数。一个正常假设似乎有点可疑（尽管我自己偶尔会将它用于delta值）。但是，在我看来，我们可以根据拟合值对方差进行建模。

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 

所以，你得到了同样的问题与使用四个参数少的模型。现在我不知道为什么物种效应在顺序方差分析中显着，如果包含方差结构，尽管相应的模型参数不显着。你可以学习Pinheiro &贝茨2000，并尝试找出自己或要求Cross Validated。