Two-way ANOVA
Exercise 1
Download the CSV file “growth” from Canvas and read it into R using the import option from the Environment tab.
Call the data frame growth_df. Read Diet and Supplement as factors.
setwd("E:\\S9510\\CAP3330")
growth_df = read.csv("growth.csv", stringsAsFactors = T, header = T)After reading the csv file, take a look at the data frame that was created:
str(growth_df)## 'data.frame': 48 obs. of 3 variables:
## $ supplement: Factor w/ 4 levels "agrimore","control",..: 3 3 3 3 2 2 2 2 4 4 ...
## $ diet : Factor w/ 3 levels "barley","oats",..: 3 3 3 3 3 3 3 3 3 3 ...
## $ gain : num 17.4 16.8 18.1 15.8 17.7 ...The data comes from a farm-scale trial of animal diets. There are two factors: diet and supplement. Diet is a factor with three levels: barley, oats and wheat. Supplement is a factor with four levels: agrimore, control, supergain, and supersupp. The outcome variable is weight gain after 6 weeks.
Analyze the data and conclude about the effects of diet and supplement on the animals’ weight gain.
Before running ANOVA, let’s answer a few questions:
Number of treatments?
3 diets
4 supplements
3 * 4 = 12 treatments (12 combinations of diets and supplements)
Number of replicates (number of animals tested per treatment)?
Total observations = # treatments * # replicates
48 observations = 12 treatments * # replicates
Number of replicates = 48/12 = 4
_Hypotheses__
For the effect of diet:
Ho: The mean weight gain is the same for all three diets (barley, oats, and wheat). Ha: There is a difference in mean weight gain for at least two diets.
Or, alternatively…
Ho: M barley = M oats = M wheat Ha: Same as above
For the effect of supplement:
H₀: The mean weight gain is the same for all four supplements (agrimore, control, supergain, and supersupp)
Ha: There is a difference in mean weight gain for at least two supplements.
Or, alternatively…
Ho: M agrimore = M control = M supergain = M supersupp Ha: Same as above
For the interaction between diet and supplement:
Ho: No interaction effect (= The effect of diet on weight gain does NOT depend on the type of supplement)
Ha: There is interaction effect (= The effect of diet on weight gain DEPENDS on the type of supplement)
Or, alternatively
Ho: No interaction effect (= The effect of supplement on weight gain does NOT depend on the type of diet)
Ha: There is interaction effect (= The effect of supplement on weight gain DEPENDS on the type of diet)
Run the ANOVA test
twoway_growth = aov(gain ~ diet * supplement, data=growth_df)
summary(twoway_growth)## Df Sum Sq Mean Sq F value Pr(>F)
## diet 2 287.17 143.59 83.52 3.00e-14 ***
## supplement 3 91.88 30.63 17.82 2.95e-07 ***
## diet:supplement 6 3.41 0.57 0.33 0.917
## Residuals 36 61.89 1.72
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpret the results of ANOVA:
The diet effect interpretation: The PV for diet is really tiny (3.00e-14 = 3*10^-14); therefore, it is less than alpha 0.05. We reject Ho and support Ha. We can claim that at least two diets yield different average weight gain.
The supplement effect interpretation: The PV for supplement is really tiny (2.95e-07 = 2.95*10^-7); therefore, it is less than alpha 0.05. We reject Ho and support Ha. We can claim that at least two supplements yield different average weight gain.
The interaction effect interpretation: The PV for the interaction is 0.917, much bigger than alpha 0.05. Therefore, we do not reject Ho and cannot support Ha in case of the interaction effect. The data DO NOT give us evidence to claim that the effect of diet on weight gain depends (is affected by) on the type of supplement.
The data DO NOT give us evidence to claim that there is an interaction between Diet and Supplement.
How to proceed based on the results of this exercise?
Since the interaction is NOT significant, but both main effects are, we should conduct a post-hoc test individually for each factor, diet and supplement.
Side comment: How do we proceed in an exercise where the interaction effect is significant?
When there is interaction, understanding the effect of the two factors on the DV becomes more challenging. It is NOT as simple as doing a post-hoc test individually for each factor and selecting the best combination. Discuss more about this in exercise 2.
End of Side comment
Post-hoc test for each factor
Before we apply the post-hoc test, lets get the sample means for the levels of each factor:
# The sample mean weight gain for the three diets
sort(tapply(growth_df$gain, growth_df$diet, mean))## wheat oats barley
## 18.43134 21.32882 24.42164# The sample mean weight gain for the four supplements
sort(tapply(growth_df$gain, growth_df$supplement, mean))## supergain control supersupp agrimore
## 19.71385 20.39861 22.36796 23.09531Now we apply Tukey to see which of these sample means are statistically different.
# Tukey for diet
TukeyHSD(twoway_growth, which = "diet")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = gain ~ diet * supplement, data = growth_df)
##
## $diet
## diff lwr upr p adj
## oats-barley -3.092817 -4.225918 -1.959715 3e-07
## wheat-barley -5.990298 -7.123399 -4.857196 0e+00
## wheat-oats -2.897481 -4.030582 -1.764379 1e-06Analysis of Tukey for diets
All PVs are less than alpha (0.05). Therefore, all diets are significantly different from each other with respect to average weight gain.
Optional discussion
sort(tapply(growth_df$gain, growth_df$diet, mean))## wheat oats barley
## 18.43134 21.32882 24.42164This is an easy case if you want to choose the best diet because we have a clear hierarchy. What’s the best diet? What’s the second best diet? What’s the worst diet?
The best diet is barley
The second best diet is oats
The worst diet is wheat
End of optional discussion
# Tukey for supplement
TukeyHSD(twoway_growth, which = "supplement")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = gain ~ diet * supplement, data = growth_df)
##
## $supplement
## diff lwr upr p adj
## control-agrimore -2.6967005 -4.138342 -1.2550592 0.0000764
## supergain-agrimore -3.3814586 -4.823100 -1.9398173 0.0000015
## supersupp-agrimore -0.7273521 -2.168993 0.7142892 0.5326710
## supergain-control -0.6847581 -2.126399 0.7568832 0.5817637
## supersupp-control 1.9693484 0.527707 3.4109897 0.0040534
## supersupp-supergain 2.6541065 1.212465 4.0957478 0.0000972Analysis of Tukey for supplements
In this case, not all PVs are less than alpha; therefore, some differences are not significant.
Optional discussion
sort(tapply(growth_df$gain, growth_df$supplement, mean))## supergain control supersupp agrimore
## 19.71385 20.39861 22.36796 23.09531In these situations where some differences are not significant, I like to start by comparing the largest sample mean to the rest and see which means are statistically different from the largest one.
The largest sample mean is obtained with agrimore. Agrimore is statistically higher than supergain and control, but has no significant difference with supersupp.
Then, I repeat this analysis with the second largest mean. The second largest mean is supersupp. Supersupp is statistically higher than control and supergain, but has no significant difference with agrimore.
I can also see that control and supergain are not statistically different from each other.
At this point, I can make a reasonable suggestion about what is the most reasonable supplement (s) to propose. I would propose either agrimore or supersupp because:
- They are better than the other two (supergain and control)
- We did not find a difference between them.
Overall suggestion: Feed the animals with barley and either agrimore or supersupp.
End of optional discussion
I will show you how to do it for this exercise. VERY IMPORTANT: We DO NOT have to do an interaction plot in this example because the interaction is not significant. I am just going to do it to illustrate how such a plot looks like in R.
interaction.plot (growth_df$diet, growth_df$supplement, growth_df$gain, xlab="Diet", trace.label= "Supplement", ylab= "Weight Gain")
END THE EXERCISE BY SHOWING THEM HOW TO GET THE INTERACTION PLOT.
EVEN WHEN DOING THE INTERACTION PLOT IS NOT NEEDED IN THIS CASE (BC THE INTERACTION EFFECT IS NOT SIGNIFICANT)
Exercise 2 An engineer is designing a battery for a device that could be deployed at different temperatures. Three choices of material are available to build the battery. Three temperatures levels (15, 70, and 125 F) are considered because they are consistent with the product end-use environment. Four batteries are tested at each combination of material and temperature, and the resulting observed battery lifetime is recorded. Analyze the data using R and conclude about the influence of the two factors on the DV.
Download the CSV file “BatterriesTemperature” from Canvas and read it as a data frame in R.
To read the csv file, I WILL RECOMMEND YOU to run the following code:
YOU MUST MOVE THE CSV TO THE SAME FOLDER WHERE THIS NOTEBOOK IS !!!!
Batteries_temp_mat = read.csv ("BatteriesTemperature.csv",colClasses = c("numeric","factor","factor"))
Batteries_temp_matstr(Batteries_temp_mat)## 'data.frame': 36 obs. of 3 variables:
## $ lifetime : num 130 155 74 180 34 40 80 75 20 70 ...
## $ material : Factor w/ 3 levels "1","2","3": 1 1 1 1 1 1 1 1 1 1 ...
## $ temperature: Factor w/ 3 levels "125","15","70": 2 2 2 2 3 3 3 3 1 1 ...Before running ANOVA, answer a few questions:
- Outcome (dependent) variable: Battery lifetime
Factors and their levels:
Factor A: Type of material. It has three levels (mat 1, mat 2, & mat 3)
Factor B: Deployment temperature. It has three levels (15, 70, & 125 F)
9 treatments (3 materials * 3 temperatures)
Total observations = # treatments * # replicates
36= 9 * # replicates
Number replicates= 36/9= 4
- State the hypotheses for the possible three effects: For Material:
Ho: The mean battery lifetime is the same for all three materials (material 1, material 2, and material 3). Ha: There is a difference in mean battery lifetime for at least two materials.
Or, alternatively…
Ho: M mat1 = M mat2 = M mat3 Ha: Same as above
For Temperature:
Ho: … (Same ideas as Material)
Ha: … (Same ideas as Material)
For the Interaction between Temp and Material
Ho: No interaction effect (= The effect of material on battery lifetime does NOT depend on the temperature)
Ha: There is interaction effect (= The effect of material on battery lifetime DEPENDS on the temperature)
twoway_temp_mat= aov(lifetime ~ material * temperature, data= Batteries_temp_mat)
summary (twoway_temp_mat)## Df Sum Sq Mean Sq F value Pr(>F)
## material 2 10684 5342 7.911 0.00198 **
## temperature 2 39119 19559 28.968 1.91e-07 ***
## material:temperature 4 9614 2403 3.560 0.01861 *
## Residuals 27 18231 675
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Interpretation of the three effects: The p values for both factors, Material and Temperature, are less than alpha (0.05). Therefore, we can reject Ho in both cases and conclude: - At least two materials yield different mean lifetime - At least two temperatures yield different mean lifetime In this case, the interaction effect is also significant because the PV for the interaction is less than alpha. This suggests that the effect of the temperature on the lifetime DEPENDS on the type of material. In other words, the temperature affects the lifetime differently depending on the material.
How do we proceed in an exercise where the interaction effect is significant?
In this class, when the interaction is significant, you will ONLY be asked to do and interpret the interaction plot.
But, what do we do in a real life analysis when the interaction is significant? Read optional discussion #### Optional discussion
In a real life analysis, the analyst is usually trying to find the best combination between the factors (the combination that leads to the most desirable average for the DV).
When there is interaction, understanding the effect of the two factors on the DV becomes more challenging. It is NOT as simple as doing a post-hoc test individually for each factor and selecting the best combination.
The first thing that we usually do when there is interaction is the interaction plot. This plot gives us a visual of how one factor affects the impact of the other factor on the DV.
The second thing we could do is to run a more complex post-hoc test where we fix one level of factor A and then we compare the sample means for all the levels of factor B. The process is repeated for each level of factor A. You can do this more complex post-hoc test using a function in R called “emmeans()”.
End of optional discussion
#Change the order of the levels for the Temperature factor
Batteries_temp_mat$temperature = factor(Batteries_temp_mat$temperature, levels= c("15", "70", "125"))#plot
interaction.plot(Batteries_temp_mat$temperature, Batteries_temp_mat$material, Batteries_temp_mat$lifetime, xlab="temperature", trace.label= "material", ylab= "lifetime")
We can see in the plot that the effect of temperature on the battery
lifetime depends on the material (which confirms the existence of a
significant interaction between temperature and material). For example,
when we change from temp = 15 to 70, the lifetime clearly decreases if
materials 1 or 2 are used, but the lifetime barely changes if material 3
is used. Also, when we change from temp = 70 to 125, the lifetime
clearly decreases if materials 2 or 3 are used, but the lifetime barely
changes if material 1 is used.
Exercise 3
Consider a clinical trial in which three competing treatments for joint pain are compared in terms of their mean time to pain relief in patients with osteoarthritis.
Because investigators hypothesize that there may be a difference in time to pain relief in men versus women, they randomly assign 15 participating men to one of the three competing treatments and randomly assign 15 participating women to one of the three competing treatments.
Participants are instructed to take the assigned medication when they experience joint pain and to record the time, in minutes, until the pain subsides.
Run the following code chunks to create the data frame with the data:
time_to_relief= c (12, 15, 16, 17, 14, 14, 17, 19, 20, 17, 25, 27, 29, 24, 22, 21, 19, 18, 24, 25, 21, 20, 23, 27, 25, 37, 34, 36, 26, 29)
patient_sex= as.factor(rep(c("Male", "Female"), each= 15))
pain_treatment= as.factor(rep(c("A", "B", "C"), each= 5, times= 2))pain_treat_df= data.frame(time_to_relief, patient_sex, pain_treatment)
pain_treat_dfstr(pain_treat_df)## 'data.frame': 30 obs. of 3 variables:
## $ time_to_relief: num 12 15 16 17 14 14 17 19 20 17 ...
## $ patient_sex : Factor w/ 2 levels "Female","Male": 2 2 2 2 2 2 2 2 2 2 ...
## $ pain_treatment: Factor w/ 3 levels "A","B","C": 1 1 1 1 1 2 2 2 2 2 ...Outcome variable: Time to pain relief (in minutes)
Factor (s) and levels:
Factor 1: The pain treatment
Levels of factor 1: Treatment A, treatment B, and treatment C.
Factor 2: Patient’s sex
Levels of factor 2: Male and Female.
Number of treatments?
2 sexes * 3 treatments= 6 treatments (or level combinations)
Number of replicates (number of patients tested per treatment)?
Number of patients= # experimental treatments * replicates
30 = 6 * # replicates
replicates= 5 (patients tested per treatment)
Hypotheses:
For the pain treatment:
Ho: M treatment A = M treatment B = treatment C Ha: There is a difference in mean time to pain relief for at least 2 treatments
or…
Ho: All treatments have the same mean time to pain relief Ha: There is a difference in mean time to pain relief for at least 2 treatments
For the sex:
Ho: M Male = M Female Ha: There is a difference in mean time to pain relief for the sexes
or…
Ho: Both sexes have the same mean time to pain relief Ha: There is a difference in mean time to pain relief for the sexes
For the interaction between Pain Treatment and Sex:
Ho: There is NO interaction (The effect of the pain treatments on time to pain relief DOES NOT depend on the patient’s sex)
Ha: There is interaction (The effect of the pain treatments on time to pain relief depends on the patient’s sex. The pain treatments effect on time to pain relief changes according to the patient’s sex)
Run ANOVA and conclude:
twoway_pain_treat= aov(time_to_relief ~ patient_sex * pain_treatment, data= pain_treat_df)
summary (twoway_pain_treat)## Df Sum Sq Mean Sq F value Pr(>F)
## patient_sex 1 313.6 313.6 33.54 5.7e-06 ***
## pain_treatment 2 651.5 325.7 34.84 8.0e-08 ***
## patient_sex:pain_treatment 2 1.9 0.9 0.10 0.905
## Residuals 24 224.4 9.4
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Since the interaction is not significant but the effect of each factor is significant, we should run a post-hoc test individually for each factor.
Before running the post-hoc test, obtain the interaction plot to visually confirm that the interaction is not significant. THIS IS OPTIONAL. WHEN THERE IS NOT INTERACTION, YOU DO NOT NEED TO DO THIS PLOT! Obtaining this plot is actually not something that you should do when the interaction is NOT significant. I am asking you to do it just to practice.
interaction.plot (pain_treat_df$pain_treatment, pain_treat_df$patient_sex, pain_treat_df$time_to_relief, xlab ="Treatment", trace.label= "Gender", ylab= "Time to relief")
Post hoc test for Treatment:
sort(tapply(pain_treat_df$time_to_relief, pain_treat_df$pain_treatment, mean))## A B C
## 18.1 20.3 28.9TukeyHSD(twoway_pain_treat, which = "pain_treatment")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = time_to_relief ~ patient_sex * pain_treatment, data = pain_treat_df)
##
## $pain_treatment
## diff lwr upr p adj
## B-A 2.2 -1.214989 5.614989 0.2613818
## C-A 10.8 7.385011 14.214989 0.0000001
## C-B 8.6 5.185011 12.014989 0.0000049Post hoc test for Gender:
sort(tapply(pain_treat_df$time_to_relief, pain_treat_df$patient_sex, mean))## Male Female
## 19.20000 25.66667TukeyHSD(twoway_pain_treat, which = "patient_sex")## Tukey multiple comparisons of means
## 95% family-wise confidence level
##
## Fit: aov(formula = time_to_relief ~ patient_sex * pain_treatment, data = pain_treat_df)
##
## $patient_sex
## diff lwr upr p adj
## Male-Female -6.466667 -8.771096 -4.162237 5.7e-06