Remember you should
This practice reviews the More ANOVAs lecuture.
Following the memory example from class, read in and check data
memory <- read.table("http://www.statsci.org/data/general/eysenck.txt", header = T,
stringsAsFactors = T)
str(memory)'data.frame': 100 obs. of 3 variables:
$ Age : Factor w/ 2 levels "Older","Younger": 2 2 2 2 2 2 2 2 2 2 ...
$ Process: Factor w/ 5 levels "Adjective","Counting",..: 2 2 2 2 2 2 2 2 2 2 ...
$ Words : num 8 6 4 6 7 6 5 7 9 7 ...Let’s put younger level first
library(plyr)
memory$Age <- relevel(memory$Age, "Younger")and graph
library(Rmisc)Loading required package: latticefunction_output <- summarySE(memory, measurevar="Words", groupvars =
c("Age", "Process"), na.rm = T)
library(ggplot2)
ggplot(function_output, aes(x=Age, y=Words,color=Process,
shape = Process)) +
geom_line(aes(group=Process, linetype = Process), size=2) +
geom_point(size = 5) +
ylab("Words remembered")+
xlab("Age") +
ggtitle("Process type interacts with \n age to impact memory")+
theme(axis.title.x = element_text(face="bold", size=28),
axis.title.y = element_text(face="bold", size=28),
axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20),
legend.text =element_text(size=20),
legend.title = element_text(size=20, face="bold"),
plot.title = element_text(hjust = 0.5, face="bold", size=32))Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
ℹ Please use `linewidth` instead.
There appears to be some interactions. Let’ build a model
memory_interactions <- lm(Words ~ Age * Process, memory)and check assumptions.
par(mfrow=c(2,2))
plot(memory_interactions)
These appear to be met, so look at output
library(car)Warning: package 'car' was built under R version 4.4.1Loading required package: carDataAnova(memory_interactions, type = "III")Anova Table (Type III tests)
Response: Words
Sum Sq Df F value Pr(>F)
(Intercept) 2190.4 1 272.9281 < 2.2e-16 ***
Age 72.2 1 8.9963 0.0034984 **
Process 1353.7 4 42.1690 < 2.2e-16 ***
Age:Process 190.3 4 5.9279 0.0002793 ***
Residuals 722.3 90
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since interaction is significant, analyze subsets. For example,
memory_interactions_young <- lm(Words ~ Process, memory[memory$Age == "Younger",])
plot(memory_interactions_young)
Anova(memory_interactions_young, type = "III")Anova Table (Type III tests)
Response: Words
Sum Sq Df F value Pr(>F)
(Intercept) 2190.4 1 343.442 < 2.2e-16 ***
Process 1353.7 4 53.064 < 2.2e-16 ***
Residuals 287.0 45
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1There is a significant difference in words recalled based on process, but why? Investigate with post-hoc tests.
library(multcomp)Loading required package: mvtnormLoading required package: survivalLoading required package: TH.dataLoading required package: MASS
Attaching package: 'TH.data'The following object is masked from 'package:MASS':
geysercomp_young <- glht(memory_interactions_young, linfct = mcp(Process = "Tukey"))
summary(comp_young)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Words ~ Process, data = memory[memory$Age == "Younger",
])
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Counting - Adjective == 0 -8.300 1.129 -7.349 < 1e-04 ***
Imagery - Adjective == 0 2.800 1.129 2.479 0.11370
Intentional - Adjective == 0 4.500 1.129 3.984 0.00222 **
Rhyming - Adjective == 0 -7.200 1.129 -6.375 < 1e-04 ***
Imagery - Counting == 0 11.100 1.129 9.828 < 1e-04 ***
Intentional - Counting == 0 12.800 1.129 11.333 < 1e-04 ***
Rhyming - Counting == 0 1.100 1.129 0.974 0.86545
Intentional - Imagery == 0 1.700 1.129 1.505 0.56457
Rhyming - Imagery == 0 -10.000 1.129 -8.854 < 1e-04 ***
Rhyming - Intentional == 0 -11.700 1.129 -10.359 < 1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)Following feather color example from class:
# more than 2? ####
feather <- read.csv("https://raw.githubusercontent.com/jsgosnell/CUNY-BioStats/master/datasets/wiebe_2002_example.csv", stringsAsFactors = T)
str(feather)'data.frame': 32 obs. of 3 variables:
$ Bird : Factor w/ 16 levels "A","B","C","D",..: 1 2 3 4 5 6 7 8 9 10 ...
$ Feather : Factor w/ 2 levels "Odd","Typical": 2 2 2 2 2 2 2 2 2 2 ...
$ Color_index: num -0.255 -0.213 -0.19 -0.185 -0.045 -0.025 -0.015 0.003 0.015 0.02 ...set.seed(25)
special <- data.frame(Bird = LETTERS[1:16], Feather = "Special",
Color_index= feather[feather$Feather == "Typical", "Color_index"] +
.3 +runif(16,1,1)*.01)
feather <- merge(feather, special, all = T)
Anova(lm(Color_index ~ Feather + Bird, data=feather), type= "III")Anova Table (Type III tests)
Response: Color_index
Sum Sq Df F value Pr(>F)
(Intercept) 0.36392 1 59.9538 1.224e-08 ***
Feather 1.67906 2 138.3093 7.208e-16 ***
Bird 0.34649 15 3.8055 0.0008969 ***
Residuals 0.18210 30
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1library(multcomp)
compare <- glht(lm(Color_index ~ Feather + Bird, data=feather), linfct = mcp("Feather" = "Tukey"))
summary(compare)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Color_index ~ Feather + Bird, data = feather)
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
Typical - Odd == 0 0.13712 0.02755 4.978 <1e-04 ***
Special - Odd == 0 0.44712 0.02755 16.232 <1e-04 ***
Special - Typical == 0 0.31000 0.02755 11.254 <1e-04 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)#note comparison doesn't work
Anova(lm(Color_index ~ Feather * Bird, data=feather), type= "III")Error in Anova.lm(lm(Color_index ~ Feather * Bird, data = feather), type = "III"): residual df = 0Swirl is an R package that provides guided lessons to help you learn and review material. These lessons should serve as a bridge between all the code provided in the slides and background reading and the key functions and concepts from each lesson. A full course lesson (all lessons combined) can also be downloaded using the following instructions.
THIS IS ONE OF THE FEW TIMES I RECOMMEND WORKING DIRECTLY IN THE CONSOLE! THERE IS NO NEED TO DEVELOP A SCRIPT FOR THESE INTERACTIVE SESSIONS, THOUGH YOU CAN!
install the “swirl” package
run the following code once on the computer to install a new course
library(swirl)
install_course_github("jsgosnell", "JSG_swirl_lessons")start swirl!
swirl()then follow the on-screen prompts to select the JSG_swirl_lessons course and the lessons you want
TIP: If you are seeing duplicate courses (or odd versions of each), you can clear all courses and then re-download the courses by
exiting swirl using escape key or bye() function
bye()uninstalling and reinstalling courses
uninstall_all_courses()
install_course_github("jsgosnell", "JSG_swirl_lessons")when you restart swirl with swirl(), you may need to select
A survey was conducted to see if athletes and non-athletes deal with anger in the same way. Data is @
angry <- read.csv(“https://docs.google.com/spreadsheets/d/e/2PACX-1vSaawG37o1ZUEs1B4keIJpZAY2c5tuljf29dWnzqQ0tHNCzfbz85AlWobYzBQ3nPPXJBLP-FWe4BNZB/pub?gid=1784556512&single=true&output=csv”, stringsAsFactors = T)
and more information is at
http://onlinestatbook.com/case_studies/angry_moods.html.
Focus on the following variables:
Sports 1 = athletes, 2 = non-athletes Gender 1 = males, 2 = females Expression (AE) index of general anger expression: (Anger-Out) + (Anger-In) - (Control-Out) - (Control-In) + 48
Is there any evidence that gender or athlete status impact how anger is expressed?
angry <- read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vSaawG37o1ZUEs1B4keIJpZAY2c5tuljf29dWnzqQ0tHNCzfbz85AlWobYzBQ3nPPXJBLP-FWe4BNZB/pub?gid=1784556512&single=true&output=csv", stringsAsFactors = T)
str(angry)'data.frame': 78 obs. of 7 variables:
$ Gender : int 2 2 2 2 1 1 1 2 2 2 ...
$ Sports : int 1 1 1 1 1 1 1 1 1 1 ...
$ Anger.Out : int 18 14 13 17 16 16 12 13 16 12 ...
$ Anger.In : int 13 17 14 24 17 22 12 16 16 16 ...
$ Control.Out : int 23 25 28 23 26 25 31 22 22 29 ...
$ Control.In : int 20 24 28 23 28 23 27 31 24 29 ...
$ Anger_Expression: int 36 30 19 43 27 38 14 24 34 18 ...angry$Gender <- as.factor(angry$Gender)
library(plyr)
angry$Gender <- revalue(angry$Gender, c("1" = "male",
"2" = "female"))
angry$Sports <- as.factor(angry$Sports)
angry$Sports <- revalue(angry$Sports, c("1" = "athlete",
"2" = "non-athlete"))
library(Rmisc)
anger_summary <- summarySE(angry, measurevar="Anger_Expression", groupvars =
c("Sports", "Gender"), na.rm = T)
library(ggplot2)
ggplot(anger_summary, aes(x=Gender, y=Anger_Expression, color=Sports,
shape = Sports)) +
geom_point(size = 3) +
geom_line(aes(group=Sports, linetype =Sports), size=2) +
geom_errorbar(aes(ymin=Anger_Expression-ci, ymax=Anger_Expression+ci), size=1.5) +
ylab("Anger level")+
xlab("Gender") +
scale_shape_discrete(guide=FALSE)+
scale_linetype_discrete(guide=FALSE)+
ggtitle("Anger level among groups")+
theme(axis.title.x = element_text(face="bold", size=28),
axis.title.y = element_text(face="bold", size=28),
axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20),
legend.text =element_text(size=20),
legend.title = element_text(size=20, face="bold"),
plot.title = element_text(hjust = 0.5, face="bold", size=32))Warning: The `guide` argument in `scale_*()` cannot be `FALSE`. This was deprecated in
ggplot2 3.3.4.
ℹ Please use "none" instead.
I first read in and recoded some data for ease and plotting. I then produced a plot to consider the null hypotheses that
angry_gender <- lm(Anger_Expression ~ Sports * Gender, angry)
plot(angry_gender)
library(car)
Anova(angry_gender, type = "III")Anova Table (Type III tests)
Response: Anger_Expression
Sum Sq Df F value Pr(>F)
(Intercept) 11200.1 1 71.8690 1.617e-12 ***
Sports 480.1 1 3.0807 0.08336 .
Gender 17.7 1 0.1135 0.73711
Sports:Gender 5.2 1 0.0336 0.85505
Residuals 11532.2 74
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#remove interaction since not significant
angry_gender <- lm(Anger_Expression ~ Sports + Gender, angry)
plot(angry_gender)
Anova(angry_gender, type = "III") #only differs among those who play sportsAnova Table (Type III tests)
Response: Anger_Expression
Sum Sq Df F value Pr(>F)
(Intercept) 17367.1 1 112.8964 < 2.2e-16 ***
Sports 1357.2 1 8.8227 0.003995 **
Gender 16.3 1 0.1061 0.745501
Residuals 11537.4 75
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1I then analyzed the data using factorial ANOVA. The outcome is continuous and both explanatory variables are categorical. The design is also fully randomized. Resildual plots indicated all assumptions were met (there is no pattern in the residuals and they are normally distributed). Analysis shows an insignificant interaction (F1,74=.04, p=.855) between sport and gender, so I removed the interaction term. The reduced model showed anger levels differed among athletes and non-athletes but not by gender. There was no need for post-hoc tests (only 2 levels/groups for each categorial variable.)
A professor carried out a long-term study to see how various factors impacted pulse rate before and after exercise. Data can be found at http://www.statsci.org/data/oz/ms212.txt With more info at http://www.statsci.org/data/oz/ms212.html. Is there evidence that frequency of exercise (Exercise column) and gender impact change in pulse rate for students who ran (Ran column = 1)?
pulse <- read.table("http://www.statsci.org/data/oz/ms212.txt", header = T,
stringsAsFactors = T)
pulse$Exercise <- factor(pulse$Exercise)
library(plyr)
pulse$Exercise <- revalue(pulse$Exercise, c("1" = "high",
"2" = "moderate",
"3" = "low"))
pulse$Gender <- factor(pulse$Gender)
pulse$Gender <- revalue (pulse$Gender, c("1" = "male", "2" = "female"))
pulse$change <- pulse$Pulse2 - pulse$Pulse1
change_summary <- summarySE(pulse[pulse$Ran == 1, ], measurevar="change", groupvars =
c("Exercise", "Gender"), na.rm = T)Warning in qt(conf.interval/2 + 0.5, datac$N - 1): NaNs producedggplot(change_summary, aes(x=Gender, shape = Exercise, color = Exercise,
y=change)) +
geom_point(size = 3) +
geom_line(aes(group=Exercise, linetype =Exercise), size=2) +
geom_errorbar(aes(ymin=change-ci, ymax=change+ci), size=1.5) +
ylab("Change in pulse \n (beats per minute)") +
scale_color_discrete(name = "Exercise level")+
scale_shape_discrete(guide=FALSE)+
scale_linetype_discrete(guide=FALSE)+
ggtitle("Change in pulse does \n not differ among groups") +
theme(axis.title.x = element_text(face="bold", size=28),
axis.title.y = element_text(face="bold", size=28),
axis.text.y = element_text(size=20),
axis.text.x = element_text(size=20),
legend.text =element_text(size=20),
legend.title = element_text(size=20, face="bold"),
plot.title = element_text(hjust = 0.5, face="bold", size=32))
I first read in and recoded some data for ease and plotting. I then produced a plot to consider the null hypotheses that
exercise <- lm(change ~ Gender * Exercise, pulse[pulse$Ran == 1, ])
summary(exercise)
Call:
lm(formula = change ~ Gender * Exercise, data = pulse[pulse$Ran ==
1, ])
Residuals:
Min 1Q Median 3Q Max
-44.231 -11.300 1.769 10.083 48.444
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 45.000 8.746 5.145 7.44e-06 ***
Genderfemale 35.000 23.139 1.513 0.1383
Exercisemoderate 9.231 10.573 0.873 0.3879
Exerciselow 12.400 12.972 0.956 0.3449
Genderfemale:Exercisemoderate -38.231 24.678 -1.549 0.1292
Genderfemale:Exerciselow -46.844 26.043 -1.799 0.0796 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.42 on 40 degrees of freedom
Multiple R-squared: 0.0828, Adjusted R-squared: -0.03185
F-statistic: 0.7222 on 5 and 40 DF, p-value: 0.6107Anova(exercise, type = "III")Anova Table (Type III tests)
Response: change
Sum Sq Df F value Pr(>F)
(Intercept) 12150.0 1 26.4739 7.444e-06 ***
Gender 1050.0 1 2.2879 0.1383
Exercise 496.3 2 0.5407 0.5865
Gender:Exercise 1488.8 2 1.6220 0.2102
Residuals 18357.7 40
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1#rerun without interaction
exercise <- lm(change ~ Gender + Exercise, pulse[pulse$Ran == 1, ])
summary(exercise)
Call:
lm(formula = change ~ Gender + Exercise, data = pulse[pulse$Ran ==
1, ])
Residuals:
Min 1Q Median 3Q Max
-43.99 -14.89 2.23 11.91 45.19
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 50.391 8.273 6.091 2.94e-07 ***
Genderfemale -2.738 6.770 -0.404 0.688
Exercisemoderate 3.603 9.572 0.376 0.708
Exerciselow 1.155 10.617 0.109 0.914
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 21.74 on 42 degrees of freedom
Multiple R-squared: 0.008416, Adjusted R-squared: -0.06241
F-statistic: 0.1188 on 3 and 42 DF, p-value: 0.9485Anova(exercise, type = "III") #no significanceAnova Table (Type III tests)
Response: change
Sum Sq Df F value Pr(>F)
(Intercept) 17532.0 1 37.1019 2.937e-07 ***
Gender 77.3 1 0.1636 0.6879
Exercise 97.1 2 0.1028 0.9025
Residuals 19846.5 42
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1I then analyzed the data using factorial ANOVA. The outcome is continuous and both explanatory variables are categorical. The design is also fully randomized. Residual plots indicated all assumptions were met (there is no pattern in the residuals and they are normally distributed). Analysis shows an insignificant interaction (F2,40=1.2, p=.21) between exercise level and gender, so I removed the interaction term. The reduced model showed neither gender (F1,42=.16, p =.69) or exercise level (F2,42=.1, p=.90) influenced change in pulse rate, so I failed to reject the related null hypotheses.
Data from Valdez et al 2023 is available @ https://docs.google.com/spreadsheets/d/e/2PACX-1vT2gaLu6pyRMlcbzarn3ej4bFmT_iHvrlNWJYSdrsLdUWIjcJi7rU11-ipvYpGnqD9qLDnbhNd2sDUW/pub?gid=1707080634&single=true&output=csv.
Import it into to R and
valdez_2023 <- read.csv("https://docs.google.com/spreadsheets/d/e/2PACX-1vT2gaLu6pyRMlcbzarn3ej4bFmT_iHvrlNWJYSdrsLdUWIjcJi7rU11-ipvYpGnqD9qLDnbhNd2sDUW/pub?gid=1707080634&single=true&output=csv", stringsAsFactors = T)
shoot_model <-lm( Shoot.density..m2~Snail.Level + Nitrogen.level + Snail.Level:Nitrogen.level, valdez_2023[valdez_2023$Snail.Level != "uncaged",])
plot(shoot_model)
Anova(shoot_model, type="III")Anova Table (Type III tests)
Response: Shoot.density..m2
Sum Sq Df F value Pr(>F)
(Intercept) 196608 1 271.0588 1.334e-09 ***
Snail.Level 36238 2 24.9804 5.277e-05 ***
Nitrogen.level 4267 1 5.8824 0.0320 *
Snail.Level:Nitrogen.level 3100 2 2.1373 0.1607
Residuals 8704 12
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1No interaction, so I can remove or not.
shoot_model_reduced <- update(shoot_model, .~.- Snail.Level:Nitrogen.level)
Anova(shoot_model_reduced)Anova Table (Type II tests)
Response: Shoot.density..m2
Sum Sq Df F value Pr(>F)
Snail.Level 45596 2 27.039 1.556e-05 ***
Nitrogen.level 11150 1 13.224 0.002696 **
Residuals 11804 14
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Both have impacts.
snail_post_hoc <- glht(shoot_model_reduced, linfct = mcp(Snail.Level = "Tukey"))
summary(snail_post_hoc)
Simultaneous Tests for General Linear Hypotheses
Multiple Comparisons of Means: Tukey Contrasts
Fit: lm(formula = Shoot.density..m2 ~ Snail.Level + Nitrogen.level,
data = valdez_2023[valdez_2023$Snail.Level != "uncaged",
])
Linear Hypotheses:
Estimate Std. Error t value Pr(>|t|)
removal - control snails == 0 50.67 16.76 3.022 0.02306 *
snail addition - control snails == 0 -72.00 16.76 -4.295 0.00201 **
snail addition - removal == 0 -122.67 16.76 -7.317 < 0.001 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Adjusted p values reported -- single-step method)significant difference among all snails.
shoot_summary <- summarySE(valdez_2023[valdez_2023$Snail.Level != "uncaged",], measurevar = "Shoot.density..m2", groupvars = c("Snail.Level", "Nitrogen.level"))
shoot_summary Snail.Level Nitrogen.level N Shoot.density..m2 sd se
1 control snails Fertilized 3 256.0000 27.712813 16.000000
2 control snails without 3 202.6667 9.237604 5.333333
3 removal Fertilized 3 320.0000 16.000000 9.237604
4 removal without 3 240.0000 27.712813 16.000000
5 snail addition Fertilized 3 165.3333 36.950417 21.333333
6 snail addition without 3 149.3333 33.306656 19.229607
ci
1 68.84244
2 22.94748
3 39.74620
4 68.84244
5 91.78992
6 82.73832shoot_summary$Snail.Level <- relevel(shoot_summary$Snail.Level, "removal")
shoot_summary$Snail.Level <- relevel(shoot_summary$Snail.Level, "uncaged")
ggplot(shoot_summary, aes(x=Snail.Level,
y=Shoot.density..m2,
fill=Nitrogen.level)) +
geom_col(color="black", position=position_dodge()) +
geom_errorbar(aes(ymin=Shoot.density..m2, ymax=Shoot.density..m2+ci), position = position_dodge()) +
labs(title="Grazing impacts depend on nitrogen levels",
x= "Grazing level",
y= expression(paste("# of shoots/" ,m^{-2})))
Find an example of a factorial ANOVA from a paper that is related to your work. Make sure you understand the connections between the methods, results, and graphs. Briefly answer the following questions