Фиксированный эффект, выделенный для анализа с использованием анализа повторных измерений 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-value для всех сезонов.

Мне интересно, если это я не понимаю, как работает lme (я ОЧЕНЬ новичок в R), или если есть что-то, что мне нужно включить в мою формулу/применить к моему набору данных, чтобы анализ был завершен на уровне всех сезонов.

Я прочитал несколько советов об интерпретации вывода lme, но я не могу понять, как именно я буду интерпретировать вывод, который я сейчас получаю, потому что сезоны разделены.

Я также пытаюсь найти подходящий пост-hoc-тест.

Любая помощь будет очень признательна, спасибо!

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 по моделям с этим параметром и без него, или вы можете наложить параметрическую форму на эффект сезона и выполнить стандартный тест гипотезы по единому модельному термину, который влечет за собой. Но я думаю, что это действительно вопрос для Cross Validated, а не для переполнения стека, потому что речь идет о методе, а не о коде.