One-Way ANOVA

In this post I want to exemplify how to implement a simple One-Way ANOVA in R (the R code is based on Field, Miles, and Field 2012).

An ANOVA (Analysis of Variance) is used when we have a numeric dependent variable and at least one categorical independent variable – if we only have a nominal independent variable, we would use a t-test. Thus, an ANOVA is in essence an extension of a t-test (which is used to test if two means differ significantly or not) but it can also be understood as a linear regression model (mathematically, ANOVAs and linear regressions are almost identical by the way).

As this post is concerned only with the implementation you should read up on the underlying reasoning and workings of an ANOVA elsewhere – I can highly recommend chapter 10 in Field, Miles, and Field (2012)! However, ANOVAs are commonly used to evaluate experimental data where different groups are compared with respect to some outcome (e.g. how long it takes to recover from an illness for a placebo, and two treatment groups – one of which received a high dose and the other a low dose of a drug). The ANOVA then tests if the residual deviance explained by the model is greater than the residual deviance not explained by the model. The residual deviance is the difference between the model prediction and the observed values (see the upper panel of the figure below).

The F-ratio is pretty much just the ratio of the residual deviance we get, if we only use the overall mean as a predictor, divided by the residual deviance we obtain, if we use the group means as predictors (It is a little more complex than this but I think it gives you a good general understanding of what ANOVAs do).

Of course, there are more complex variants of ANOVAs but this is basically the core of an ANOVA. And we will consider more complex scenarios below, but let us start with a rather simple and straight forward example to get started.

One-Way ANOVA

The example we are will be dealing with here is a prototypical and simple experimental design: we are investigating, if alcohol consumption increases speech errors and if more alcohol leads to more mistakes. Participants were assigned to one of three groups and all participants had to recite a tongue twister. Our dependent variable is the number of mistakes they made during the recital. One group did not have any alcohol, the second group drank 400ml beer half an hour before the recital, and the third group drank 400ml of tequila half an hour before the recital.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
###############################################################
### One-way ANOVA
###############################################################
### --- Prepare data
# Remove all lists from the current workspace
rm(list=ls(all=T))
 
# Install packages we need or which may be useful
# (to activate just delete the #)
#install.packages("granova")
#install.packages("car")
#install.packages("pastecs")
#install.packages("multcomp")
#install.packages("compute.es")
#install.packages("sciplot")
# To install the WS package is a little more hustle
# first: install dependent packages
#install.packages(c("MASS", "akima", "robustbase"))
# second: install suggested packages
#install.packages(c("cobs", "robust", "mgcv", "scatterplot3d", "quantreg", "rrcov", "lars", "pwr", "trimcluster", "parallel", "mc2d", "psych", "Rfit"))
# third: install WRS
#install.packages("WRS", repos="http://R-Forge.R-project.org", type="source")
 
# Initiate the packages
library(ggplot2)
library(granova)
library(car)
library(pastecs)
library(multcomp)
library(compute.es)
library(WRS)
library(sciplot)
# load functions for robust alternatives to ANOVA
source("http://www-rcf.usc.edu/~rwilcox/Rallfun-v13") 
###############################################################
### Load and manipulate data
###############################################################
# WARNING: To use this script you need to set our own paths! 
# Your path should be the path to the corpus on your own computer! 
# Remember to use double backslash instead of single backslash, if 
# you use Windows on your machine. 
 
# Read in data
anova.data <- read.delim("C:\\MyProjects\\ANOVA/anova-data-1.txt", header = TRUE)
attach(anova.data) # we attach the data so we don't need to specify the path all the time
anova.data$group <- factor(anova.data$group,  levels = c(1:3), labels = c("None", "Beer", "Tequila"))
 
###############################################################
### Visualize and eyeball data
###############################################################
par(mfrow = c(2, 3)) # plot two plots in one window
 
# Plot residuals from overall mean
plot(anova.data$errors, ylab = "Number of Errors")
abline(mean(anova.data$errors),0)
for(i in 1:length(anova.data$errors)) lines(c(i,i), c(mean(anova.data$errors), anova.data$errors[i]))
###############################################################
 
 
# Plot group means
plot(anova.data$errors, ylab = "Number of Errors")
lines(1:5, rep(mean(anova.data$errors[anova.data$group == "None"]), 5))
lines(6:10, rep(mean(anova.data$errors[anova.data$group == "Beer"]), 5))
lines(11:15, rep(mean(anova.data$errors[anova.data$group == "Tequila"]), 5))
 
# Plot residuals from group means
for(i in 1:length(anova.data$errors[anova.data$group == "None"])) lines(c(i, i), c(mean(anova.data$errors[anova.data$group == "None"]), anova.data$errors[anova.data$group == "None"] [i]))
for(i in 1:length(anova.data$errors[anova.data$group == "Beer"])) lines(c(i+5, i+5), c(mean(anova.data$errors[anova.data$group == "Beer"]), anova.data$errors[anova.data$group == "Beer"] [i]))
for(i in 1:length(anova.data$errors[anova.data$group == "Tequila"])) lines(c(i+10, i+10), c(mean(anova.data$errors[anova.data$group == "Tequila"]), anova.data$errors[anova.data$group == "Tequila"] [i]))
 
###############################################################
# Plot group means
plot(anova.data$errors, ylab = "Number of Errors")
abline(mean(anova.data$errors),0)
lines(1:5, rep(mean(anova.data$errors[anova.data$group == "None"]), 5))
lines(6:10, rep(mean(anova.data$errors[anova.data$group == "Beer"]), 5))
lines(11:15, rep(mean(anova.data$errors[anova.data$group == "Tequila"]), 5))
 
# Plot residuals from group means
for(i in 1:length(anova.data$errors[anova.data$group == "None"])) lines(c(i, i), c(mean(anova.data$errors[anova.data$group == "None"]), rep(mean(anova.data$errors), 5) [i]))
 
for(i in 1:length(anova.data$errors[anova.data$group == "Beer"])) lines(c(i+5, i+5), c(mean(anova.data$errors[anova.data$group == "Beer"]), rep(mean(anova.data$errors), 5) [i]))
 
for(i in 1:length(anova.data$errors[anova.data$group == "Tequila"])) lines(c(i+10, i+10), c(mean(anova.data$errors[anova.data$group == "Tequila"]), rep(mean(anova.data$errors), 5) [i]))
 
 
###############################################################
# Linegraph of means
plot(tapply(anova.data$errors, anova.data$group, mean),
  type = "b",
  ylim = c(1, 7),
  xlim = c(0.5, 3.5),
  ann = F,
  axes = F)
box()
axis(side = 1, at = c(1:3), labels = F)
axis(side = 2, at = c(1:7))
mtext(text = c("None", "Beer", "Tequila"), side = 1, line = 1, at = c(1:3), cex = .8)
mtext(text = "Number of Errors", side = 2, line = 3, at = 4, cex = .8)
arrows(1, mean(anova.data$errors[anova.data$group == "None"]), 1,
  mean(anova.data$errors[anova.data$group == "None"] + se(anova.data$errors[anova.data$group == "None"])),
  length = .1, angle = 90)
arrows(1, mean(anova.data$errors[anova.data$group == "None"]), 1,
  mean(anova.data$errors[anova.data$group == "None"] - se(anova.data$errors[anova.data$group == "None"])),
  length = .1, angle = 90)
arrows(2, mean(anova.data$errors[anova.data$group == "Beer"]), 2,
  mean(anova.data$errors[anova.data$group == "Beer"] + se(anova.data$errors[anova.data$group == "Beer"])),
  length = .1, angle = 90)
arrows(2, mean(anova.data$errors[anova.data$group == "Beer"]), 2,
  mean(anova.data$errors[anova.data$group == "Beer"] - se(anova.data$errors[anova.data$group == "Beer"])),
  length = .1, angle = 90)
arrows(3, mean(anova.data$errors[anova.data$group == "Tequila"]), 3,
  mean(anova.data$errors[anova.data$group == "Tequila"] + se(anova.data$errors[anova.data$group == "Tequila"])),
  length = .1, angle = 90)
arrows(3, mean(anova.data$errors[anova.data$group == "Tequila"]), 3,
  mean(anova.data$errors[anova.data$group == "Tequila"] - se(anova.data$errors[anova.data$group == "Tequila"])),
  length = .1, angle = 90)
grid()
 
###############################################################
# Barplot with CI
bargraph.CI(group, errors, data = anova.data, xlab = "",
  ylab = "Number of Errors", ylim = c(1, 7),
  cex.lab = 1, x.leg = 1,
  col = gray.colors(3), angle = 45, cex.names = 1.25)
box()
grid()
 
###############################################################
# Boxplots of errors as a function of groups
boxplot(anova.data$errors ~ anova.data$group, notch = T, col = "lightgrey", ylab = "Number of Errors")
grid ()
 
par(mfrow = c(1, 1)) # restore original plotting parameters
Figure 1: Plots showing teh data structure.

Figure 1: Plots showing teh data structure.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
###############################################################
###############################################################
# Retrieve descriptive statistics
by(anova.data$errors, anova.data$group, stat.desc)
 
#anova.data$group: None
#     nbr.val     nbr.null       nbr.na          min          max        range 
#   5.0000000    0.0000000    0.0000000    1.0000000    4.0000000    3.0000000 
#         sum       median         mean      SE.mean CI.mean.0.95          var 
#  11.0000000    2.0000000    2.2000000    0.5830952    1.6189318    1.7000000 
#     std.dev     coef.var 
#   1.3038405    0.5926548 
#------------------------------------------------------------ 
#anova.data$group: Beer
#     nbr.val     nbr.null       nbr.na          min          max        range 
#   5.0000000    0.0000000    0.0000000    2.0000000    5.0000000    3.0000000 
#         sum       median         mean      SE.mean CI.mean.0.95          var 
#  16.0000000    3.0000000    3.2000000    0.5830952    1.6189318    1.7000000 
#     std.dev     coef.var 
#   1.3038405    0.4074502 
#------------------------------------------------------------ 
#anova.data$group: Tequila
#     nbr.val     nbr.null       nbr.na          min          max        range 
#   5.0000000    0.0000000    0.0000000    3.0000000    7.0000000    4.0000000 
#         sum       median         mean      SE.mean CI.mean.0.95          var 
#  25.0000000    5.0000000    5.0000000    0.7071068    1.9632432    2.5000000 
#     std.dev     coef.var 
#   1.5811388    0.3162278 
 
 
###############################################################
### Check conditions
###############################################################
# Perform Levene's test to check if variance homogeneiety is given
# (ok if p > .05)
leveneTest(errors, group, center = median) 
 
#Levene's Test for Homogeneity of Variance (center = median)
#      Df F value Pr(>F)
#group  2  0.1176   0.89
#      12               
#Warnmeldung:
#In leveneTest.default(errors, group, center = median) :
#  group coerced to factor.
 
###############################################################
### Perform ANOVA
###############################################################
anova.aov <- aov(errors ~ group, data = anova.data)
 
# anova.lm <- lm(errors ~ group, data = anova.data)
 
summary(anova.aov)
 
#            Df Sum Sq Mean Sq F value Pr(>F)  
#group        2  20.13  10.067   5.119 0.0247 *
#Residuals   12  23.60   1.967                 
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 
summary.lm(anova.aov)
 
#Call:
#aov(formula = errors ~ group, data = anova.data)
#
#Residuals:
#   Min     1Q Median     3Q    Max 
#  -2.0   -1.2   -0.2    0.9    2.0 
#
#Coefficients:
#             Estimate Std. Error t value Pr(>|t|)   
#(Intercept)    2.2000     0.6272   3.508  0.00432 **
#groupBeer      1.0000     0.8869   1.127  0.28158   
#groupTequila   2.8000     0.8869   3.157  0.00827 **
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#
#Residual standard error: 1.402 on 12 degrees of freedom
#Multiple R-squared:  0.4604,    Adjusted R-squared:  0.3704 
#F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469
 
par(mfrow = c(2,2))
plot(anova.aov)
par(mfrow = c(1,1))
Figure 2: Diagnostic plots for the ANOVA.

Figure 2: Diagnostic plots for the ANOVA.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
 
###############################################################
### Robust alternatives to ANOVA
###############################################################
# If the conditions for ANOVA are not met, we need to use
# robust alternatives which are not affected by e.g. unequal variances.
# In the following we will go through some of these alternatives.
 
### --- Welch Test
# If Lavene's test is significant, i.e. in case we are dealing
# with unequal variances, we can use the Welch Test which combat 
# problems arising from violations of the homogeneity of
# variance assumption (Field, Miles & Field 2012: 414)
oneway.test(errors ~ group, data = anova.data)
 
#        One-way analysis of means (not assuming equal variances)
#
#data:  errors and group
#F = 4.3205, num df = 2.000, denom df = 7.943, p-value = 0.05374
 
###############################################################
# We can also take an approach which trims the means but before we can use
# this approach, we need to re-organize our data:
 
anova.wide <- unstack(anova.data, errors ~ group) # Prepare data
granova.1w(anova.wide)
 
#$grandsum
#    Grandmean        df.bet       df.with        MS.bet       MS.with 
#         3.47          2.00         12.00         10.07          1.97 
#       F.stat        F.prob SS.bet/SS.tot 
#         5.12          0.02          0.46 
#
#$stats
#        Size Contrast Coef Wt'd Mean Mean Trim'd Mean Var. St. Dev.
#None       5         -1.27       2.2  2.2           2  1.7     1.30
#Beer       5         -0.27       3.2  3.2           3  1.7     1.30
#Tequila    5          1.53       5.0  5.0           5  2.5     1.58
Figure 3: Graphical summary of the newly arranged data.

Figure 3: Graphical summary of the newly arranged data.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
 
# We can now continue with the analysis on trimmed means
t1way(anova.wide, tr = .1)
 
#$TEST
#[1] 4.320451
#
#$nu1
#[1] 2
#
#$nu2
#[1] 7.943375
#
#$n
#$n[[1]]
#[1] 5
#
#$n[[2]]
#[1] 5
#
#$n[[3]]
#[1] 5
#
#$p.value
#[1] 0.05373847
 
###############################################################
# Instead of trimmed means, we can also consider differences between medians instead of means
med1way(anova.wide)
 
#[1] "NOTE: This function was modified in Dec. 2004"
#[1] "A new approximate critical value is used if crit=NA"
#[1] "This might improve type I error probabilities substantially"
#[1] "For discrete data with ties, this function is NOT recommended."
#[1] "Use the function medpb; it is best for general use"
#[1] "WARNING: tied values detected."
#[1] "Estimate of standard error might be highly inaccurate, even with n large"
#[1] "WARNING: tied values detected."
#[1] "Estimate of standard error might be highly inaccurate, even with n large"
#$TEST
#[1] 4.782879
#
#$crit.val
#[1] 5.472958
#
#$p.value
#[1] 0.07
 
###############################################################
# Alternatively, we can use bootstrapping.
t1waybt(anova.wide, tr = .2, nboot = 599)
 
#[1] "Taking bootstrap samples. Please wait."
#[1] "Working on group  1"
#[1] "Working on group  2"
#[1] "Working on group  3"
#[1] "Some bootstrap estimates of the test statistic could not be computed"
#[1] "Effective number of bootstrap samples was"
#[1] 395
#$test
#[1] 3
#
#$p.value
#[1] 0.08860759
 
 
###############################################################
### Include Planned Contrasts
###############################################################
# Contrasts are a handy way to inspect if there are differences
# not only between the groups, but e.g. between those that have
# not had any alcohol and those who had. 
 
summary.lm(anova.aov)
 
contrast1 <- c(-2, 1, 1)
contrast2 <- c(0, -1, 1)
 
contrasts(anova.data$group) <- cbind(contrast1, contrast2)
 
contrasts(anova.data$group) <- cbind(c(-2, 1, 1), c(0, -1, 1))
anova.aov.cont <- aov(errors ~ group, data = anova.data)
summary.lm(anova.aov.cont)
 
#Call:
#aov(formula = errors ~ group, data = anova.data)
#Residuals:
#   Min     1Q Median     3Q    Max 
 # -2.0   -1.2   -0.2    0.9    2.0 
 
#Coefficients:
 #           Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
#group1        0.6333     0.2560   2.474   0.0293 *  
#group2        0.9000     0.4435   2.029   0.0652 .  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 1.402 on 12 degrees of freedom
#Multiple R-squared:  0.4604,    Adjusted R-squared:  0.3704 
#F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469
 
 
contrasts(anova.data$group) <- contr.poly(3)
anova.aov.trend <- aov(errors ~ group, data = anova.data)
summary.lm(anova.aov.trend)
 
#Call:
#aov(formula = errors ~ group, data = anova.data)
#Residuals:
#   Min     1Q Median     3Q    Max 
#  -2.0   -1.2   -0.2    0.9    2.0 
#Coefficients:
 #           Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
#group.L       1.9799     0.6272   3.157  0.00827 ** 
#group.Q       0.3266     0.6272   0.521  0.61201    
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 1.402 on 12 degrees of freedom
#Multiple R-squared:  0.4604,    Adjusted R-squared:  0.3704 
#F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469
 
 
contrasts(anova.data$group) <- contr.helmert(3)
anova.aov.helmert <- aov(errors ~ group, data = anova.data)
summary.lm(anova.aov.helmert)
 
# Call:
#aov(formula = errors ~ group, data = anova.data)
#Residuals:
#   Min     1Q Median     3Q    Max 
#  -2.0   -1.2   -0.2    0.9    2.0 
#Coefficients:
#            Estimate Std. Error t value Pr(>|t|)    
#(Intercept)   3.4667     0.3621   9.574 5.72e-07 ***
#group1        0.5000     0.4435   1.127   0.2816    
#group2        0.7667     0.2560   2.994   0.0112 *  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#Residual standard error: 1.402 on 12 degrees of freedom
#Multiple R-squared:  0.4604,    Adjusted R-squared:  0.3704 
#F-statistic: 5.119 on 2 and 12 DF,  p-value: 0.02469
 
 
# An alternative to setting contrasts consists of performing
# post-hoc tests. We will now have a look at how to implement
# such post hoc test in R.
 
pairwise.t.test(anova.data$errors, anova.data$group, p.adjust.method = "bonferroni")
 
#        Pairwise comparisons using t tests with pooled SD 
#data:  anova.data$errors and anova.data$group 
#        None  Beer 
#Beer    0.845 -    
#Tequila 0.025 0.196
#P value adjustment method: bonferroni
 
pairwise.t.test(anova.data$errors, anova.data$group, p.adjust.method = "BH")
 
#        Pairwise comparisons using t tests with pooled SD 
#data:  anova.data$errors and anova.data$group 
#        None  Beer 
#Beer    0.282 -    
#Tequila 0.025 0.098
#P value adjustment method: BH
 
postHocs <- glht(anova.aov, linfct = mcp(group = "Tukey"))
summary(postHocs)
 
#         Simultaneous Tests for General Linear Hypotheses
#Multiple Comparisons of Means: Tukey Contrasts
#Fit: aov(formula = errors ~ group, data = anova.data)
#Linear Hypotheses:
 #                   Estimate Std. Error t value Pr(>|t|)  
#Beer - None == 0      1.0000     0.8869   1.127   0.5164  
#Tequila - None == 0   2.8000     0.8869   3.157   0.0208 *
#Tequila - Beer == 0   1.8000     0.8869   2.029   0.1475  
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#(Adjusted p values reported -- single-step method)
 
confint(postHocs)
 
#         Simultaneous Confidence Intervals
#Multiple Comparisons of Means: Tukey Contrasts
#Fit: aov(formula = errors ~ group, data = anova.data)
#Quantile = 2.6718
#95% family-wise confidence level
 #Linear Hypotheses:
 #                  Estimate lwr     upr    
#Beer - None == 0     1.0000  -1.3698  3.3698
#Tequila - None == 0  2.8000   0.4302  5.1698
#Tequila - Beer == 0  1.8000  -0.5698  4.1698
 
 
postHocs <- glht(anova.aov, linfct = mcp(group = "Dunnett"), base = 1)
summary(postHocs)
 
#         Simultaneous Tests for General Linear Hypotheses
#Multiple Comparisons of Means: Dunnett Contrasts
#Fit: aov(formula = errors ~ group, data = anova.data)
#Linear Hypotheses:
#                    Estimate Std. Error t value Pr(>|t|)  
#Beer - None == 0      1.0000     0.8869   1.127   0.4459  
#Tequila - None == 0   2.8000     0.8869   3.157   0.0152 *
#---
#Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#(Adjusted p values reported -- single-step method)
 
confint(postHocs)
#         Simultaneous Confidence Intervals
#Multiple Comparisons of Means: Dunnett Contrasts
#Fit: aov(formula = errors ~ group, data = anova.data)
#Quantile = 2.5023
#95% family-wise confidence level
#Linear Hypotheses:
#                    Estimate lwr     upr    
#Beer - None == 0     1.0000  -1.2194  3.2194
#Tequila - None == 0  2.8000   0.5806  5.0194
 
lincon(anova.wide, tr = .1)
 
#[1] "Note: confidence intervals are adjusted to control FWE"
#[1] "But p-values are not adjusted to control FWE"
#$n
#[1] 5 5 5
#
#$test
#     Group Group     test     crit        se       df
#[1,]     1     2 1.212678 2.960000 0.8246211 8.000000
#[2,]     1     3 3.055050 2.985397 0.9165151 7.719912
#[3,]     2     3 1.963961 2.985397 0.9165151 7.719912
#
#$psihat
#     Group Group psihat  ci.lower    ci.upper    p.value
#[1,]     1     2   -1.0 -3.440879  1.44087853 0.25984505
#[2,]     1     3   -2.8 -5.536161 -0.06383861 0.01638290
#[3,]     2     3   -1.8 -4.536161  0.93616139 0.08644456
 
 
mcppb20(anova.wide, tr = .2, nboot = 2000)
 
#[1] "Taking bootstrap samples. Please wait."
#$psihat
#     con.num psihat       se  ci.lower   ci.upper p-value
#[1,]       1     -1 1.154701 -3.333333  1.3333333  0.3250
#[2,]       2     -3 1.154701 -5.333333 -0.3333333  0.0055
#[3,]       3     -2 1.154701 -4.333333  0.6666667  0.0840
#
#$crit.p.value
#[1] 0.017
#
#$con
#     [,1] [,2] [,3]
#[1,]    1    1    0
#[2,]   -1    0    1
#[3,]    0   -1   -1
 
lincon(anova.wide)
 
#[1] "Note: confidence intervals are adjusted to control FWE"
#[1] "But p-values are not adjusted to control FWE"
#$n
#[1] 5 5 5
#
#$test
#     Group Group      test crit       se df
#[1,]     1     2 0.8660254 3.74 1.154701  4
#[2,]     1     3 2.5980762 3.74 1.154701  4
#[3,]     2     3 1.7320508 3.74 1.154701  4
#
#$psihat
#     Group Group psihat ci.lower ci.upper    p.value
#[1,]     1     2     -1 -5.31858  3.31858 0.43533094
#[2,]     1     3     -3 -7.31858  1.31858 0.06016985
#[3,]     2     3     -2 -6.31858  2.31858 0.15830242
 
mcppb20(anova.wide)
 
#[1] "Taking bootstrap samples. Please wait."
#$psihat
#     con.num psihat       se  ci.lower   ci.upper p-value
#[1,]       1     -1 1.154701 -3.333333  1.3333333  0.3250
#[2,]       2     -3 1.154701 -5.333333 -0.3333333  0.0055
#[3,]       3     -2 1.154701 -4.333333  0.6666667  0.0840
#
#$crit.p.value
#[1] 0.017
#
#$con
#     [,1] [,2] [,3]
#[1,]    1    1    0
#[2,]   -1    0    1
#[3,]    0   -1   -1
 
# Last but not least, we will calculate the effect sizes.
# You can do this in a more efficient manner by writing
# a function that you can call to do the work for you. 
rcontrast <- function(t, df) {
  r <- sqrt(t^2/(t^2 + df))
  print(paste("r = ", r))  }
 
rcontrast(2.474, 12)
# [1] "r =  0.581182458413787"
 
rcontrast(2.029, 12)
# [1] "r =  0.505407970122564"
 
omega <- function(SSm, SSr, dfm, MSr) {
  SSt = SSm + SSr
  omega = (SSm-(dfm*MSr))/(SSt+MSr)
  print(paste("Omega-Squared: ", omega))  }
 
omega(20.133, 23.600, 2, 1.9667)
# [1] "Omega-Squared:  0.35447935106795"
 
omega(450.66,38.09, 5, 0.334)
# [1] "Omega-Squared:  0.918022262024519"
 
omega(4180.6, 4356.6, 3, 167.56)
# [1] "Omega-Squared:  0.422518254380362"

So that was it on simple independent ANOVA in R.

References

Field, Andy, Jeremy Miles & Zoe Field. 22012. Discovering statistics using R. London, Thousand oaks, CA, New Delhi, Singapore: SAGE.

Schreibe einen Kommentar

Deine E-Mail-Adresse wird nicht veröffentlicht. Erforderliche Felder sind mit * markiert.