固定效應分離使用nlme重複測量分析R

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值的所有季節。

我想知道這是不是我誤解了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() 

您也有通過觀察QQ圖得到這樣的結論：

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

Answer 2

由於您指定的模型將季節作爲四個級別的因子對待，因此您會得到每個季節的效果估計值。爲了評估季節是否相關，你可以比較有和沒有這個參數的模型的AIC，或者你可以對季節的影響強加一個參數形式，並對單個模型術語進行標準假設檢驗。但我認爲這是真正的交叉驗證問題而不是堆棧溢出問題，因爲它涉及的是方法而不是代碼。