The commands below apply to the freeware statistical environment called R (R Development Core Team 2010). Each set of commands can be copy-pasted directly into R. Example datasets can be copy-pasted into .txt files from Examples of Analysis of Variance and Covariance (Doncaster & Davey 2007). For a given design and dataset in the format of the linked example, the commands will work for any number of factor levels and observations per level.
Each numbered section below shows commands for a design to be analyzed either by the function 'aov' which assumes residuals have a normal distribution, or by the function 'glm' which takes any named error structure. An example of a dataset of proportions analysed with a binomial error structure is given for split-plot model 5.9.
Model descriptors in titles and comments use the notation of Doncaster & Davey (2007), with '|' signifying 'cross-factored with', and '(' signifying 'nested in'. Factor and variable names take uppercase A, B, C, S, P, Q, with a prime identifying a factor as random. The numbers of factor levels take lowercase a, b, c, s, p, q (with no replication of p and q in split-plots). The symbol 'ε' signifies full replication. The R commands for model formulae use notation 'A:B' to request analysis of the interaction of A with B; 'A*B' to request main effects and interaction (i.e., 'A + B + A:B' with 'a' levels of A and 'b' levels of B); 'A/B' to request the main effect A, and B nested in A (i.e., 'A + A:B' with 'a' levels of A and 'b' levels of B per level of A).
Doncaster, C. P. & Davey, A. J. H. (2007) Analysis of Variance and Covariance: How to Choose and Construct Models for the Life Sciences. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511611377. Example datasets: http://southampton.likn.co/~cpd/anovas/datasets/
R Development Core Team (2010). R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org.
- R commands for analysis of ANOVA and ANCOVA datasets
3 Fully replicated factorial designs
- Analysis of datasets for figures in Doncaster & Davey (2007)
- Power calculations for ANOVA designs
In the following command lines, built-in functions are highlighted in dark blue. Clicking on a function will link to explanatory help on its uses. These and other built-in functions and programming controls are summarised in the Short R Reference Card, and R Reference Card 2.0.
# this hash-tag symbol on a command line means that everything following it on the same line is comment and not instruction
help.start() # general introduction to the R language
help(anova) # help on anova
?anova # help on anova
??anova # help on help files for anova
?read.table # help on importing data
?par # help on graphical parameters
setwd("C:/R/ANOVA datasets in R")
getwd() # list the working directory
dir() # list files in the working directory
E <- pi*2/3 # create a vector E which gets the single-element value of pi*2/3
F <- c(1, -0.5, 5.2, 3) # create a vector F of numeric elements
G <- c("1", "one", ".") # create a vector G of character elements
H <- 1:5 # create a vector H of integers 1 through 5
J <- gl(4, 3) # create vector J of 3 replicates of each of factor levels "1" through "4"
K <- c(rep(max(as.numeric(B)), length(B))) # create a vector K containing the maximum value in vector B, repeated as many times as the length of B
YM <- as.vector(tapply(Y,list(A, B), mean)) # create a vector YM which gets the mean value of vector Y at each combination of levels of vectors A and B
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean))) # create a vector AM which gets the levels of factor A that correspond to each element of YM
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B), mean))) # create a vector BM which gets the levels of factor B that correspond to each element of YM
aovdatam <- data.frame(cbind(AM, BM, YM)) # bind together the three vectors AM, BM, YM, into a data frame object, here called 'aovdatam'
Hold the data for an analysis of variance in a 'data frame', comprising equal-length columns (the 'vectors'). The response measurements go in one column, one measurement per row, and each explanatory factor or variable takes a column, showing its level for the corresponding response measurement. Use NA for any missing data. See model 1.1 and other models for examples of databases suitable for reading into a data frame.
aovdata <- read.table("Model1_1.txt", header = T) # the named data frame (here called 'aovdata') gets (the two-character '<-' symbol) data read from the named text-file table, with column headers
attach(aovdata) # make the data frame accessible by name within the R session
aovdata # see the data frame of column vectors
aovdata <- read.table("Model1_1.csv", sep = ",", col.names = c("A", "Y"), colClasses = c("factor", "numeric"))
attach(aovdata)
ls() # list all objects
search() # list all attached data frames and packages
names(aovdata) # list the vector names in the data frame object
aovdata # view the contents of the named data frame or vector object
Y[as.numeric(A) > 1] # list the values of vector Y for which vector A takes values above 1
class(Y) # show the class of vector Y
length(A) # show the length of the vector A
levels(A) # list the levels of the factor A
nlevels(A) # number of levels of factor A
table(A) ; table(A, B) # tabulate the number of data rows at each level of factor A, or at each combination of levels of factors A and B
Y <- as.numeric(Y) ; A <- as.factor(A) # convert vector Y to numeric and vector A to a factor
Y1 <- log(Y) # create a vector Y1 that is the natural log of the vector Y
Y2 <- log10(Y) # create a vector Y2 that is the log to base 10 of the vector Y
Y3 <- asin(sqrt(Y))) # create a vector Y3 that is the arcsin-root transformation of the vector Y
summary(aovdata) # mean, median, quartiles, min-max of all vectors in the data frame
mean(Y) ; median(Y) # arithmetic mean and median of the vector
min(Y) ; max(Y) # minimum and maximum values in the vector
var(Y) # sample variance for values in the vector
tapply(Y, A, mean) # mean of vector Y at each level of factor A
tapply(Y, list(A, B), mean) # mean of vector Y at each combination of levels of factors A and B
tapply(Y, A, var)^0.5 # sample standard deviation of vector Y at each level of factor A
tapply(Y, A, length) # number of observations of Y at each level of A
coef(aov(Y ~ X)) # intercept and slope of linear regression of vector Y on numeric vector X
plot(aov(Y ~ A)) # plot residuals to an ANOVA in four consecutive graphs: Residuals vs Fitted, Normal Q-Q, Scale-Location, Constant Leverage
plot(fitted(aov(Y ~ A)), resid(aov(Y ~ A)), xlab = "Fitted values", ylab = "Residuals", main = "residuals vs Fitted") # plot a diagnostic check for heteroscedasticity in the residuals from the response Y to a factor or numeric A
qqnorm(resid(aov(Y ~ A)), main = "Rankit plot of residuals") # plot normal scores to check for skewness (convex bow for right skew, concave bow for left skew), kurtosis ('S' shape for flattened, inverted 'S' for peaked) and outliers in the residuals from the response Y to a factor or numeric A
hist(resid(aov(Y ~ A))) # histogram of residuals
bartlett.test(Y ~ A) # Bartlett test of homogeneity of variances amongst levels of factor A
bartlett.test(Y ~ interaction(A, B)) # Bartlett test of homogeneity of variances amongst cross-factors A and B
shapiro.test(resid(aov(Y ~ A))) # Shapiro-Wilk normality test of residuals
plot(A, Y) # box and whisker plot of the response Y to a factor A, showing for each level of A its median straddled by the box of 25-50% and 50-75% quartiles, which together make up the interquartile range, and whiskers of minimum and maximum values lying within 1.5 interquartile ranges below the bottom and above the top of the box, and locating any outliers beyond the whiskers
plot(A:B, Y, las = 1, xlab = "cross-factors A:B", ylab = "Response Y") # box and whisker plot of the response Y to cross-factored A:B, with reading direction of y-axis numbers aligned left to right, and specified axis labels
abline(h = 0.0) # add a reference line to the plot at Y = 0.0
interaction.plot(A, B, Y, las = 1, xlab = "Factor A", ylab = "Response Y", trace.label = "Factor B", xtick = TRUE) # interaction plot of cross-factored means
The suite of commands below will create a panel of b plots for the means of numeric variable Y at each level of factor A given factor B, from a data frame containing equal-length vectors Y, A with a levels, and B with b levels.
library(lattice)
xyplot(Y ~ A|B,
panel = function(x, y)
{
panel.xyplot(x, y) # plot Y against A at each level of B
panel.xyplot(x, y, type = "a") # join the means of consecutive levels of A at each level of B
})
Click here for a user-defined function that can be saved for future use as a script (e.g., use the RStudio menus: File - New script, copy-pasting it in and saving it as, for example, "PlotMeans.R"). Having run the script, a plot is obtained on each call of the function plot_means(factor, response, xlabel, ylabel, bars) for specified vectors of a factor and response, x-axis label, y-axis label, and type of error bars: SD or SE or CI. For example, plot the means +/-1 standard error given by the data provided for factorial model 1.1 by attaching the data frame and then calling:
source(file="http://southampton.likn.co/~cpd/anovas/datasets/PlotMeans.R") # run the script, here stored in the script file PlotMeans.R on the web
plot_means(A, Y, "Factor A", "Response Y", "SE") # plot means and s.e.
Means +/-95% confidence intervals given by the data provided for factorial model 3.1 are plotted by attaching the data frame and then calling:
source(file="http://southampton.likn.co/~cpd/anovas/datasets/PlotMeans.R")
plot_means(A:B, Y, "Cross-factors A:B", "Response Y", "CI")
The suite of commands below will plot means and 95% confidence intervals of Y at each level of factor A
XM <- tapply(as.numeric(A), A, mean)
YM <- tapply(Y, A, mean)
SE <- (tapply(Y, A, var)/tapply(Y, A, length))^0.5
CI <- SE*qt(0.975, (tapply(Y, A, length) - 1))
plot(XM, YM, type = "p", lwd = 6, xaxt = "n", xlim = c(min(XM), max(XM)), ylim = c(min(YM-CI), max(YM+CI)), las = 1, xlab = "A", ylab = "Y means and 95% C.I.") # plot means
for (i in 1:length(XM)) {lines(c(XM[i], XM[i]), c(YM[i]-CI[i], YM[i]+CI[i]))} # add 95% confidence intervals
axis(1, at = XM, labels = names(XM)) # plot x-axis labels at each level of A
Use the plot commands one at a time, because a new plot will overwrite any previous plot.
plot(X, Y, las = 1, xlab = "covariate X", ylab = "Response Y") # scatter plot of the response Y to a numeric covariate X
abline(coef(aov(Y ~ X))) # add the least-squares regression line to the scatter plot
library(nlme) ; plot(groupedData(Y ~ X|A), outer = ~1) # scatter plot of Y against covariate X with coloured lines joining consecutive points within each level of A
The suite of commands below will plot a graph with linear regressions of numeric variable Y against numeric variable A for each level of factor B, from a data frame containing equal-length vectors Y, A, and B with b levels
plot(A, Y, type = "n", las = 1, xlab = "A", ylab = "Response Y") # plot the graph axes
SA <- split(A, B) # create new vector SA that splits A by levels of B
SY <- split(Y, B) # create new vector SY that splits Y by levels of B
for (i in 1:nlevels(B)) points(SA[[i]], SY[[i]], pch = i - 1) # add symbols for data points at each level of B (squares, circles, etc)
for (i in 1:nlevels(B)) abline(lm(SY[[i]] ~ SA[[i]]), lty = i) # add regression lines for each level of B (continuous line, broken line, etc)
The suite of commands below will create a panel of plots for the linear regression of numeric variable Y against numeric variable A given factor B, from a data frame containing equal-length vectors Y, A, and B with b levels.
library(lattice)
xyplot(Y ~ A|B,
panel = function(x, y)
{
panel.xyplot(x, y) # plot Y against A at each level of B
panel.abline(coef(aov(y ~ x))) # add the linear regressions of Y against A at each level of B
})
library() # list all available packages
library(nlme) # open the package of linear and non-linear mixed effects models
library(help = nlme) # help on the 'nlme' package
detach(aovdata) # detach the named data frame
detach() # detach the most recently attached data frame
rm(aovdata) # remove the named object (e.g. a data frame or a vector)
rm(list = ls()) # remove all objects
<ctrl>L
Click here for a zip file containing all of the datasets named below. Copy-paste your own data into a .txt file with the same structure of tab-delimited columns with headers.
The commands below use data file 'Model1_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model1_1.txt", header = T)
attach(aovdata)
A <- factor(A)
model1_1i <- aov(Y ~ A) # Create an object (here called 'model1_1i') that gets an analysis of variance with the model: numeric vector Y is explained by factor A
summary(model1_1i) # report the results of the analysis
anova(glm(Y ~ A, family = gaussian(link = identity)), test = "F") # deviances equal SS
plot(A, Y, las = 1, xlab = "A", ylab = "Y")# box and whisker plot of the response Y to factor A
round(tapply(Y, A, mean),2) # report the mean of Y at each level of factor A, rounded to 2 decimal places
A <- as.numeric(A)
model1_1i <- aov(Y ~ A) # Create an object (here called 'model1_1i') that gets an analysis of variance with the model: numeric vector Y is explained by numeric vector A
summary(model1_1i) # report the results of the regression
A <- as.numeric(A)
model1_1i <- lm(Y ~ A)
summary(model1_1i)
A <- as.numeric(A)
anova(glm(Y ~ A, family = gaussian(link = identity)), test = "F")
plot(A, Y, las = 1, xlab = "A", ylab = "Y") # scatter plot of the response Y to numeric covariate A
abline(coef(model1_1i)) # add regression line to scatter plot
coef(model1_1i) # report the intercept and slope of the linear regression
The commands below use data file 'Model1_1contrasts3.txt' on the web for an example analysis.
aovcont <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model1_1contrasts3.txt", header = T)
attach(aovcont)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model1_1c <- aov(Y ~ B + C)
summary(model1_1c)
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model1_1i <- aov(Y ~ A)
summary.lm(model1_1i)
options(contrasts = c("contr.helmert", "contr.helmert"))
model1_1i <- lm(Y ~ A)
summary(model1_1i)
anova(glm(Y ~ B + C, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model1_1i <- glm(Y ~ A, family = gaussian(link = identity), contrasts = A)
summary(model1_1i)
The commands below use data file 'Model1_1contrasts5.txt' on the web for an example analysis.
aovcont <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model1_1contrasts5.txt", header = T)
attach(aovcont)
A <- factor(A) ; B <- factor(B) ; C <- factor(C) ; D <- factor(D) ; E <- factor(E)
model1_1c <- aov(Y ~ B + C + D + E)
summary(model1_1c)
contrasts(A) <- cbind(c(4,-1,-1,-1,-1), c(0,1,1,-1,-1), c(0,1,-1,0,0), c(0,0,0,1,-1))
model1_1i <- aov(Y ~ A)
summary.lm(model1_1i)
anova(glm(Y ~ B + C + D + E, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(4,-1,-1,-1,-1), c(0,1,1,-1,-1), c(0,1,-1,0,0), c(0,0,0,1,-1))
model1_1i <- glm(Y ~ A, family = gaussian(link = identity), contrasts = A)
summary(model1_1i)
The commands below use data file 'Model2_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model2_1.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
A_effect <- aov(Y ~ A + Error(A:B)) # test A against B nested in A
summary(A_effect)
anova(aov(Y ~ A), aov(Y ~ A/B), test = "F") # test A:B
# Step 1. Between-B ANOVA for A tested against B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-subjects ANOVA for B
between_B <- glm(Y ~ A, family = gaussian(link = identity))
between_S <- glm(Y ~ A/B, family = gaussian(link = identity))
anova(between_B, between_S, test = "F")
# Step 1. Convert A to numeric, and B from b factor levels nested in each level of A into a*b factor levels
A <- as.numeric(A)
b <- max(as.numeric(B)) ; B1 <- as.factor(as.numeric(B) + b*(A - 1))
# Step 2. Between-B1 ANOVA for covariate A tested against B1
A_effect <- aov(Y ~ A + Error(A:B1))
summary(A_effect)
# Step 3. Between-subjects ANOVA for B1
anova(aov(Y ~ A), aov(Y ~ A/B1), test = "F")
# Step 1. Convert A to numeric, and B from b levels nested in each level of A into a*b levels
A <- as.numeric(A)
b <- max(as.numeric(B)) ; B1 <- as.factor(as.numeric(B) + b*(A - 1))
# Step 2. Between-B1 ANOVA for covariate A tested against B1
YM <- as.vector(tapply(Y,list(A, B1), mean))
AM <- as.vector(tapply(as.numeric(A),list(A, B1), mean))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 3. Between-subjects ANOVA for B1
between_B <- glm(Y ~ A, family = gaussian(link = identity))
between_S <- glm(Y ~ A/B1, family = gaussian(link = identity))
anova(between_B, between_S, test = "F")
The commands below use data file 'Model2_2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model2_2.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
A_effect <- aov(Y ~ A + Error(A:B)) # test A
summary(A_effect)
B_effect <- aov(Y ~ A/B + Error(A:B:C)) # test A:B
summary(B_effect) # ignore the F-value for A
anova(aov(Y ~ A/B), aov(Y ~ A/B/C), test = "F") # test A:B:C
# Step 1. Between-B ANOVA for A tested against B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-C ANOVA for B tested against C
YM <- as.vector(tapply(Y,list(A, B, C), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, C), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, C), mean)))
between_B <- glm(YM ~ AM, family = gaussian(link = identity))
between_C <- glm(YM ~ AM/BM, family = gaussian(link = identity))
anova(between_B, between_C, test = "F")
# Step 3. Between-subjects ANOVA for C
between_C <- glm(Y ~ A/B, family = gaussian(link = identity))
between_S <- glm(Y ~ A/B/C, family = gaussian(link = identity))
anova(between_C, between_S, test = "F")
The commands below use data file 'Model3_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_1.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model3_1i <- aov(Y ~ A*B) # A*B requests A + B + A:B
summary(model3_1i)
anova(glm(Y ~ A*B, family = gaussian(link = identity)), test = "F")
interaction.plot(A, B, Y, las = 1, xlab = "Factor A", ylab = "Response Y", trace.label = "Factor B", xtick = TRUE)
A <- as.numeric(A) ; B <- as.factor(B)
model3_1i <- aov(Y ~ A*B)
summary(model3_1i)
plot(A, Y, type = "n", las = 1, xlab = "A", ylab = "Response Y") # plot the graph axes
SA <- split(A, B) # create new vector SA that splits A by levels of B
SY <- split(Y, B) # create new vector SY that splits Y by levels of B
for (i in 1:nlevels(B)) points(SA[[i]], SY[[i]], pch = i - 1) # add symbols for data points at each level of B (squares and circles)
for (i in 1:nlevels(B)) abline(lm(SY[[i]] ~ SA[[i]]), lty = i) # add regression lines for each level of B (continuous and broken)
library(lattice)
xyplot(Y ~ A|B,
panel = function(x, y)
{
panel.xyplot(x, y) # plot Y against A at each level of B
panel.abline(coef(aov(y ~ x))) # add the linear regressions of Y against A at each level of B
})
A <- as.numeric(A) ; B <- as.numeric(B)
model3_1i <- aov(Y ~ A*B)
summary(model3_1i)
The unreplicated version of this design is analyzed by model 4.1.
The commands below use data file 'Model3_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_1.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model3_1ii <- aov(Y ~ A + B + Error(A:B)) # unrestricted model for B
summary(model3_1ii)
anova(aov(Y ~ A), aov(Y ~ A + B), aov(Y ~ A*B), test = "F") # restricted alternative for B, and A:B under either model
# Step 1. Main effects A + B averaged across subjects and tested against A:B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B), mean)))
anova(glm(YM ~ AM + BM, family = gaussian(link = identity)), test = "F")
# Step 2. Interaction A:B
main_effects <- glm(Y ~ A + B, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B, family = gaussian(link = identity))
anova(main_effects, full_model, test = "F")
The commands below use data file 'Model3_1contrasts3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_1contrasts3.txt", header = T)
attach(aovdata)
A <- factor(A) ; C <- factor(C) ; D <- factor(D) ; B <- factor(B)
model3_1c <- aov(Y ~ C/D + B*C/D)
summary(model3_1c)
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model3_1i <- aov(Y ~ A*B)
summary.lm(model3_1i)
anova(glm(Y ~ C/D + B*C/D, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model3_1i <- glm(Y ~ A*B, family = gaussian(link = identity), contrasts = A)
summary(model3_1i)
The commands below use data file 'Model3_1contrasts3x4.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_1contrasts3x4.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C) ; D <- factor(D) ; E <- factor(E) ; F <- factor(F) ; G <- factor(G)
model3_1c <- aov(Y ~ C + D + E + F + G + C*E + C*F + C*G + D*E + D*F + D*G)
summary(model3_1c)
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
contrasts(B) <- cbind(c(1,1,-1,-1), c(1,-1,0,0), c(0,0,1,-1))
model3_1i <- aov(Y ~ A*B)
summary.lm(model3_1i)
anova(glm(Y ~ C + D + E + F + G + C*E + C*F + C*G + D*E + D*F + D*G, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
contrasts(B) <- cbind(c(1,1,-1,-1), c(1,-1,0,0), c(0,0,1,-1))
model3_1i <- glm(Y ~ A*B, family = gaussian(link = identity), contrasts = A, B)
summary(model3_1i)
The commands below use data file 'Model3_2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_2.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model3_2i <- aov(Y ~ A*B*C)
summary(model3_2i)
anova(glm(Y ~ A*B*C, family = gaussian(link = identity)), test = "F")
The unreplicated version of this design is analyzed by model 4.2.
The commands below use data file 'Model3_2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_2.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
# Step 1. Main effects A + B (restricted) + C and interaction A:C
model3_2ii <- aov(Y ~ A*B*C - A:B - C:B - A:B:C + Error(A:B + B:C + A:B:C))
summary(model3_2ii)
# Step 2. Add restricted-model analysis of interactions A:B, B:C, and A:B:C
all_but_B_interactions <- aov(Y ~ A*B*C - A:B - B:C - A:B:C)
add_AB <- aov(Y ~ A*B*C - B:C - A:B:C)
add_BC <- aov(Y ~ A*B*C - A:B:C)
full_model <- aov(Y ~ A*B*C)
anova(all_but_B_interactions, add_AB, add_BC, full_model, test = "F")
# Step 1. Main effects A + B averaged across C and tested against A:B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B), mean)))
anova(glm(YM ~ AM + BM, family = gaussian(link = identity)), test = "F")
# Step 2. Add C averaged across A, and tested against B:C
YM <- as.vector(tapply(Y,list(B, C), mean))
BM <- factor(as.vector(tapply(as.numeric(B),list(B, C), mean)))
CM <- factor(as.vector(tapply(as.numeric(C),list(B, C), mean)))
add_B <- glm(YM ~ BM, family = gaussian(link = identity))
add_C <- glm(YM ~ BM + CM, family = gaussian(link = identity))
anova(add_B, add_C, test = "F")
# Step 3. Add A:C averaged across subjects and tested against A:B:C
YM <- as.vector(tapply(Y,list(A, B, C), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, C), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, C), mean)))
CM <- factor(as.vector(tapply(as.numeric(C),list(A, B, C), mean)))
all_but_AC <- glm(YM ~ AM*BM + BM*CM, family = gaussian(link = identity))
full_model <- glm(YM ~ AM*BM*CM - AM:BM:CM, family = gaussian(link = identity))
anova(all_but_AC, full_model, test = "F")
# Step 4. Add restricted-model analysis of interactions A:B, B:C, and A:B:C
all_but_B_interactions <- glm(Y ~ A*B*C - A:B - B:C - A:B:C, family = gaussian(link = identity))
add_AB <- glm(Y ~ A*B*C - B:C - A:B:C, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*B*C - A:B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B*C, family = gaussian(link = identity))
anova(all_but_B_interactions, add_AB, add_BC, full_model, test = "F")
The unreplicated version of this design is analyzed by model 5.6.
The commands below use data file 'Model3_3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_3.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model3_3i <- aov(Y ~ A*C + Error(A:B + A:B:C)) # test A*C
summary(model3_3i)
anova(aov(Y ~ A/B + A*C), aov(Y ~ A/B*C), test = "F") # test C:A/B (and calculation of F-value for A:B by hand from residual error)
# Step 1. Between-B ANOVA for A averaged across C and tested against B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-B ANOVA by model comparisons for sequentially added terms C, A:C, averaged across subjects and tested against A:B:C with (c-1)*(b-a)*a error d.f.
YM <- as.vector(tapply(Y,list(A, B, C), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, C), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, C), mean)))
CM <- factor(as.vector(tapply(as.numeric(C),list(A, B, C), mean)))
between_B <- glm(YM ~ AM/BM, family = gaussian(link = identity))
add_C <- glm(YM ~ AM/BM + CM, family = gaussian(link = identity))
add_AC <- glm(YM ~ AM/BM + CM/AM, family = gaussian(link = identity))
anova(between_B, add_C, add_AC, test = "F")
# Step 3. Within-B ANOVA by model comparison for sequentially added term A:B:C (and calculation of F-value for A:B by hand from residual error)
within_B <- glm(Y ~ A/B + C/A, family = gaussian(link = identity))
full_model <- glm(Y ~ A/B + C/A + A:B:C, family = gaussian(link = identity))
anova(within_B, full_model, test = "F")
A <- as.factor(A) ; B <- as.factor(B) ; C <- as.numeric(C)
model3_3v <- aov(Y ~ A*C + Error(A:B + A:B:C)) # test A*C
summary(model3_3v)
anova(aov(Y ~ A/B + A*C), aov(Y ~ A/B*C), test = "F") # test C:A/B (and calculation of F-value for A:B by hand from residual error)
The commands below use data file 'Model3_4.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model3_4.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model3_4i <-aov(Y ~ A*B + Error(A:B:C)) # test A*B
summary(model3_4i)
anova(aov(Y ~ A*B), aov(Y ~ A*B*C), test = "F") # test C
# Step 1. Between-C ANOVA for A|B averaged across subjects and tested against C
YM <- as.vector(tapply(Y,list(A, B, C), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, C), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, C), mean)))
anova(glm(YM ~ AM*BM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-subjects ANOVA for C
within_C <- glm(Y ~ A*B, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B + A:B:C, family = gaussian(link = identity))
anova(within_C, full_model, test = "F")
The fully replicated version of this design is analyzed by the random-factor version of model 3.1.
The commands below use data file 'Model4_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_1.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A)
model4_1i <- aov(Y ~ S + A)
summary(model4_1i) # Model-2 output for S
library(nlme)
anova(lme(Y ~ A, random = ~ 1|S), test = "F")
anova(glm(Y ~ S + A, family = gaussian(link = identity)), test = "F")
The commands below use data file 'Model4_1contrasts.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_1contrasts.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model4_1c <- aov(Y ~ S + B/C + Error(S + S:B + S:B:C))
summary(model4_1c)
model4_1c <- aov(Y ~ S + B/C)
summary(model4_1c)
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model4_1i <- aov(Y ~ S + A)
summary.lm(model4_1i)
anova(glm(Y ~ S + B/C, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model4_1i <- glm(Y ~ S + A, family = gaussian(link = identity), contrasts = A)
summary(model4_1i)
The commands below use data file 'Model4_1BIB.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_1BIB.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A)
library(nlme)
anova(lme(Y ~ S + A, random = ~ 1|S), test = "F")
# Step 1. S adjusted for A
A_only <- glm(Y ~ A, family = gaussian(link = identity))
full_model <- glm(Y ~ A + S, family = gaussian(link = identity))
anova(A_only, full_model, test = "F")
# Step 2. A adjusted for S
S_only <- glm(Y ~ S, family = gaussian(link = identity))
anova(S_only, full_model, test = "F")
The commands below use data file 'Model4_1LS.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_1LS.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model4_1LS <- aov(Y ~ A + B + C)
summary(model4_1LS)
anova(glm(Y ~ A + B + C, family = gaussian(link = identity)), test = "F")
The commands below use data file 'Model4_1YS.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_1YS.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B) ; C <- factor(C)
anova(aov(Y ~ B + C), aov(Y ~ A + B + C)) # A, adjusted for B and C
anova(aov(Y ~ A + C), aov(Y ~ A + B + C)) # B, adjusted for A and C
anova(aov(Y ~ A + B), aov(Y ~ A + B + C)) # C, adjusted for A and B
# Step 1. A adjusted for B and C
BC_only <- glm(Y ~ B + C, family = gaussian(link = identity))
full_model <- glm(Y ~ A + B + C, family = gaussian(link = identity))
anova(BC_only, full_model, test = "F")
# Step 2. B adjusted for A and C
AC_only <- glm(Y ~ A + C, family = gaussian(link = identity))
anova(AC_only, full_model, test = "F")
# Step 3. C adjusted for A and B
AB_only <- glm(Y ~ A + B, family = gaussian(link = identity))
anova(AB_only, full_model, test = "F")
The fully replicated version of this design is analyzed by the random-factor version of model 3.2.
The commands below use data file 'Model4_2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_2.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model4_2i <- aov(Y ~ S + A*B + Error(S + S:A + S:B))
summary(model4_2i)
model4_2i <- aov(Y ~ S + A*B)
summary(model4_2i)
library(nlme)
anova(lme(Y ~ S + A*B, random = ~ 1|S), test = "F")
anova(glm(Y ~ S + A*B, family = gaussian(link = identity)), test = "F")
The commands below use data file 'Model4_2contrasts.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_2contrasts.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; C <- factor(C) ; D <- factor(D); B <- factor(B)
model4_2c <- aov(Y ~ S + C/D + B*C/D - S:B:C:D + Error(S + S:B + S:C + S:C:D + S:B:C + S:B:C/D))
summary(model4_2c)
model4_2i <- aov(Y ~ S + C/D + B*C/D)
summary(model4_2i)
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model4_2i <- aov(Y ~ S + A*B)
summary.lm(model4_2i)
anova(glm(Y ~ S + C/D + B*C/D, family = gaussian(link = identity)), test = "F")
contrasts(A) <- cbind(c(2,-1,-1), c(0,1,-1))
model4_2i <- glm(Y ~ S + A*B, family = gaussian(link = identity), contrasts = A)
summary(model4_2i)
The commands below use data file 'Model4_3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model4_3.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model4_3i <- aov(Y ~ S + A*B*C + Error(S + S:A + S:B + S:C + S:A:B + S:A:C + S:B:C))
summary(model4_3i)
model4_3i <- aov(Y ~ S + A*B*C)
summary(model4_3i)
library(nlme)
anova(lme(Y ~ S + A*B*C, random = ~ 1|S), test = "F")
anova(glm(Y ~ S + A*B*C, family = gaussian(link = identity)), test = "F")
The commands below use data file 'Model5_1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_1.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model5_1i <- aov(Y ~ A*B + Error(S/A))
summary(model5_1i)
# Step 1. Between-P ANOVA for A averaged across B and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, S), mean)))
anova(glm(YM ~ SM + AM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-P ANOVA by model comparisons for sequentially added terms B, A:B, both tested against pooled interactions with S with a*(b-1)*(s-1) error d.f.
between_P <- glm(Y ~ A*S, family = gaussian(link = identity))
add_B <- glm(Y ~ A*S + B, family = gaussian(link = identity))
add_AB <- glm(Y ~ A*S + B/A, family = gaussian(link = identity))
anova(between_P, add_B, add_AB, test = "F")
The commands below use data file 'Model5_2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_2.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_2i <- aov(Y ~ A*B*C + Error(S/A:B))
summary(model5_2i)
# Step 1. Between-P ANOVA for A averaged across C and tested against S
YM <- as.vector(tapply(Y,list(A, B, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, S), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, B, S), mean)))
anova(glm(YM ~ SM + AM*BM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-P ANOVA by model comparisons for sequentially added terms C, A:C, B:C, A:B:C, all tested against pooled interactions with S with (a*b-1)*(c-1)*(s-1) error d.f.
between_P <- glm(Y ~ A*B*S, family = gaussian(link = identity))
add_C <- glm(Y ~ A*B*S + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A*B*S + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*B*S + C/A + B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B*S + A*B*C, family = gaussian(link = identity))
anova(between_P, add_C, add_AC, add_BC, full_model, test = "F")
The commands below use data file 'Model5_3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_3.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_3i <- aov(Y ~ A*B*C + Error(S/A))
summary(model5_3i)
# Step 1. Between-P ANOVA for A averaged across B & C and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, S), mean)))
anova(glm(YM ~ SM + AM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-P ANOVA by model comparisons for sequentially added terms B, A:B, C, A:C, B:C, A:B:C, all tested against pooled interactions with S with a*(b*c-1)*(s-1) error d.f.
between_P <- glm(Y ~ A*S, family = gaussian(link = identity))
add_B <- glm(Y ~ A*S + B, family = gaussian(link = identity))
add_AB <- glm(Y ~ A*S + B/A, family = gaussian(link = identity))
add_C <- glm(Y ~ A*S + B/A + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A*S + B/A + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*S + A*B*C - A:B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*S + A*B*C, family = gaussian(link = identity))
anova(between_P, add_B, add_AB, add_C, add_AC, add_BC, full_model, test = "F")
The commands below use data file 'Model5_4.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_4.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_4i <- aov(Y ~ A*B*C + Error(S/A/B))
summary(model5_4i)
# Step 1. Between-P ANOVA for A averaged across B & C and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, S), mean)))
anova(glm(YM ~ SM + AM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-Q ANOVA for B and A:B averaged across C and tested against pooled interactions with S, both with a*(b-1)*(s-1) d.f.
YM <- as.vector(tapply(Y,list(A, B, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, S), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, B, S), mean)))
between_P <- glm(YM ~ AM*SM, family = gaussian(link = identity))
add_B <- glm(YM ~ AM*SM + BM, family = gaussian(link = identity))
add_AB <- glm(YM ~ AM*SM + BM/AM, family = gaussian(link = identity))
anova(between_P, add_B, add_AB, test = "F")
# Step 3. Within-Q ANOVA by model comparisons for sequentially added terms C, A:C, B:C, A:B:C, all tested against pooled interactions with S with a*b*(c-1)*(s-1) error d.f.
between_Q <- glm(Y ~ A*B*S, family = gaussian(link = identity))
add_C <- glm(Y ~ A*B*S + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A*B*S + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*B*S + C/A + B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B*S + A*B*C, family = gaussian(link = identity))
anova(between_Q, add_C, add_AC, add_BC, full_model, test = "F")
The commands below use data file 'Model5_5.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_5.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_5i <- aov(Y ~ A*B*C + Error(S:A + S:A:B))
summary(model5_5i)
# Step 1. Between-P ANOVA for A averaged across B & C and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, S), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-P ANOVA for B and A:B averaged across C and tested against the interaction with S, both with a*(b-1)*(s-1) d.f.
YM <- as.vector(tapply(Y,list(A, B, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, S), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, S), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(A, B, S), mean)))
between_S <- glm(YM ~ AM/SM, family = gaussian(link = identity))
add_B <- glm(YM ~ AM/SM + BM, family = gaussian(link = identity))
add_AB <- glm(YM ~ AM/SM + BM/AM, family = gaussian(link = identity))
anova(between_S, add_B, add_AB, test = "F")
# Step 3. Within-P ANOVA by model comparisons for sequentially added terms C, A:C, B:C, A:B:C, all tested against pooled interactions with S with a*b*(c-1)*(s-1) error d.f.
between_P <- glm(Y ~ A*B*S, family = gaussian(link = identity))
add_C <- glm(Y ~ A*B/S + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A*B/S + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*B/S + C/A + B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B/S + A*B*C, family = gaussian(link = identity))
anova(between_P, add_C, add_AC, add_BC, full_model, test = "F")
The fully replicated version of this design is analyzed by model 3.3.
The commands below use data file 'Model5_6.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_6.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model5_6i <- aov(Y ~ A*B + Error(S:A))
summary(model5_6i)
# Step 1. Between-S ANOVA for A averaged across B and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-S ANOVA by model comparisons for sequentially added terms B, A:B, both tested against A:B:S with a*(b-1)*(s-1) error d.f.
between_S <- glm(Y ~ A/S, family = gaussian(link = identity))
add_B <- glm(Y ~ A/S + B, family = gaussian(link = identity))
add_AB <- glm(Y ~ A/S + B/A, family = gaussian(link = identity))
anova(between_S, add_B, add_AB, test = "F")
See model 6.5 for a Model-1 analysis of this design. The commands below use data file 'Model5_7.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_7.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_7i <- aov(Y ~ A*B*C + Error(S:A))
summary(model5_7i)
# Step 1. Between-S ANOVA for A averaged across B & C, and tested against S
YM <- as.vector(tapply(Y,list(A, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, S), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-S ANOVA by model comparisons for sequentially added terms B, C, A:B, A:C, A:B:C, all tested against pooled interactions with S, with a*(b*C-1)*(s-1) error d.f.
between_S <- glm(Y ~ A/S, family = gaussian(link = identity))
add_B <- glm(Y ~ A/S + B, family = gaussian(link = identity))
add_C <- glm(Y ~ A/S + B + C, family = gaussian(link = identity))
add_AB <- glm(Y ~ A/S + B/A + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A/S + B/A + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A/S + A*B*C - A:B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A/S + A*B*C, family = gaussian(link = identity))
anova(between_S, add_B, add_C, add_AB, add_AC, add_BC, full_model, test = "F")
The commands below use data file 'Model5_7contrasts.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_7contrasts.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C) ; D <- factor(D) ; E <- factor(E)
model5_7iv_c <- aov(Y ~ A*E*C + A*E*C:D + Error(A:S + A:E + A:C + A:C:D + A:E:C + A:E:C:D))
summary(model5_7iv_c) # then calculate F-values by hand
The commands below use data file 'Model5_8.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_8.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
# To avoid calculating F-values by hand, use the suite of commands in the next section for analysis by glm
model5_8i <- aov(Y ~ A/B/S*C + Error(A:B + S:A:B))
summary(model5_8i) # then calculate F-values by hand
# Step 1. Between-B ANOVA for A averaged across S and C and tested against B
YM <- as.vector(tapply(Y,list(A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B), mean)))
anova(glm(YM ~ AM, family = gaussian(link = identity)), test = "F")
# Step 2. Between-subjects ANOVA for B averaged across C and tested against S
YM <- as.vector(tapply(Y,list(A, B, S), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, S), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, S), mean)))
between_B <- glm(YM ~ AM, family = gaussian(link = identity))
between_S <- glm(YM ~ AM/BM, family = gaussian(link = identity))
anova(between_B, between_S, test = "F")
# Step 3. Within-subjects ANOVA by model comparisons for sequentially added terms C, A:C, averaged across S and both tested against A:B:C with (c-1)*(b-a)*a error d.f.
YM <- as.vector(tapply(Y,list(A, B, C), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(A, B, C), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(A, B, C), mean)))
CM <- factor(as.vector(tapply(as.numeric(C),list(A, B, C), mean)))
between_S <- glm(YM ~ AM/BM, family = gaussian(link = identity))
add_C <- glm(YM ~ AM/BM + CM, family = gaussian(link = identity))
add_AC <- glm(YM ~ AM/BM + CM/AM, family = gaussian(link = identity))
anova(between_S, add_C, add_AC, test = "F")
# Step 4. Within-subjects ANOVA by model comparison for sequentially added term A:B:C
between_S <- glm(Y ~ A/B/S + C/A, family = gaussian(link = identity))
full_model <- glm(Y ~ A/B/S + C/A + A:B:C, family = gaussian(link = identity))
anova(between_S, full_model, test = "F")
The commands below use data file 'Model5_9.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_9.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model5_9i <- aov(Y ~ A*B*C + Error(S:A:B))
summary(model5_9i)
# Step 1. Between-subjects ANOVA for A*B averaged across C, all tested against A:B:S with error d.f. = (s-1)*a*b
YM <- as.vector(tapply(Y,list(S, A, B), mean))
AM <- factor(as.vector(tapply(as.numeric(A),list(S, A, B), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(S, A, B), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(S, A, B), mean)))
anova(glm(YM ~ AM*BM/SM - AM:BM:SM, family = gaussian(link = identity)), test = "F")
# Step 2. Within-subjects ANOVA by model comparisons for sequentially added terms C, A:C, B:C, A:B:C all tested against A:B:C:S with error d.f. = a*b*(c-1)*(s-1)
between_S <- glm(Y ~ A*B + A:B/S, family = gaussian(link = identity))
add_C <- glm(Y ~ A*B + A:B/S + C, family = gaussian(link = identity))
add_AC <- glm(Y ~ A*B + A:B/S + C/A, family = gaussian(link = identity))
add_BC <- glm(Y ~ A*B + A:B/S + C/A + B:C, family = gaussian(link = identity))
full_model <- glm(Y ~ A*B*C + A:B/S, family = gaussian(link = identity))
anova(between_S, add_C, add_AC, add_BC, full_model, test = "F")
The commands below use data file 'Model5_9binomial.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model5_9binomial.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
Apply an arcsine-root transformation to the response, and then use commands as in the previous subsection for aov or glm.
Y <- asin(sqrt(Su/(Su + Fa)))
# Step 1. Between-subjects ANOVA for A*B averaged across C, all tested against A:B:S with error d.f. = (s-1)*a*b
SuM <- as.integer(tapply(Su,list(S, A, B), mean))
FaM <- as.integer(tapply(Fa,list(S, A, B), mean))
YM <- cbind(SuM,FaM)
AM <- factor(as.vector(tapply(as.numeric(A),list(S, A, B), mean)))
BM <- factor(as.vector(tapply(as.numeric(B),list(S, A, B), mean)))
SM <- factor(as.vector(tapply(as.numeric(S),list(S, A, B), mean)))
anova(glm(YM ~ AM*BM/SM - AM:BM:SM , family = binomial), test = "Chisq")
# Step 2. Within-subjects ANOVA by model comparisons for sequentially added terms C, A:C, B:C, A:B:C all tested against A:B:C:S with error d.f. = a*b*(c-1)*(s-1)
Y <- cbind(Su,Fa)
between_S <- glm(Y ~ A*B + A:B/S, family = binomial)
add_C <- glm(Y ~ A*B + A:B/S + C, family = binomial)
add_AC <- glm(Y ~ A*B + A:B/S + C/A, family = binomial)
add_BC <- glm(Y ~ A*B + A:B/S + C/A + B:C, family = binomial)
full_model <- glm(Y ~ A*B*C + A:B/S, family = binomial)
anova(between_S, add_C, add_AC, add_BC, full_model, test = "Chisq")
# If the residual deviance >> residual d.f., then significances from chi-square tests will be inflated by "overdispersion". One way to compensate for overdispersion is to treat the deviance for each term as a Sum of Squares, and calculate the Mean Squares and F by hand (e.g., M. J. Crawley 2002 Statistical Computing, Wiley, pp 526-532). This needs doing manually for each term, obtaining its Error Mean Square by dividing the residual deviance of the last-entered term in a step by its residual d.f. (requesting: test = "F" will produce an F-value calculated by an alternative method that does not compensate for overdispersion, and may differ little in significance from the chi-square).
Analysis of this model is exactly as for model 4.1.
The fully replicated version of this design is analyzed by the random-factor version of model 3.1.
Analysis of this model is exactly as for model 4.2.
Analysis of this model is exactly as for model 5.6.
The fully replicated version of this design is analyzed by model 3.3.
The commands below use data file 'Model6_4.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model6_4.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model6_4i <- aov(aov(Y ~ A/S + B/C + A:B + A/S:B + A:B/C + A/S:B/C))
summary(model6_4i) # then calculate F-values by hand
See model 5.7 for a Model-2 analysis of this design.
The commands below use data file 'Model6_5.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Model6_5.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B) ; C <- factor(C)
model6_5i <- aov(Y ~ A*B*C + Error(S:A + B:S:A + C:S:A))
summary(model6_5i)
Analysis of this model is exactly as for model 5.8.
Analysis of this model is exactly as for model 5.9.
The commands below use data file 'Fig1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig1.txt", header = T)
attach(aovdata)
A <- factor(A)
model1_1i <- aov(Y ~ A)
summary(model1_1i)
The commands below use data file 'Fig2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig2.txt", header = T)
attach(aovdata)
model1_1reg <- aov(Y ~ X)
summary(model1_1reg)
model1_1reg <- lm(Y ~ X)
summary(model1_1reg)
The commands below use data file 'Fig3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig3.txt", header = T)
attach(aovdata)
model1_1reg <- aov(Y ~ X)
summary(model1_1reg)
The commands below use data file 'Fig4.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig4.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model2_1i <- aov(Y ~ A + Error(A:B)) # test A
summary(model2_1i)
anova(aov(Y ~ A), aov(Y ~ A/B), test = "F") # test A:B
The commands below use data file 'Fig5.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig5.txt", header = T)
attach(aovdata)
anova(glm(Y ~ A*B), test = "F")
The commands below use data file 'Fig6a.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig6a.txt", header = T)
attach(aovdata)
model3_1 <- aov(Y ~ A*B)
summary(model3_1)
The commands below use data file 'Fig6b.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig6b.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model4_2i <- aov(Y ~ S + A*B + Error(S + S:A + S:B))
summary(model4_2i)
model4_2i <- aov(Y ~ S + A*B)
summary(model4_2i)
The commands below use data file 'Fig7a.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig7a.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model5_1i <- aov(Y ~ A*B + Error(S/A))
summary(model5_1i)
The commands below use data file 'Fig7b.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig7b.txt", header = T)
attach(aovdata)
S <- factor(S) ; A <- factor(A) ; B <- factor(B)
model5_6i <- aov(Y ~ A*B + Error(S:A))
summary(model5_6i)
The commands below use data file 'Fig8.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig8.txt", header = T)
attach(aovdata)
A <- as.factor(A) ; X <- as.numeric(X)
model3_1 <- aov(Y ~ X*A)
summary(model3_1)
The commands below use data file 'Fig9.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig9.txt", header = T)
attach(aovdata)
model1_1reg <- aov(Y ~ X)
summary(model1_1reg)
The commands below use data file 'Fig10.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig10.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model3_1 <- aov(Y1 ~ A*B)
summary(model3_1)
The commands below use data file 'Fig11a.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig11a.txt", header = T)
attach(aovdata)
t <- as.numeric(t) ; tplus1 <- as.numeric(tplus1)
anova(aov(r ~ tplus1), aov(r ~ t + tplus1)) # t, adjusted for tplus1
anova(aov(r ~ t), aov(r ~ t + tplus1)) # tplus1, adjusted for t
The commands below use data file 'Fig11b.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Fig11b.txt", header = T)
attach(aovdata)
t <- as.numeric(t) ; tplus1 <- as.numeric(tplus1)
anova(aov(r ~ tplus1), aov(r ~ t + tplus1)) # t, adjusted for tplus1
anova(aov(r ~ t), aov(r ~ t + tplus1)) # tplus1, adjusted for t
The commands below use data file 'Worked1.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Worked1.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model2_1i <- aov(Y ~ A + Error(A:B)) # test A
summary(model2_1i)
anova(aov(Y ~ A), aov(Y ~ A/B), test = "F") # test A:B
The commands below use data file 'Worked2.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Worked2.txt", header = T)
attach(aovdata)
A <- factor(A) ; B <- factor(B)
model3_1 <- aov(Y ~ A*B)
summary(model3_1)
The commands below use data file 'Worked3.txt' on the web for an example analysis.
aovdata <- read.table("http://southampton.likn.co/~cpd/anovas/datasets/R/Worked3.txt", header = T)
attach(aovdata)
Shore <- factor(Shore) ; Recruit <- factor(Recruit) ; Trtmnt <- factor(Trtmnt)
model3_3i <- aov(Density ~ Recruit*Trtmnt + Error(Recruit:Shore + Recruit:Shore:Trtmnt))
summary(model3_3i)
anova(aov(Density ~ Recruit/Shore + Recruit*Trtmnt), aov(Density ~ Recruit/Shore*Trtmnt), test = "F") # test Trtmnt:Recruit/Shore (and calculation of F-value for Recruit:Shore by hand from residual error)
Shore <- as.factor(Shore) ; Recruit <- as.factor(Recruit) ; Trtmnt <- as.numeric(Trtmnt)
model3_3i <- aov(Density ~ Recruit*Trtmnt + Error(Recruit:Shore + Recruit:Shore:Trtmnt))
summary(model3_3i)
anova(aov(Density ~ Recruit/Shore + Recruit*Trtmnt), aov(Density ~ Recruit/Shore*Trtmnt), test = "F") # test Trtmnt:Recruit/Shore (and calculation of F-value for Recruit:Shore by hand from residual error)
1 - pf(q = 5.98, df1 = 1, df2 = 6)
[1] 0.05010252
2*(1 - pt(q = sqrt(5.98), df = 6)) # two-tailed p at t = sqrt(F) equals p at F
[1] 0.05010252
qf(p = 1-0.05, df1 = 1, df2 = 6)
[1] 5.987378
qt(p = 1-0.05/2, df = 6)^2 # t^2 at two-tailed alpha equals F at alpha
[1] 5.987378
Power calculations for a term in any ANOVA design require values of all but any one of these parameters:
- Critical probability of Type-I error, α
- Test degrees of freedom, df1
- Error degrees of freedom, df2
- Effective sample size, n
- Effect size variance in the sampled population, θ2
- Error variance in the sampled population, σ2
- Power of the design to detect an effect at α, given θ2/σ2
If the term is fixed, its non-centrality parameter, ncp = df1*n*θ2/σ2
If the term is random, its expected value of F, evf = 1 + n*θ2/σ2
# Declare quantities
alpha <- 0.05 ; df1 <- 2 ; df2 <- 9 ; n <- 24 ; theta2 <- 0.25 ; sigma2 <- 1
# Calculate power of fixed term
ncp <- df1*n*theta2/sigma2
power_fixed <- 1 - pf(q = qf(p = 1-alpha, df1, df2), df1, df2, ncp)
power_fixed
[1] 0.743147
# Declare quantities
alpha <- 0.05 ; df1 <- 2 ; df2 <- 9 ; n <- 24 ; theta2 <- 0.25 ; sigma2 <- 1
# Calculate power of random term
evf <- 1 + n*theta2/sigma2
power_random <- 1 - pf(q = qf(p = 1-alpha, df1, df2)/evf, df1, df2)
power_random
[1] 0.5653277
power.anova.test(groups = 2, n = NULL, within.var = 1, between.var = 1.6, power = 0.9, sig.level = 0.05)
# yields n = 7.668471
power.anova.test(groups = 2, n = 6, within.var = 1, between.var = NULL, power = 0.9, sig.level = 0.05)
# yields θ^2 = 2.164521
power.t.test(n = NULL, sd = 45, delta = 26, power = 0.9, sig.level = 0.05, type = "two.sample", alternative = "two.sided")
# yields n = 63.92708
power.t.test(n = 10, sd = 45, delta = NULL, power = 0.9, sig.level = 0.05, type = "one.sample", alternative = "one.sided")
# yields delta = 45.21825
Doncaster, C. P. & Davey, A. J. H. (2007) Analysis of Variance and Covariance: How to Choose and Construct Models for the Life Sciences. Cambridge: Cambridge University Press. https://doi.org/10.1017/CBO9780511611377. Example datasets: http://southampton.likn.co/~cpd/anovas/datasets/
Code formatted by Pretty R at inside-R.org ... RIP. Pages maintained by C. Patrick Doncaster, 2 July 2024