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

FervikagreiningBakgrunnur

Hér munun við sýna hvernig má framkvæma einhliða fervika greiningu fyrir endurteknar mælingar. Fyrir gögn sem ekki eru með endurteknum mælingum þá má lesa um einhliða fervikagreiningu (e. one-way ANOVA = “analysis of variance”) hér og sjá myndband hér. R sýnikóða fyrir einhliða fervikagreiningu má finna hér og fyrir tvíhliða (e. two-way) fervikagreiningu  hér. Fervikagreining er notuð til að bera saman meðaltöl milli hópa.

Við getum notað einhliða (e. one-way) ANOVA fyrir endurteknar mælingar þegar við mælum alla einstaklinga við öll gildi á einni flokkabreytu, t.d. þegar við mælum alla einstaklinga á öllum tímapunktum eða eftir að hafa fengið allar gerðir af lyfjum. Þar sem hver einstaklingur hefur mælingar á öllum tímapunktum, þá kallast tímapunktar “cross-over factor”.

Þegar framkvæmd er einhliða ANOVA fyrir endurteknar mælingar þá er mikilvægt að allir einstaklingar hafi mælingar á öllum tímapunktum (í því tilfelli gefa allar aðferðirnar fyrir neðan sömu niðurstöðu). Ef einhver einstaklingur hefur ekki allar mælingar þá þarf annað hvort að nota aðrar tölfræðiaðferðir, sleppa öllum mælingum á þessum einstakling eða átta sig á nákvæmlega hvað hver aðferð prófar og nota það sem svarar rannsóknaspurningunni best.

Fyrir neðan sýnum við nokkrar leiðir til að framkvæma einhliða fervikagreiningu fyrir enduteknar mælingar í R. Síðar sýnum við hvernig má framkvæma tvíhlið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.

demo2 <- read.csv("http://www.ats.ucla.edu/stat/data/demo2.csv")
demo2 <- within(demo2, {
  group <- factor(group)
  time <- factor(time)
  id <- factor(id)
})
summary(demo2)
##        id    group      pulse       time 
##  1      :3   1:12   Min.   :10.00   1:8  
##  2      :3   2:12   1st Qu.:16.00   2:8  
##  3      :3          Median :22.50   3:8  
##  4      :3          Mean   :22.38        
##  5      :3          3rd Qu.:26.75        
##  6      :3          Max.   :35.00        
##  (Other):6
str(demo2)
## '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  14 19 29 15 25 26 16 16 31 12 ...
##  $ time : Factor w/ 3 levels "1","2","3": 1 2 3 1 2 3 1 2 3 1 ...

1. Notum aov fallið

demo2.aov <- aov(pulse ~ time  + Error(id), data = demo2)
demo2.aov<-summary(demo2.aov)
demo2.aov$`Error: Within`
##           Df Sum Sq Mean Sq F value  Pr(>F)    
## time       2  978.2   489.1   62.02 1.1e-07 ***
## Residuals 14  110.4     7.9                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Hér að ofan sjáum við að p=1.1e-7 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.

2. Notum Anova fallið í car pakkanum

library(car)
library(reshape2)

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

# prepare arguments for Anova function
# model
mod <- lm(cbind(t1,t2,t3)~1,demo2.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
demo2.Anova <- Anova(mod, idata=time.frame,idesign=~time.factor, type=2)
## Note: model has only an intercept; equivalent type-III tests substituted.
demo2.Anova <- summary(demo2.Anova)
demo2.Anova$univariate.tests
##                  SS num Df Error SS den Df       F    Pr(>F)    
## (Intercept) 12015.4      1   122.96      7 684.034 3.057e-08 ***
## time.factor   978.2      2   110.42     14  62.017 1.104e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Aftur fáum við sömu niðurstöðu.

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 jafnmargi einstaklingar með mælingar á öllum tímapunktum.

library(lme4)
form <- as.formula("pulse~time + (1|id)")
demo2.mixed <- lmer(form,data=demo2,REML = T)
x<-summary(demo2.mixed)
y<-anova(demo2.mixed)
y
## Analysis of Variance Table
##      Df Sum Sq Mean Sq F value
## time  2 978.25  489.12  62.017

Við fáum sama F gildi. Það er með ráðum gert að frí gráðurnar eru ekki reiknaðar (sjá tæknilegar útskýringar hér). Við getum hins vegar borið F=62.0173561 saman við F dreifingu með frígráðurnar 2 og 14, vegna þess að við gögnin eru “balanced”, þ.e. allir 8 einstaklingar eru með mælingar á öllum 3 tímapunktum. Þá fáum við sömu niðurstöðu.

Athugasemdir