 |
| See Video ⮞ ☝ |
AGRON Info-Tech
June 20, 2020
This tutorial introduces a two-factor analysis with repeated observations. A repeated-measures ANOVA is appropriate when the same experimental units are observed under each level of the within-subject factor or factors. If one factor identifies groups to which different subjects belong, the design is mixed rather than fully repeated.
A two-factor repeated-measures analysis is also described as a two-way within-subject ANOVA when both factors are repeated on the same experimental units.
- Two factor repeated measures anova
- Two factor or two way anova with repeated measures
- Within within subjects anova
Common applications include measurements taken from the same participants across time points and conditions. The unit that is repeatedly measured must be identified, and the dependence among its observations must be modeled.
- A dependent variable over two or more time points
- When subjects have undergone two or more conditions
A mixed ANOVA has at least one between-subject factor and one within-subject factor. In the example, if male and female participants belong to only one group each and are measured at several times, group is between-subject and time is within-subject. The analysis must reflect this structure.
Description of data
The data include a subject identifier, group, time, and response. Confirm whether each subject occurs in every group-by-time combination. If each subject belongs to only one group, the code’s Error(subID/(group * time)) term treats group as repeated and does not match that design; use an appropriate mixed model or mixed-design ANOVA instead.
Clear R environment
Clear objects only when needed and preserve any work that should remain in the session. Closing graphics devices and clearing the console are optional housekeeping tasks.
rm(list = ls(all = TRUE))
graphics.off()
shell("cls")Importing data
Import the CSV file with the correct path and header setting. Check the first rows, dimensions, missing values, and whether subject IDs map to the intended repeated observations.
Convert group and time to factors when they represent categories. Keep the subject identifier consistent within each participant and verify that each subject has the expected number of observations.
data <- read.csv(file = "data.csv",
header = TRUE)
head(data)
data$group = as.factor(data$group)
data$time = as.factor(data$time)
attach(data)# subID group time resp.var
# 1 1 male day1 26
# 2 2 male day1 27
# 3 3 male day1 24
# 4 4 female day1 18
# 5 5 female day1 17
# 6 6 female day1 16Fitting ANOVA model
The displayed aov() syntax specifies within-subject error strata for both group and time. Use it only when both factors are repeated within each subject. When group is between-subject and time is repeated, use a mixed-design model with subjects nested in group and the repeated-time covariance handled appropriately.
model.aov <- aov(resp.var ~
group * time +
Error(subID/(group*time)))
#OR
model.aov <- aov(resp.var ~
group * time +
Error(subID/(group + time + group:time)))
summary(model.aov)#
# Error: subID
# Df Sum Sq Mean Sq
# group 1 408 408
#
# Error: subID:group
# Df Sum Sq Mean Sq
# group 1 67.25 67.25
#
# Error: subID:time
# Df Sum Sq Mean Sq
# time 3 356.7 118.9
#
# Error: subID:group:time
# Df Sum Sq Mean Sq
# time 3 49.94 16.65
#
# Error: Within
# Df Sum Sq Mean Sq F value Pr(>F)
# group 1 13.86 13.863 6.105 0.03866 *
# time 3 52.79 17.596 7.749 0.00943 **
# group:time 3 2.60 0.868 0.382 0.76887
# Residuals 8 18.17 2.271
# ---
# Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1The Error() term partitions variation according to the repeated-measures structure. Inspect the resulting strata and identify the correct denominator mean square and degrees of freedom for each effect; a single residual term is not generally appropriate for all tests.
- Subject ID
- Subject ID and group interaction
- Subject ID and time interaction
- Three way interaction of subject ID, group and time
Interpret main effects and interactions from the strata that correspond to the design. A non-significant interaction does not prove identical effects; report estimates and uncertainty, and consider the scientific importance of the interaction.
Mean comparison test
Follow-up comparisons should be tied to the fitted model and the effects supported by the analysis. For interactions, examine simple effects or prespecified contrasts rather than interpreting marginal means in isolation.
SNK test
The Student–Newman–Keuls procedure is a stepwise range test and is less protective against family-wise false positives than procedures such as Tukey’s method. Select a multiplicity adjustment that matches the inferential goal and report it clearly.
Edf <- df.residual(model.aov$Within)
Edf# [1] 8EMS <- deviance(model.aov$Within)/Edf
EMS# [1] 2.270833Before using SNK.test(), identify the error degrees of freedom and mean square for the particular effect being compared. The code shown extracts one residual stratum for multiple factors; verify that this is the correct denominator for each contrast before relying on the results.
- Y argument that specify the response variable
- TRT argument that represents the treatment factor variable applied to each experimental unit as group and time in this example.
- DF error and MS error represents degree of freedom and mean square for error term as determined earlier
- Alpha shows the significance level. Default value is 5 percent.
- Group argument specify the formation of treatment groups. True value for group will result in treatment groups with lettering. False value for group will provide information regarding comparisons between treatments.
library(agricolae)
SNK.test1 <- SNK.test(y = resp.var,
trt = group,
DFerror = Edf,
MSerror = EMS,
alpha = 0.05,
group = TRUE)
print(SNK.test1)# $statistics
# MSerror Df Mean CV
# 2.270833 8 14.33333 10.51345
#
# $parameters
# test name.t ntr alpha
# SNK group 2 0.05
#
# $snk
# Table CriticalRange
# 2 3.261182 1.418656
#
# $means
# resp.var std r Min Max Q25 Q50 Q75
# female 9.833333 4.687184 12 4 18 6.75 9.0 11.50
# male 18.833333 4.687184 12 13 27 15.75 17.5 21.75
#
# $comparison
# NULL
#
# $groups
# resp.var groups
# male 18.833333 a
# female 9.833333 b
#
# attr(,"class")
# [1] "group"For group comparisons, use the group-specific model estimates and the correct between-subject error term when group is between-subject. Repeated-measures software output should not be replaced by an unrelated residual mean square.
SNK.test2 <- SNK.test(y = resp.var,
trt = time,
DFerror = Edf,
MSerror = EMS,
alpha = 0.05,
group = TRUE)
print(SNK.test2)# $statistics
# MSerror Df Mean CV
# 2.270833 8 14.33333 10.51345
#
# $parameters
# test name.t ntr alpha
# SNK time 4 0.05
#
# $snk
# Table CriticalRange
# 2 3.261182 2.006283
# 3 4.041036 2.486050
# 4 4.528810 2.786128
#
# $means
# resp.var std r Min Max Q25 Q50 Q75
# day1 21.333333 4.885352 6 16 27 17.25 21.0 25.50
# day2 13.833333 6.145459 6 7 21 9.00 13.5 18.75
# day3 12.500000 3.937004 6 8 17 9.25 12.5 15.75
# day4 9.666667 5.240865 6 4 16 5.25 9.5 13.75
#
# $comparison
# NULL
#
# $groups
# resp.var groups
# day1 21.333333 a
# day2 13.833333 b
# day3 12.500000 b
# day4 9.666667 c
#
# attr(,"class")
# [1] "group"Visualizing results
Barplot for time variable
For time comparisons, account for repeated observations within subjects and use contrasts from the fitted repeated-measures model. For a significant interaction, compare times within groups or other prespecified simple effects with multiplicity control.
library(ggplot2)
SE1 = SNK.test1$means$std/sqrt(SNK.test1$means$r)
ggplot(data=SNK.test1$means,
aes(x= rownames(SNK.test1$means),
y= resp.var,
fill=rownames(SNK.test1$means))) +
geom_bar(stat="identity") +
geom_errorbar(aes(ymax = resp.var + std/sqrt(r),
ymin = resp.var - std/sqrt(r)),
position = position_dodge(width=0.9),
width = 0.25) +
labs(title = "",
x = "group",
y = "Response variable",
fill = "group") +
geom_text(data = SNK.test1$groups,
aes(x = rownames(SNK.test1$groups),
y = resp.var + SE1,
label = as.matrix(groups)),
position = position_dodge(width = 0.9),
vjust = -(0.5))The plots show means with standard errors and letters from the selected comparison procedure. Ensure the summaries and letters are based on the same model estimands and factor levels.
Barplot for group variable
SE2 = SNK.test2$means$std/sqrt(SNK.test2$means$r)
ggplot(data=SNK.test2$means,
aes(x= rownames(SNK.test2$means),
y= resp.var,
fill=rownames(SNK.test2$means))) +
geom_bar(stat="identity") +
geom_errorbar(aes(ymax = resp.var + std/sqrt(r),
ymin = resp.var - std/sqrt(r)),
position = position_dodge(width=0.9),
width = 0.25) +
labs(title = "",
x = "time",
y = "Response variable",
fill = "time") +
geom_text(data = SNK.test2$groups,
aes(x = rownames(SNK.test2$groups),
y = resp.var + SE2,
label = as.matrix(groups)),
position = position_dodge(width = 0.9),
vjust = -(0.5))For each factor, check the standard-error calculation and group sizes. In repeated-measures data, a conventional SE across raw observations may not represent uncertainty for a within-subject contrast; model-based estimates and confidence intervals can be more informative.
Please comment below if you have any questions.
Download data file — Click_here
Download Rscript — Click_here
Download R program —
Click_hereDownload R studio —
Click_here
For some reason, the "error within" section of the summary(model.aov) is not showing up in my R console. Is there a common reason why?
ReplyDelete