lme with varIdent 用于重复测量

Question

在使用 lme 和重复测量和 varIdent 时，我遇到了奇怪的结果。对此的任何帮助将不胜感激！

我正在测试两个物种（A 和 B）在时间序列上的 13C 信号是否不同。我基本上对物种之间的整体差异感兴趣，而不是特定时间点。

这是我的数据集：

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，而在其他一些时间点则相反。

> anova(bulk.var)
                         numDF denDF  F-value p-value
(Intercept)                  1    15 13.25772  0.0024
SpeciesCode                  1     2 67.08281  0.0146
SamplingTime                 4    15  4.42320  0.0147
SpeciesCode:SamplingTime     4    15  1.27659  0.3227

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

我想知道问题是否可能是每个物种在每个采样时间 (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。请记住，delta 13 C 是（转换后的）分数 13C/12C。一个正态假设似乎有点可疑（尽管我自己偶尔将它用于增量值）。但是，在我看来，我们可以根据拟合值对方差进行建模。

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 & Bates 2000 并尝试找出自己或在Cross Validated 上提问。