Chapter 1 Treatment Structures and Comparisons
Copyright 2018 Emilio A. Laca
If you use this, please cite it.
1.1 Learning Objectives for Chapter
- Identify treatments and factor levels in a factorial experiment.
- State the null and alternative hypotheses of a factorial experiment.
- Define and compare simple and main effects in a factorial.
- Given a graph of factorial experiment results, determine if there are main effects and/or interaction effects.
- Calculate PLSD confidence intervals for a pairwise comparison of means and determine if two means are significantly different based on the CIs.
- Design factorial treatment combinations to answer a research question.
- Transform a real-world question into a null hypothesis about means.
- Create linear combinations of means that address practical real-world questions.
- Define statistical interaction and explain it graphically.
- Test hypotheses involving more than one and more than 2 means.
- Interpret result of ANOVA with factorial treatments.
- Interpret a graph of a factorial experiment, including description of L, Q and C trends.
- Given a set of treatments, formulate and test useful contrasts.
- Determine if a linear combination of means is a contrast.
- Create contrasts to estimate the value and variance of interactions.
- Estimate the variance of a linear combination of means.
- Make confidence intervals for linear combination of estimated means.
- Calculate and use the Least Significant Difference (LSD) to test for differences among treatments.
- Apply the Protected LSD and explain how it ameliorates the inflation of alpha?
- Create linear combinations of estimated treatment means.
- Create contrasts to answer research questions such as effects of treatments, linear and quadratic effects of continuous treatment variables.
1.2 Pairwise Comparisons
ANOVA performs one test for the whole experiment, and one for each effect or line item in the model, including treatments. The single F test for treatments is used to test the null hypothesis that there are no treatment effects. If the F-value is greater than the critical value “from table” we reject Ho and accept the alternative, but we still do not know which treatments or factor levels are significantly different from each other. Estimated means need to be compared to determine which ones are different.
1.2.1 Least Significant Difference
Means are compared by subtracting one from the other and testing if the difference is significantly different form zero.
\[H_0: \mu_1 = \mu_2 \\ H_0: \mu_1 - \mu_2 = 0\]
Because the means are not known we use the estimates. Estimates of means are statistics, and because they are calculated as linear functions of normally distributed random variables, they are themselves normally distributed random variables.
- Recall: Linear combinations of normal random variables are normal random variables.
- The mean of a linear combination is the linear combination of the means.
- The variance of the difference between two RV’s is the sum of the individual variances plus twice their covariance.
- Estimated means of balanced experiments are independent random variables (i.e. their covariance is zero).
The variance of an estimated mean depends on the number of observations used to estimate the mean (number of observations averaged) and on the variance of individual observations, which comes fro the random errors. The relationship between the variance of an average and the variance of individual observations is given by the PRIME DIRECTIVE of this course.
1.2.2 Experiment-wise and Family-wise Error Rates
When multiple means are compared sequentially, a new problem arises. To understand the solution, first we need to fully understand the problem. In testing hypotheses we are particularly interested in avoiding problems of “false positives” where two averages appear to be significantly different while they really are not. This is called “Type I” error. The nominal probability of making a Type I error in a test is \(\alpha\). This is referred to as the “test-wise” or “comparison-wise” error rate, because it is the error rate in each single comparison or test.
The probability of making at least one error Type I in a set of tests comparing treatments or testing other contrasts is called “family-wise” or “experiment-wise” error rate. Family-wise error rate tends to increase as the number of tests in the family increases, because an error in any of the tests is an error for the family. If tests are independent and test-wise error rate is kept constant at \(\alpha\), the number of errors in a family has a Binomial distribution. The probability of not making any errors in n comparisons, assuming independence, is \((1-\alpha)^n\). Thus, the probability of making at least one error is \(1-(1-\alpha)^n\).
As an example, consider an experiment where there are 6 treatments. The total number of pairwise comparisons is
\[\left( \begin{array}{c} 6 \\ 2 \end{array} \right) = \frac{6!}{2!(6-2)!} = 15\]
With 15 comparisons, the probability of making at least one error in the set is 0.5367087698. Obviously, this error is much greater than the nominal rate for each comparison. Several methods are available to make sure that the family-wise error rate remains close to the nominal of 5%.
1.3 Factorials
Factorials are experiments where treatments have a special structure resulting from combining levels of different factors. The treatment structure allows us to assess the effects of several treatment factors as well as the interaction between factors. The concept of interaction is fundamental in statistics and in life in general and it is treated in detail below. Let’s concentrate on the factorial structure for now.
For a factorial set of treatments we need at least two factors or variables, each of which will take at least two values. Suppose you are trying to invite someone to have a meal with you and you want to maximize the chances of being accepted. Two main factors come into play: the type of cuisine you propose and the meal. In your town you have Mexican, Thai, Chinese, American and Japanese restaurants. You can have lunch or dinner. Therefore, you have 5 x 2 options for the invitation. For another example, imagine that you want to determine the best combination of speed and tire pressure to maximize your fuel efficiency while commuting. The car manufacturer gives a recommendation for tire pressure, but that is not specifically adjusted for the roads in your commute, or for your typical range of speeds. You can choose several tire pressures around the manufacturers recommendation, and several speeds that are reasonable around your average. Most likely, the effect of tire pressure will depend on speed because the relative importance of tire resistance and bouncing depend on speed. Say you choose 32, 34 and 36 psi for tire pressure and 65, 70 and 75 mph for speed. You will have to test 3 x 3 combinations. These are called complete factorials because each level of each factor appears with each level of the other factors. Sometimes we say that the factors are crossed as in a table where rows are levels of one factor and columns are levels of the other factor (Table 1.1)
| 32 psi | 34 psi | 36 psi | |
|---|---|---|---|
| 65 mph | 65-32 | 65-34 | 65-36 |
| 70 mph | 70-32 | 70-34 | 70-36 |
| 75 mph | 75-32 | 75-34 | 75-36 |
1.3.1 Interactions
Statistically, a two way interaction (or simply an interaction) is present when the effect of one factor depends on the level of another. The combined effects of changes in two factors is larger (synergistic or positive interaction) or smaller (antagonistic or negative interaction) than the sum of individual effects obtained when each factors is changed at a time. Interactions are common and we are familiar with them in daily life, but we typically do not take the time to think of them as interactions. Take for example the effect of putting on a nice warm coat on your level of comfort. It is a desirable effect if it is cold and not so good if it is really hot. There is an interaction between clothing and weather on comfort.
Three variables are involved in a two-way interaction, two predictor variables or treatment factors and one response variable. We say that predictor variable 1 (clothing) interacts with predictor variable 2 (ambient temperature) to determine the value of the response variable (level of comfort). Typical examples of interactions are:
- environment by genotype interaction on phenotype;
- drug interaction of alcohol with opiods on alertness;
- nitrogen and water supply interaction on crop yield.
Interactions are scale dependent, which is a specific effect of the fact that interactions depend on the response variable. An interaction that is significant in one scale may become non-significant after a logarithmic transformation of the response variable.
Hamilton and Saito [@HamiltonSaito2014] describe an interaction between level of education and political affiliation on belief in anthropogenic climate change (1.1). The effect of increasing education on the proportion of people that believe that climate change is mostly anthropogenic depends on their political affiliation: while for Democrats belief increases with increasing education, for Republicans there is no effect and for Tea Party members it.
Figure 1.1: Example of an interaction graph where political party and level of education interact in determining the probability that a person believes that climate change is mostly anthropogenic. See citation in text.
Olson et al. [@OlsonEtAl1991] reported a typical genotype by environment on calf birth weight showing that different cattle breed have specific adaptation to environmental conditions. Brahman is a breed of Bos indicus that evolved in Asia and has been selected for tropical conditions, whereas Angus and Hereford are breeds that were selected for colder climates with higher quality forages. The effect of location on calf birth weight (kg) depends on the breed, therefore, there is an interaction.
Figure 1.2: Example of genotype by environment interaction on calf birth weight. British crosses (Angus-Hereford) perform better in Nebraska, whereas Brahman crosses do better in the semitropical climate of Florida. Vertical bars are standard errors of estimated means
Factorial treatment combinations can result from combining levels of more than two factors. There is no theoretical limit to the number of factors that can be included, but the number of treatments grows very fast with number of treatments and levels per factor. For example, a factorial of three factors, each with three levels has 27 different treatments.
When there are more than two factors, new levels of interactions appear in the partition of variation and each observation. With three factors, say A, B and C there are three 2-way interactions (AB, AC and BC) and one 3-way interaction (ABC). A 3-way interaction means that the two way interaction between two of the three factors is different in the two levels of the third factor.
Photosynthesis and plant growth are known to require nitrogen (mostly for enzymes, nucleic acids and structural proteins), light and carbon dioxide. Growth rate typically responds nonlinearly to each of these factors and is limited by all of them. Figure 1.3 shows a representation of a typical pattern of response.
Figure 1.3: Representation of effects of light intensity, nitrogen and carbon dioxide concentration on photosynthesis rate by a well watered plant. \(N_0\) and \(N_1\) represents high and low nitrogen fertilization. \(CO_{2_0}\) and \(CO_{2_1}\) are low and high concentrations of carbon dioxide. Effects of increasing \(CO_2\) concentration with high light intensity are labeled \(e_1\) and \(e_2\) for high and low N. Effects of increasing \(CO_2\) concentration with low light intensity are labeled \(e_3\) and \(e_4\) for low and high N.
To work our way up to the 3-way interaction we start with “effects” of increasing \(CO_2\) concentration. The effect of increasing \(CO_2\) concentration in with high light and N is calculated on the red curve in the top graph. The four curves are separated into different panels to facilitate the visualization of effects (1.4)
Figure 1.4: Same information as in previous figure, but presented in one separate set of axes for each combination of light intensity and N fertilizer addition. Top row has high N, and right colunm has high light intensity. The effect of increasing \(CO_2\) concentration is depiced by a vertical bar and labeled \(e_1\) to \(e_4\).
To calculate the effect of increasing \(CO_2\) concentration in with high light and N, we move vertically on \(CO_{2_0}\) and \(CO_{2_1}\) until we reach the red curve. At each of the two points reached on the red curve we move horizontally to the left until we reach the Y axis, were we read the photosynthesis rate (PR) for low and high \(CO_2\). The difference between the two PR’s is the effect of increasing \(CO_2\) for high light intensity and high level of N fertilization, and it is labeled \(e_1\). The process is repeated for each one of the graphs to get the effect of increasing \(CO_2\) in each of the other three combinations of light and N, which are labeled \(e_2\) through \(e_4\).
First, consider only the two graphs with high light intensity (right column of Figure 1.4. We have a large effect of \(CO_2\) with high N (\(e_1\), red vertical bar) and a smaller effect of \(CO_2\) with low N (\(e_2\), yellow vertical bar). Therefore, the effect of \(CO_2\) depends on the level of N, which fits the definition of a significant 2-way interaction. The value of this interaction is the difference in length between the red and the yellow bars or \(e_1 - e_2\), which is represented in Figure 1.5 with a green vertical bar.
Now consider only the two graphs with low light intensity (left column of Figure 1.4). By repeating the same procedure applied in the right column we get two new effects of increasing \(CO_2\) (\(e_3\) with high and \(e_4\) with low N) and their difference, which is the value of the 2-way interaction between \(CO_2\) and N with low light (dark red or maroon vertical bar).
Finally, just like the difference in effects of one factor between levels or a second one is the 2-way interaction, the difference between 2-way interactions between levels of a third factor is the 3-way interaction, which is represented by the dark pink vertical bar.
Figure 1.5: Two-way interactions between \(CO_2\) and nitrogen are calculated as differences of effects of \(CO_2\) between levels of N. There is a large difference in the high light condition and a smaller difference in the low light condifion. The difference between the 2-way \(CO_2\) x N interactions in the two light conditions constitute the 3-way light x \(CO_2\) x N interaction.
Four factors produce six 2-way interactions (AB, AC, AD, BC, BD, CD), four 3-way interactions (ABC, ABD, ACD, BCD) and one 4-way interaction (ABCD), and so on. Two-way interactions are common and easy to understand intuitively; 3-way interactions may be common and are a bit harder to understand intuitively; 4-way interactions are almost impossible to understand in general terms for most people. However, I would guess that in the “real world” there are all sorts of complex interactions. Consider for example that you study plant growth rate and you use a factorial of four factors, each with 2 levels: with or without water, with or without carbon dioxide, with or without light and at -20 and 20 C. The plant will only grow in one of the 16 treatments, the one that includes water, carbon dioxide, light and a positive temperature. That is a 4-way interaction.
1.3.2 Interactions and Main Effects
When a factorial has a significant interaction, the main effects may have little meaning, because they average over conditions that have different consequences for the factor in question. Take for example a typical (realistic but fictitious data) experiment in which a strong ecotype by environment is present. Two ecotypes (Montanus and Desertorum) are grown in two environments (desert and mountain) in what is called a reciprocal transplant experiment, where each ecotype is grown in its native environment and in the other’s environment, with replication. Mature size of the organisms (could be plants or animals or anything else that is alive and has a mature size) is measured and analyzed.
Figure 1.6: Interaction plot showing the effects of ecotpye and environment on mature size. Although there are obsious responses to the treatments, main effects show no response for either factor.
Each ecotype does better than the other when grown in its original environment (Panel A, Figure 1.6). However, on average over ecotypes and over environments, there are no differences in mature size (Panels B and C, Figure 1.6). The main effects hide the clear pattern of response exhibited by the individual treatment means. It would not make sense to say that either factor has no effect, however it is true that averaged over the levels of ecotype, environment had not effect, and vice versa. Thus, whenever interactions are present and significant, interpretation of results should start with an inspection of interaction graphs.
1.3.3 Model and Calculations in Factorials
Factorial treatments can be used with different experimental designs, because the treatment structure is independent of the experimental design that determines how treatments are assigned to experimental units. Consider an experiment in which two factors (say nitrogen fertilization and amount of irrigation) are used to create a factorial treatment structure by combining all levels of both factors. Each combination is a treatment, and we can (but some times we may not want to) ignore the treatment structure when assigning treatments to experimental units. We start with the simplest experimental design, the completely randomized design (CRD) where treatments are assigned to experimental units completely at random.
To make the equations general, we will call
nitrogen factor A with levels \(a_1, \dots, a_{k_a}\) and
irrigation will be factor B with levels \(b_1, \dots, b_{k_b}\).
For example, the third nitrogen fertilization level is called \(a_3\), and the first level of irrigation is \(b_1\), thus, the treatment with level 3 of nitrogen and level 1 of irrigation is called \(a_3b_1\). The symbol \(k_a\) represent the number of levels of factor A, and \(k_b\) is the number of levels of factor B. The statistical model for this experiment is:
\[\begin{equation} Y_{ijh} = \mu_{ij.} + \epsilon_{ijh} \\[20pt] = \mu_{...} + \tau_{ij} + \epsilon_{ijh} \\[20pt] = \mu_{...} + \alpha_i + \beta_j + \alpha\beta_{ij} + \epsilon_{ijh} \\[20pt] \epsilon_{ijh} \sim N(0,\sigma) \\[20pt] \text{where: } \quad \mu_{...} \quad \text{is the overall mean} \\[20pt] \mu_{ij.} = \mu_{...} + \alpha_i + \beta_j + \alpha\beta_{ij} \quad \text{is the mean for treatment} \ a_ib_j\\[20pt] \tau_{ij} = \alpha_i + \beta_j + \alpha\beta_{ij} \quad \text{is the effect of treatment} \ a_ib_j\\[20pt] \alpha_i \quad \text{is the effect of level } a_i \text{ of factor } A\\[20pt] \beta_j \quad \text{is the effect of level } b_j \text{ of factor } B \\[20pt] \alpha\beta_{ij} \quad \text{is the interaction between } a_i \text{ and } b_j \tag{1.1} \end{equation}\]
The model implies that each observation is partitioned into the five components present in the model: overall mean, effect of A, effect of B, interaction AB and random error. As usual, the model has a counterpart that uses estimates instead of the true but unknown parameter values. When the design is complete and balanced such that all treatments have the same number of replications,
\[\begin{equation} Y_{ijh} = \hat{\mu}_{ij.} + \hat{\epsilon}_{ijh} \\[20pt] = \hat{\mu}_{...} + \hat{\tau}_{ij} + \hat{\epsilon}_{ijh} \\[20pt] = \hat{\mu}_{...} + \hat{\alpha}_i + \hat{\beta}_j + \hat{\alpha\beta}_{ij} + e_{ijh} \\[20pt] \text{where: } \quad \hat{\mu}_{...} = \bar{Y}_{...} \quad \text{is the overall average} \\[20pt] \hat{\tau}_{ij} = \bar{Y}_{ij.} - \bar{Y}_{...} = \hat{\alpha}_i + \hat{\beta}_j + \hat{\alpha\beta}_{ij} \quad \text{is the estimated effect of treatment} \ a_ib_j\\[20pt] \hat{\alpha}_i = \bar{Y}_{i..} - \bar{Y}_{...} \quad \text{is the estimated effect of level } a_i \text{ of factor } A\\[20pt] \hat{\beta}_j = \bar{Y}_{.j.} - \bar{Y}_{...} \quad \text{is the estimated effect of level } b_j \text{ of factor } B \\[20pt] \hat{\alpha\beta}_{ij} = \bar{Y}_{ij.} - \bar{Y}_{i..} - \bar{Y}_{.j.} + \bar{Y}_{...} \quad \text{is the estimated interaction between } a_i \text{ and } b_j \tag{1.2} \end{equation}\]
If the experiment is laid out in a randomized complete block design, we simply add a term for blocks to the model (Equation (1.3)). Block effects are calculated as the difference of block averages and overall average, as usual
\[\begin{equation} Y_{ijh} = \mu_{ij.} + B_h + \epsilon_{ijh} \\[20pt] = \mu_{...} + \tau_{ij} + B_h + \epsilon_{ijh} \\[20pt] \text{where: } \quad B_h \quad \text{is the effect of block h} \\[20pt] B_h = \mu_{..j} - \mu_{...} \\[20pt] \text{which is estimated as} \\[20pt] \hat{B}_h = \bar{Y}_{..j} - \bar{Y}_{...} \\[20pt] \tag{1.3} \end{equation}\]
1.3.3.1 R Code for Factorial in RCBD
The equations and calculations are illustrated with a numerical example. We conducted an experiment to determine the competitive ability of two naturalized grasses of the California Annual Grassland (Laca, unpublished data), Lolium multiflorum (Lam.) and Bromus hordeaceus (L.) which we will call Lomu and Brho for short. The experiment was laid out in an RCBD with 10 blocks and 6 treatments. Blocks were used because the field had fertility gradients caused by differences in soil types. Treatment structure was a 2 x 3 factorial of target plant (the species of the plant measured, which was Lomu or Brho) and competitor (self, both and none), which was the species surrounding the target individual. “Self” competitor was a competitor of the same species as the target (intraspecific competition), “both” was a mix of both species, and “none” meant that the target plant was grown in an area that was kept clear of any other species within a radius of 10 cm of the target plant. All plots were circular and had a radius of 10 cm around the target plant (Figure 1.7). Areas around plots were kept clear of vegetation. Target plant mass and height were measured at maturity. Plant mass was transformed to meet the assumption of homogeneity of variance, but the range of values is comparable to the original units in mg per plant. We analyze plant mass as a function of blocks and the factorial set of treatments.
Goal and hypothesis The goal of the study (for this section) was to determine if intra and interspecific competition were different for Lomu and Brho. The null hypothesis is that the effect of adding competitors depends on whether the competitor is of the same species or not. For logistic reasons, instead of using treatments with only interspecific competitors, we used a mixture of both species. Therefore, the main reason for the experiment was to determine if there was a significant statistical interaction between target and competitor species. We also had more specific hypotheses about which species was the superior competitor and the relative effects of intra and interspecific competition. Those are addressed in the Section about Contrasts.
Figure 1.7: Treatments used in the experiment to study competition between ryegrass (Lomu) and softchess (Brho). Rows are for target plants measured and columns are for competitor. Plots were circles of 20 cm diameter established in a field.
|
|
First we analyze the data using high-level R functions for linear models. Then, we will recreate the calculations “by hand” using basic functions and applying the formulas directly derived from the model. Prior to looking at results we check for homogeneity of residuals and normality using graphics. Residuals should be distributed symmetrically about zero, and should not exhibit any patterns when plotted against the predicted values for each observation. The quantile-quantile scatterplot should be closely grouped around a straight line.
lb.m1 <- lm(mass.mg ~ block + species * competitor, data = lb.compDat)
plot(residuals(lb.m1) ~ fitted(lb.m1),
xlab = "Fitted values (Y hats)",
ylab = "Residuals (e's)",
pch = 16,
col = "blue"
)
abline(h = 0, lty = 2, col = "red")
qqp(residuals(lb.m1), pch = 16, col = "blue")## [1] 12 15
Figure 1.8: Residual plots for factorial of target and competitor in RCBD.
Observations 12 and 15 are a little out of the assumed distribution (Figure 1.8) because they have extremely large residuals compared to the rest of the observations. In a real analysis, we would explore the reasons why those observations have large residuals by double checking data entry and field notes, as well as integrity of lab samples. We note but ignore the outliers in the rest of the explanation, as the extreme residuals do not change the way calculations are done.
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| block | 9 | 147639 | 16404 | 3.08 | 0.00586 |
| species | 1 | 34225 | 34225 | 6.42 | 0.01486 |
| competitor | 2 | 365809 | 182905 | 34.29 | 0.00000 |
| species:competitor | 2 | 57388 | 28694 | 5.38 | 0.00804 |
| Residuals | 45 | 240016 | 5334 | NA | NA |
All effects in the model were significant. The line items in the ANOVA table (Table 1.3) correspond to specific hypotheses, and the fact that the probabilities, listed in the last column, are all smaller than 0.05 means that all null hypotheses are rejected. In general, hypotheses tested are that effects are all zero, but the origin of the hypotheses is a set of questions about intra and interspecific competition. Our main question whether intra and interspecific competition is the same in both species. This is convenient because it focuses first on the interaction. In general, the interactions should be tested first in factorials. If they are significant, then specific comparisons of cell means should be performed and interpreted, because marginal means may no be interpretable (see Section about Contrasts).
Interaction between species and competitor: Is the effect of competition the same in both species? If competitive effects are equal, the effect of adding competitors will be the same in all cases. This effect is the difference between the control with no competition and each of the treatments with competitors, within species. To say that competitive strengths differs between species is the same as stating that the effects of competition (factor B) differ between levels of species (factor A), which is the definition of interaction. So our first interesting hypothesis refers directly to the interaction.
\[\text{ Null hypothesis: competition effects are the same for Lomu and Brho } \\[15pt] H_0: \alpha\beta_{ij} = 0 \quad \text{for all pairs } \ i,j \\[15pt] \text{which is equivalent to: } \\[15pt] \begin{equation} H_0: \quad \begin{cases} \mu_{LomuNone} - \mu_{LomuSelf} = \mu_{BrhoNone} - \mu_{BrhoSelf} & \quad \text{and} \\ \mu_{LomuNone} - \mu_{LomuBoth} = \mu_{BrhoNone} - \mu_{BrhoBoth} & \quad \text{and} \\ \mu_{LomuBoth} - \mu_{LomuSelf} = \mu_{BrhoBoth} - \mu_{BrhoSelf} \end{cases} \\[45pt] \text{Alternative hypothesis: at least one of the three pairs of effects is not equal} \\[15pt] \end{equation}\]
The hypothesis about interactions is tested in the line of the ANOVA labeled species:competitor. The calculated F value is 5.38, which is rather large. The probability of observing a calculated F of such value or larger when \(H_0\) is true (called “p-value”) is 0.008, which is smaller than the \(\alpha = 0.05\) required for rejecting the null hypothesis. Thus, we reject \(H_0\) and accept \(H_a\), which means that competition effects differ between the species. Because the interaction is “significant,” we can proceed to more detailed comparisons of cell means to answer the question more specifically: which species is the superior competitor, and how does competition against both competitors differ from competing against self? Those specific hypotheses are addressed below (see Section about Contrasts).
Marginal competition effect: is there an effect of competition on mature mass? We pose the null hypothesis that leads to a strong inference by its rejection: on average over species, competition does not affect mature size. If there were no competition, then plant mass at maturity for each species should not be affected by the competitor treatment. Lomu and Brho may differ in mean size, but their size would not be affected by the competitor. Keep in mind that a significant interaction was detected, so the interpretation of marginal results may be misleading and have to be considered with caution. The line corresponding to competitor in the ANOVA table (Table 1.3) shows a p-value of 0, which is significant. The null hypothesis that competition has no effect, or that mature size is the same, regardless of competition, is rejected.
\[\text{ Null hypothesis: there is no competition effect on mature size } \\[15pt] H_0: \mu_{.None} = \mu_{.Self} = \mu_{.Both} \\[15pt] \text{Alternative hypothesis: at least one of the three pairs is not equal} \\[15pt] \begin{equation} H_a: \quad \begin{cases} \mu_{.None} \neq \mu_{.Self} & \quad \text{or} \\ \mu_{.None} \neq \mu_{.Both} & \quad \text{or} \\ \mu_{.Self} \neq \mu_{.Both} \end{cases} \end{equation}\]
Note that the hypotheses were posed as means or effects being equal or not equal. It would be equivalent to state that the differences between the terms on each side of the “=” sign are zero. For example:
\[H_0: \mu_{.None} = \mu_{.Self} = \mu_{.Both} \\[15pt]\] is equivalent to
\[H_0:\mu_{.None} - \mu_{.Self} = 0 \quad \text{AND} \quad \mu_{.None} - \mu_{.Both} = 0 \quad \text{AND} \quad \mu_{.Self} - \mu_{.Both} = 0 \]
Marginal species effect: do species differ in mature size, averaged over competition treatments? Because of the significant interaction detected, the same note about caution applies here. The interaction means that the difference in mature size between species depends on competition. Ignoring the interaction for now, we test the marginal effect of species. Based on the line for species in the ANOVA table, we reject the null hypothesis that mature plant mass is the same for both species.
\[\text{ Null hypothesis: there is no species effect on mature size } \\[15pt] H_0: \mu_{Lomu.} = \mu_{Brho.} \\[15pt] \text{Alternative hypothesis: mature size differs between species} \\[15pt] \begin{equation} H_a: \quad \mu_{Lomu.} \neq \mu_{Brho.} \end{equation}\]
Results can be presented using a table or a figure. A figure facilitates the interpretation of interactions. A table also includes standard errors and marginal means, which is more complete but harder to interpret quickly. Both are presented below for illustration.
Figure 1.9: Effects of species and competitor on mature plant size. Vertical bars are 95% confidence intervals. Plant mass is in transformed units, but values are commensurate with mass in mg.
cmm <- as.data.frame(emmeans(lb.m1, "competitor"))
smm <- as.data.frame(emmeans(lb.m1, "species"))
cspmm <- as.data.frame(emmeans(lb.m1, "species", "competitor"))| Competitor | ||||
|---|---|---|---|---|
| Species | Self | Both | None | Spp means |
| Brho | 275.2 | 186.8 | 335.8 | 265.9 |
| Lomu | 245.6 | 237.9 | 457.6 | 313.7 |
| Comp. means | 260.4 | 212.3 | 396.7 |
R Code Summary
The analysis of the factorial experiment, without the code for making the pretty figures and table was very simple. It only required 5 short lines of code:
# Fit model to data. lb.compDat is the data frame with data
lb.m1 <- lm(mass.mg ~ block + species * competitor, data = lb.compDat)
anova(lb.m1) # Perform ANOVA based on model## Analysis of Variance Table
##
## Response: mass.mg
## Df Sum Sq Mean Sq F value Pr(>F)
## block 9 147638.82 16404.313 3.07561 0.0058621 **
## species 1 34224.82 34224.817 6.41673 0.0148598 *
## competitor 2 365809.43 182904.717 34.29239 8.9656e-10 ***
## species:competitor 2 57388.23 28694.117 5.37980 0.0080369 **
## Residuals 45 240015.68 5333.682
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1plot(residuals(lb.m1) ~ fitted(lb.m1)) # Check homogeneity of variance
qqp(residuals(lb.m1)) # Check normality
## [1] 12 15# Get marginal means, S.E. and CI's for competitor
emmeans(lb.m1, "competitor")## competitor emmean SE df lower.CL upper.CL
## 1.Self 260 16.3 45 228 293
## 2.Both 212 16.3 45 179 245
## 3.None 397 16.3 45 364 430
##
## Results are averaged over the levels of: block, species
## Confidence level used: 0.95# Get marginal means, S.E. and CI's for species
emmeans(lb.m1, "species")## species emmean SE df lower.CL upper.CL
## Brho 266 13.3 45 239 293
## Lomu 314 13.3 45 287 341
##
## Results are averaged over the levels of: block, competitor
## Confidence level used: 0.95# Get cell means, S.E. and CI's for treatments
emmeans(lb.m1, "species", "competitor")## competitor = 1.Self:
## species emmean SE df lower.CL upper.CL
## Brho 275 23.1 45 229 322
## Lomu 246 23.1 45 199 292
##
## competitor = 2.Both:
## species emmean SE df lower.CL upper.CL
## Brho 187 23.1 45 140 233
## Lomu 238 23.1 45 191 284
##
## competitor = 3.None:
## species emmean SE df lower.CL upper.CL
## Brho 336 23.1 45 289 382
## Lomu 458 23.1 45 411 504
##
## Results are averaged over the levels of: block
## Confidence level used: 0.951.3.3.2 Detailed Calculations for Factorial in RCBD
In this section we show the basic calculations to obtain parameter estimates and complete the ANOVA table for a factorial in a Randomized Complete Block Design. First, we obtain the parameter estimates by using overall and marginal averages. Then we partition each observation using the parameters and calculate the corresponding sums of squares. Finally, degrees of freedom are computed and the ANOVA table is constructed. Treatment averages:
\[\hat{\mu}_{ij.} = \bar{Y}_{ij.} = \frac{1}{r}\sum^{h=r}_{h=1}Y{ijh} \\[25pt] \text{For example: } \\[25pt] \hat{\mu}_{Brho,1.Self} = \frac{203 + 268 + 413 + 346 + 136 + 212 + 205 + 242 + 401 + 326}{10} \\[25pt] = 275.2\]
We can calculate treatment or cell averages using the aggregate function in R. Averages are merged with the data set so we can calculate treatment effects and later partition them into main effects and interactions.
Writing code in R can be greatly facilitated by carefully planning the names of objects. Note how the names below are organized such that it is easy to modify the code for species to get the code for competitors. The first chunk for cell averages was copied and quickly edited to get the chunks for marginal and block averages.
(cell.avgs <- aggregate(mass.mg ~ species + competitor,
data = lb.compDat,
FUN = mean
))## species competitor mass.mg
## 1 Brho 1.Self 275.2
## 2 Lomu 1.Self 245.6
## 3 Brho 2.Both 186.8
## 4 Lomu 2.Both 237.9
## 5 Brho 3.None 335.8
## 6 Lomu 3.None 457.6names(cell.avgs)[3] <- "cell.avg"
lb.compDat <- merge(lb.compDat,
cell.avgs,
all = TRUE
)Marginal averages are calculated as the averages of the corresponding cell averages.
\[\hat{\mu}_{i..} = \bar{Y}_{i..} = \frac{1}{k_b}\sum^{j=k_b}_{j=1}\bar{Y}_{ij.} = \frac{1}{k_b}\sum^{j=k_b}_{j=1}\hat{\mu}_{ij.} \\[25pt] \text{where } k_b = 3 \text{ is the number of levels of factor B (competitor). For example: } \\[25pt] \hat{\mu}_{Brho} = \frac{275.2 + 186.8 + 342.9}{3} = 265.9333333333\]
(s.avgs <- aggregate(cell.avg ~ species, # get species averages
data = cell.avgs,
FUN = mean
))## species cell.avg
## 1 Brho 265.9333333
## 2 Lomu 313.7000000names(s.avgs)[2] <- "sp.avg" # change column name to facilitate merge
(comp.avgs <- aggregate(cell.avg ~ competitor, # get competitor averages
data = cell.avgs,
FUN = mean
))## competitor cell.avg
## 1 1.Self 260.40
## 2 2.Both 212.35
## 3 3.None 396.70names(comp.avgs)[2] <- "comp.avg" # change column name to facilitate merge\[\hat{\mu}_{.j.} = \bar{Y}_{.j.} = \frac{1}{k_a}\sum^{i=k_a}_{i=1}\bar{Y}_{ij.} = \frac{1}{k_a}\sum^{i=k_a}_{i=1}\hat{\mu}_{ij.} \\[25pt] \text{where } k_a = 2 \text{ is the number of levels of factor A (target species). For example: } \\[25pt] \hat{\mu}_{Self} = \frac{275.2 + 245.6}{2} = 260.4\\[30pt]\]
We follow an analogous procedure to estimate overall2 and block means with block averages. As an exercise, write down the equations to calculate block averages for this data set before read any further. The equations are shown further down (Equation (1.4)).
# merge data set with marginal species and competitor averages
lb.compDat <- merge(lb.compDat,
s.avgs,
all = TRUE
)
lb.compDat <- merge(lb.compDat,
comp.avgs,
all = TRUE
)
# calculate block averages and merge with data
b.avgs <- aggregate(mass.mg ~ block,
data = lb.compDat,
FUN = mean
)
names(b.avgs)[2] <- "block.avg" # change column name to facilitate merge
lb.compDat <- merge(lb.compDat,
b.avgs,
all = TRUE
)
print(lb.compDat[, -4], digits = 5) # data with averages## block competitor species mass.mg cell.avg sp.avg comp.avg block.avg
## 1 1 1.Self Brho 203 275.2 265.93 260.40 241.00
## 2 1 3.None Lomu 391 457.6 313.70 396.70 241.00
## 3 1 2.Both Brho 133 186.8 265.93 212.35 241.00
## 4 1 1.Self Lomu 205 245.6 313.70 260.40 241.00
## 5 1 2.Both Lomu 245 237.9 313.70 212.35 241.00
## 6 1 3.None Brho 269 335.8 265.93 396.70 241.00
## 7 10 1.Self Brho 205 275.2 265.93 260.40 246.67
## 8 10 2.Both Brho 243 186.8 265.93 212.35 246.67
## 9 10 2.Both Lomu 143 237.9 313.70 212.35 246.67
## 10 10 1.Self Lomu 166 245.6 313.70 260.40 246.67
## 11 10 3.None Lomu 410 457.6 313.70 396.70 246.67
## 12 10 3.None Brho 313 335.8 265.93 396.70 246.67
## 13 11 3.None Brho 345 335.8 265.93 396.70 273.67
## 14 11 2.Both Brho 129 186.8 265.93 212.35 273.67
## 15 11 3.None Lomu 563 457.6 313.70 396.70 273.67
## 16 11 1.Self Lomu 186 245.6 313.70 260.40 273.67
## 17 11 1.Self Brho 242 275.2 265.93 260.40 273.67
## 18 11 2.Both Lomu 177 237.9 313.70 212.35 273.67
## 19 13 3.None Brho 394 335.8 265.93 396.70 331.17
## 20 13 1.Self Brho 401 275.2 265.93 260.40 331.17
## 21 13 3.None Lomu 508 457.6 313.70 396.70 331.17
## 22 13 2.Both Lomu 320 237.9 313.70 212.35 331.17
## 23 13 1.Self Lomu 229 245.6 313.70 260.40 331.17
## 24 13 2.Both Brho 135 186.8 265.93 212.35 331.17
## 25 14 3.None Brho 395 335.8 265.93 396.70 305.83
## 26 14 2.Both Brho 162 186.8 265.93 212.35 305.83
## 27 14 3.None Lomu 540 457.6 313.70 396.70 305.83
## 28 14 2.Both Lomu 222 237.9 313.70 212.35 305.83
## 29 14 1.Self Brho 326 275.2 265.93 260.40 305.83
## 30 14 1.Self Lomu 190 245.6 313.70 260.40 305.83
## 31 2 3.None Brho 340 335.8 265.93 396.70 282.00
## 32 2 2.Both Brho 221 186.8 265.93 212.35 282.00
## 33 2 1.Self Brho 268 275.2 265.93 260.40 282.00
## 34 2 2.Both Lomu 226 237.9 313.70 212.35 282.00
## 35 2 3.None Lomu 268 457.6 313.70 396.70 282.00
## 36 2 1.Self Lomu 369 245.6 313.70 260.40 282.00
## 37 4 1.Self Brho 413 275.2 265.93 260.40 412.67
## 38 4 2.Both Brho 317 186.8 265.93 212.35 412.67
## 39 4 3.None Lomu 649 457.6 313.70 396.70 412.67
## 40 4 2.Both Lomu 453 237.9 313.70 212.35 412.67
## 41 4 3.None Brho 279 335.8 265.93 396.70 412.67
## 42 4 1.Self Lomu 365 245.6 313.70 260.40 412.67
## 43 5 1.Self Brho 346 275.2 265.93 260.40 300.17
## 44 5 3.None Brho 248 335.8 265.93 396.70 300.17
## 45 5 1.Self Lomu 240 245.6 313.70 260.40 300.17
## 46 5 2.Both Brho 287 186.8 265.93 212.35 300.17
## 47 5 2.Both Lomu 205 237.9 313.70 212.35 300.17
## 48 5 3.None Lomu 475 457.6 313.70 396.70 300.17
## 49 7 1.Self Brho 136 275.2 265.93 260.40 259.17
## 50 7 2.Both Lomu 181 237.9 313.70 212.35 259.17
## 51 7 3.None Brho 427 335.8 265.93 396.70 259.17
## 52 7 3.None Lomu 425 457.6 313.70 396.70 259.17
## 53 7 2.Both Brho 132 186.8 265.93 212.35 259.17
## 54 7 1.Self Lomu 254 245.6 313.70 260.40 259.17
## 55 8 3.None Brho 348 335.8 265.93 396.70 245.83
## 56 8 1.Self Brho 212 275.2 265.93 260.40 245.83
## 57 8 2.Both Lomu 207 237.9 313.70 212.35 245.83
## 58 8 1.Self Lomu 252 245.6 313.70 260.40 245.83
## 59 8 2.Both Brho 109 186.8 265.93 212.35 245.83
## 60 8 3.None Lomu 347 457.6 313.70 396.70 245.83\[\begin{equation} \hat{\mu}_{..h} = \bar{Y}_{..h} = \frac{1}{k_a k_b} \sum^{i=k_a}_{i=1}\sum^{j=k_b}_{j=1}Y_{ijh} \tag{1.4} \end{equation}\]
\[\begin{equation} \hat{\mu}_{...} = \frac{1}{k_a k_b} \sum^{i=k_a}_{i=1}\sum^{j=k_b}_{j=1}\hat{\mu}_{ij.} \tag{1.5} \end{equation}\]
Sums of Squares Now that we have all averages we can calculate sums of squares. Sums of squares (SS) correspond to the effects in the model. Sums of squares are
TSS or total sum of squares that is the sum of squares of total differences between observed value (\(Y_{ijh}\)) and overall average (\(\hat{\mu}_{...}\)).
SST or treatment sum of squares is the sum of squared differences between treatment (cell) averages\(\hat{\mu}_{ij.}\) and overall average (\(\hat{\mu}_{...}\)).
SSB or block sum of squares is the sum of squared differences between block average (\(\hat{\mu}_{..h}\)) and overall average (\(\hat{\mu}_{...}\)).
SSE or residual (error) sum of squares of differences between observed value (\(Y_{ijh}\)) and predicted value (\(\hat{Y}_{ijh}\)), where predicted value is the sum of overall average and all effects.
Because the treatments have a factorial structure and each treatment effect \(\tau_{ij}\) is partitioned into its components \(\alpha_i\), \(\beta_j\) and \(\alpha\beta_{ij}\).Accordingly, SST is partitioned into SSA, SSB and SSAxB.
Total sum of squares
We use R to implement the equations that define each of the effects and corresponding sum of squares. First we calculate the overall average and use it to compute TSS. In the code below all.avg = \(\mu_{...}\) and the column lb.compDat$mass.mg is the set of all \(Y_{ijh}\). Thus, sum((lb.compDat\$mass.mg - all.avg) ^ 2)= \(\sum_{i=1}^{k_a}\sum_{j=1}^{k_b}\sum_{h=1}^{r}(Y_{ijh} - \hat{\mu}_{...})^2\)
(all.avg <- mean(cell.avgs$cell.avg)) # get overall average## [1] 289.8166667# Total sum of squares
(TSS <- sum((lb.compDat$mass.mg - all.avg)^2))## [1] 845076.9833Sum of squares of Blocks
Now we compute SSB. In the code below, lb.compDat\$block.avg is \(\hat{\mu}_{..h}\), thus, sum((lb.compDat\$block.avg - all.avg) ^ 2) = \(\sum_{i=1}^{k_a}\sum_{j=1}^{k_b}\sum_{h=1}^{r}(\hat{\mu}_{..h} - \hat{\mu}_{...})^2 = k_a k_b \sum_{h=1}^{r}(\hat{\mu}_{..h} - \hat{\mu}_{...})^2\)
# Sum of squares of blocks
(SSB <- sum((lb.compDat$block.avg - all.avg)^2))## [1] 147638.8167Sum of squares of Treatments
For the SST we use the column lb.compDat$cell.avg, which contains the estimated mean for each combination of target species and competitor assigned to each observation.
Keep in mind that each observation is partitioned into an estimate of the overall mean plus all effects: block, factor A, factor B, interactions AxB and residual. All sums for SS are over all observations, although some estimates are repeated. For example, the estimated effect of Brho is the same for all plots where Brho was the target species.
1.3.4 Advantages of factorials
Consider that you have to study the effects of irrigation and fertilization with N on corn yield. You could take a simple approach and do one experiment for irrigation and another for nitrogen, or you could study both factor at the same time in a factorial experiment. Note that “factorial” refers to the way treatments are organized and it says nothing about the experimental design in terms of how treatments are assigned to experimental units.
Here we demonstrate two main advantages of factorials:
- ability to detect interactions between factors, and
- much greater efficiency than separate experiments.
To show the advantage of a factorial, we first simulate two separate experiments and then a single combined experiment with both factors. By comparing the results we will be able to see that the factorial approach has advantages.
1.3.4.1 Effects of Irrigation alone
Using the simple approach, the experiment for irrigation includes two treatments: one or two irrigations. For this experiment, a moderate level of nitrogen of 160 lb/ac is selected. Treatments are repeated three times in each of two blocks (Table 1.5). In this case we include replicates of treatments within blocks to have a total number of observations that will be comparable to case when we test both factors at once. The original data set has 20 observations, so we will simulate 12 observations 3 for the irrigation experiment and 10 observations for the nitrogen experiment.
| Block | Irrigation | Nitrogen | rep | Yield |
|---|---|---|---|---|
| 1 | 1 | 160 | 1 | 56.82849058 |
| 2 | 1 | 160 | 1 | 58.28407460 |
| 1 | 2 | 160 | 1 | 71.85949159 |
| 2 | 2 | 160 | 1 | 68.05398749 |
| 1 | 1 | 160 | 2 | 58.82404173 |
| 2 | 1 | 160 | 2 | 63.91184633 |
| 1 | 2 | 160 | 2 | 68.68347448 |
| 2 | 2 | 160 | 2 | 64.19005738 |
| 1 | 1 | 160 | 3 | 56.46286929 |
| 2 | 1 | 160 | 3 | 57.66065748 |
| 1 | 2 | 160 | 3 | 70.89138521 |
| 2 | 2 | 160 | 3 | 68.89097567 |
We analyze these data using the following model (Equation (1.6)), which is known to have the correct form because the data were simulated based on the same model. Note that the model does not mention the effect of nitrogen because the same dose of nitrogen was applied to all plots.
\[\begin{equation} Y_{ijh} = \mu_{...} + \tau_i + B_j + \epsilon_{ijh} \\[20pt] \tau_i = \text{ effect of irrigation treatment with level i} \\ B_j = \text{ effect of block j} \\ \epsilon_{ijh} = \text{ random error in treatment i, block j and replicate h} \tag{1.6} \end{equation}\]
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Irrigation | 1 | 306.0 | 306.0 | 38.28 | 0.00016 |
| Block | 1 | 0.5 | 0.5 | 0.07 | 0.79982 |
| Residuals | 9 | 71.9 | 8.0 | NA | NA |
| Irrigation | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| 1 | 58.7 | 1.15 | 9 | 56.05 | 61.27 |
| 2 | 68.8 | 1.15 | 9 | 66.15 | 71.37 |
The ANOVA table shows that there is a significant difference in yield between treatments. Table 1.7 contains the estimated means for both irrigations and shows that yield is significantly greater with two than with a single irrigation. The precision of the estimates is reflected in both the MSE in the ANOVA table and in the width of the confidence intervals of the table of means, which is about 5.22 bushels/ac (1 bushel of wheat is equivalent to 60 lb). The width of the confidence interval is calculated as twice the product of the critical t value, which depends on the degrees of freedom of the residuals (dfe = 9), and the standard error, which is the square root of the MSE divided by the square root of the number of experimental units averaged to get each estimated mean. There were 12 experimental units equally divided between the two treatments in the experiment, so six experimental units were used for the average of each treatment.
\[\begin{equation} \text{CI width } = 2 \ t_{dfe, \alpha} \ \sqrt{\frac{MSE}{r}} \\[20pt] \text{where the number of plots averaged per treatment is } r = 6 \\[20pt] dfe = 9 \\[20pt] t_{dfe, \alpha} = 2.2622 \\[20pt] \text{therefore, CI width } = 2 \times 2.2622 \times 1.154 = 5.222 \end{equation}\]
1.3.4.2 Effects of Nitrogen alone
Following with the simple approach of one factor at a time, we do a second experiment to determine the effect of nitrogen fertilization on wheat yield. Five equally-spaced doses of nitrogen (0, 80, 160, 240 and 320 lb/ac) are selected. For this experiment, one irrigation is used for all plots. The five treatments are replicated in each of two blocks, such that each treatment appears in a single experimental unit per block.
| Block | Irrigation | Nitrogen | Yield |
|---|---|---|---|
| 1 | 1 | 0 | 33.52849058 |
| 2 | 1 | 0 | 34.98407460 |
| 1 | 1 | 80 | 53.55949159 |
| 2 | 1 | 80 | 49.75398749 |
| 1 | 1 | 160 | 58.82404173 |
| 2 | 1 | 160 | 63.91184633 |
| 1 | 1 | 240 | 59.43347448 |
| 2 | 1 | 240 | 54.94005738 |
| 1 | 1 | 320 | 51.36286929 |
| 2 | 1 | 320 | 52.56065748 |
We analyze these data using the following model (Equation (1.7)), which is known to have the correct form because the data were simulated based on the same model. Note that the model does not mention the effect of irrigation because the same single irrigation was applied to all plots.
\[\begin{equation} \begin{aligned} Y_{ij} &= \mu_{...} + \tau_i + B_j + \epsilon_{ij} \\[20pt] \tau_i &= \text{ effect of nitrogen treatment with level i} \\ B_j &= \text{ effect of block j} \\ \epsilon_{ij} &= \text{ random error in treatment i, block j} \end{aligned} \tag{1.7} \end{equation}\]
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Nitrogen | 4 | 854.1 | 213.5 | 26.67 | 0.00382 |
| Block | 1 | 0.0 | 0.0 | 0.00 | 0.95329 |
| Residuals | 4 | 32.0 | 8.0 | NA | NA |
| Nitrogen | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| 0 | 34.3 | 2 | 4 | 28.70 | 39.81 |
| 80 | 51.7 | 2 | 4 | 46.10 | 57.21 |
| 160 | 61.4 | 2 | 4 | 55.81 | 66.92 |
| 240 | 57.2 | 2 | 4 | 51.63 | 62.74 |
| 320 | 52.0 | 2 | 4 | 46.41 | 57.52 |
The ANOVA table for the nitrogen experiment shows that there is a significant difference in yield between treatments. Table 1.10 contains the estimated means for nitrogen doses and shows that yield is highest when a dose of 160 lb N/ac is applied. Note that the response to nitrogen in not a straight line, and that the effect of adding 80 lb N/ac depends on how much was added already. An important part of this result is that more is not always better when it comes to nitrogen and plants. Too much nitrogen has a negative effect on plant growth, and in wheat it is known to be detrimental to the metabolism necessary for grain filling. Too much nitrogen results in reduced antioxidant ability, which lead to increased lipid peroxidation [KongEtAl2017]. Moreover, wheat is unable to absorb all the applied nitrogen, which leads to nitrogen percolating into deeper soil layers and water contamination [DongEtAl2011]. A graph of the response to nitrogen levels reveals the non-linearity of the response and the reduction of yield for high N doses (Figure 1.10).
Figure 1.10: Effects of nitrogen fertilization level on wheat yield. Vertical error bars are confidence intervals.
Precision of mean estimates is reflected in both the MSE in the ANOVA table and in the width of the confidence intervals of the table of means, which is about 11.11 bushels/ac. Again, confidence interval width is calculated as twice the product of the critical t value, which depends on the degrees of freedom of the residuals (in this case dfe = 4), and the standard error, which is the square root of the MSE divided by the square root of the number of experimental units averaged to get each estimated mean. There were 10 experimental units equally divided between the 5 treatments in the experiment, so 2 experimental units were used for the average of each nitrogen treatment.
\[\begin{equation} \text{CI width } = 2 \ t_{dfe, \alpha} \ \sqrt{\frac{MSE}{r}} \\[20pt] \text{where the number of plots averaged per treatment is } r = 2 \\[20pt] dfe = 4 \\[20pt] t_{dfe, \alpha} = 2.7764 \\[20pt] \text{therefore, CI width } = 2 \times 2.7764 \times 2.001 = 11.11 \end{equation}\]
If we base our irrigation and fertilization decisions in these two experiments we would decide to irrigate twice and apply 160 lb N/ac. As we will see below, it turns out that this is not the best decision possible, and that only a factorial experiment can reveal the optimal combination of irrigation and nitrogen because those two factors interact biologically and statistically.
1.3.4.3 Combined Effects of Irrigation and Nitrogen
Instead of two experiments, we could do a single one where all levels of irrigation are combined with all levels of nitrogen in a factorial treatment structure. Each treatment is replicated once in each of two blocks. The experiment has 5 levels of nitrogen x 2 levels of irrigation x 2 blocks = 20 plots or experimental units. This is 2 fewer than the two individual experiments together.
| Irrig | Nitro | Block 1 | Block 2 |
|---|---|---|---|
| 1 | 0 | 31.7 | 39.1 |
| 1 | 80 | 50.0 | 48.6 |
| 1 | 160 | 56.8 | 60.6 |
| 1 | 240 | 57.3 | 59.4 |
| 1 | 320 | 53.1 | 54.1 |
| 2 | 0 | 38.9 | 40.0 |
| 2 | 80 | 61.4 | 53.1 |
| 2 | 160 | 68.5 | 66.7 |
| 2 | 240 | 73.5 | 72.9 |
| 2 | 320 | 70.6 | 72.3 |
Data from the factorial experiment are analyzed using the model in Equation (1.8). The equation shows that each treatment effect is partitioned into nitrogen effect, irrigation effect and nitrogen by irrigation interaction. The interaction reflects the part of the response that is not explained by simply adding the individual effects of each factor.
\[\begin{equation} \begin{aligned} Y_{ijh} &= \mu_{...} + \tau_{ij} + B_j + \epsilon_{ij} \\[20pt] \tau_{ij} &= \text{ effect of treatment with levels i of nitrogen and j of irrigation} \\ &= N_i + I_j + NI_{ij} \\ N_i &= \text{ effect of level i of nitrogen} \\ I_j &= \text{ effect of level j of irrigation} \\ NI_{ij} &= \text{ interaction effect of levels i and j} \\ B_j &= \text{ effect of block j} \\ \epsilon_{ij} &= \text{ random error in treatment i, block j} \end{aligned} \tag{1.8} \end{equation}\]
| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Irrigation | 1 | 574.6 | 574.6 | 68.65 | 0.00002 |
| Nitrogen | 4 | 2163.1 | 540.8 | 64.61 | 0.00000 |
| Block | 1 | 1.2 | 1.2 | 0.15 | 0.70814 |
| Irrigation:Nitrogen | 4 | 123.4 | 30.8 | 3.68 | 0.04826 |
| Residuals | 9 | 75.3 | 8.4 | NA | NA |
Because we did the experiment using a factorial arrangement of treatments, we have three types of comparisons we can do:
- compare average yield for levels of irrigation averaged over levels of nitrogen,
- compare average yield for levels of nitrogen averaged over levels of irrigation, and
- compare averages of individual treatments resulting from any two combinations of irrigation and nitrogen.
The estimated means for irrigation levels and the estimated means for nitrogen levels are called marginal means because they are at the margins of the table and are the averages over all levels of the other factors. The estimated means for individual treatments or combinations of the factors are frequently called cell means.
Compare tables 1.10 vs. 1.13 and 1.7 vs. 1.14 to see the advantages of the factorial. In the case of nitrogen, the factorial experiment has larger dfe than the nitrogen experiment. In both cases the standard error (SE) is smaller in the factorial, mostly because estimated means come from averaging more observations than in the separate experiments, although the separate experiments had more observations in total. By combining the treatments in a factorial we increase the degrees of freedom in the residuals, which results in a smaller critical t value, and we increase the total number of observations averaged to get the estimated marginal means.
| Nitrogen | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| 0 | 37.4 | 1.45 | 9 | 34.15 | 40.70 |
| 80 | 53.3 | 1.45 | 9 | 50.00 | 56.55 |
| 160 | 63.2 | 1.45 | 9 | 59.88 | 66.42 |
| 240 | 65.8 | 1.45 | 9 | 62.50 | 69.05 |
| 320 | 62.5 | 1.45 | 9 | 59.25 | 65.80 |
| Irrigation | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|
| 1 | 51.1 | 0.91 | 9 | 49.00 | 53.14 |
| 2 | 61.8 | 0.91 | 9 | 59.72 | 63.86 |
| Nitrogen | Irrigation | emmean | SE | df | lower.CL | upper.CL |
|---|---|---|---|---|---|---|
| 0 | 1 | 35.4 | 2.05 | 9 | 30.8 | 40.0 |
| 80 | 1 | 49.3 | 2.05 | 9 | 44.7 | 53.9 |
| 160 | 1 | 58.7 | 2.05 | 9 | 54.1 | 63.3 |
| 240 | 1 | 58.4 | 2.05 | 9 | 53.7 | 63.0 |
| 320 | 1 | 53.6 | 2.05 | 9 | 49.0 | 58.2 |
| 0 | 2 | 39.5 | 2.05 | 9 | 34.8 | 44.1 |
| 80 | 2 | 57.3 | 2.05 | 9 | 52.6 | 61.9 |
| 160 | 2 | 67.6 | 2.05 | 9 | 63.0 | 72.2 |
| 240 | 2 | 73.2 | 2.05 | 9 | 68.6 | 77.8 |
| 320 | 2 | 71.4 | 2.05 | 9 | 66.8 | 76.1 |
| Irrigation levels | |||
|---|---|---|---|
| **N levels* | 1 Irrigation | 2 Irrigations | N means |
| 0 | 35.4 | 39.5 | 37.4 |
| 80 | 49.3 | 57.3 | 53.3 |
| 160 | 58.7 | 67.6 | 63.2 |
| 240 | 58.4 | 73.2 | 65.8 |
| 320 | 53.6 | 71.4 | 62.5 |
| Irrig. means | 51.1 | 61.8 |
High precision, reflected in narrow confidence intervals and small LSD’s are desirable because it means we have more certainty about the results of the experiment. Both the LSD and the width of confidence intervals have a factor that is the product of the critical t and the square root of the MSE divided by the number of observations averaged. Both the LSD and the CI width become smaller due of the joint effects of increased dfe and increased number of observations averaged for the marginal means.
\[LSD = t_{dfe, \alpha} \ \sqrt{\frac{2 \ MSE}{n_{obs}}}\]
\[CI_{width} = 2 \ t_{dfe, \alpha}\sqrt{\frac{MSE}{n_{obs}}}\]
A second and perhaps more important advantage of the factorial approach is that it reveals the potential interaction between factors. If there is an interaction, the effect of one factor depends on the level of the other factor. In the case of nitrogen and irrigation, there is a typical positive interaction, for example, because the effect of adding an extra 80 lb N/ac, from 160 to 240 lb N/ac, is greater with two than with a single irrigation (Figure 1.11). The value of the difference in effects can be calculated using a linear combination of treatment averages that is also a contrast:
Figure 1.11: Effects of nitrogen fertilization and irrigation on wheat yield. Vertical error bars are confidence intervals.
1.4 Contrasts
The Lomu Brho experiment (Figure 1.7) was done because we wanted to answer specific questions about how plants compete with each other. Based on previous observations, it looked like Lomu was a stronger competitor than Brho. This should result in an asymmetry of the interaction between species. First, we wanted to check that competition within 10 cm of a plant did have an impact on plant size, so a “control” treatment where plants were grown without competitors was included. We also wanted to test the hypothesis that intraspecific competition was greater in Lomu than in Brho, and that having half of the surrounding plants being of the other species (“2.Both” level of competitor) was sufficient for Lomu to experience reduced and Brho to experience increased competition relative to having only con-specifics. If indeed Brho is the weaker competitor, then Brho-1.Self should be larger than Brho-2.Both, and Lomu-1.Self should be smaller than Lomu-2.Both. These questions translate directly in contrasts that are special linear combinations of treatment means.
Recall that a linear combination of means (“ingredients”) is like a recipe: each ingredient is multiplied by a known constant (coefficient) and added together. Contrasts are linear combinations of parameters in which the sum of the coefficients is zero.
\[Linear \ combination \ of \ estimated \ means = \sum_{i=1}^k c_i \ \hat{\mu}_{i} \\[20pt] \text{where } c_i \text{ are known constants}\]
- A contrast is a linear combination of estimated parameters whose coefficients add up to zero.
- Two contrasts, \(C_1\) and \(C_2\), are orthogonal if the sum of the products of corresponding coefficients equals zero.
\[Contrast \ C \ = \sum_{i=1}^k c_i \ \hat{\mu}_{i} \\[20pt] \text{where } c_i \text{ are known constants and }\sum_{i=1}^k c_i = 0\]
\[\begin{equation} C_1 = \sum_{i=1}^kc_{1i}\hat{\mu}_{i} \ \text{ and } \ C_2 = \sum_{i=1}^k c_{2i}\hat{\mu}_{i} \ \text{ are orthogonal if }\\[20pt] \sum_{i=1}^kc_{1i}c_{2i} = 0 \tag{1.9} \end{equation}\]
Contrasts are useful to test hypotheses because they have an expected value of 0 if the corresponding null hypothesis is true. Therefore, a contrast whose value is significantly different from zero leads to rejection of the null hypothesis.
Each contrast has an associated sum of squares and one degree of freedom.
Sums of squares of treatments, factors and interactions can be partitioned into as many orthogonal contrasts as their degrees of freedom.
If the null hypothesis that the contrast of true parameter values is zero is true, the mean square of a contrast has an F distribution with 1 degree of freedom in the numerator and dfe degrees of freedom in the denominator.
If the null hypothesis that the contrast of true parameter values is zero is true, the quotient of the estimated contrast and its standard error has a t distribution with dfe degrees of freedom.
- All testable hypothesis can be expressed as linear combinations of estimated parameters.
Contrasts can be classified as Class comparisons or Trends.
Class comparisons are comparisons among treatments that have categorical factor levels. For example, the competition experiment involved class comparisons.
Trends are comparisons among treatments that have quantitative levels. For example, in the irrigation by nitrogen experiment we can determine if the shape of the response to N was linear, quadratic, etc.
Class Comparisons
Treatments can have structures other than factorials. A pair of typical classes of treatment is what is called “control” and “placebo.” The control usually is a “treatment” for which experimental units receive no manipulation at all, whereas the placebo receives the same manipulation required to deliver an active principle, treatment of medicine, without actually including the compound or procedure that is being tested. This “controls” for the potential effects of the procedure that are not caused by the compound or procedure of interest. The use of placebos in medical experimentation is well established and it is generally understood. Patients can respond not only to the chemical compound in a medications but also to the idea of being cared for and to the process of ingesting a pill that might help them.
In agricultural sciences, controls are also necessary in some experiments. In order to study the effect of a dietary additive on animal production, we would want to have a control treatment that adds the same amount of energy and nutrients, but without the specific chemical compounds whose effect we want to test. To test the effect of a chemical in a pesticide application, a control treatment should involve the application of water and other additives as in the main treatment, but without the active principle to be tested.
One interesting example of necessary control is in the study of conditioned flavor aversions in animals. To test if sheep can be conditioned to avoid eating grape leaves to use them to graze between rows in vineyards, we dosed sheep with lithium chloride (Laca, unpublished work), which causes digestive discomfort. Sheep were given grape leaves and then dosed with LiCl delivered directly to the rumen with a tube. To determine if the LiCL was the cause of the aversion instead of the intubation procedure used to deliver it, we had a treatment in which sheep were intubated and given an equivalent amount of tap water without LiCl.
The design of proper controls is an art beyond statistics. It involves concepts and interpretations that can be highly creative. Correct experimentation and interpretation of results requires much more than formally correct statistical methods.
In addition to controls and placebos, treatments may have natural classes that can be compared as such. For example, we may want to test the effects of several methods to control weeds on crop yield (fictitious example). We are not interested in isolating the effect of the chemicals’ active principles, because those are already known to exist. The seven treatments are:
- Control: Undisturbed control
- chemA: Spraying chemical A
- chemB: Spraying chemical B
- Flame: Flaming
- Steam: Covered steaming
- Hoe: Rotary hoe
- Blast: Abrasive blasting
The effects of the chemical sprays can be determined with the following contrasts:
- Effect of spraying chemical A: chemA - Control
- Effect of spraying chemical B: chemB - Control
- Difference between effects of spraying A or B = (chemA - Control) - (chemB - Control) = chemA - chemB
Weed control methods tested can also be seen as being in two classes: chemical and non-chemical. Furthermore, non-chemical treatments can be classified as thermal or non-thermal. These classes offer a natural (but not unique) way to design contrasts. The difference between chemical and non-chemical treatments is calculated as the following contrast:
- Difference between chemical and non-chemical treatments = 2 (chemA + chemB) - (Flame + Steam + Hoe + Blast)
If we are interested in the difference between thermal and non-thermal treatments, we can use the following contrast:
- Difference between thermal and non-thermal methods = (Flame + Steam) - (Hoe + Blast)
Finally, we want to test if Flame method is different from Blast4
However, these comparisons (contrasts) may not be all orthogonal to each other.
In order to check for orthogonality we need to specify all coefficients for all treatment averages, including the coefficients that are zero. Coefficients for treatment averages are organized in a vector with 7 columns to make a table of coefficients with 6 rows (contrasts) and 7 columns, one for each treatment mean. For the tests of hypotheses, the corresponding contrast is calculated by multiplying the row of the table times the vector of treatment averages:
\[\begin{equation} \text{Vector with 7 columns, one for each estimated treatment mean} \\[10pt] \hat{\vec{\mu}} = (\hat{\mu}_{Control}, \ \hat{\mu}_{chemA}, \ \hat{\mu}_{chemB}, \ \hat{\mu}_{Flame}, \ \hat{\mu}_{Steam}, \ \hat{\mu}_{Hoe}, \ \hat{\mu}_{Blast}) \\[20pt] \\ \text{Vector of contrast coefficients to estimate} \\ \text{the effect of chemical A relative to the control} \\[10pt] \vec{c}_{chemA} = (-1, \ 1, \ 0, \ 0, \ 0, \ 0, \ 0) \\[20pt] \\ \text{Contrast chemA is obtained by summing the products of each} \\ \text{coefficient with the corresponding estimated mean: }\\[10pt] chemA = \vec{c}_{chemA} . \hat{\vec{\mu}} = \\[10pt] = -1 \ \hat{\mu}_{Control} + 1 \ \hat{\mu}_{chemA} + 0 \ \hat{\mu}_{chemB} + 0 \ \hat{\mu}_{Flame} + 0 \ \hat{\mu}_{Steam} + 0 \ \hat{\mu}_{Hoe} + 0 \ \hat{\mu}_{Blast} \\[10pt] = \hat{\mu}_{chemA} - \hat{\mu}_{Control}\\[20pt] \text{Vector of contrast coefficients to estimate} \\ \text{the effect of chemical B relative to the control} \\[10pt] \vec{c}_{chemB} = (-1, \ 0, \ 1, \ 0, \ 0, \ 0, \ 0) \\[20pt] \\ \vec{c}_{chemA} \text{ and } \vec{c}_{chemB} \text{ are not orthogonal because } \ \vec{c}_{chemA} . \vec{c}_{chemB} = 1 \ne 0 \end{equation}\]
In the following R chunk we us R functions to facilitate the calculations and determine which pairs of contrasts are orthogonal and winch ones are not. First we create a vector of coefficients for each contrast and give it a name that helps us remember what it is. Then, the vectors are joined into a matrix using rbind and the columns are given names that indicate the average that is multiplied by each coefficient. A few pairs are tested individually and finally the matrix is used to calculate the sum of products for all possible pairs at once.
# Name the estimated treatment means
avg.names <- c("Control", "chemA", "chemB", "Flame", "Steam", "Hoe", "Blast")
chemAC <- c(-1, 1, 0, 0, 0, 0, 0)
chemBC <- c(-1, 0, 1, 0, 0, 0, 0)
chemAB <- c(0, 1, -1, 0, 0, 0, 0)
ch.nch <- c(0, 2, 2, -1, -1, -1, -1)
th.nth <- c(0, 0, 0, 1, -1, 1, -1)
FlameBlast <- c(0, 0, 0, 1, 0, 0, -1)
contr.mat <- rbind(
chemAC,
chemBC,
chemAB,
ch.nch,
th.nth,
FlameBlast
)
colnames(contr.mat) <- avg.names
# Are chemA and chemB orthogonal? NO
sum(chemAC * chemBC)## [1] 1# Are chemA and chemAB orthogonal? NO
sum(chemBC * chemAB)## [1] -1# Are A and ch.nch orthogonal? NO
sum(chemAC * ch.nch)## [1] 2# Are A and th.nth orthogonal? YES
sum(chemAC * th.nth)## [1] 0# A table of all possible pairs is easily constructed by using matrix multiplication:
orth.mat <- contr.mat %*% t(contr.mat)| Control | chemA | chemB | Flame | Steam | Hoe | Blast | |
|---|---|---|---|---|---|---|---|
| chemAC | -1 | 1 | 0 | 0 | 0 | 0 | 0 |
| chemBC | -1 | 0 | 1 | 0 | 0 | 0 | 0 |
| chemAB | 0 | 1 | -1 | 0 | 0 | 0 | 0 |
| ch.nch | 0 | 2 | 2 | -1 | -1 | -1 | -1 |
| th.nth | 0 | 0 | 0 | 1 | -1 | 1 | -1 |
| FlameBlast | 0 | 0 | 0 | 1 | 0 | 0 | -1 |
| chemAC | chemBC | chemAB | ch.nch | th.nth | FlameBlast | |
|---|---|---|---|---|---|---|
| chemAC | 2 | 1 | 1 | 2 | 0 | 0 |
| chemBC | 1 | 2 | -1 | 2 | 0 | 0 |
| chemAB | 1 | -1 | 2 | 0 | 0 | 0 |
| ch.nch | 2 | 2 | 0 | 12 | 0 | 0 |
| th.nth | 0 | 0 | 0 | 0 | 4 | 2 |
| FlameBlast | 0 | 0 | 0 | 0 | 2 | 2 |
Table 1.17 shows that 6 out of the 15 pairs of contrasts are not orthogonal. For example, the number 2 in row 1 (chemA) and column 4 (ch.nch) of the table in the right results from the sum of the products of coefficients in rows 1 and 4 of the table on the left:
\[\begin{align} 2 = (-1)&\times0 \\ + 1 &\times 2 \\ + 0 &\times 2 \\ + 0 &\times (-1) \\ + 0 &\times (-1) \\ + 0 &\times (-1) \\ + 0 &\times (-1) \\ = \vec{c}_{chemA} &. \vec{c}_{ch.nch} \end{align}\]
Although not mandatory, orthogonal contrasts are generally considered desirable particularly if they are posed before the data are collected. These so-called a priori sets of hypothesis are independent form each other and account for all the variation due to experimental treatments. Because they are a priori, it is accepted that no corrections are needed to account for inflation of experiment wise error rate.
A complete set of orthogonal contrasts for the weed control example can be created as follows by taking a class-comparisons approach. The set has to have 6 contrasts because there are 7 treatments, and there is one orthogonal contrast per degree of freedom. This is not the only set of orthogonal contrasts possible; in fact there are 15 different possible sets of orthogonal contrasts for 7 treatment levels. The number of sets gets large very quickly as the number of levels increases.
TvsC <- c(-6, 1, 1, 1, 1, 1, 1) # Average treatment effect.
CvsNC <- c(0, -2, -2, 1, 1, 1, 1) # Chemical vs. non-chemical
AvsB <- c(0, 1, -1, 0, 0, 0, 0) # Do chemical treatments differ?
THvsNTH <- c(0, 0, 0, 1, 1, -1, -1) # Thermal vs. non-thermal
CSvsFF <- c(0, 0, 0, -1, 1, 0, 0) # Does Steam differ from Flame?
RHvsAB <- c(0, 0, 0, 0, 0, -1, 1) # Does Hoe differ from Blast?
Cmat <- rbind(
TvsC,
CvsNC,
AvsB,
THvsNTH,
CSvsFF,
RHvsAB
)
(orth.mat <- Cmat %*% t(Cmat)) # All othtogonal!!## TvsC CvsNC AvsB THvsNTH CSvsFF RHvsAB
## TvsC 42 0 0 0 0 0
## CvsNC 0 12 0 0 0 0
## AvsB 0 0 2 0 0 0
## THvsNTH 0 0 0 4 0 0
## CSvsFF 0 0 0 0 2 0
## RHvsAB 0 0 0 0 0 2In the first contrast we compare all treatments against the control. In the second one we compare the chemical vs. non chemical treatments. The rationale behind the contrasts is that we compare the average of one class vs. the average of another class. Although not strictly necessary, the coefficients are multiplied by the largest denominator to make them all integers. For example, the average of chemical treatments is \[\frac{1}{2}\hat{\mu}_{chemA} + \frac{1}{2}\hat{\mu}_{chemB}\] and the average of non-chemical treatments is \[\frac{1}{4}\hat{\mu}_{Flame} + \frac{1}{4}\hat{\mu}_{Steam} + \frac{1}{4}\hat{\mu}_{Hoe} + \frac{1}{4}\hat{\mu}_{Blast}\] so we multiply all coefficients times 4 to obtain the set of integers that make the contrast coefficients.
Now that we have a meaningful set of orthogonal contrasts we can proceed with data analysis. Because contrasts are a priori, no correction for multiple comparisons will be applied.
Caution: - Some authors recommend applying corrections, even if comparisons were posed a priori.
# Give treatments the same order used in the text:
weed.ctl.dat$trt <- factor(weed.ctl.dat$trt,
levels = c(
"Control",
"chemA",
"chemB",
"Flame",
"Steam",
"Hoe",
"Blast"
)
)
# Create a column for treatments with defined contrasts.
weed.ctl.dat$Ctrt <- weed.ctl.dat$trt
# Assign contrast structure to treatments.
# Cmat has to be transposed because R takes the contrasts as columns.
contrasts(weed.ctl.dat$Ctrt) <- t(Cmat)
# Second, run the model with orthogonal contrasts.
wc.m2 <- lm(weed.mass ~ Ctrt, data = weed.ctl.dat)
# The summary of an lm model gives us the t-tests for the contrasts.
# Where Ho is that the contrast equals zero.
summary(wc.m2)##
## Call:
## lm(formula = weed.mass ~ Ctrt, data = weed.ctl.dat)
##
## Residuals:
## Min 1Q Median 3Q Max
## -318.48372 -99.72823 44.47685 104.27446 226.27340
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1470.05059 28.64963 51.31133 < 2.22e-16 ***
## CtrtTvsC -106.47088 11.69616 -9.10306 7.3447e-10 ***
## CtrtCvsNC 73.31724 21.88152 3.35065 0.00231915 **
## CtrtAvsB 213.67098 53.59855 3.98651 0.00043568 ***
## CtrtTHvsNTH 157.43191 37.89990 4.15389 0.00027778 ***
## CtrtCSvsFF -42.39589 53.59855 -0.79099 0.43559938
## CtrtRHvsAB 33.80725 53.59855 0.63075 0.53331950
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 169.4935 on 28 degrees of freedom
## Multiple R-squared: 0.820815, Adjusted R-squared: 0.7824182
## F-statistic: 21.37718 on 6 and 28 DF, p-value: 2.927676e-09When we run an lm it uses the contrasts specified for the treatment column, assuming that the contrast coefficients were given in the order of the levels of the factor column. The contrasts of the contrast matrix had names, and the program is able to read those and prints the tests with the corresponding names. The output of summary(wc.m2) shows the model call at the top, followed by summary statistics for the residuals. By definition, the average of residuals is zero. The median is greater than zero, indicating a bit of skewness. Under the heading coefficients we find a list of contrasts, their estimated value, standard error and calculated t (t value) with the corresponding probabilities for two-tailed tests; t value = Estimate / Std. Error. The estimated value of each contrast is simply the sum of the products of each mean and the corresponding contrast coefficient, for example,
\(AvsB = -1 \ \hat{\mu}_{Control} + 1 \ \hat{\mu}_{chemA} + 0 \ \hat{\mu}_{chemB} + 0 \ \hat{\mu}_{Flame} + 0 \ \hat{\mu}_{Steam} + 0 \ \hat{\mu}_{Hoe} + 0 \ \hat{\mu}_{Blast} = \hat{\mu}_{chemA} - \hat{\mu}_{Control}\)
The first row of effects is an estimate of the overall mean, called (Intercept), which is calculated based on a linear combination that is not a contrast. This line is usually not relevant for the analyses, as it tests that the overall mean is zero. Below the overall average there are lines for each contrast. On average, treatments reduced weed mass by r round(summary(wc.m2)$coefficients[2, 1], 1) kg/ha, a value that was significantly different from 0. Non-chemical methods were slightly but significantly worse than chemical ones because the had r round(summary(wc.m2)$coefficients[3, 1], 1) kg/ha more of weed mass. Interpretation of the rest of the contrasts is left as an exercise. Keep in mind that the sign of the estimated contrast and its meaning depends on which class was given the negative coefficients, which is arbitrary (i.e., we can choose).
# The summary of the analysis of variance shows the
# complete partition of total SS into contrasts.
(wc.aov.tbl <- summary(aov(weed.mass ~ Ctrt, data = weed.ctl.dat),
split = list(Ctrt = list(
"TvsC" = 1,
"CvsNC" = 2,
"AvsB" = 3,
"THvsNTH" = 4,
"CSvsFF" = 5,
"RHvsAB" = 6
))
))## Df Sum Sq Mean Sq F value Pr(>F)
## Ctrt 6 3684748 614124.6 21.37718 2.9277e-09 ***
## Ctrt: TvsC 1 2380570 2380570.3 82.86572 7.3447e-10 ***
## Ctrt: CvsNC 1 322525 322525.0 11.22683 0.00231915 **
## Ctrt: AvsB 1 456553 456552.9 15.89224 0.00043568 ***
## Ctrt: THvsNTH 1 495696 495696.1 17.25478 0.00027778 ***
## Ctrt: CSvsFF 1 17974 17974.1 0.62566 0.43559938
## Ctrt: RHvsAB 1 11429 11429.3 0.39784 0.53331950
## Residuals 28 804385 28728.0
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1kable(wc.aov.tbl[[1]],
caption = "Partition of SS into contrasts for weed control example. Meaning of contrasts are defined in R object Cmat above."
) %>% kable_styling(full_width = FALSE)| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Ctrt | 6 | 3684747.77831 | 614124.62972 | 21.3771799647 | 0.0000000029 |
| Ctrt: TvsC | 1 | 2380570.33975 | 2380570.33975 | 82.8657215630 | 0.0000000007 |
| Ctrt: CvsNC | 1 | 322525.03346 | 322525.03346 | 11.2268346681 | 0.0023191459 |
| Ctrt: AvsB | 1 | 456552.88868 | 456552.88868 | 15.8922355373 | 0.0004356827 |
| Ctrt: THvsNTH | 1 | 495696.10324 | 495696.10324 | 17.2547790692 | 0.0002777760 |
| Ctrt: CSvsFF | 1 | 17974.11409 | 17974.11409 | 0.6256643245 | 0.4355993825 |
| Ctrt: RHvsAB | 1 | 11429.29909 | 11429.29909 | 0.3978446258 | 0.5333194985 |
| Residuals | 28 | 804385.31465 | 28728.04695 | NA | NA |
# Show treatment means and confidence intervals.
emmeans(wc.m2, "Ctrt")## Ctrt emmean SE df lower.CL upper.CL
## Control 2109 75.8 28 1954 2264
## chemA 1431 75.8 28 1275 1586
## chemB 1003 75.8 28 848 1159
## Flame 1637 75.8 28 1481 1792
## Steam 1552 75.8 28 1397 1707
## Hoe 1246 75.8 28 1090 1401
## Blast 1313 75.8 28 1158 1469
##
## Confidence level used: 0.95Each contrast has 1 df and an associated SS contrast. The SS associated with each contrast in Table 1.18 is calculated as follows:
\[\begin{align} &\text{let } C_h \text{ be a row of coefficients for contrast } h \\[15pt] &\text{and } \mathbf{C} = \begin{Bmatrix}c_{hi}\end{Bmatrix} \text{be a }(k-1) \times k \text{ matrix of orthogonal contrasts} \\[15pt] &\text{for a set of k treatments in a balanced experiment with r replicates;} \\[15pt] &\text{the sum of squares associated with each row of } \mathbf{C} \text{ is: } \\[15pt] &SSC_h = \frac{r \sum_{i=1}^k \left( c_{hi}\bar{Y}_{i.} \right)^2 }{\sum_{i=1}^kc_{hi}^2} \tag{1.10} \end{align}\]
Let’s calculate the SSC for the first two contrasts in the weed control experiment “by hand” using basic R functions. Note that \(\mathbf{C}\) in Equation (1.10) is the matrix called Cmat in the previous code. Recall that the first set of contrast coefficients, which is the first row of Cmat, is called TvsC and compares all weed control methods against the undisturbed control.
# Make a vector with the estimated means for each treatment in the weed control experiment.
# We extract them from the estimated marginal means to avoid problems if the
# experiment were unbalanced.
mu.hat <- summary(emmeans(wc.m2, "Ctrt"))$emmean
names(mu.hat) <- summary(emmeans(wc.m2, "Ctrt"))$Ctrt
cbind(names(mu.hat), mu.hat)## mu.hat
## Control "Control" "2108.87589826802"
## chemA "chemA" "1430.61621605275"
## chemB "chemB" "1003.27425055837"
## Flame "Flame" "1636.72473944646"
## Steam "Steam" "1551.9329613307"
## Hoe "Hoe" "1245.65779113132"
## Blast "Blast" "1313.27228392747"# In this balanced case, aggregate also works:
aggregate(weed.mass ~ Ctrt, data = weed.ctl.dat, FUN = mean)## Ctrt weed.mass
## 1 Control 2108.875898
## 2 chemA 1430.616216
## 3 chemB 1003.274251
## 4 Flame 1636.724739
## 5 Steam 1551.932961
## 6 Hoe 1245.657791
## 7 Blast 1313.272284r <- max(weed.ctl.dat$rep)
(SS.TvsC <- r * (sum(TvsC * mu.hat))^2 / sum(TvsC^2))## [1] 2380570.34# If we use a small number of digits in the 'hand' calculations
# we obtain a similar result, but not exactly the same due to rounding.
# Let this be a warning for the corroboration of formula results
# by actually typing in rounded numbers. (It does not always work!)
5 * ((-6) * 2108.876 +
1 * 1430.616 +
1 * 1003.274 +
1 * 1636.725 +
1 * 1551.933 +
1 * 1245.658 +
1 * 1313.272)^2 /
((-6)^2 +
1^2 +
1^2 +
1^2 +
1^2 +
1^2 +
1^2)## [1] 2380571.248(SS.CvsNC <- r * (sum(CvsNC * mu.hat))^2 / sum(CvsNC^2))## [1] 322525.03355 * (0 * 2108.876 +
(-2) * 1430.616 +
(-2) * 1003.274 +
1 * 1636.725 +
1 * 1551.933 +
1 * 1245.658 +
1 * 1313.272)^2 /
(0^2 +
(-2)^2 +
(-2)^2 +
1^2 +
1^2 +
1^2 +
1^2)## [1] 322525.882Treatment SS (first row in Table 1.18) is the sum of all 6 contrast SS, because we used a complete set of orthogonal contrasts. Residual SS can be calculated by subtracting the treatment average from each observation, squaring and adding all results. MSE is SSE/dfe, where dfe = number of observations minus number of treatments, or we can get it directly from the anova in R using anova(wc.m2)$'Mean Sq'[2].
(MSE <- anova(wc.m2)$"Mean Sq"[2])## [1] 28728.04695# Calculated F for treatments vs control.
(Fcalc.TvsC <- SS.TvsC / MSE)## [1] 82.86572156df.num <- 1
(df.den <- length(weed.ctl.dat$weed.mass) - length(levels(weed.ctl.dat$Ctrt)))## [1] 28# Probability of obtaining an F this large or larger is Ho is true.
(Prob.TvsC <- pf(Fcalc.TvsC,
df1 = df.num,
df2 = df.den,
lower.tail = FALSE
))## [1] 7.34473127e-10-
Output of analyses in R typically are lists or other structures with components that can be extracted or even modified. Use the
strfunction to see what parts are there and to figure out how to “bring them out.”
1.5 Trends
In the nitrogen by irrigation fertilization experiment (Table 1.12) there are 9 df for treatments that are partitioned into 1 df for irrigation, 4 df for nitrogen and 4 df for the irrigation x nitrogen interaction. Thus, the corresponding SS can be partitioned into a total of 9 orthogonal contrasts: 1 for irrigation, 4 for nitrogen and 4 for the interaction. The specific contrasts are not unique, and we can design them to answer specific biological questions. A natural way to partition SS and df in the nitrogen by irrigation factorial is as follows:
- Does one extra irrigation increase wheat yield, averaging across levels of nitrogen?
- What is the simplest curve we can use to summarize the effect of nitrogen on wheat yield? Is the response linear, quadratic, cubic or quartic?
- Does the shape of the response to nitrogen depend on the number of irrigations? How does it differ?
Because the levels of nitrogen are quantitative, we can create contrasts that test for the significance of increasingly complex shapes of the relationship between yield and nitrogen. The fact that N levels are equally spaced greatly facilitates the calculations because we can use the following standard formulas for contrast coefficients (Equation (1.11)):
\[\begin{equation} Linear: C_1 = \gamma_1 \left[ \frac{X-\bar{X}}{\Delta}\right] \\[25pt] Quadratic: C_2 = \gamma_2 \left[ \left(\frac{X-\bar{X}}{\Delta}\right)^2-\left( \frac{k_f^2-1}{12} \right)\right] \\[25pt] Cubic: C_3 = \gamma_3 \left[ \left(\frac{X-\bar{X}}{\Delta}\right)^3- \left(\frac{X-\bar{X}}{\Delta}\right)\left( \frac{3k_f^2-7}{20} \right)\right] \\[25pt] Quartic: C_4 = \\[10pt] =\gamma_4 \left[ \left( \frac{X - \bar{X}}{\Delta} \right)^4 - \left( \frac{X - \bar{X}}{\Delta} \right)^2 \left( \frac{3k_f^2 - 13}{14} \right)+ \frac{3 \left(k_f^2 - 1\right) \left(k_f^2 - 9 \right)}{560} \right] \\[25pt] \tag{1.11} \end{equation}\]
where \(X\) are the levels of the quantitative factor, \(\bar{X}\) is the average of all factor levels, \(\Delta\) is the distance between levels, \(k_f\) is the number of levels in the factor and \(\gamma_i\) are arbitrary values used to make all contrast coefficients into integers for simplicity. In this section we will use these equations to create a complete set of orthogonal contrasts to test for trends in the irrigation by nitrogen experiment. Then we apply the contrasts and interpret the results.
There are 5 levels of nitrogen in the experiment: 0, 80, 160, 240, 320. The distance between levels is 80, and the average of levels is 160.
- Linear coefficients for N:
\[\begin{align} X &= 0\\[15pt] \bar{X} &= 160\\[15pt] \Delta &= 80\\[15pt] c_{L,0} &= \gamma_1 \left[ \frac{0-160}{80}\right] = -2 \gamma_1 \\[15pt] c_{L,80} &= \gamma_1 \left[ \frac{80-160}{80}\right] = -1 \gamma_1 \\[15pt] c_{L,160} &= \gamma_1 \left[ \frac{160-160}{80}\right] = 0 \gamma_1 \\[15pt] c_{L,240} &= \gamma_1 \left[ \frac{240-160}{80}\right] = 1 \gamma_1 \\[15pt] c_{L,320} &= \gamma_1 \left[ \frac{320-160}{80}\right] = 2 \gamma_1 \\[15pt] \text{Choose } \gamma_1 &= 1 \\[15pt] C_L &= (-2, -1, 0, 1, 2) \tag{1.12} \end{align}\]
- Quadratic coefficients for N:
\[\begin{align} k_f &= 5 \\[15pt] \left(\frac{k_f^2-1}{12} \right) &= \left( \frac{5^2-1}{12} \right) = 2 \\[15pt] c_{Q,0} &= \gamma_2 \left[ \left( \frac{0-160}{80} \right) ^2-2 \right] = 2 \gamma_2 \\[15pt] c_{Q,80} &= \gamma_2 \left[ \left( \frac{80-160}{80} \right) ^2-2 \right] = -1 \gamma_2 \\[15pt] c_{Q,160} &= \gamma_2 \left[ \left( \frac{160-160}{80} \right) ^2-2 \right] = -2 \gamma_2 \\[15pt] c_{Q,240} &= \gamma_2 \left[ \left( \frac{240-160}{80} \right) ^2-2 \right] = -1 \gamma_2 \\[15pt] c_{Q,320} &= \gamma_2 \left[ \left( \frac{320-160}{80} \right) ^2-2 \right] = 2 \gamma_2 \\[15pt] \text{Choose } \gamma_2 &= 1 \\[15pt] C_Q &= (2, -1, -2, -1, 2) \tag{1.13} \end{align}\]
- Cubic coefficients for N:
\[\begin{align} k_f &= 5 \\[15pt] \left(\frac{3 \times k_f^2-7}{20} \right) &= \left(\frac{3 \times 5^2-7}{20} \right) = 3.4\\[15pt] c_{C,0} &= \gamma_3 \left[ \left(\frac{0-160}{80}\right)^3- \left(\frac{0-160}{80}\right) \times 3.4\right] = -1.2 \gamma_3 \\[15pt] c_{C,80} &= \gamma_3 \left[ \left(\frac{80-160}{80}\right)^3- \left(\frac{80-160}{80}\right) \times 3.4\right] = 2.4 \gamma_3 \\[15pt] c_{C,160} &= \gamma_3 \left[ \left(\frac{160-160}{80}\right)^3- \left(\frac{160-160}{80}\right) \times 3.4\right] = 0 \gamma_3 \\[15pt] c_{C,240} &= \gamma_3 \left[ \left(\frac{240-160}{80}\right)^3- \left(\frac{240-160}{80}\right) \times 3.4\right] = -2.4 \gamma_3 \\[15pt] c_{C,320} &= \gamma_3 \left[ \left(\frac{380-160}{80}\right)^3- \left(\frac{380-160}{80}\right) \times 3.4\right] = 1.2 \gamma_3 \\[15pt] \text{Choose } \gamma_3 &= 5/6 \\[15pt] C_C &= (-1, 2, 0, -2, 1) \tag{1.14} \end{align}\]
- Quartic coefficients for N:
\[\begin{align} k_f &= 5 \\[15pt] \left( \frac{3k_f^2 - 13}{14} \right) = \left( \frac{3 \times 5^2 - 13}{14} \right) &= 31/7 \\[15pt] \frac{3 \left(k_f^2 - 1\right) \left(k_f^2 - 9 \right)}{560} = \frac{3 \left(5^2 - 1\right) \left(5^2 - 9 \right)}{560} &= 72/35\\[15pt] c_{Qt,0} = \gamma_4 \left[ \left( \frac{X - 160}{80} \right)^4 - \left( \frac{X - 160}{80} \right)^2 \frac{31}{7} + \frac{72}{35} \right] &= \frac{12}{35} \gamma_4 \\[15pt] c_{Qt,0}= \gamma_4 \left[ \left( \frac{X - 160}{80} \right)^4 - \left( \frac{X - 160}{80} \right)^2 \frac{31}{7} + \frac{72}{35} \right] &= \frac{-48}{35} \gamma_4 \\[15pt] c_{Qt,0} = \gamma_4 \left[ \left( \frac{X - 160}{80} \right)^4 - \left( \frac{X - 160}{80} \right)^2 \frac{31}{7} + \frac{72}{35} \right] &= \frac{70}{35} \gamma_4 \\[15pt] c_{Qt,0} = \gamma_4 \left[ \left( \frac{X - 160}{80} \right)^4 - \left( \frac{X - 160}{80} \right)^2 \frac{31}{7} + \frac{72}{35} \right] &= \frac{-48}{35} \gamma_4 \\[15pt] c_{Qt,0} = \gamma_4 \left[ \left( \frac{X - 160}{80} \right)^4 - \left( \frac{X - 160}{80} \right)^2 \frac{31}{7} + \frac{72}{35} \right] &= \frac{12}{35} \gamma_4 \\[15pt] \text{Choose } \gamma_4 &= 35/12 \\[15pt] C_Qt = (1, -4, 6, -4, 1) \tag{1.15} \end{align}\]
1.6 Linear Combinations
From a very practical point of view we may also interested in determining if linear combinations of parameters are different from known fixed values. These linear combinations do not have to be contrasts. For example, a producer who needs to choose the level of fertilization of wheat may be interested in the probability that gross income per ac of applying two irrigations and 240 lb N/ac is at least $100 more than the gross income form applying 1 irrigation and 160 lb N/ac, because the difference in costs is $100/ac. The difference in question is a linear combination of wheat yield in the two treatments, where the coefficients are determined by the price of wheat.
\[ I_{ij} \ \text{ is gross income from level i of nitrogen and level j of irrigation}\\[20pt] p \ \text{ is the price of wheat}\\[20pt] I_{240,2} = p \ \hat{\mu}_{240,2}\\[20pt] I_{160,1} = p \ \hat{\mu}_{160,1}\\[20pt] P((I_{240,2} - I_{160,1}) \ge 100) = P((p \ \hat{\mu}_{240,2} - p \ \hat{\mu}_{160,1}) \ge 100)\\[20pt] = P((p \ (\hat{\mu}_{240,2} - \ \hat{\mu}_{160,1}) - 100) \ge 0)\]
If residuals of ANOVA are normally distributed as is assumed, theory indicates that the linear combination of estimated parameters also has a normal distribution, and the ratio of the linear combination over its standard error has a Student’s t distribution. These facts allow us to use the distribution to calculate the probability sought.
The big picture: Anova tells you whether there is at least one contrast that is significantly different from zero, but it does not tell what it is. Without the restriction of orthogonality, there is an large number of contrasts that could be tested, so we need to choose carefully. When we do not choose a few contrasts a priori and start calculating comparisons (contrasts) looking for something significant, the probability of at least one Error Type I increases very fast and can reach 100% very quickly, so something needs to be done. One approach is to make orthogonal contrasts before collecting the data, during the planning of the experiment. The restriction of orthogonality will keep the number of contrasts limited to one fewer than the treatment degrees of freedom, and the a priori condition makes it acceptable not to have to correct for multiple comparisons. Another approach is to make a few post-hoc (after seeing the data) contrasts and applying some method to contain the inflation of experiment wise error rate, for example the Protected LSD, Bonferroni’s correction or Tukey’s HSD. In the case of quantitative treatments such as doses, rates, times, age, etc., contrasts can be used to naturally partition the shape of the response into incremental orthogonal components that increase the “wiggliness” of the response in an orderly way. Non-significant complexity can be ignored, which simplifies comparisons and explanations.
1.7 Exercises and Solutions
Use the data below and basic R functions to calculate Estimate, Std. Error, t value, Pr(>|t|), Sum Sq, Mean Sq, F value and Pr(>F) for the last 3 contrasts in the weed control example (Tables 1.18 and 1.19).
weed.ctl.dat## trt rep wc.fx weed.mass Ctrt
## 1 Blast 1 -100 1292.590817 Blast
## 2 Blast 4 -100 1033.306997 Blast
## 3 Blast 3 -100 1413.875790 Blast
## 4 Blast 2 -100 1421.217694 Blast
## 5 Blast 5 -100 1405.370122 Blast
## 6 chemA 1 -200 1640.896760 chemA
## 7 chemA 4 -200 1572.566873 chemA
## 8 chemA 3 -200 1199.940497 chemA
## 9 chemA 2 -200 1621.172691 chemA
## 10 chemA 5 -200 1118.504260 chemA
## 11 chemB 1 -200 711.372525 chemB
## 12 chemB 4 -200 1047.751096 chemB
## 13 chemB 3 -200 1145.176056 chemB
## 14 chemB 2 -200 1053.694085 chemB
## 15 chemB 5 -200 1058.377490 chemB
## 16 Control 1 600 1790.392174 Control
## 17 Control 4 600 2052.653449 Control
## 18 Control 3 600 2158.150912 Control
## 19 Control 2 600 2313.133287 Control
## 20 Control 5 600 2230.049669 Control
## 21 Flame 1 100 1643.910088 Flame
## 22 Flame 4 100 1515.366131 Flame
## 23 Flame 3 100 1662.122917 Flame
## 24 Flame 2 100 1862.998144 Flame
## 25 Flame 5 100 1499.226417 Flame
## 26 Hoe 1 -300 1309.283161 Hoe
## 27 Hoe 4 -300 1324.548520 Hoe
## 28 Hoe 3 -300 1117.377184 Hoe
## 29 Hoe 2 -300 1272.584375 Hoe
## 30 Hoe 5 -300 1204.495715 Hoe
## 31 Steam 1 100 1403.964871 Steam
## 32 Steam 4 100 1652.322679 Steam
## 33 Steam 3 100 1555.416272 Steam
## 34 Steam 2 100 1473.835118 Steam
## 35 Steam 5 100 1674.125867 Steamkable(summary(wc.m2)$coefficients, caption = "Summary of t-tests for contrasts in the weed control example with simulated data.") %>% kable_styling(full_width = FALSE)| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 1470.05059153 | 28.64963074 | 51.3113277071 | 0.0000000000 |
| CtrtTvsC | -106.47088446 | 11.69616277 | -9.1030611095 | 0.0000000007 |
| CtrtCvsNC | 73.31723688 | 21.88151692 | 3.3506469029 | 0.0023191459 |
| CtrtAvsB | 213.67098275 | 53.59855124 | 3.9865066835 | 0.0004356827 |
| CtrtTHvsNTH | 157.43190643 | 37.89989904 | 4.1538872239 | 0.0002777760 |
| CtrtCSvsFF | -42.39588906 | 53.59855124 | -0.7909894592 | 0.4355993825 |
| CtrtRHvsAB | 33.80724640 | 53.59855124 | 0.6307492575 | 0.5333194985 |
Interpret the results for the last 3 contrasts in the weed control example. Include statements about the effectiveness of treatments in the real world and your certainty about them.
1.8 Homework
1.9 Laboratory Exercises
1.9.1 Plant Sciences
Prepare an .Rmd document starting with the following text, where you substitute the corresponding information for author name and date.
---
title: "Lab07 Factorial ANOVA and LSD"
author: "YourFirstName YourLastName"
date: "enter date here"
output: html_document
---1.9.1.1 Instructions
The same instructions for completion and submission of work used in previous labs applies here. Refer to previous labs for the details.
For this lab you have to be familiar with the following sections of the textbook:
?? RCBD without subsamples 1 Factorial Experiments ?? Fisher’s Protected LSD
In this exercise we analyze the data resulting from an experiment at the UC Davis Plant Sciences research fields. These are based on the same experiment used in Lab 06, but now we will consider that seeding and fertilizer are two factors that create the set of 6 treatments in a factorial treatment structure. We ignore the factor water for simplicity and because it did not cause any effects. The 6 treatments were applied in a Completely Randomized Block Design in which each of 6 treatments were applied randomly to 2 plots in each of 3 blocks for a total of 36 experimental units or plots.
Simulated block effects of 0.50, 0.25, and -0.75 were added to logMass observations in blocks 1, 2 and 3 to show block effects. Blocks were defined areas that were more or less homogeneous in soil and plant properties and that were contiguous. We will assume that the variances of the errors or residuals are the same in all treatments. Each observation represents the average mass of seeds produced by 4-5 medusahead plants (logMass). Data were log transformed as \(logMass = log(SeedMass + 0.2)\) to achieve normality and equality of variances.
1.9.1.2 Part 1. Factorial ANOVA using R functions [45 points]
1.9.1.2.1 1a. Random assignment of treatments to plots [15 points]
In a RCBD the treatments are assigned randomly to plots within blocks such that each treatment appears at least once in each block. Because we ignore the water factor, treatments will have 2 replicate plots in each block, one for each original level of the water factor. There are a total of 6 treatments that result from the combination each level of the factor “nitrogen” (n and N) with each of the three levels of the factor “seeding” (s, S and E). The treatments are identified by the two letters together. For example. NE is the treatment with addition of nitrogen and at the edge between naturalized annual and seeded native perennial grasses.
First, we pretend that the experiment is about to start and we need to assign treatments to experimental units. The random assignment of treatments to plots within blocks can be done using the ’*‘sample’ R function. To preform the random assignment we start by creating a dataframe with a row for each treatment and a column for each factor. A new column called “Treatment” is created by pasting together the two letters that represent the corresponding level of the factor.
# Read and understand the code in this part but do not type or modify anything.
# Create a data frame with all combinations of levels of the 2 factors. First, create treatment names for all 6 treatments.
treat <- expand.grid(
nitro = c("n", "N"),
seeding = c("s", "S", "E")
)
# Create a column containing the treatment name in the data frame
treat$Treatment <- paste(
treat$nitro,
treat$seeding,
sep = ""
)In the second step, plots in each block are labeled with the numbers 1 through 12. We make a character vector with two copies of each treatment label, shuffle the order of treatment with ‘sample’ and column-bind the results with the plot labels. At this point, the objects called ‘block1’, ‘block2’ and ‘block3’ list the plot number and the treatment it should receive.
# Number the 12 plots in each block 1 through 12, then shuffle the treatment names and assign to plots of each block.
# The sample() function is used to create a vector that has each treatment randomly assigned to each block twice, because there are two plots for each treatment in each block.
# Plots 1-12 of each block get:
(block1 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "ns"
## [2,] "2" "nE"
## [3,] "3" "nE"
## [4,] "4" "ns"
## [5,] "5" "NS"
## [6,] "6" "NE"
## [7,] "7" "nS"
## [8,] "8" "NE"
## [9,] "9" "nS"
## [10,] "10" "Ns"
## [11,] "11" "NS"
## [12,] "12" "Ns"(block2 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "nE"
## [2,] "2" "NE"
## [3,] "3" "NE"
## [4,] "4" "nS"
## [5,] "5" "ns"
## [6,] "6" "Ns"
## [7,] "7" "nS"
## [8,] "8" "NS"
## [9,] "9" "Ns"
## [10,] "10" "nE"
## [11,] "11" "NS"
## [12,] "12" "ns"(block3 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "nE"
## [2,] "2" "ns"
## [3,] "3" "NE"
## [4,] "4" "ns"
## [5,] "5" "nS"
## [6,] "6" "nE"
## [7,] "7" "Ns"
## [8,] "8" "NS"
## [9,] "9" "NE"
## [10,] "10" "Ns"
## [11,] "11" "nS"
## [12,] "12" "NS"1.9.1.2.2 1.b Add columns with factor levels to the data [15 points]
Treatments were applied, plants grown and measured. Data were collected and organized, and we enter them below. Because the data have the labels for the original treatments, we need to split the labels into the columns that correspond to nitrogen and seeding treatment. In order to analyze the experiment as a factorial of nitrogen x seeding, each factor has to be in a separate column.
# Read and understand the code in this part but do not type or modify anything.
# Read data (seedF for "factorial")
# seedF <- read.csv(file = "Lab07SeedMassData.txt", header = TRUE)
seedF <- read.table(header = TRUE, text = "
Treatment block logMass
wnE Block1 -0.210836005431437
wNE Block1 1.55288938986704
WnE Block1 0.055235139345024
WNE Block1 0.835167804254832
wns Block1 1.24627323200302
wnS Block1 0.540066519204629
wNs Block1 1.57199447282598
wNS Block1 1.23870355136554
Wns Block1 0.986995928669204
WnS Block1 -0.170749893385349
WNs Block1 0.664136397253749
WNS Block1 0.736888706368722
wnE Block2 0.041449411568171
wNE Block2 1.428708818124
WnE Block2 0.319992367391131
WNE Block2 1.03098409320956
wns Block2 0.643797294741899
wnS Block2 -0.511211922892043
wNs Block2 1.84080033644501
wNS Block2 0.813650468511317
Wns Block2 1.69986757761409
WnS Block2 -0.154140967169852
WNs Block2 1.423732115253
WNS Block2 0.2968644855256
wnE Block3 -0.655926955677754
wNE Block3 -0.27934035269992
WnE Block3 -0.5542030003591
WNE Block3 -1.12888399571028
wns Block3 0.45627765803375
wnS Block3 -1.04875645623447
wNs Block3 0.49759250163951
wNS Block3 -0.097792008066854
Wns Block3 -0.513584856362292
WnS Block3 -1.36796547632004
WNs Block3 1.02290365359422
WNS Block3 -0.527897010482484
")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 contain labels for all 12 treatments. We create the column by pasting the letters for nitrogen and seeding levels, make it into a factor, and add it to the data frame. Recall that a factor is a special type of column in data frames that behaves as a categorical variable with a certain set of possible values. The same column of strings could be set to be a character vector, in which case it would not have recognizable levels. You can explore this by running ‘levels(seedF$NStreat)’ and then ‘levels(as.character(seedF$NStreat))’
# Add columns $nitro and $seeding to the data, to show the levels for each factor.
seedF$nitro <- factor(substr(seedF$Treatment, 2, 2))
seedF$seeding <- factor(substr(seedF$Treatment, 3, 3))
# Add column NStreat for treatments resulting from nitro and seeding combinations.
seedF$NStreat <- factor(paste(seedF$nitro, seedF$seeding, sep = ""))1.9.1.2.3 1.c Perform a test of the Ho that all treatment means are equal [15 points]
For this part you need to test the null hypotheses that: - there are no interactions between nitrogen and seeding, - there are no “main effects” of nitro or seeding.
“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 simple effect of S for level n is:
\[\bar{Y}_{nS} - \bar{Y}_{n.}\]
We do these tests with an ANOVA that partitions all variation (sum of squares) in the response variable (logMass) into Block, Treatment and Error, and then partition the Treatment variation into main effect of nitro, main effect of seeding, and the interaction between nitro and seeding (nitro x seeding).
# First do the ANOVA ignoring the treatment factorial structure.
# This yields the first partition of SS
rcbd.mod <- lm(formula = logMass ~ NStreat + block, data = seedF)
pander(anova(rcbd.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| NStreat | 5 | 10.76 | 2.153 | 11.05 | 6.049e-06 |
| block | 2 | 9.62 | 4.81 | 24.68 | 6.617e-07 |
| Residuals | 28 | 5.457 | 0.1949 | NA | NA |
# Now, analyze the same data, this time partitioning the Treatment SS according to a factorial treatment structure.
rcbdF.mod <- lm(formula = logMass ~ nitro * seeding + block, data = seedF)
# complete the code below to make an ANOVA table for the same data partitioning the Treatment SS according to factorial, as above:
pander(anova(rcbdF.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| nitro | 1 | 4.079 | 4.079 | 20.93 | 8.851e-05 |
| seeding | 2 | 6.367 | 3.183 | 16.34 | 1.988e-05 |
| block | 2 | 9.62 | 4.81 | 24.68 | 6.617e-07 |
| nitro:seeding | 2 | 0.3178 | 0.1589 | 0.8155 | 0.4527 |
| Residuals | 28 | 5.457 | 0.1949 | NA | NA |
Questions:
Report MSE in the RCBD and RCBDF models.
State the null hypothesis for the main effects of N in the RCBDF model.
State the null hypothesis for the main effects of seeding in in the RCBDF model.
Do we reject or fail to reject the following null hypothesis for the RCBDF model?
Ho: There is no interaction between nitrogen and seeding affecting medusahead seed mass?
1.9.1.3 Part 2. Factorial ANOVA detailed calculations [30 points]
Create a new column in the data for each of the following: - overall average - nitro treatment average - nitro effect = nitro average - overall average - seeding treatment average - seeding effect = seeding average - overall average - overall average + nitro effect + seeding effect - nitro x seeding average - nitro x seeding interaction = = nitro x seeding average - (overall average + nitro effect + seeding effect) - block average -block effect = block average - overall average - residual
Then calculate the sum of squares and degrees of freedom for nitro effect, seeding effect, interaction and residuals. Make sure that the numbers match the results from the previous section.
# First calculate the average medusahead seed mass per factor and add to a data frame
# Type the column name of the response variable after the dollar sign:
seedF$overall.avg <- mean(seedF$logMass)
nitro.avg <- aggregate(
formula = logMass ~ nitro,
data = seedF,
FUN = mean
)
names(nitro.avg)[2] <- "nitro.avg"
seedF <- merge(seedF, nitro.avg)
# Type the function for calculating averages:
seed.avg <- aggregate(
formula = logMass ~ seeding,
data = seedF,
FUN = mean
)
names(seed.avg)[2] <- "seed.avg"
seedF <- merge(seedF, seed.avg)
nitro.seed.avg <- aggregate(
formula = logMass ~ nitro + seeding,
data = seedF,
FUN = mean
)
names(nitro.seed.avg)[3] <- "nitro.seed.avg"
seedF <- merge(seedF, nitro.seed.avg)
# Create a vector called block.avg to aggregate the response variable by block and calculate the means of each block
block.avg <- aggregate(
formula = logMass ~ block,
data = seedF,
FUN = "mean"
)
names(block.avg)[2] <- "block.avg"
seedF <- merge(seedF, block.avg)
# Second, calculate the effects and residuals
seedF$nitro.fx <- seedF$nitro.avg - seedF$overall.avg
# Type the vector name for the average of the seeded/unseeded treatment:
seedF$seed.fx <- seedF$seed.avg - seedF$overall.avg
seedF$nitro.seed.fx <- with(seedF, nitro.seed.avg - (overall.avg + nitro.fx + seed.fx))
seedF$block.fx <- seedF$block.avg - seedF$overall.avg
seedF$resdl <- with(
seedF,
logMass -
overall.avg -
nitro.fx -
seed.fx -
nitro.seed.fx -
block.fx
)
# Finally, type your equations to calculate the sum of squares.
(SSnitro <- sum(seedF$nitro.fx^2))## [1] 4.079406108(SSseed <- sum(seedF$seed.fx^2))## [1] 6.366998066(SSnitro.seed <- sum(seedF$nitro.seed.fx^2))## [1] 0.3178363507(SSblock <- sum(seedF$block.fx^2))## [1] 9.620035671(SSresdl <- sum(seedF$resdl^2))## [1] 5.456590694# HINT: check your answers of the SS with the anova table in Part1. It should be the same if your calculations are correct:
pander(anova(rcbdF.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| nitro | 1 | 4.079 | 4.079 | 20.93 | 8.851e-05 |
| seeding | 2 | 6.367 | 3.183 | 16.34 | 1.988e-05 |
| block | 2 | 9.62 | 4.81 | 24.68 | 6.617e-07 |
| nitro:seeding | 2 | 0.3178 | 0.1589 | 0.8155 | 0.4527 |
| Residuals | 28 | 5.457 | 0.1949 | NA | NA |
1.9.1.4 Part 3. Multiple comparison of means using Fisher’s PLSD [25 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 nitro, seeding and nitro x seeding levels. For each one, construct a table that shows the means that are NOT significantly different connected by a common letter.
# First, use the LSD.test() function of the agricolae package.
LSD.test(rcbdF.mod, "nitro", group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ "nitro"
##
## LSD t Test for logMass
##
## Mean Square Error: 0.1948782391
##
## nitro, means and individual ( 95 %) CI
##
## logMass std r LCL UCL Min
## n 0.04458775526 0.7946030192 18 -0.1685506451 0.2577261557 -1.367965476
## N 0.71783907929 0.8054141606 18 0.5047006789 0.9309774797 -1.128883996
## Max
## n 1.699867578
## N 1.840800336
##
## Alpha: 0.05 ; DF Error: 28
## Critical Value of t: 2.048407142
##
## least Significant Difference: 0.3014232165
##
## Treatments with the same letter are not significantly different.
##
## logMass groups
## N 0.71783907929 a
## n 0.04458775526 bLSD.test(rcbdF.mod, "seeding", group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ "seeding"
##
## LSD t Test for logMass
##
## Mean Square Error: 0.1948782391
##
## seeding, means and individual ( 95 %) CI
##
## logMass std r LCL UCL Min
## E 0.20293639282 0.8506851147 12 -0.05810376998 0.4639765556 -1.1288839957
## s 0.96173219264 0.6623600110 12 0.70069202984 1.2227723554 -0.5135848564
## S -0.02102833363 0.7797215931 12 -0.28206849643 0.2400118292 -1.3679654763
## Max
## E 1.552889390
## s 1.840800336
## S 1.238703551
##
## Alpha: 0.05 ; DF Error: 28
## Critical Value of t: 2.048407142
##
## least Significant Difference: 0.3691665386
##
## Treatments with the same letter are not significantly different.
##
## logMass groups
## s 0.96173219264 a
## E 0.20293639282 b
## S -0.02102833363 bLSD.test(rcbdF.mod, c("nitro", "seeding"), group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ c("nitro", "seeding")
##
## LSD t Test for logMass
##
## Mean Square Error: 0.1948782391
##
## nitro:seeding, means and individual ( 95 %) CI
##
## logMass std r LCL UCL Min
## n:E -0.1673815072 0.3797097928 6 -0.53654804575 0.20178503136 -0.6559269557
## N:E 0.5732542928 1.0577140454 6 0.20408775428 0.94242083140 -1.1288839957
## n:s 0.7532711391 0.7614245756 6 0.38410460056 1.12243767767 -0.5135848564
## n:S -0.4521263661 0.6856580171 6 -0.82129290469 -0.08295982758 -1.3679654763
## N:s 1.1701932462 0.5302111982 6 0.80102670761 1.53935978473 0.4975925016
## N:S 0.4100696989 0.6491290541 6 0.04090316031 0.77923623743 -0.5278970105
## Max
## n:E 0.3199923674
## N:E 1.5528893899
## n:s 1.6998675776
## n:S 0.5400665192
## N:s 1.8408003364
## N:S 1.2387035514
##
## Alpha: 0.05 ; DF Error: 28
## Critical Value of t: 2.048407142
##
## least Significant Difference: 0.5220803256
##
## Treatments with the same letter are not significantly different.
##
## logMass groups
## N:s 1.1701932462 a
## n:s 0.7532711391 ab
## N:E 0.5732542928 b
## N:S 0.4100696989 b
## n:E -0.1673815072 c
## n:S -0.4521263661 c# Second, calculate the LSD's and CI's "by hand."
# Hint: check your answers of LSD with the lsd test table. It should be the same if your calculations are correct.
(dfe.rcbd.F <- 36 - 1 - 2 - 2 - 2 - 1)## [1] 28(t.critical <- qt(0.975, dfe.rcbd.F))## [1] 2.048407142(MSE.rcbd.F <- SSresdl / dfe.rcbd.F)## [1] 0.1948782391# nitrogen
(se.nitro <- sqrt(MSE.rcbd.F / (length(seedF$logMass) / length(levels(seedF$nitro)))))## [1] 0.1040507993(LSD.nitro <- t.critical * sqrt(2) * se.nitro)## [1] 0.3014232165# seeding
(se.seed <- sqrt(MSE.rcbd.F / (length(seedF$logMass) / length(levels(seedF$seeding)))))## [1] 0.1274356828# Type your equation and calculate the LSD for seeding.
(LSD.seed <- t.critical * sqrt(2) * se.seed)## [1] 0.3691665386# nitrogen x seeding
(se.nitro.seed <- sqrt(MSE.rcbd.F /
(length(seedF$logMass) /
(length(levels(seedF$nitro)) *
length(levels(seedF$seeding))))))## [1] 0.180221271# Type your equation and calculate the LSD for nitro * seeding
(LSD.nitro.seed <- t.critical * sqrt(2) * se.nitro.seed)## [1] 0.5220803256Question Read the last lsd table and answer following questions:.
Do treatments “n:s” and “N:s” differ?
Do treatments “n:s” and “N:E” differ?
Do treatments “n:s” and “n:E’ differ?
1.9.2 Animal Sciences
---
title: "Lab07 Factorial ANOVA and LSD"
author: "YourFirstName YourLastName"
date: "enter date here"
output: html_document
---1.9.2.1 Instructions
The same instructions for completion and submission of work used in previous labs applies here. Refer to previous labs for the details.
For this lab you have to be familiar with the following sections of the textbook:
8.1 RCBD without subsamples 8.5 Factorial Experiments 9.1 Fisher’s Protected LSD
In this exercise we analyze the milk protein data modified from the Milk R dataset. These are based on the same experiment used in Lab 06, but now we will consider that lactation period (A, B, C, D) and diet (barley, lupin, barley+lupine) are two factors that create the set of 12 treatments in a factorial treatment structure. We ignore the factor water for simplicity and because it did not cause any effects. The 12 treatments were applied in a Completely Randomized Block Design in which each of 12 treatments were applied randomly to cows blocked by body weight.
Simulated block effects of 0.10, 0.20, and -0.30 were added to logProt observations in blocks 1, 2 and 3. We will assume that the variances of the errors or residuals are the same in all treatments.
1.9.2.2 Part 1. Factorial ANOVA using R functions [ points]
1.9.2.2.1 1a. Random assignment of treatments to plots
In a RCBD the treatments are assigned randomly to cows within blocks such that each treatment appears at least once in each block. The random assignment of treatments to cows within blocks can be done using the sample() R function.
# Read and understand the code in this part but do not type or modify anything.
# Create a data frame with all combinations of levels of the 2 factors. First, create treatment names for all 6 treatments.
treat <- expand.grid(
diet = c("bar", "lup", "bar+lup"),
period = c("A", "B", "C", "D")
)
# Create a column containing the treatment name in the data frame
treat$Treatment <- paste(
treat$diet,
treat$period,
sep = ""
)
# Number the 12 plots in each block 1 through 12, then shuffle the treatment names and assign to plots of each block.
# The sample() function is used to create a vector that has each treatment randomly assigned to each block twice, because there are two plots for each treatment in each block.
# Plots 1-12 of each block get:
(block1 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "bar+lupA"
## [2,] "2" "bar+lupD"
## [3,] "3" "lupD"
## [4,] "4" "lupA"
## [5,] "5" "lupB"
## [6,] "6" "bar+lupB"
## [7,] "7" "bar+lupA"
## [8,] "8" "lupB"
## [9,] "9" "barC"
## [10,] "10" "lupC"
## [11,] "11" "barB"
## [12,] "12" "lupA"
## [13,] "1" "barC"
## [14,] "2" "bar+lupC"
## [15,] "3" "bar+lupC"
## [16,] "4" "barD"
## [17,] "5" "barB"
## [18,] "6" "bar+lupB"
## [19,] "7" "lupC"
## [20,] "8" "lupD"
## [21,] "9" "barA"
## [22,] "10" "barA"
## [23,] "11" "bar+lupD"
## [24,] "12" "barD"(block2 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "lupA"
## [2,] "2" "barD"
## [3,] "3" "lupB"
## [4,] "4" "bar+lupC"
## [5,] "5" "barC"
## [6,] "6" "barA"
## [7,] "7" "bar+lupB"
## [8,] "8" "bar+lupA"
## [9,] "9" "lupD"
## [10,] "10" "barB"
## [11,] "11" "bar+lupC"
## [12,] "12" "lupD"
## [13,] "1" "lupC"
## [14,] "2" "barC"
## [15,] "3" "lupC"
## [16,] "4" "bar+lupD"
## [17,] "5" "barB"
## [18,] "6" "barD"
## [19,] "7" "lupB"
## [20,] "8" "lupA"
## [21,] "9" "bar+lupA"
## [22,] "10" "bar+lupD"
## [23,] "11" "barA"
## [24,] "12" "bar+lupB"(block3 <- cbind(1:12, sample(c(treat$Treatment, treat$Treatment))))## [,1] [,2]
## [1,] "1" "lupA"
## [2,] "2" "lupB"
## [3,] "3" "barA"
## [4,] "4" "lupC"
## [5,] "5" "bar+lupC"
## [6,] "6" "bar+lupD"
## [7,] "7" "barD"
## [8,] "8" "bar+lupD"
## [9,] "9" "barC"
## [10,] "10" "lupD"
## [11,] "11" "lupB"
## [12,] "12" "barC"
## [13,] "1" "bar+lupB"
## [14,] "2" "barA"
## [15,] "3" "barB"
## [16,] "4" "barD"
## [17,] "5" "lupA"
## [18,] "6" "bar+lupA"
## [19,] "7" "lupC"
## [20,] "8" "bar+lupA"
## [21,] "9" "lupD"
## [22,] "10" "barB"
## [23,] "11" "bar+lupC"
## [24,] "12" "bar+lupB"1.9.2.2.2 1.b Add columns with factor levels to the data
# Read and understand the code in this part but do not type or modify anything.
# Read data (periodF for "factorial")
milk.prot <- read.csv(file = "Datasets/Lab06MilkProteinData.csv", header = TRUE)
# Add columns $diet and $period to the data, to show the levels for each factor.
milk.prot$diet <- factor(gsub("[A-Z]", "", milk.prot$Treatment)) # eliminate capital letters
milk.prot$period <- factor(substr(milk.prot$Treatment, 4, 4))
milk.prot$period <- factor(substr(milk.prot$Treatment, nchar(as.character(milk.prot$Treatment)), nchar(as.character(milk.prot$Treatment)))) # take just the last character1.9.2.2.3 1.c Perform a test of the Ho that all treatment means are equal.
For this part you need to test the null hypotheses that: - there are no interactions between diet and period, - there are no “main effects” of diet or period.
“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.
We do these tests with an ANOVA that partitions all variation (sum of squares) in the response variable (logProt) into Block, Treatment and Error, and then partition the Treatment variation into main effect of diet, main effect of period, and the interaction between diet and period (diet x period).
# First do the ANOVA ignoring the treatment factorial structure.
# This yields the first partition of SS
rcbd.mod <- lm(formula = logProt ~ Treatment + block, data = milk.prot)
pander(anova(rcbd.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| Treatment | 11 | 0.9113 | 0.08284 | 10.67 | 1.862e-06 |
| block | 2 | 1.712 | 0.856 | 110.3 | 3.414e-12 |
| Residuals | 22 | 0.1708 | 0.007762 | NA | NA |
# Now, analyze the same data, this time partitioning the Treatment SS according to a factorial treatment structure.
rcbdF.mod <- lm(formula = logProt ~ diet * period + block, data = milk.prot)
# complete the code below to make an ANOVA table for the same data partitioning the Treatment SS according to factorial, as above:
pander(anova(rcbdF.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| diet | 2 | 0.6591 | 0.3295 | 42.46 | 2.802e-08 |
| period | 3 | 0.1553 | 0.05178 | 6.671 | 0.002263 |
| block | 2 | 1.712 | 0.856 | 110.3 | 3.414e-12 |
| diet:period | 6 | 0.09688 | 0.01615 | 2.08 | 0.0972 |
| Residuals | 22 | 0.1708 | 0.007762 | NA | NA |
Questions; a) Report MSE in the RCBD and RCBDF models. b) State the null hypothesis for the main effects of N in the RCBDF model c) State the null hypothesis for the main effects of period in in the RCBDF model. c) Do we reject or fail to reject the following null hypothesis for the RCBDF model? Ho: There is no interaction between diet and period affecting milk production?
1.9.2.3 Part 2. Factorial ANOVA detailed calculations [ points]
Create a new column in the data for each of the following: - overall average - diet treatment average - diet effect = diet average - overall average - period treatment average - period effect = period average - overall average - overall average + diet effect + period effect - diet x period average - diet x period interaction = = diet x period average - (overall average + diet effect + period effect) - block average -block effect = block average - overall average - residual
Then calculate the sum of squares and degrees of freedom for diet effect, period effect, interaction and residuals. Make sure that the numbers match the results from the previous section.
# First calculate the average milk produced and add to a data frame
# type the column name of the response variable after the dollar sign:
milk.prot$overall.avg <- mean(milk.prot$logProt)
diet.avg <- aggregate(
formula = logProt ~ diet,
data = milk.prot,
FUN = mean
)
names(diet.avg)[2] <- "diet.avg"
milk.prot <- merge(milk.prot, diet.avg)
# type the function for calculating averages:
period.avg <- aggregate(
formula = logProt ~ period,
data = milk.prot,
FUN = mean
)
names(period.avg)[2] <- "period.avg"
milk.prot <- merge(milk.prot, period.avg)
diet.period.avg <- aggregate(
formula = logProt ~ diet + period,
data = milk.prot,
FUN = mean
)
names(diet.period.avg)[3] <- "diet.period.avg"
milk.prot <- merge(milk.prot, diet.period.avg)
# Create a vector called block.avg to aggregate the response variable by block and calculate the means of each block
block.avg <- aggregate(
formula = logProt ~ block,
data = milk.prot,
FUN = "mean"
)
names(block.avg)[2] <- "block.avg"
milk.prot <- merge(milk.prot, block.avg)
# Second, calculate the effects and residuals
milk.prot$diet.fx <- milk.prot$diet.avg - milk.prot$overall.avg
# type the vector name for the average of the perioded/unperioded treatment:
milk.prot$period.fx <- milk.prot$period.avg - milk.prot$overall.avg
milk.prot$diet.period.fx <- with(milk.prot, diet.period.avg - (overall.avg + diet.fx + period.fx))
milk.prot$block.fx <- milk.prot$block.avg - milk.prot$overall.avg
milk.prot$resdl <- with(milk.prot, logProt - overall.avg - diet.fx - period.fx - diet.period.fx - block.fx)
# Finally, type your equations to calculate the sum of squares.
(SSdiet <- sum(milk.prot$diet.fx^2))## [1] 0.6590641499(SSperiod <- sum(milk.prot$period.fx^2))## [1] 0.1553369656(SSdiet.period <- sum(milk.prot$diet.period.fx^2))## [1] 0.09688021954(SSblock <- sum(milk.prot$block.fx^2))## [1] 1.712091214(SSresdl <- sum(milk.prot$resdl^2))## [1] 0.1707607133# HINT: check your answers of the SS with the anova table in Part1. It should be the same if your calculations are correct:
pander(anova(rcbdF.mod))| Df | Sum Sq | Mean Sq | F value | Pr(>F) | |
|---|---|---|---|---|---|
| diet | 2 | 0.6591 | 0.3295 | 42.46 | 2.802e-08 |
| period | 3 | 0.1553 | 0.05178 | 6.671 | 0.002263 |
| block | 2 | 1.712 | 0.856 | 110.3 | 3.414e-12 |
| diet:period | 6 | 0.09688 | 0.01615 | 2.08 | 0.0972 |
| Residuals | 22 | 0.1708 | 0.007762 | NA | NA |
1.9.2.4 Part 3. Multiple comparison of means using Fisher’s PLSD [ 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 diet, period and diet x period levels. For each one, construct a table that shows the means that are NOT significantly different connected by a common letter.
# First, use the LSD.test() function of the agricolae package.
LSD.test(rcbdF.mod, "diet", group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ "diet"
##
## LSD t Test for logProt
##
## Mean Square Error: 0.007761850607
##
## diet, means and individual ( 95 %) CI
##
## logProt std r LCL UCL
## bar+lup -0.2777315400 0.2529537598 12 -0.3304756767 -0.22498740327
## barley -0.1520857663 0.2486223980 12 -0.2048299030 -0.09934162959
## lupins -0.4805077371 0.2613397992 12 -0.5332518738 -0.42776360041
## Min Max
## bar+lup -0.7764267356 0.04643648580
## barley -0.5587781208 0.16287801246
## lupins -0.9981033424 -0.07277649498
##
## Alpha: 0.05 ; DF Error: 22
## Critical Value of t: 2.073873068
##
## least Significant Difference: 0.07459147349
##
## Treatments with the same letter are not significantly different.
##
## logProt groups
## barley -0.1520857663 a
## bar+lup -0.2777315400 b
## lupins -0.4805077371 cLSD.test(rcbdF.mod, "period", group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ "period"
##
## LSD t Test for logProt
##
## Mean Square Error: 0.007761850607
##
## period, means and individual ( 95 %) CI
##
## logProt std r LCL UCL Min
## A -0.1973961870 0.2209948551 9 -0.2582998700 -0.1364925039 -0.5298895083
## B -0.3318335871 0.2319143397 9 -0.3927372701 -0.2709299040 -0.6841721634
## C -0.3753809789 0.3029214556 9 -0.4362846620 -0.3144772959 -0.8458581716
## D -0.3091559716 0.3680558951 9 -0.3700596547 -0.2482522885 -0.9981033424
## Max
## A 0.09367342884
## B -0.02434783968
## C -0.05673715804
## D 0.16287801246
##
## Alpha: 0.05 ; DF Error: 22
## Critical Value of t: 2.073873068
##
## least Significant Difference: 0.08613081459
##
## Treatments with the same letter are not significantly different.
##
## logProt groups
## A -0.1973961870 a
## D -0.3091559716 b
## B -0.3318335871 b
## C -0.3753809789 bLSD.test(rcbdF.mod, c("diet", "period"), group = TRUE, console = TRUE)##
## Study: rcbdF.mod ~ c("diet", "period")
##
## LSD t Test for logProt
##
## Mean Square Error: 0.007761850607
##
## diet:period, means and individual ( 95 %) CI
##
## logProt std r LCL UCL
## bar+lup:A -0.21366941489 0.2218517024 3 -0.3191576883 -0.10818114145
## bar+lup:B -0.30915393911 0.1831156463 3 -0.4146422126 -0.20366566567
## bar+lup:C -0.34781054154 0.3790472142 3 -0.4532988150 -0.24232226810
## bar+lup:D -0.24029226443 0.3292239112 3 -0.3457805379 -0.13480399099
## barley:A -0.09010314859 0.2516073774 3 -0.1955914220 0.01538512485
## barley:B -0.20671148464 0.3049643590 3 -0.3121997581 -0.10122321120
## barley:C -0.22388679243 0.2384502655 3 -0.3293750659 -0.11839851899
## barley:D -0.08764163959 0.3203164064 3 -0.1931299130 0.01784663385
## lupins:A -0.28841599741 0.2296489141 3 -0.3939042708 -0.18292772397
## lupins:B -0.47963533743 0.1777063506 3 -0.5851236109 -0.37414706399
## lupins:C -0.55444560287 0.2878335936 3 -0.6599338763 -0.44895732943
## lupins:D -0.59953401080 0.3516735150 3 -0.7050222842 -0.49404573736
## Min Max
## bar+lup:A -0.4694223183 -0.07310371007
## bar+lup:B -0.5204241566 -0.19610186835
## bar+lup:C -0.7764267356 -0.05673715804
## bar+lup:D -0.5998240551 0.04643648580
## barley:A -0.3768654981 0.09367342884
## barley:B -0.5587781208 -0.02434783968
## barley:C -0.4984340597 -0.06854723285
## barley:D -0.4485548725 0.16287801246
## lupins:A -0.5298895083 -0.07277649498
## lupins:B -0.6841721634 -0.36311693599
## lupins:C -0.8458581716 -0.27032973184
## lupins:D -0.9981033424 -0.33293569797
##
## Alpha: 0.05 ; DF Error: 22
## Critical Value of t: 2.073873068
##
## least Significant Difference: 0.149182947
##
## Treatments with the same letter are not significantly different.
##
## logProt groups
## barley:D -0.08764163959 a
## barley:A -0.09010314859 a
## barley:B -0.20671148464 ab
## bar+lup:A -0.21366941489 ab
## barley:C -0.22388679243 ab
## bar+lup:D -0.24029226443 b
## lupins:A -0.28841599741 b
## bar+lup:B -0.30915393911 b
## bar+lup:C -0.34781054154 bc
## lupins:B -0.47963533743 cd
## lupins:C -0.55444560287 d
## lupins:D -0.59953401080 d# Second, calculate the LSD's and CI's "by hand."
# Hint: check your answers of LSD with the lsd test table. It should be the same if your calculations are correct.
(dfe.rcbd.F <- 36 - 1 - 2 - 3 - 2 - 6)## [1] 22(t.critical <- qt(0.975, dfe.rcbd.F))## [1] 2.073873068(MSE.rcbd.F <- SSresdl / dfe.rcbd.F)## [1] 0.007761850607# dietgen
(se.diet <- sqrt(MSE.rcbd.F / (length(milk.prot$logProt) / length(levels(milk.prot$diet)))))## [1] 0.02543267355(LSD.diet <- t.critical * sqrt(2) * se.diet)## [1] 0.07459147349# period
(se.period <- sqrt(MSE.rcbd.F / (length(milk.prot$logProt) / length(levels(milk.prot$period)))))## [1] 0.02936712184(LSD.period <- t.critical * sqrt(2) * se.period) # type your equation and calculate the LSD for period## [1] 0.08613081459# dietgen x period
(se.diet.period <- sqrt(MSE.rcbd.F /
(length(milk.prot$logProt) /
(length(levels(milk.prot$diet)) *
length(levels(milk.prot$period))))))## [1] 0.0508653471(LSD.diet.period <- t.critical * se.diet.period) # type your equation and calculate the LSD for diet * period## [1] 0.1054882734Question Read the last lsd table and answer following questions:
Question Read the last lsd table and answer following questions:. a) Do treatments barleyD and barleyB differ? b) Do treatments barley+lupineD and lupineB differ? c) Do treatments barleyA and lupineA differ?
\(Y = 1000 \ \frac{[4 + \log(PlantMass.g + 0.03)]} {6}\)↩︎
In this experiment, the overall average is the same as the estimated overall mean because the data are balanced (all treatments appear once in each block). When the number of observations differs among treatments, the estimated overall mean is not the overall average, but the average of all estimated treatment means.↩︎
We use 12 instead of 10 for simplicity, to have the same number of replicates equal in all treatments. This is not absolutely necessary, but different sample sizes per treatment makes things a bit more complicated.↩︎
Note that for simplicity, we use the language in a lax manner here and in other parts of the book. When writing “if Flame method is different from Blast,” we actually mean “if the effect of the Flame method on response variable Y is the same as the effect of the Blast method on Y.” Of course, the methods themselves ARE different as one involves burning weeds and in the other, weeds are blasted with an abrasive agent.↩︎