Anova Examples
Daniel
January 13, 2016
This document lists codes that are used to run ANOVA (Analysis of Variance) procedures
The examples include:
- One-way ANOVA
- One-way ANCOVA
- Two-way Factorial Design
library(knitr)
opts_chunk$set(message=FALSE, warning=FALSE, comment=NA)Load Libraries
library(gplots)
library(multcomp)
library(car)One-way ANOVA
One-way ANOVA is used to test:
- Equal population means among more than two groups
- Equal population variance between two groups
Regular One-Way ANOVA
attach(cholesterol)
aggregate(response, by=list(trt),mean) Group.1 x
1 1time 5.78197
2 2times 9.22497
3 4times 12.37478
4 drugD 15.36117
5 drugE 20.94752aggregate(response, by=list(trt),sd) Group.1 x
1 1time 2.878113
2 2times 3.483054
3 4times 2.923119
4 drugD 3.454636
5 drugE 3.345003fit <- aov(response~trt)
plotmeans(response~trt,
xlab = 'Treatment',
ylab = 'Response',
main = 'Mean Plot\n with 95% CI')
detach(cholesterol)Multiple Comparison
ANOVA can be used to tell whether there are differences among more than two treatment groups. But it does not tell which groups differ.
A multiple Comparison Procedure, like Tukey’s Method, will make pairwise comparisons between levels.
TukeyHSD(fit) Tukey multiple comparisons of means
95% family-wise confidence level
Fit: aov(formula = response ~ trt)
$trt
diff lwr upr p adj
2times-1time 3.44300 -0.6582817 7.544282 0.1380949
4times-1time 6.59281 2.4915283 10.694092 0.0003542
drugD-1time 9.57920 5.4779183 13.680482 0.0000003
drugE-1time 15.16555 11.0642683 19.266832 0.0000000
4times-2times 3.14981 -0.9514717 7.251092 0.2050382
drugD-2times 6.13620 2.0349183 10.237482 0.0009611
drugE-2times 11.72255 7.6212683 15.823832 0.0000000
drugD-4times 2.98639 -1.1148917 7.087672 0.2512446
drugE-4times 8.57274 4.4714583 12.674022 0.0000037
drugE-drugD 5.58635 1.4850683 9.687632 0.0030633par(las=2)
par(mar=c(5,8,4,2))
plot(TukeyHSD(fit))
Checking Assumptions
Anova has the following assumptions:
- The dependent variable is normally distributed
- The dependent variable has equal variances among different treatment groups
- ANOVA is sensitive to outliers, so we need to check outlier presence
QQ plot testing Normality
qqPlot(lm(response~trt,data=cholesterol),
simulate=TRUE, main='Q-Q plot',labels=FALSE)
Bartlett Test for equal variance
bartlett.test(response ~ trt, data=cholesterol)
Bartlett test of homogeneity of variances
data: response by trt
Bartlett's K-squared = 0.5797, df = 4, p-value = 0.9653The test result p-value is 0.97, which is much lower than most threshold. We don’t reject the Null hypothesis that The population variances are equal
Outlier Test
outlierTest(fit)
No Studentized residuals with Bonferonni p < 0.05
Largest |rstudent|:
rstudent unadjusted p-value Bonferonni p
19 2.251149 0.029422 NASince p-value is low, we don’t reject the Null that there are no outliers
One-way ANCOVA
One-way Analysis of Covariance (ANCOVA) is a modified version of ANOVA. The difference is that ANCOVA has included one or more quantitative covariates, which we are not interested but affects the dependent
ANCOVA procedure
data(litter,package='multcomp')
attach(litter)
aggregate(weight,
by=list(dose),
FUN=mean) Group.1 x
1 0 32.30850
2 5 29.30842
3 50 29.86611
4 500 29.64647fit <- aov(weight~gesttime+dose)
summary(fit) Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.049 0.00597 **
dose 3 137.1 45.71 2.739 0.04988 *
Residuals 69 1151.3 16.69
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Assessing Assumptions
ANCOVA makes the following assumptions about the data:
- The covariate’s effects are the same in all different treatment groups
- Normality of dependent variable
- Equal variances
Equal variance and Normality can be tested using the same method as ANOVA
Homogeneity of Covariate effects can be tested via interaction term. A significant interaction implie the effects of treatements are different at different covariate levels
fit2 <- aov(weight~gesttime*dose, data=litter)
summary(fit2) Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.3 134.30 8.289 0.00537 **
dose 3 137.1 45.71 2.821 0.04556 *
gesttime:dose 3 81.9 27.29 1.684 0.17889
Residuals 66 1069.4 16.20
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since the p-value is high, we don’t reject the hypothesis that there is no effect of the interaction.
Visualizing the results
library(HH)
ancova(weight~gesttime+dose, data=litter)Analysis of Variance Table
Response: weight
Df Sum Sq Mean Sq F value Pr(>F)
gesttime 1 134.30 134.304 8.0493 0.005971 **
dose 3 137.12 45.708 2.7394 0.049883 *
Residuals 69 1151.27 16.685
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
The slopes from the graph, for different groups, are close or similar
Two-way Factorial ANOVA
In two way factorial ANOVA, subjects are assigned to two-way cross groups (for instance, 2 factors, each with 2 levels will produce 2x2 groups)
Not like ANCOVA, we are interested in both of the two factors, and both of them are categorical variables
Balanced design is the design that each factor-level combination has same sample size
attach(ToothGrowth)
table(supp,dose) dose
supp 0.5 1 2
OJ 10 10 10
VC 10 10 10aggregate(len,by=list(supp,dose),FUN=mean) Group.1 Group.2 x
1 OJ 0.5 13.23
2 VC 0.5 7.98
3 OJ 1.0 22.70
4 VC 1.0 16.77
5 OJ 2.0 26.06
6 VC 2.0 26.14aggregate(len, by=list(supp, dose), FUN=sd) Group.1 Group.2 x
1 OJ 0.5 4.459709
2 VC 0.5 2.746634
3 OJ 1.0 3.910953
4 VC 1.0 2.515309
5 OJ 2.0 2.655058
6 VC 2.0 4.797731dose <- factor(dose)
fit <- aov(len~supp*dose)
summary(fit) Df Sum Sq Mean Sq F value Pr(>F)
supp 1 205.4 205.4 15.572 0.000231 ***
dose 2 2426.4 1213.2 92.000 < 2e-16 ***
supp:dose 2 108.3 54.2 4.107 0.021860 *
Residuals 54 712.1 13.2
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The summary indicates that both main effects of supp and dose and the interaction effect are significant
Visualize interactions
interaction.plot(dose, supp,
len, type='b',
col=c('red','blue'), pch=c(16,18),
main="Interaction between Dose and Supplement Types")
Repeated Measure ANOVA
Repeated measure ANOVA is a modified version of regular ANOVA, in which a subjects’ measurement is taken at multiple time points or under multiple conditions.
In regular ANOVA, the F-statistic is calculated as the ratio between ‘Between Group Variance’ and ‘Within Group Variance’. In repeated measure ANOVA, the ‘Within Group Variance’ is decomposed into ‘Among different subjects’ and ‘In a subject itself’. We will use ‘Among different subjects’ as the denominator, which will increase the F-statistic and thus make it **easier* to reject Null
data(CO2)
CO2$conc <- factor(CO2$conc)
w1b1 <- subset(CO2, Treatment=='chilled')
w1b1 Plant Type Treatment conc uptake
22 Qc1 Quebec chilled 95 14.2
23 Qc1 Quebec chilled 175 24.1
24 Qc1 Quebec chilled 250 30.3
25 Qc1 Quebec chilled 350 34.6
26 Qc1 Quebec chilled 500 32.5
27 Qc1 Quebec chilled 675 35.4
28 Qc1 Quebec chilled 1000 38.7
29 Qc2 Quebec chilled 95 9.3
30 Qc2 Quebec chilled 175 27.3
31 Qc2 Quebec chilled 250 35.0
32 Qc2 Quebec chilled 350 38.8
33 Qc2 Quebec chilled 500 38.6
34 Qc2 Quebec chilled 675 37.5
35 Qc2 Quebec chilled 1000 42.4
36 Qc3 Quebec chilled 95 15.1
37 Qc3 Quebec chilled 175 21.0
38 Qc3 Quebec chilled 250 38.1
39 Qc3 Quebec chilled 350 34.0
40 Qc3 Quebec chilled 500 38.9
41 Qc3 Quebec chilled 675 39.6
42 Qc3 Quebec chilled 1000 41.4
64 Mc1 Mississippi chilled 95 10.5
65 Mc1 Mississippi chilled 175 14.9
66 Mc1 Mississippi chilled 250 18.1
67 Mc1 Mississippi chilled 350 18.9
68 Mc1 Mississippi chilled 500 19.5
69 Mc1 Mississippi chilled 675 22.2
70 Mc1 Mississippi chilled 1000 21.9
71 Mc2 Mississippi chilled 95 7.7
72 Mc2 Mississippi chilled 175 11.4
73 Mc2 Mississippi chilled 250 12.3
74 Mc2 Mississippi chilled 350 13.0
75 Mc2 Mississippi chilled 500 12.5
76 Mc2 Mississippi chilled 675 13.7
77 Mc2 Mississippi chilled 1000 14.4
78 Mc3 Mississippi chilled 95 10.6
79 Mc3 Mississippi chilled 175 18.0
80 Mc3 Mississippi chilled 250 17.9
81 Mc3 Mississippi chilled 350 17.9
82 Mc3 Mississippi chilled 500 17.9
83 Mc3 Mississippi chilled 675 18.9
84 Mc3 Mississippi chilled 1000 19.9fit <- aov(uptake~conc*Type + Error(Plant/(conc)),data=w1b1)
summary(fit)
Error: Plant
Df Sum Sq Mean Sq F value Pr(>F)
Type 1 2667.2 2667.2 60.41 0.00148 **
Residuals 4 176.6 44.1
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Error: Plant:conc
Df Sum Sq Mean Sq F value Pr(>F)
conc 6 1472.4 245.40 52.52 1.26e-12 ***
conc:Type 6 428.8 71.47 15.30 3.75e-07 ***
Residuals 24 112.1 4.67
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Error() accounts for the effects of repeated measure