Tvíhliða fervikagreining (ANOVA) fyrir endurteknar mælingar í R

Bakgrunnur

Hér munun við sýna hvernig má framkvæma tvíhliða fervikagreiningu fyrir endurteknar mælingar. Fyrir gögn sem ekki eru með endurteknum mælingum þá má sjá R sýnikóða fyrir tvíhliða fervikagreiningu hér. Fervikagreining er notuð til að bera saman meðaltöl milli hópa.

Við getum notað tvíhliða (e. two-way) ANOVA fyrir endurteknar mælingar þegar við mælum alla einstaklinga við öll gildi á einni flokkabreytu (t.d. tímapunktum) og skiptum einstaklingunum í hópa miðað við gildi á annarri flokkabreytu (t.d. þá sem fá lyf og ekki lyf). Þar sem hver einstakligur hefur mælingar á öllum tímapunktum, þá eru tímapunktar “cross-over factor” og þar sem ekki allir einstaklingar eru mældir fyrir öll gildi á seinni flokkabreytunni (lyfjabreytunni) þá er hún “nest factor”.

Þegar framkvæmd er tvíhliða ANOVA fyrir endurteknar mælingar þá er mikilvægt að hafa jafnmargar mælingar á öllum samsetningum af báðum flokkabreytum, t.d. jafnmargir sem fá lyf og ekki lyf og allir mældir á öllum tímapunktum. Ef þetta gildir er sagt að gögnin séu “balanced”.

Fyrir neðan sýnum við nokkrar leiðir til að framkvæma tvíhliða fervikagreiningu fyrir enduteknar mælingar í R. Áður sýndum við hvernig mátti framkvæma einhliða fervikagreinigu fyrir endurteknar mælingar.

Gagnasett

Við lesum inn gagnasett héðan.

Við þurfum að skilgreina flokkabreytur sem “factor” þannig R líti ekki á þær sem talnabreytur.

demo3 <- read.csv("http://www.ats.ucla.edu/stat/data/demo3.csv")
demo3 <- within(demo3, {
  group <- factor(group)
  time <- factor(time)
  id <- factor(id)
})
summary(demo3)
##        id    group      pulse       time 
##  1      :3   1:12   Min.   :12.00   1:8  
##  2      :3   2:12   1st Qu.:22.00   2:8  
##  3      :3          Median :29.00   3:8  
##  4      :3          Mean   :32.79        
##  5      :3          3rd Qu.:46.00        
##  6      :3          Max.   :60.00        
##  (Other):6
str(demo3)
## 'data.frame':    24 obs. of  4 variables:
##  $ id   : Factor w/ 8 levels "1","2","3","4",..: 1 1 1 2 2 2 3 3 3 4 ...
##  $ group: Factor w/ 2 levels "1","2": 1 1 1 1 1 1 1 1 1 1 ...
##  $ pulse: int  35 25 16 32 23 12 36 22 14 34 ...
##  $ time : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...

1. Notum aov fallið

demo3.aov <- aov(pulse ~ group + time + group:time + Error(id), data = demo3)
demo3.aov<-summary(demo3.aov)
demo3.aov
## 
## Error: id
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## group      1 2035.0  2035.0   343.1 1.6e-06 ***
## Residuals  6   35.6     5.9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Error: Within
##            Df Sum Sq Mean Sq F value   Pr(>F)    
## time        2 2830.3  1415.2   553.8 1.52e-12 ***
## group:time  2  200.3   100.2    39.2 5.47e-06 ***
## Residuals  12   30.7     2.6                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hér að ofan sjáum við að p=1.5e-12 fyrir tengsl milli tímapunkta (time) og pulse, þ.e. pulse virðist mærktæt mismunandi milli tímapunktanna þriggja að teknu tilliti til þess að hver einstaklingur er með pulse mælingar á hverjum tímapuntki. Að auki eru áhrifin af group marktæk og einnig víxlhrifin milli group og time.

2. Notum Anova fallið í car pakkanum

library(car)
library(reshape2)

# reformat to wide format
demo3.wide <- dcast(demo3,id+group~time,value.var="pulse")
names(demo3.wide)[3:5] <- c("t1","t2","t3")

# prepare arguments for Anova function
# model
mod <- lm(cbind(t1,t2,t3)~1+group,demo3.wide)
# prepare 3-level factor for 3 time-points
time.level <- c(1,2,3)
time.factor <- as.factor(time.level)
time.frame <- data.frame(time.factor)
# run Anova
demo3.Anova <- Anova(mod, idata=time.frame,idesign=~time.factor, type=2)
demo3.Anova <- summary(demo3.Anova)
demo3.Anova$univariate.tests
##                        SS num Df Error SS den Df        F    Pr(>F)    
## (Intercept)       25807.0      1   35.583      6 4351.539 8.162e-10 ***
## group              2035.0      1   35.583      6  343.145 1.596e-06 ***
## time.factor        2830.3      2   30.667     12  553.761 1.517e-12 ***
## group:time.factor   200.3      2   30.667     12   39.196 5.474e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aftur fáum við sömu niðurstöðu. Takið eftir að hér settum við ekki víxlhrifin sérstaklega inn í modelið, Anova fallið sér um það.

3. Notum línuleg slembiþáttalíkön (linear mixed model)

Einnig er hægt að nota línuleg slembiþáttalíkön. Athugið að hér er mikilvægt að nota sviga utan um “1|id” í formúlunni. Slembiþáttalíkön eru sérstaklega gagnleg þegar gögnin er ekki “balanced”, þ.e. til dæmis þegar vantar einhverjar mælingar þannig ekki eru jafnmargir einstaklingar með mælingar á öllum tímapunktum.

library(lme4)
form <- as.formula("pulse~ group + time + group:time + (1|id)")
demo3.mixed <- lmer(form,data=demo3,REML = T)
x<-summary(demo3.mixed)
y<-anova(demo3.mixed)
y
## Analysis of Variance Table
##            Df  Sum Sq Mean Sq F value
## group       1  876.93  876.93 343.145
## time        2 2830.33 1415.17 553.761
## group:time  2  200.33  100.17  39.196

Við fáum sömu F gildi. Áður bentum við á af hverju p-gildunum er ekki skilað en það eru tæknilegar skýringar á því.

Athugasemdir