R中nlme重复测量分析中的固定效应分离分析

Question

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

> 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-value来表示所有的季节。

我想知道这是不是我误解了lme的工作原理(我对R非常陌生)，或者我需要在我的公式/应用到我的数据集中包含一些东西，以便在所有季节的水平上完成分析。

我已经阅读了许多关于解释lme输出的公告板，但我不知道如何准确地解释我目前获得的输出，因为季节是分开的。

我也在尝试寻找一个合适的后自组织测试。

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

Answer 1

要获得季节效果的全面测试，只需使用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

Answer 2

您将获得每个季节的影响估计，因为您指定的模型将季节视为一个具有四个级别的因素。要评估季节是否相关，您可以比较具有和不具有该参数的模型之间的AIC，或者您可以对季节的影响强加一个参数形式，并对所需的单个模型项进行标准假设检验。但我认为这对于Cross Validated来说是一个真正的问题，而不是Stack Overflow，因为它是关于方法而不是代码。