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í.
