# packages
pacman::p_load(car, sjstats, jtools, agricolae, tidyverse, dplyr, tidyr, readxl, ggplot2)
ath = read_excel('DataSets.xlsx',
sheet = 'anova',
range = 'Q15:T123')
str(ath)
## tibble [108 × 4] (S3: tbl_df/tbl/data.frame)
## $ sr. : num [1:108] 1 2 3 4 5 6 7 8 9 10 ...
## $ Treatment: chr [1:108] "None" "None" "None" "None" ...
## $ Exercise : chr [1:108] "Rest" "Rest" "Rest" "Rest" ...
## $ Rate : num [1:108] 53.5 88.5 76.7 44.9 66 ...
ath$Treatment = as.factor(ath$Treatment)
ath$Exercise = as.factor(ath$Exercise)
str(ath)
## tibble [108 × 4] (S3: tbl_df/tbl/data.frame)
## $ sr. : num [1:108] 1 2 3 4 5 6 7 8 9 10 ...
## $ Treatment: Factor w/ 3 levels "Drug-lesion",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ Exercise : Factor w/ 2 levels "During","Rest": 2 2 2 2 2 2 2 2 2 2 ...
## $ Rate : num [1:108] 53.5 88.5 76.7 44.9 66 ...
# two-way anova: additive
mod2 = aov(Rate ~ Treatment + Exercise, ath)
summary(mod2)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 2 2525 1262 10.43 7.47e-05 ***
## Exercise 1 154431 154431 1275.35 < 2e-16 ***
## Residuals 104 12593 121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Interactive model
modint = aov(Rate ~ Treatment + Exercise + Treatment:Exercise, ath)
summary(modint)
## Df Sum Sq Mean Sq F value Pr(>F)
## Treatment 2 2525 1262 10.401 7.75e-05 ***
## Exercise 1 154431 154431 1272.299 < 2e-16 ***
## Treatment:Exercise 2 213 106 0.876 0.42
## Residuals 102 12381 121
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(ath$Treatment)
## Drug-lesion Lesion None
## 36 36 36
summary(ath$Exercise)
## During Rest
## 54 54
# interaction plot
with(ath, interaction.plot(Treatment, Exercise, Rate))
interaction.plot(ath$Treatment, ath$Exercise, ath$Rate)
# interaction is not significant, so proceed with additive model
# post hoc
tukey = HSD.test(mod2, trt = c('Treatment', 'Exercise'))
a = tukey$mean %>% rownames_to_column("Variable")
b = tukey$groups %>% rownames_to_column("Variable")
figdata = left_join(a[c(1:4)], b[c(1,3)], by = 'Variable')
figdata
## Variable Rate std r groups
## 1 Drug-lesion:During 151.84744 15.368802 18 a
## 2 Drug-lesion:Rest 72.25709 8.452339 18 c
## 3 Lesion:During 150.13009 9.284707 18 a
## 4 Lesion:Rest 76.28621 9.080305 18 c
## 5 None:During 139.14789 10.968733 18 b
## 6 None:Rest 65.69673 11.474478 18 c
figdata = figdata %>% mutate(
se = std/sqrt(r),
UL = Rate + se*1.96,
LL = Rate - se*1.96
)
figdata = figdata %>%
separate(Variable, sep = ":",
into = c('Lesion', 'Exercise'))
figdata
## Lesion Exercise Rate std r groups se UL
## 1 Drug-lesion During 151.84744 15.368802 18 a 3.622461 158.94747
## 2 Drug-lesion Rest 72.25709 8.452339 18 c 1.992235 76.16187
## 3 Lesion During 150.13009 9.284707 18 a 2.188426 154.41940
## 4 Lesion Rest 76.28621 9.080305 18 c 2.140248 80.48110
## 5 None During 139.14789 10.968733 18 b 2.585355 144.21519
## 6 None Rest 65.69673 11.474478 18 c 2.704560 70.99767
## LL
## 1 144.74742
## 2 68.35231
## 3 145.84077
## 4 72.09133
## 5 134.08060
## 6 60.39579
p = position_dodge(width = 0.9)
ggplot(figdata) +
aes(x = Lesion, y = Rate, color = Lesion, fill = Exercise) +
geom_col(position = p) +
geom_errorbar(aes(ymin = LL, ymax = UL), width = 0.2, color = 'black', position = p) +
geom_text(aes(y = UL+10, label = groups), color = 'grey20', position = p)+
theme_bw() +
theme(legend.position = 'top')
