使用 R 中的 nlme 重复测量分析分离固定效应以进行分析答案

【问题标题】：Fixed Effect Being Separated for Analysis Using nlme Repeated Measures Analysis in R使用 R 中的 nlme 重复测量分析分离固定效应以进行分析
【发布时间】：2015-09-30 23:26:45
【问题描述】：

我有一组水质数据，我试图从中确定是否存在季节性趋势。我按季节对数据进行了分类，我的数据如下所示：

> Cond
  Site      Cond Season Watershed  logCond
1  BICO201  41.86667 Spring      BICO 1.621868
2  BICO301  53.16000 Spring      BICO 1.725585
3  MIDO301  42.63333 Spring      MIDO 1.629749
4  MIDO601  52.10000 Spring      MIDO 1.716838
5  MIDO704  82.70000 Spring      MIDO 1.917506
6  MIDO801  74.36667 Spring      MIDO 1.871378
7  MIDO802  73.43333 Spring      MIDO 1.865893
8  MIDO803  85.72000 Spring      MIDO 1.933082
9  NORO401  43.30000 Spring      NORO 1.636488
10 NORO502 132.05000 Spring      NORO 2.120738
11 NORO503  61.36667 Spring      NORO 1.787933
12 NORO517 142.40000 Spring      NORO 2.153510
13 NORO520  95.20000 Spring      NORO 1.978637
14 NORO527  81.08000 Spring      NORO 1.908914
15 NORO601 479.75000 Spring      NORO 2.681015
16 BICO201  47.73333 Summer      BICO 1.678822
17 BICO301  58.46667 Summer      BICO 1.766908
18 MIDO301  45.75000 Summer      MIDO 1.660391
19 MIDO601  51.80000 Summer      MIDO 1.714330
20 MIDO704 112.30000 Summer      MIDO 2.050380
21 MIDO801  90.10000 Summer      MIDO 1.954725
22 MIDO802  74.58000 Summer      MIDO 1.872622
23 MIDO803 112.70000 Summer      MIDO 2.051924
24 NORO401  71.40000 Summer      NORO 1.853698
25 NORO502 192.88000 Summer      NORO 2.285287
26 NORO503  80.42500 Summer      NORO 1.905391
27 NORO517 156.50000 Summer      NORO 2.194514
28 NORO520 114.22500 Summer      NORO 2.057761
29 NORO527 109.00000 Summer      NORO 2.037426
30 NORO601 420.00000 Summer      NORO 2.623249
31 BICO201  46.85000   Fall      BICO 1.670710
32 BICO301  55.43333   Fall      BICO 1.743771
33 MIDO301  42.52500   Fall      MIDO 1.628644
34 MIDO601  69.26667   Fall      MIDO 1.840524
35 MIDO704 102.40000   Fall      MIDO 2.010300
36 MIDO801  81.67500   Fall      MIDO 1.912089
37 MIDO802  62.05000   Fall      MIDO 1.792742
38 MIDO803  86.90000   Fall      MIDO 1.939020
39 NORO401  62.85000   Fall      NORO 1.798305
40 NORO502 149.60000   Fall      NORO 2.174932
41 NORO503  57.90000   Fall      NORO 1.762679
42 NORO517  92.90000   Fall      NORO 1.968016
43 NORO520 118.31667   Fall      NORO 2.073046
44 NORO527 123.15000   Fall      NORO 2.090434
45 NORO601 522.33333   Fall      NORO 2.717948
46 BICO201 101.96000 Winter      BICO 2.008430
47 BICO301  69.47500 Winter      BICO 1.841829
48 MIDO301  43.58333 Winter      MIDO 1.639320
49 MIDO601  49.78000 Winter      MIDO 1.697055
50 MIDO704  94.73333 Winter      MIDO 1.976503
51 MIDO801  76.28000 Winter      MIDO 1.882411
52 MIDO802  65.86667 Winter      MIDO 1.818666
53 MIDO803 119.13333 Winter      MIDO 2.076033
54 NORO401  54.20000 Winter      NORO 1.733999
55 NORO502 171.76000 Winter      NORO 2.234922
56 NORO503  83.76667 Winter      NORO 1.923071
57 NORO517 191.07500 Winter      NORO 2.281204
58 NORO520 118.31667 Winter      NORO 2.073046
59 NORO527 123.15000 Winter      NORO 2.090434
60 NORO601 576.00000 Winter      NORO 2.760422

我正在尝试使用混合效应来运行重复测量分析，其中季节作为我的固定效应，站点作为随机效应。我正在使用 nlme 包，我的代码如下所示：

> mod.1.2<-lme(Cond~Season, random=~1|Site,data=Cond)

然后我运行我的模型的摘要并获得以下输出：

> summary(mod.1.2)
 Linear mixed-effects model fit by REML
  Data: Cond 
        AIC      BIC    logLik
   595.4271 607.5792 -291.7136

Random effects:
 Formula: ~1 | Site
        (Intercept) Residual
StdDev:    111.1618 22.68229

Fixed effects: Cond ~ Season 
                 Value Std.Error DF   t-value p-value
(Intercept)  111.61000 29.293255 42  3.810092  0.0004
SeasonSpring  -8.86822  8.282401 42 -1.070731  0.2904
SeasonSummer   4.24733  8.282401 42  0.512814  0.6108
SeasonWinter  17.66200  8.282401 42  2.132474  0.0389
 Correlation: 
             (Intr) SsnSpr SsnSmm
SeasonSpring -0.141              
SeasonSummer -0.141  0.500       
SeasonWinter -0.141  0.500  0.500

Standardized Within-Group Residuals:
       Min         Q1        Med         Q3        Max 
-3.3746755 -0.3431503 -0.0313137  0.3702357  2.9115215 

Number of Observations: 60
Number of Groups: 15

我很困惑，因为 R 将我的固定因子分解为不同的季节，但我期望我的输出只是给我一个值/StdDev/DF/p 值用于所有季节。

我想知道这是否是我误解了 lme 的工作原理（我对 R 非常陌生），或者我是否需要在我的公式中包含一些内容/应用于我的数据集以完成分析所有季节。

我已经阅读了许多关于解释 lme 输出的板，但我无法弄清楚我将如何准确地解释我目前得到的输出，因为季节是分开的。

我也在尝试寻找合适的事后测试。

任何帮助将不胜感激，谢谢！

【问题讨论】：

标签： r mixed-models

【解决方案1】：

要全面测试Season的效果，只需使用anova()：

anova(mod.1.2)
##             numDF denDF   F-value p-value
## (Intercept)     1    42 15.852534  0.0003
## Season          3    42  3.558053  0.0221

顺便说一下，查看您的数据，我建议您进行日志转换你的响应变量：

library(ggplot2); theme_set(theme_bw())
ggplot(Cond,aes(Season,Cond,group=Site))+
    geom_line(colour="gray")+
        geom_point()+
            scale_y_log10()

您也可以通过查看 Q-Q 图得出这个结论：

qqnorm(mod.1.2)
mod.1.3 <- update(mod.1.2,log(Cond)~.)
qqnorm(mod.1.3)  ## better

【讨论】：

非常感谢！实际上，我确实对其进行了日志转换以进行实际分析，为简单起见，我只是使用了初始公式！

【解决方案2】：

因为您指定的模型将季节视为具有四个水平的因素，所以您得到了每个季节的效果估计值。要评估季节是否相关，您可以比较具有和不具有该参数的模型之间的 AIC，或者您可以对季节的影响施加参数形式，并对包含的单个模型项进行标准假设检验。但我认为这确实是 Cross Validated 而不是 Stack Overflow 的问题，因为它是关于方法而不是代码。

【讨论】：

像anova() 所做的那样，仅对参数的综合影响进行 F 检验有什么问题... ??
我认为你不能用广义线性混合模型的参数来做到这一点，因为分布假设不成立。我不确定这是否正确，但同样，这确实是一个方法问题，而不是编码问题。
对于具有平衡设计且没有交叉随机效应的线性混合模型，它确实成立。