The analysis of variance (Part III)

stripchart(consumptionDifferenceFromControl ~ ppmCaffeine, data = strungOutBees, vertical = TRUE, method = "jitter", xlab="Caffeine (ppm)")

Discuss: State the null and alternative hypotheses appropriate for this question.

Answer: \[ \begin{eqnarray*} H_{0} & : & \mu_{50} = \mu_{100} = \mu_{150} = \mu_{200} \\ H_{A} & : & \mathrm{At \ least \ one \ of \ the \ means \ is \ different} \end{eqnarray*} \]

caffResults <- lm(consumptionDifferenceFromControl ~ ppmCaffeine, data=strungOutBees)
anova(caffResults)

Analysis of Variance Table

Response: consumptionDifferenceFromControl
            Df Sum Sq Mean Sq F value  Pr(>F)  
ppmCaffeine  3 1.1344 0.37814  4.1779 0.02308 *
Residuals   16 1.4482 0.09051                  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

So there is a difference between means amongst all groups. Now what? I want to know which means are different from one another!!!

Discuss: Which groups do you think have significantly different means?

Definition: A planned comparison is a comparison between means planned during the design of the study, identified before the data are examined.

A planned comparison must have a strong a priori justification, such as an expectation from theory or a prior study.

Only one or a small number of planned comparisons is allowed, to minimize inflating the Type I error rate.

Definition: An unplanned comparison is one of multiple comparisons, such as between all pairs of means, carried out to help determine where differences between means lie.

Unplanned comparisons are a form of data dredging, so we need to minimize the rising Type I errors that we get from performing many tests.

A planned comparison is essentially a 2-sample \( t \)-test, except you use the \( \mathrm{MS}_{\mathrm{error}} \) in place of the pooled sample variance when computing the standard error, i.e.

\[ \mathrm{SE} = \sqrt{\mathrm{MS}_{\mathrm{error}}\left(\frac{1}{n_{1}} + \frac{1}{n_{2}}\right)} \]

In R, we use the multcomp package to do this (could do it by hand, of course, but ain't nobody got time for that!)

library(multcomp)
caffPlanned <- glht(caffResults, linfct = mcp(ppmCaffeine = c("ppm200 - ppm50 = 0")))
confint(caffPlanned)

ppm200 - ppm50 = 0 is called a contrast. In this case, it simply means “Test if the means between the two groups ppm50 and ppm200 are the same.”

     Simultaneous Confidence Intervals

Multiple Comparisons of Means: User-defined Contrasts


Fit: lm(formula = consumptionDifferenceFromControl ~ ppmCaffeine, 
    data = strungOutBees)

Quantile = 2.1199
95% family-wise confidence level


Linear Hypotheses:
                    Estimate lwr      upr     
ppm200 - ppm50 == 0  0.37000 -0.03336  0.77336

Definition: With the Tukey-Kramer method, the probability of making at least one Type I error throughout the course of testing all pairs of means is no greater than the significance level \( \alpha \).

library(multcomp)
tukeyResults <- glht(caffResults, linfct = mcp(ppmCaffeine = "Tukey"))

Groups in the figure are assigned the same symbol if their means are not significantly different.

Groups in the figure are assigned the same symbol if their means are not significantly different.

• Same assumptions as ANOVA.
• \( P \)-value for T-K is exact for balanced designs.
• \( P \)-value for T-K is conservative for unbalanced designs.
• Conservative means real probability of making at least one Type I error is smaller than \( \alpha \), which makes it harder to reject \( H_{0} \).

Definition: Fixed-effects ANOVA is ANOVA on fixed groups, i.e. when different categories of the explanatory variable are
     - predetermined,
     - of direct interest,
     - and repeatable.

Any conclusion reached about differences among fixed groups apply only to those fixed groups.

Definition: Random-effects ANOVA is ANOVA on randomly chosen groups, which are groups sampled from a much larger “population” of groups.

Conclusions reached about differences among randomly chosen groups can be generalized to the whole population of groups.

Examples:

• Geographical site for an observational field study
• Subject or individual, in a study involving repeated measures on an individual.

Definition: An explanatory variable is called a fixed effect if the groups are predefined and are of direct interest. An explanatory variable is called a random effect if the groups are randomly sampled from a population of possible groups.

Important: The main use of random-effects ANOVA is to estimate variance components, i.e. the amount of the variance in the data that is among random groups and the amount that is within groups.

Examples:

• Explanatory variables: Genes, Environment
• Response variable: Phenotypic trait

What explains more of the variance in the phenotype - genes or environment?

Important: The main use of random-effects ANOVA is to estimate variance components, i.e. the amount of the variance in the data that is among random groups and the amount that is within groups.

Examples:

• Explanatory variables: Genes, Environment
• Response variable: Phenotypic trait

Question: What explains more of the variance in the phenotype - genes or environment?

Single-factor ANOVA with random effects has two levels of random variation in the response variable \( Y \) - the error and the groups.

Definition: The variance within groups in the population is written \( \sigma^2 \), with
\[ \sigma^2 \approx \mathrm{MS}_{\mathrm{error}}. \]

Definition: The variance between groups (second level of random variation in random-effects ANOVA) is the variance among the group means in the population of groups and is denoted \( \sigma_{A}^{2} \). The grand mean \( \mu_{A} \) is the mean of group means.

Random-effects ANOVA assumes that the group means are normally distributed, i.e.

\[ \mu_{G} \sim N(\mu_{A},\sigma^{2}_{A}). \]

Definition: The parameters \( \sigma^{2} \) and \( \sigma_{A}^{2} \) are called variance components. They describe all the variance in the response variable \( Y \).

For a balanced design, \( \sigma_{A}^{2} \) can be estimated by

\[ s_{A}^{2} = \frac{\mathrm{MS}_{\mathrm{groups}} - \mathrm{MS}_{\mathrm{error}}}{n}. \]

Definition: Repeatability is the fraction of the summed variance that is present among groups:

\[ \mathrm{Repeatability} = \frac{s_{A}^2}{s_{A}^{2} + \mathrm{MS}_{\mathrm{error}}}. \]

A repeatability near zero indicates that most of the variation is within groups.

A repeatability near one indicates that most of the variation is between groups.

Practice Problem #11

One way to assess whether a trait in males has a genetic basis is to determine how similar the measurements of that trait are among his offspring born to different, randomly chosen females. In a lab experiment, Kotiaho et al. (2001) randomly sampled 12 male dung beetles, Onthophagus taurus, and mated each of them to three different virgin females. The average body-condition score of offspring born to each of the three females was measured. ANOVA was used to test whether males differed in the mean condition of their offspring using the three measurements for each male.

Question: Is this a random-effects or fixed-effects ANOVA?

Let's load the data:

dungData <- read.csv("http://whitlockschluter.zoology.ubc.ca/wp-content/data/chapter15/chap15q11DungBeetleCondition.csv")
str(dungData)

'data.frame':   36 obs. of  2 variables:
 $ male              : int  1 1 1 2 2 2 3 3 3 4 ...
 $ offspringCondition: num  0.82 0.44 0.92 0.35 0.19 1.39 0.12 0.84 0.16 0.49 ...

Let's visualize the data (see Chap. 15 R-Pubs):

Let's get the variance components. We will use a new package called nlme:

library(nlme)
dungBeetleAnova <- lme(fixed = offspringCondition ~ 1, random = ~ 1|male, data = dungData)
tmp <- VarCorr(dungBeetleAnova)
tmp

male = pdLogChol(1) 
            Variance  StdDev   
(Intercept) 0.2361843 0.4859879
Residual    0.1950916 0.4416917

male = pdLogChol(1) 
            Variance  StdDev   
(Intercept) 0.2361843 0.4859879
Residual    0.1950916 0.4416917

(Intercept) corresponds to \( s_{A}^{2} \), which estimates \( \sigma_{A}^{2} \) (variance between groups)

Residual corresponds to \( \mathrm{MS}_{\mathrm{error}} \), which estimates \( \sigma^{2} \) (variance within groups)

Definition: The heritability of a trait is the fraction of variation in the trait in the population that is genetic rather than environmental.

In our example, differences among males indicate a genetic component, because they were randomly mated and their offspring were raised in a common (lab) environment (if properly designed, of course).

In our case, heritability is the same as the repeatability defined earlier, i.e.

\[ \mathrm{Heritability} = \frac{s_{A}^2}{s_{A}^{2} + \mathrm{MS}_{\mathrm{error}}}. \]

Let's calculate heritability:

varAmong  <- as.numeric( tmp[1,1] )
varWithin <- as.numeric( tmp[2,1] )
heritability <- varAmong / (varAmong + varWithin)
heritability

[1] 0.5476408

So an estimated 55% of the variation in offspring body-condition scores was due to genetic differences among males.