The methodology comes from https://gribblelab.wordpress.com/2009/03/09/repeated-measures-anova-using-r/
Load the various packages.
# packages for analysis of repeated measures
library(nlme)
library(car)
# document definition and construction
library(knitr)
opts_knit$set(eval.after = "fig.cap")
library(rmarkdown)
#reading and writing files of various types (in this case, csv and text files)
library(readr)
# data manipulation
library(tidyr)
library(dplyr)
library(tibble)
# pipelining commands
library(magrittr)
# plots
library(ggplot2)Read the data from the CSV file (which was copied and pasted from the source Excel spreadsheet). Sort the data by subject so things will line up later on. Remove the unneeded columns (not critical to the process).
obs <- read_csv("mep_observations.csv")Parsed with column specification:
cols(
Subject = col_integer(),
TestNumber = col_integer(),
MEPHR = col_double(),
Weight = col_double(),
MaxRelVO2 = col_double(),
FatBIA = col_double(),
FatSF = col_double(),
Avgkcals = col_integer(),
AvgFat = col_integer(),
AvgCHO = col_integer(),
AvgPro = col_integer()
)# these models work with factors
obs$Subject <- as.factor(obs$Subject)
obs$TestNumber <- as.factor(obs$TestNumber)
obs <- as_tibble(obs)
obsMEPHR.Subject.TestNumber.aov.model <- aov(MEPHR ~ TestNumber +
Error(Subject/TestNumber),
data=obs)
summary(aov.model)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 888.2 222
Error: Subject:TestNumber
Df Sum Sq Mean Sq F value Pr(>F)
TestNumber 3 505.9 168.6 2.148 0.147
Residuals 12 942.0 78.5 This method fits a linear mixed-effects model.
lme.model <- lme(MEPHR ~ TestNumber,
random = ~1|Subject/TestNumber,
data=obs)
summary(lme.model)Linear mixed-effects model fit by REML
Data: obs
AIC BIC logLik
139.8128 145.2209 -62.9064
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 5.990346
Formula: ~1 | TestNumber %in% Subject
(Intercept) Residual
StdDev: 8.220278 3.306014
Fixed effects: MEPHR ~ TestNumber
Value Std.Error DF t-value p-value
(Intercept) 149.000 4.783031 12 31.151792 0.0000
TestNumber1 14.200 5.603666 12 2.534055 0.0262
TestNumber2 6.666 5.603666 12 1.189578 0.2572
TestNumber3 6.400 5.603666 12 1.142109 0.2757
Correlation:
(Intr) TstNm1 TstNm2
TestNumber1 -0.586
TestNumber2 -0.586 0.500
TestNumber3 -0.586 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.503809254 -0.240710015 -0.002141907 0.152251725 0.650097511
Number of Observations: 20
Number of Groups:
Subject TestNumber %in% Subject
5 20 Display F- and p-values for the factor effect.
anova(lme.model) numDF denDF F-value p-value
(Intercept) 1 12 2186.8866 <.0001
TestNumber 3 12 2.1481 0.1474obs.matrix <- with(obs, # data from obs data frame
cbind(
MEPHR[TestNumber=="0"], # subject observations baseline
MEPHR[TestNumber=="1"], # subject observations test 1
MEPHR[TestNumber=="2"], # etc
MEPHR[TestNumber=="3"]
) # cbind column bind, rows are flipped to columns
)
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$TestNumber)
# display the matrix
obs.matrix 0 1 2 3
112121 144 180 160.33 157
258996 161 151 150.00 140
456121 142 165 158.00 164
522210 163 168 165.00 163
563751 135 152 145.00 153mlm.model <- lm(obs.matrix ~ 1)
mlm.model
Call:
lm(formula = obs.matrix ~ 1)
Coefficients:
0 1 2 3
(Intercept) 149.0 163.2 155.7 155.4The coefficients are equal to the means across subjects at each test.
rfactor <- as.factor(levels(obs$TestNumber))
rfactor[1] 0 1 2 3
Levels: 0 1 2 3mlm.model.aov.mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
idata = data.frame(rfactor),
idesign = ~rfactor,
type = "III")
summary(mlm.model.aov,
multivariate = FALSE) # don't show multivariate tests
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 485576 1 888.16 4 2186.8866 1.251e-06 ***
rfactor 506 3 942.03 12 2.1481 0.1474
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
rfactor 0.013187 0.049816
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
rfactor 0.49191 0.1985
HF eps Pr(>F[HF])
rfactor 0.7102564 0.1746291Weight.Subject.TestNumber.aov.model <- aov(Weight ~ TestNumber + Error(Subject/TestNumber), data=obs)
summary(aov.model)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 4609 1152
Error: Subject:TestNumber
Df Sum Sq Mean Sq F value Pr(>F)
TestNumber 3 24.01 8.002 3.005 0.0725 .
Residuals 12 31.96 2.663
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1This method fits a linear mixed-effects model.
lme.model <- lme(Weight ~ TestNumber,
random = ~1|Subject/TestNumber,
data=obs)
summary(lme.model)Linear mixed-effects model fit by REML
Data: obs
AIC BIC logLik
105.7959 111.204 -45.89794
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 16.9536
Formula: ~1 | TestNumber %in% Subject
(Intercept) Residual
StdDev: 1.213168 1.09147
Fixed effects: Weight ~ TestNumber
Value Std.Error DF t-value p-value
(Intercept) 141.20 7.616922 12 18.537672 0.0000
TestNumber1 -0.58 1.032101 12 -0.561960 0.5845
TestNumber2 -2.86 1.032101 12 -2.771046 0.0169
TestNumber3 -1.70 1.032101 12 -1.647125 0.1254
Correlation:
(Intr) TstNm1 TstNm2
TestNumber1 -0.068
TestNumber2 -0.068 0.500
TestNumber3 -0.068 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.7318710 -0.3871376 -0.1167672 0.4124592 1.1006606
Number of Observations: 20
Number of Groups:
Subject TestNumber %in% Subject
5 20 Display F- and p-values for the factor effect.
anova(lme.model) numDF denDF F-value p-value
(Intercept) 1 12 339.7583 <.0001
TestNumber 3 12 3.0047 0.0725obs.matrix <- with(obs, # data from obs data frame
cbind(
Weight[TestNumber=="0"], # subject observations baseline
Weight[TestNumber=="1"], # subject observations test 1
Weight[TestNumber=="2"], # etc
Weight[TestNumber=="3"]
) # cbind column bind, rows are flipped to columns
)
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$TestNumber)
# display the matrix
obs.matrix 0 1 2 3
112121 147 150.0 144.5 145.0
258996 128 128.0 126.2 128.5
456121 165 162.0 159.0 159.0
522210 148 145.0 144.0 146.0
563751 118 118.1 118.0 119.0mlm.model <- lm(obs.matrix ~ 1)
mlm.model
Call:
lm(formula = obs.matrix ~ 1)
Coefficients:
0 1 2 3
(Intercept) 141.2 140.6 138.3 139.5The coefficients are equal to the means across subjects at each test.
rfactor <- as.factor(levels(obs$TestNumber))
rfactor[1] 0 1 2 3
Levels: 0 1 2 3mlm.model.aov.mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
idata = data.frame(rfactor),
idesign = ~rfactor,
type = "III")
summary(mlm.model.aov,
multivariate = FALSE) # don't show multivariate testsHF eps > 1 treated as 1
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 391524 1 4609.4 4 339.7583 5.097e-05 ***
rfactor 24 3 32.0 12 3.0047 0.0725 .
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
rfactor 0.17448 0.47272
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
rfactor 0.6815 0.1045
HF eps Pr(>F[HF])
rfactor 1.40161 0.07249535MaxRelVO2.Subject.TestNumber.aov.model <- aov(MaxRelVO2 ~ TestNumber + Error(Subject/TestNumber), data=obs)
summary(aov.model)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 1401 350.1
Error: Subject:TestNumber
Df Sum Sq Mean Sq F value Pr(>F)
TestNumber 3 48.72 16.239 1.837 0.194
Residuals 12 106.10 8.841 This method fits a linear mixed-effects model.
lme.model <- lme(MaxRelVO2 ~ TestNumber,
random = ~1|Subject/TestNumber,
data=obs)
summary(lme.model)Linear mixed-effects model fit by REML
Data: obs
AIC BIC logLik
115.4304 120.8385 -50.71519
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 9.237129
Formula: ~1 | TestNumber %in% Subject
(Intercept) Residual
StdDev: 2.539497 1.546705
Fixed effects: MaxRelVO2 ~ TestNumber
Value Std.Error DF t-value p-value
(Intercept) 55.152 4.339721 12 12.708650 0.0000
TestNumber1 1.008 1.880568 12 0.536008 0.6017
TestNumber2 2.150 1.880568 12 1.143271 0.2752
TestNumber3 4.200 1.880568 12 2.233367 0.0453
Correlation:
(Intr) TstNm1 TstNm2
TestNumber1 -0.217
TestNumber2 -0.217 0.500
TestNumber3 -0.217 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.83507455 -0.21908727 -0.01385418 0.31417074 0.62269071
Number of Observations: 20
Number of Groups:
Subject TestNumber %in% Subject
5 20 Display F- and p-values for the factor effect.
anova(lme.model) numDF denDF F-value p-value
(Intercept) 1 12 185.52779 <.0001
TestNumber 3 12 1.83674 0.1942obs.matrix <- with(obs, # data from obs data frame
cbind(
MaxRelVO2[TestNumber=="0"], # subject observations baseline
MaxRelVO2[TestNumber=="1"], # subject observations test 1
MaxRelVO2[TestNumber=="2"], # etc
MaxRelVO2[TestNumber=="3"]
) # cbind column bind, rows are flipped to columns
)
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$TestNumber)
# display the matrix
obs.matrix 0 1 2 3
112121 66.17 70.01 68.76 70.11
258996 60.41 59.63 61.71 64.89
456121 42.86 48.47 50.25 44.33
522210 55.91 58.37 60.57 63.48
563751 50.41 44.32 45.22 53.95mlm.model <- lm(obs.matrix ~ 1)
mlm.model
Call:
lm(formula = obs.matrix ~ 1)
Coefficients:
0 1 2 3
(Intercept) 55.15 56.16 57.30 59.35The coefficients are equal to the means across subjects at each test.
rfactor <- as.factor(levels(obs$TestNumber))
rfactor[1] 0 1 2 3
Levels: 0 1 2 3mlm.model.aov.mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
idata = data.frame(rfactor),
idesign = ~rfactor,
type = "III")
summary(mlm.model.aov,
multivariate = FALSE) # don't show multivariate tests
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 64961 1 1400.6 4 185.5278 0.0001682 ***
rfactor 49 3 106.1 12 1.8367 0.1942032
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
rfactor 0.071141 0.23259
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
rfactor 0.42679 0.2402
HF eps Pr(>F[HF])
rfactor 0.5395167 0.2311513FatBIA.Subject.TestNumber.aov.model <- aov(FatBIA ~ TestNumber + Error(Subject/TestNumber), data=obs)
summary(aov.model)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 1516 379
Error: Subject:TestNumber
Df Sum Sq Mean Sq F value Pr(>F)
TestNumber 3 18.90 6.30 2.258 0.134
Residuals 12 33.47 2.79 This method fits a linear mixed-effects model.
lme.model <- lme(FatBIA ~ TestNumber,
random = ~1|Subject/TestNumber,
data=obs)
summary(lme.model)Linear mixed-effects model fit by REML
Data: obs
AIC BIC logLik
101.9044 107.3126 -43.95222
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 9.697831
Formula: ~1 | TestNumber %in% Subject
(Intercept) Residual
StdDev: 1.336557 1.001598
Fixed effects: FatBIA ~ TestNumber
Value Std.Error DF t-value p-value
(Intercept) 37.2 4.400852 12 8.452908 0.0000
TestNumber1 -1.5 1.056330 12 -1.420011 0.1811
TestNumber2 -2.7 1.056330 12 -2.556019 0.0252
TestNumber3 -1.0 1.056330 12 -0.946674 0.3625
Correlation:
(Intr) TstNm1 TstNm2
TestNumber1 -0.12
TestNumber2 -0.12 0.50
TestNumber3 -0.12 0.50 0.50
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.74569729 -0.32742017 -0.03294121 0.20867144 1.05087005
Number of Observations: 20
Number of Groups:
Subject TestNumber %in% Subject
5 20 Display F- and p-values for the factor effect.
anova(lme.model) numDF denDF F-value p-value
(Intercept) 1 12 68.01445 <.0001
TestNumber 3 12 2.25840 0.134obs.matrix <- with(obs, # data from obs data frame
cbind(
FatBIA[TestNumber=="0"], # subject observations baseline
FatBIA[TestNumber=="1"], # subject observations test 1
FatBIA[TestNumber=="2"], # etc
FatBIA[TestNumber=="3"]
) # cbind column bind, rows are flipped to columns
)
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$TestNumber)
# display the matrix
obs.matrix 0 1 2 3
112121 33.0 35.5 31.0 31.5
258996 28.0 27.0 27.0 30.0
456121 55.5 53.0 52.0 51.0
522210 35.5 33.0 32.0 34.0
563751 34.0 30.0 30.5 34.5mlm.model <- lm(obs.matrix ~ 1)
mlm.model
Call:
lm(formula = obs.matrix ~ 1)
Coefficients:
0 1 2 3
(Intercept) 37.2 35.7 34.5 36.2The coefficients are equal to the means across subjects at each test.
rfactor <- as.factor(levels(obs$TestNumber))
rfactor[1] 0 1 2 3
Levels: 0 1 2 3mlm.model.aov.mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
idata = data.frame(rfactor),
idesign = ~rfactor,
type = "III")
summary(mlm.model.aov,
multivariate = FALSE) # don't show multivariate tests
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 25776.2 1 1515.93 4 68.0144 0.001179 **
rfactor 18.9 3 33.47 12 2.2584 0.133953
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
rfactor 0.12845 0.37598
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
rfactor 0.56057 0.1791
HF eps Pr(>F[HF])
rfactor 0.9214632 0.1409945FatSF.Subject.TestNumber.aov.model <- aov(FatSF ~ TestNumber + Error(Subject/TestNumber), data=obs)
summary(aov.model)
Error: Subject
Df Sum Sq Mean Sq F value Pr(>F)
Residuals 4 1224 306.1
Error: Subject:TestNumber
Df Sum Sq Mean Sq F value Pr(>F)
TestNumber 3 3.686 1.2285 1.776 0.205
Residuals 12 8.302 0.6918 This method fits a linear mixed-effects model.
lme.model <- lme(FatSF ~ TestNumber,
random = ~1|Subject/TestNumber,
data=obs)
summary(lme.model)Linear mixed-effects model fit by REML
Data: obs
AIC BIC logLik
84.31884 89.72696 -35.15942
Random effects:
Formula: ~1 | Subject
(Intercept)
StdDev: 8.73833
Formula: ~1 | TestNumber %in% Subject
(Intercept) Residual
StdDev: 0.6174158 0.5573429
Fixed effects: FatSF ~ TestNumber
Value Std.Error DF t-value p-value
(Intercept) 28.62 3.925564 12 7.290673 0.0000
TestNumber1 0.44 0.526054 12 0.836415 0.4193
TestNumber2 -0.48 0.526054 12 -0.912453 0.3795
TestNumber3 -0.66 0.526054 12 -1.254623 0.2335
Correlation:
(Intr) TstNm1 TstNm2
TestNumber1 -0.067
TestNumber2 -0.067 0.500
TestNumber3 -0.067 0.500 0.500
Standardized Within-Group Residuals:
Min Q1 Med Q3 Max
-0.66970985 -0.40163206 -0.07515769 0.24426689 1.19928863
Number of Observations: 20
Number of Groups:
Subject TestNumber %in% Subject
5 20 Display F- and p-values for the factor effect.
anova(lme.model) numDF denDF F-value p-value
(Intercept) 1 12 52.86185 <.0001
TestNumber 3 12 1.77572 0.2052obs.matrix <- with(obs, # data from obs data frame
cbind(
FatSF[TestNumber=="0"], # subject observations baseline
FatSF[TestNumber=="1"], # subject observations test 1
FatSF[TestNumber=="2"], # etc
FatSF[TestNumber=="3"]
) # cbind column bind, rows are flipped to columns
)
# set row and column identities
rownames(obs.matrix) <- levels(obs$Subject)
colnames(obs.matrix) <- levels(obs$TestNumber)
# display the matrix
obs.matrix 0 1 2 3
112121 22.4 25.0 22.0 22.6
258996 23.1 22.6 21.6 22.5
456121 43.3 43.6 44.5 42.0
522210 25.1 24.1 23.7 24.0
563751 29.2 30.0 28.9 28.7mlm.model <- lm(obs.matrix ~ 1)
mlm.model
Call:
lm(formula = obs.matrix ~ 1)
Coefficients:
0 1 2 3
(Intercept) 28.62 29.06 28.14 27.96The coefficients are equal to the means across subjects at each test.
rfactor <- as.factor(levels(obs$TestNumber))
rfactor[1] 0 1 2 3
Levels: 0 1 2 3mlm.model.aov.mlm.model.aov <- Anova(mlm.model, # the multivariate linear model
idata = data.frame(rfactor),
idesign = ~rfactor,
type = "III")
summary(mlm.model.aov,
multivariate = FALSE) # don't show multivariate testsHF eps > 1 treated as 1
Univariate Type III Repeated-Measures ANOVA Assuming Sphericity
SS num Df Error SS den Df F Pr(>F)
(Intercept) 16182.4 1 1224.5 4 52.8619 0.001901 **
rfactor 3.7 3 8.3 12 1.7757 0.205231
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Mauchly Tests for Sphericity
Test statistic p-value
rfactor 0.014122 0.053215
Greenhouse-Geisser and Huynh-Feldt Corrections
for Departure from Sphericity
GG eps Pr(>F[GG])
rfactor 0.66036 0.2305
HF eps Pr(>F[HF])
rfactor 1.305202 0.2052308