Lab07 Factorial ANOVA and LSD
Danielle Spencer
15 November 2022
Instructions
The same instructions for completion and submission of work used in previous labs applies here. Refer to previous labs for the details.
Objectives
Evaluate factorials with interactions
Interpret the results of ANOVA with factorial treatment combinations
Give an example of a statistical interaction and its biological meaning.
Determine when to compare main effects or cell means and use LSD to perform the appropriate comparisons.
Part 1. Factorial ANOVA using R functions.
These data are fictitious, but they should resemble real plant responses. A crop sensitive to cold temperatures was grown with and without irrigation and exposed to different numbers of days with cold temperatures (0, 5 and 10 days). All combinations of irrigation and cold exposure were observed, and each was replicated in two plots for a total of 12 experimental units. Irrigation is a factor with two levels and cold exposure is a factor with 3 levels. Each combination of irrigation and cold exposure is a treatment. Thus, there are 6 treatments in this experiment, and each one is replicated twice in a completely randomized design.
1A. Read and plot data. [5 points]
Treatments were assigned to plots by shuffling the order of plot numbers and then adding it to the table with the list of treatments and reps.
# Read and understand the code in this part but do not type or modify anything.
# Read data (yield.dat for "factorial")
yield.dat <- read.table(header = TRUE, text = "
irrig cold cold.days rep yield
irrig 00days 0 1 8.5
irrig 00days 0 2 9.5
irrig 05days 5 1 8.0
irrig 05days 5 2 9.0
irrig 10days 10 1 6.0
irrig 10days 10 2 7.0
no.irrig 00days 0 1 7.0
no.irrig 00days 0 2 8.0
no.irrig 05days 5 1 3.0
no.irrig 05days 5 2 2.0
no.irrig 10days 10 1 0.0
no.irrig 10days 10 2 1.0
")
# Assignment of numbers randomly to plots, done prior to measuring yields. This is need to do some calculations as described below.
set.seed(123)
yield.dat$plot.id <- sample(x = 1:12,
size = 12,
replace = FALSE)The two factor columns are sufficient to perform the analysis using preexisting R functions for factorial treatment combinations, but to do the calculations using only “basic” R functions (we refer to this as “by hand”) we need a new column that contains labels for all 12 treatments. We create the column by pasting together the names of the levels of both factors.
After adding the treatment column, we make a scatter plot of the data using the scatterplot function of the car package. This function requires a formula specifying the variables or columns that are used for the X and Y axes. The default is to add a curve and a straight line to the graph, so we need to specify the arguments smooth and regLine to be FALSE to prevent the addition of lines.
The graph shows that the response to number of cold days seems to be different between the irrig and no.irrig levels of the irrigation factor. Therefore, there appears to be an interaction between cold days and irrigation on yield. This type of graph is essential to be able to interpret interactions.
# Add treatment column combining columns "irrig" and "cold".
yield.dat$treat <- factor(paste(yield.dat$irrig, yield.dat$cold, sep = "."))
table(yield.dat$irrig, yield.dat$cold) #makes a simple table examining obs in the two treatments##
## 00days 05days 10days
## irrig 2 2 2
## no.irrig 2 2 2# Make a plot to inspect data
scatterplot(yield ~ cold.days | irrig, # Asks to make on plot per irrigation level
data = yield.dat, # Specifies the data to be used
smooth = FALSE, # Prevents the addition of a curve to plot
regLine = FALSE) # Prevent the addition of a line to plot
1.B Perform a test of the Ho that all treatment means are equal. [25 points]
For this part you need to test the null hypotheses that: - there are no interactions between irrigation and cold effects, - there are no “main effects” of irrigation or cold.
“Main effects” refers to the overall effect of a factor averaged over all levels of other factors. “Simple effects” are the effects of one factor for a given individual level of another factor. For example, the main effect of irrigation is the difference between the average of all plots with irrigation minus the average of all plots without irrigation, regardless of cold exposure.
We do these tests with an ANOVA that partitions all variation (sum of squares) of the response variable (yield) into Treatment and Error, and then, we partition the Treatment variation into main effect of irrig, main effect of cold, and interaction between irrig and cold (irrig x cold).
# First do the ANOVA ignoring the treatment factorial structure.
# This yields the first partition of SS between treatments and error.
m1 <- lm(formula = yield ~ treat, data = yield.dat) #This is a cell means model
anova(m1)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 5 119.75 23.95 47.9 9.318e-05 ***
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1# Now, analyze the same data, this time partitioning the Treatment SS according to a factorial treatment structure.
# In an R formula, irrig * cold indicates that the model includes the effect of irrig, the effect of cold and the interaction irrig by cold (irrig:cold)
m2 <- lm(formula = yield ~ irrig * cold, data = yield.dat) #This is the full model
# COMPLETE CODE BELOW to make an ANOVA table for the same data partitioning the Treatment SS according to factorial, as above:
anova(m2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## irrig 1 60.75 60.75 121.5 3.316e-05 ***
## cold 2 45.50 22.75 45.5 0.0002367 ***
## irrig:cold 2 13.50 6.75 13.5 0.0060105 **
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Questions:
- Report MSEs for both models, m1 and m2.
Your answer here: MSE m1=0.50 MSE m2=0.50
- State the null hypothesis for the main effects of irrigation in m2.
Your answer here: Irrigation has no effect on yeild.
- State the null hypothesis for the main effects of cold in m2.
Your answer here: Cold exposure has no effect on yeild.
- Do you reject or fail to reject the following null hypothesis?
Ho: There is no interaction between irrigation and cold exposure affecting yield
Your answer here: Reject null hypothesis because of very small p-value.
Part 2. Calculate SS by hand [40 points].
Type an equation with numbers or using R functions to calculate the total sum of squares.
R can be used as a calculator. Simply type the numbers and the operators. Review order of operations if necessary. Your calculations have to be accurate within 0.5% to be graded as correct. You can also use R functions if you wish.
For example, the average yield with irrigation is:
(8.5 + 9.5 + 8.0 + 9.0 + 6.0 + 7.0)/6## [1] 8mean(yield.dat$yield[yield.dat$irrig =="irrig"])## [1] 8And the average yield without irrigation is:
(7.0 + 8.0 + 3.0 + 2.0 + 0.0 + 1.0)/6## [1] 3.5mean(yield.dat$yield[yield.dat$irrig =="no.irrig"])## [1] 3.5Question: Enter the calculations for the averages of yield: - overall average = 5.75 - average for 00 days cold = 8.25 - average for 05 days cold = 5.5 - average for 10 days cold = 3.5
# overall average
(8.5+9.5+8.0+9.0+6.0+7.0+7.0+8.0+3.0+2.0+0.0+1.0)/12 # basic calculator approach## [1] 5.75mean(yield.dat$yield) # fast way## [1] 5.75# average for 00 days cold
(8.5+9.5+7.0+8.0)/4## [1] 8.25mean(yield.dat$yield[yield.dat$cold =="00days"])## [1] 8.25# average for 05 days cold
mean(yield.dat$yield[yield.dat$cold =="05days"])## [1] 5.5# average for 10 days cold
mean(yield.dat$yield[yield.dat$cold =="10days"])## [1] 3.5# average for all combination of irrigation and cold
aggregate(yield ~ treat, yield.dat, mean)## treat yield
## 1 irrig.00days 9.0
## 2 irrig.05days 8.5
## 3 irrig.10days 6.5
## 4 no.irrig.00days 7.5
## 5 no.irrig.05days 2.5
## 6 no.irrig.10days 0.5#aggregate(yield ~ treat, yield.dat, FUN = "mean") If the above aggregate() function doesn't work, try this one.Total sum of squares (TSS) can be calculated as follows:
# TSS
TSS<-(8.5 - 5.75)^2 +
(9.5 - 5.75)^2 +
(8.0 - 5.75)^2 +
(9.0 - 5.75)^2 +
(6.0 - 5.75)^2 +
(7.0 - 5.75)^2 +
(7.0 - 5.75)^2 +
(8.0 - 5.75)^2 +
(3.0 - 5.75)^2 +
(2.0 - 5.75)^2 +
(0.0 - 5.75)^2 +
(1.0 - 5.75)^2Effects of irrigation can be calculated as follows:
# effect of irrig
8 - 5.75 # repeats for each line of data with irrig (6 times)## [1] 2.25# effect of no irrig
3.5 - 5.75 # repeats for each line of data with no.irrig (6 times)## [1] -2.25The SS for the irrigation factor (SSirrig) can be calculated as:
# SS irrigation
SSI<-(8 - 5.75)^2 +
(8 - 5.75)^2 +
(8 - 5.75)^2 +
(8 - 5.75)^2 +
(8 - 5.75)^2 +
(8 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2
# check the value with the anova table
anova(m2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## irrig 1 60.75 60.75 121.5 3.316e-05 ***
## cold 2 45.50 22.75 45.5 0.0002367 ***
## irrig:cold 2 13.50 6.75 13.5 0.0060105 **
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1See how the effect of irrig and no.irrig repeat for each line with the corresponding level of irrigation factor.
**Question: follow the example above to calculate the effects of cold levels and the sum of squares for the cold factor (SScold).
# effect of 00days cold
8.25 - 5.75 #repeats 4 times## [1] 2.5# effect of 05days cold
5.5 - 5.75 #repeats 4 times## [1] -0.25# effect of 10days cold
3.5 - 5.75 #repeats 4 times## [1] -2.25aggregate(yield ~ cold, yield.dat, mean)[,2] - mean(yield.dat$yield) # a fast way to get the same thing## [1] 2.50 -0.25 -2.25# SScold
SSC<-(8.25 - 5.75)^2 +
(8.25 - 5.75)^2 +
(8.25 - 5.75)^2 +
(8.25 - 5.75)^2 +
(5.5 - 5.75)^2 +
(5.5 - 5.75)^2 +
(5.5 - 5.75)^2 +
(5.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2 +
(3.5 - 5.75)^2
# check the value with the anova table
anova (m2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## irrig 1 60.75 60.75 121.5 3.316e-05 ***
## cold 2 45.50 22.75 45.5 0.0002367 ***
## irrig:cold 2 13.50 6.75 13.5 0.0060105 **
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Following the examples above, calculate the Treatment SS (SST) where each combination of irrigation and cold is a treatment (you calculated the averages in the question about “Enter the calculations for the averages of yield” above)
# SS treatment
SST<-(9 - 5.75)^2 +
(9 - 5.75)^2 +
(8.5 - 5.75)^2 +
(8.5 - 5.75)^2 +
(6.5 - 5.75)^2 +
(6.5 - 5.75)^2 +
(7.5 - 5.75)^2 +
(7.5 -5.75)^2 +
(2.5 - 5.75)^2 +
(2.5 - 5.75)^2 +
(0.5 - 5.75)^2 +
(0.5 - 5.75)^2
# check value with the m1 anova table, which give SS treatment
anova(m1)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 5 119.75 23.95 47.9 9.318e-05 ***
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Calculate the sum of squares for the interaction between cold days and irrigation (SSCI) using sum of squares for treatment, sum of squares for irrigation and sum of squares for cold
SSCI<- SST - SSI - SSC
# check value with anova table
anova(m2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## irrig 1 60.75 60.75 121.5 3.316e-05 ***
## cold 2 45.50 22.75 45.5 0.0002367 ***
## irrig:cold 2 13.50 6.75 13.5 0.0060105 **
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Calculate the SSE using the total sum of squares and the sum of squares of the treatment
SSE<- TSS - SST
# check value with anova table
anova(m2)## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## irrig 1 60.75 60.75 121.5 3.316e-05 ***
## cold 2 45.50 22.75 45.5 0.0002367 ***
## irrig:cold 2 13.50 6.75 13.5 0.0060105 **
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1anova(m1) # both models have the same SSE## Analysis of Variance Table
##
## Response: yield
## Df Sum Sq Mean Sq F value Pr(>F)
## treat 5 119.75 23.95 47.9 9.318e-05 ***
## Residuals 6 3.00 0.50
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Part 3. Multiple comparison of means using Fisher’s PLSD [30 points]
Fisher’s Protected Least Significant Difference is the LSD that is applied only if the overall F test is significant. In this section you need to calculate the Least Significant Differences for irrig, cold and irrig x cold levels.
Estimated means (i.e., averages) that differ by less than the LSD are not significantly different.
We do the LSD for both factors separately and for the interaction to show the procedure. However, when there is a significant interaction as in this example, the meaning of testing for main effects is not uniquely defined. Do not test the main effects when there are significant interactions because the main effects may not be meaningful. When there are significant interactions, compare the means of each combination of factors.
To do comparisons by hand we need to get the critical t and the SE for each average. Critical t is the same for all comparisons because it depends only on alpha and dfe, which are constant for the whole experiment.
Standard errors depend on the number of observations used to calculate each estimate.
MSE <- 0.50 # enter the MSE from the ANOVA table
dfe <- 6 # enter the dfe from the ANOVA table
#start out by determining how many observations in each factor (cold and irrig) and cell means (treat). That way you can use the correct "n" to divide the MSE when calculating the standard for irrig, cold, and treat.
(table(yield.dat$irrig))##
## irrig no.irrig
## 6 6(table(yield.dat$cold))##
## 00days 05days 10days
## 4 4 4(table(yield.dat$treat))##
## irrig.00days irrig.05days irrig.10days no.irrig.00days no.irrig.05days
## 2 2 2 2 2
## no.irrig.10days
## 2alpha = 0.05
# Firstly, calculate the LSD's and CI's "by hand."
# Hint: check your answers of LSD with the results below. Both results should be the same if your calculations are correct.
(t.critical <- qt(p = 1 - alpha/2, df = dfe))## [1] 2.446912# irrig
(se.dif.irrig = sqrt(2*MSE/6)) #Each irrigation average involves 6 observations## [1] 0.4082483(LSD.irrig <- t.critical * se.dif.irrig)## [1] 0.9989476# cold ### COMPLETE CODE BELOW ###
(se.dif.cold <- sqrt(2*MSE/4) )# Each cold average involves 4 observations## [1] 0.5# Type your equation and calculate the LSD for cold.
### COMPLETE CODE BELOW ###
(LSD.cold <- t.critical * se.dif.cold )## [1] 1.223456# irrig x cold ### COMPLETE CODE BELOW ###
(se.dif.irrig.cold <- sqrt(2*MSE/2) ) # Each irrig.cold average involves 2 observations## [1] 0.7071068# Type your equation and calculate the LSD for irrig * cold
### COMPLETE CODE BELOW ###
(LSD.irrig.cold <- t.critical * se.dif.irrig.cold )## [1] 1.730228# we can also use the LSD.test function from the agricolae package to get a nice summary of all possible pairwise comparisons with letters to identify nsd
LSD.test(m1, "treat")$group ## yield groups
## irrig.00days 9.0 a
## irrig.05days 8.5 a
## no.irrig.00days 7.5 ab
## irrig.10days 6.5 b
## no.irrig.05days 2.5 c
## no.irrig.10days 0.5 dQuestion: Use the calculated LSD.irrig.cold and answer the following questions: yield groups irrig.00days 9.0 a irrig.05days 8.5 a no.irrig.00days 7.5 ab irrig.10days 6.5 b no.irrig.05days 2.5 c no.irrig.10days 0.5 d LSD=1.73
- Do treatments irrig:00days and irrig:10days differ?
Your answer here: Yes because 9.0-6.5=2.5 which is greater than the LSD indicating a significant difference. Also, the different letters in the r output signify a differnce.
- Do treatments irrig:00days and no.irrig:00days differ?
Your answer here: no because 9.0-7.5=1.5 which is smaller than the LSD indicating no signoficant difference. Also, sharing a letter in the r code output indicates no significant difference.
- Do treatments no.irrig:00days and no.irrig:10days differ?
Your answer here: Yes because 7.5-0.5=7 which is greater than the LSD indicating a significant difference.
- Explain in your own words how irrigation affects the effect of cold days. Hint: Look at the plot generated earlier in addition to the LSD table for irrig:cold.
Your answer here: Irrigation postively affects the yield on cold days. With greater days of cold exposure irrigation positively affects yield. The yield is higher as days increase when irrigation is done. Cold exposure reduces the yield. When irrigation is done during cold exposure it helps increase the yield.