R Notebook
Fahad Taimur
Lawson, John. Design and Analysis of Experiments with R (Chapman & Hall/CRC Texts in Statistical Science) (p. 52). CRC Press. Kindle Edition.
Load the libraries
library(daewr)
library(mvtnorm)
library(survival)
library(multcomp)
library(dplyr)
library(agricolae)- In Section 2.8.2 an experiment for determining the effect of the download site selected upon the time to download a file was discussed. In this experiment:
- Describe the experimental unit.
# The individuals using the sites to download a file are the experimental unit. - Describe the treatment factor. # The download site used is the treatment factor (A,B,C,D, etcetera).
- Describe the response. # The time it takes to download a file.
- Discuss the causes for experimental error in this experiment and why the principles of replication and randomization would be important in reaching a valid conclusion. # The causes of experimental error can include various background variables e.g wifi connection/internet speed, computer performance, web browsers being used and etcetera. There could also be other sources of experimental error such as the precision of the timer being used, human error and others. # We can utilize replication in our study by assigning multiple experimental units to a treatment factor. Using replication, we can ensure that the group means of the time to download a file are in fact due to the download site used, and not simply a random result of the experimental unit. # Using randomization, we can ensure that bias is minimized. By randomly assigning different experimental units to different sites being tested, we can accomplish that. Considering the various sources of variance and background variables, randomization allows us to perform error control and bias reduction. It allows us to conclude whether that the response is due to the treatment factor, and not due to some inherent bias in the experimental unit.
- In an experiment to study the effect of the amount of baking powder in a biscuit dough upon the rise heights of the biscuits, four levels of baking powder were tested and four replicate biscuits were made with each level in a random order. The results are shown in the table below.
.25 tsp .5 tsp .75 tsp 1 tsp
11.4 27.8 47.6 61.6
11.0 29.2 47.0 62.4
11.3 26.8 47.3 63.0
9.5 26.0 45.5 63.9
- What is the experimental unit? # Biscuit is the experimental unit.
- Perform the analysis of variance to test the hypothesis of no treatment effect.
# Create the dataset
id <- rep(c("A", "B", "C", "D"), each = 4)
fac <- factor(rep(c(0.25, 0.5, 0.75, 1), each = 4))
fac## [1] 0.25 0.25 0.25 0.25 0.5 0.5 0.5 0.5 0.75 0.75 0.75 0.75 1 1
## [15] 1 1
## Levels: 0.25 0.5 0.75 1ht <- c(11.4, 11.0, 11.3, 9.5,
27.8, 29.2, 26.8, 26.0,
47.6, 47.0, 47.3, 45.5,
61.6, 62.4, 63.0, 63.9)
experiment <- data.frame(id, treatment = fac, response = ht)ANOVA to test hypothesis
Ho: u1 = u2 = u3 = u4 Ha: Atleast one of the means is different # Based on the p-value, we can reject the null hypothesis and state that at least one of the treatment is different than the rest. At the chosen level of significance i.e 0.05, there are significant differences in rise heights based on applied baking powder factors.
mod <- aov(response ~ treatment, data = experiment)
summary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 3 6146 2048.6 1823 3.23e-16 ***
## Residuals 12 13 1.1
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Extra
- Formulate a contrast to test the hypothesis that increase in rise height is a linear function of the increase in baking powder in the dough, and test this hypothesis. # The results show there is a significant linear trend. The quadratic trend is insignificant.
contrasts(experiment$treatment) <- contr.poly(4)
contrasts(experiment$treatment)## .L .Q .C
## 0.25 -0.6708204 0.5 -0.2236068
## 0.5 -0.2236068 -0.5 0.6708204
## 0.75 0.2236068 -0.5 -0.6708204
## 1 0.6708204 0.5 0.2236068summary.lm(aov(response ~ treatment, data = experiment))##
## Call:
## aov(formula = response ~ treatment, data = experiment)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.4500 -0.7688 0.2375 0.5250 1.7500
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.9562 0.2650 139.435 <2e-16 ***
## treatment.L 39.1703 0.5301 73.894 <2e-16 ***
## treatment.Q -0.3875 0.5301 -0.731 0.4788
## treatment.C -1.4031 0.5301 -2.647 0.0213 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.06 on 12 degrees of freedom
## Multiple R-squared: 0.9978, Adjusted R-squared: 0.9973
## F-statistic: 1823 on 3 and 12 DF, p-value: 3.231e-16- Estimate the variance of the experimental error σ2.
# To measure experimental error, which is the difference between observed response for a particular experiment and the long run average of all experiments conducted at the same setting, we mutate the dataset and calculate the experimental error for each group.
experiment %>%
select(id, response) %>%
group_by(id) %>%
mutate(mean_response = mean(response),
exp_error = response - mean_response) -> var_df
var_df# Variance by treatment factor
# Take note of the variance by treatment factor (specially for treatment B = 0.5 tsp). This is a cause for concern.
var_df %>%
select(id, exp_error) %>%
group_by(id) %>%
summarise(exp_error_var = var(exp_error))# Total variance
var(mod$residuals)## [1] 0.8991667- Make a plot of residuals versus predicted values and normal plot of residuals and comment on whether the assumptions of the linear model are justified.
# The residuals vs. predicted plot shows that variance changes which violates the linear model assumptions. # The probability plot shows that residuals stray from the straight line of the theoretical quantiles which is another violation. # Conclusion: Considering the small number of observations, this violation should be taken with a grain of salt. A larger number of replicates might actually lead to a different conclusion. However, in this case the assumptions cannot be justified. Transforming the data might be useful in this case.
par(mfrow = c(1,2))
plot(mod, which = 1)
plot(mod, which = 2)
- If the dough were made in batches and the four replicate biscuit rise heights in each column (shown in the table above) were all from the same batch, would your answer to (a) be different? How could the data be analyzed if this were the case?
In such a case, we would be treating the each of the 4 batches as the experimental units. The results would be different as we would average the response of the experimental unit. It would be difficult to conclude whether the response observed is due to the effect of the applied factor, or something inherent within the dough. Simple analysis might be performed to make some loose conclusions, however it is an improper way to determine cause and effect.
- The effect of plant growth regulators and spear bud scales on spear elongation in asparagus was investigated by Yang-Gyu and Woolley (2006). Elongation rate of spears is an important factor determining final yield of asparagus in many temperate climatic conditions. Spears were harvested from 6-year-old Jersey Giant asparagus plants grown in a commercial planting at Bulls (latitude 40.2S, longitude 175.4E), New Zealand. Spears were harvested randomly and transported from field to lab for investigation. After trimming to 80mm length, spears were immersed completely for 1 h in aqueous solutions of 10 mg l-1 concentration of indole-3-acetic acid (IAA), abscisic acid (ABA), GA3, or CPPU (Sitofex EC 2.0%; SKW, Trostberg, Germany) in test tubes. Control spears were submerged in distilled water for 1 h. The experiment was a completely randomized design with five replications (spears) per treatment. The resulting data (final spear length in mm) is shown below.
Control IAA ABA GA3 CPPU
94.7 89.9 96.8 99.1 104.4
96.1 94.0 87.8 95.3 98.9
86.5 99.1 89.1 94.6 98.9
98.5 92.8 91.1 93.1 106.5
94.9 99.4 89.4 95.7 104.8
- Perform the analysis of variance to test the hypothesis of no treatment effect.
# Create the dataset
fac <- factor(rep(c("Control", "IAA", "ABA", "GA3", "CPPU"), each = 5))
fac## [1] Control Control Control Control Control IAA IAA IAA
## [9] IAA IAA ABA ABA ABA ABA ABA GA3
## [17] GA3 GA3 GA3 GA3 CPPU CPPU CPPU CPPU
## [25] CPPU
## Levels: ABA Control CPPU GA3 IAAres <- c(94.7, 96.1, 86.5, 98.5, 94.9,
89.9, 94.0, 99.1, 92.8, 99.4,
96.8, 87.8, 89.1, 91.1, 89.4,
99.1, 95.3, 94.6, 93.1, 95.7,
104.4, 98.9, 98.9, 106.5, 104.8)
res## [1] 94.7 96.1 86.5 98.5 94.9 89.9 94.0 99.1 92.8 99.4 96.8
## [12] 87.8 89.1 91.1 89.4 99.1 95.3 94.6 93.1 95.7 104.4 98.9
## [23] 98.9 106.5 104.8experiment <- data.frame(treatment = fac, response = res)
experimentANOVA
Ho: u1 = u2 = u3 = u4 = u5 Ha: Atleast one of the means is different # We can reject the null hypothesis of no treatment effect # Based on the p-value, we can reject the null hypothesis and state that at least one of the treatment is different than the rest. At the chosen level of significance i.e 0.05, there are significant differences in the effect of plant growth regulators and spear bud scales on spear elongation in asparagus.
mod <- aov(response ~ treatment, data = experiment)
mod## Call:
## aov(formula = response ~ treatment, data = experiment)
##
## Terms:
## treatment Residuals
## Sum of Squares 377.4936 270.2680
## Deg. of Freedom 4 20
##
## Residual standard error: 3.676058
## Estimated effects may be unbalancedsummary(mod)## Df Sum Sq Mean Sq F value Pr(>F)
## treatment 4 377.5 94.37 6.984 0.00109 **
## Residuals 20 270.3 13.51
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Extra
- Use the Tukey method to test all pairwise comparisons of treatment means.
# Tukey methods allows us to do parwise comparisons. The diff column shows the difference in means. To test the hypothesis that the means are equal, we can utilize the p-value. CPPU-ABA, CPPU-Control, CPPU-IAA, CPPU-GA3 show a significant difference.
tukeys <- TukeyHSD(mod, ordered = T)
tukeys## Tukey multiple comparisons of means
## 95% family-wise confidence level
## factor levels have been ordered
##
## Fit: aov(formula = response ~ treatment, data = experiment)
##
## $treatment
## diff lwr upr p adj
## Control-ABA 3.30 -3.6571003 10.2571 0.6229509
## IAA-ABA 4.20 -2.7571003 11.1571 0.3973453
## GA3-ABA 4.72 -2.2371003 11.6771 0.2881487
## CPPU-ABA 11.86 4.9028997 18.8171 0.0004696
## IAA-Control 0.90 -6.0571003 7.8571 0.9948496
## GA3-Control 1.42 -5.5371003 8.3771 0.9717058
## CPPU-Control 8.56 1.6028997 15.5171 0.0114471
## GA3-IAA 0.52 -6.4371003 7.4771 0.9993935
## CPPU-IAA 7.66 0.7028997 14.6171 0.0265612
## CPPU-GA3 7.14 0.1828997 14.0971 0.0425233Extra
- Use the Dunnett procedure to compare all treatment group means to the control mean. # Running SNK.test shows that CPPU has significant differences as compared to the rest of the groups labelled as b.
dunnet <- SNK.test(mod, "treatment", alpha = 0.05)
print(dunnet)## $statistics
## MSerror Df Mean CV
## 13.5134 20 95.656 3.842997
##
## $parameters
## test name.t ntr alpha
## SNK treatment 5 0.05
##
## $snk
## Table CriticalRange
## 2 2.949998 4.849746
## 3 3.577935 5.882064
## 4 3.958293 6.507367
## 5 4.231857 6.957100
##
## $means
## response std r Min Max Q25 Q50 Q75
## ABA 90.84 3.533129 5 87.8 96.8 89.1 89.4 91.1
## Control 94.14 4.530784 5 86.5 98.5 94.7 94.9 96.1
## CPPU 102.70 3.557387 5 98.9 106.5 98.9 104.4 104.8
## GA3 95.56 2.213143 5 93.1 99.1 94.6 95.3 95.7
## IAA 95.04 4.123469 5 89.9 99.4 92.8 94.0 99.1
##
## $comparison
## NULL
##
## $groups
## response groups
## CPPU 102.70 a
## GA3 95.56 b
## IAA 95.04 b
## Control 94.14 b
## ABA 90.84 b
##
## attr(,"class")
## [1] "group"