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

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

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

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