2.1 Acuity experiment data

An experiment is conducted that measures visual acuity in groups of infants ranging from 1 to 5 months of age. Load the data into R with the following command:

load(file=url('http://pnb.mcmaster.ca/bennett/psy710/datasets/acuityWide.rda'))

Tasks:

1. Identify the number, names, and types of variables in acuityWide.df.

2. Convert acuityWide.df into a long-format data frame named acuityLong.df. Store the five ages in a single variable named age and the acuity measures in a single variable named acuity.

3. Use summary to get a summary of the variables in the new data frame. Verify that age is a factor and acuity is numeric.

4. Calculate the mean acuity for each age.

5. Plot a box plot of acuity for each age.

3.1 coefficients

We can view the best-fitting coefficients with the coef and dummy.coef commands.

coef(aov.res.01) # coefficients for 3 groups listed here
dummy.coef(aov.res.01) # coefficients for all 4 groups listed here... why?

coef(lm.res.01) # coefficients for 3 groups listed here
dummy.coef(lm.res.01) # coefficients for all 4 groups listed here... why?

The coefficients for an lm object, but not an aov object, can be submitted easily to a t test.

summary(lm.res.01) # t test on each coefficient

Question: What null hypothesis is being evaluated by each \(t\) test?

Because we are using the sum-to-zero constraint on our alpha’s, the effect of being group j is defined as the difference between the mean for group j and the unweighted mean (i.e., the average of the group means), and the intercept corresponds to the unweighted group mean. We have equal n in each group, so the unweighted mean equals the grand mean of the scores.

Task: Calculate the group means, grand mean, the unweighted mean, and the effect for each group (i.e., difference between each group mean and the unweighted mean). Verify that these alphas are the same as those estimated by aov and/or lm.

dummy.coef(aov.res.01) 
( grpMeans <- with(theAOVdata,tapply(score,condition,mean)) ) # group means
( grandMean <- mean(theAOVdata$score) ) # equals Intercept
( unweightedMean <- mean(grpMeans)) # equal n so this equals grand mean & intercept
grpMeans - grandMean # equals alphas/coefficients

3.2 treatment definition of effects

It is possible to use alternative definitions of effects. For example we could define one of our conditions (c1) as a baseline condition, the others (i.e., c2, c3, & c4) as so-called treatment conditions, and the effect of a treatment as the difference between the treatment group mean and the baseline group mean. Using the treatment definition of effects, which is the default definition in R, will change the values of the best-fitting coefficients estimated by ANOVA: the intercept equals the mean of the baseline group and the three coefficients correspond to the difference between each treatment group and the baseline group.

options(contrasts=c("contr.treatment","contr.poly")) # important!
aov.res.02 <- aov(score~condition,data=theAOVdata) # recalculate least-squares coefficients
coef(aov.res.02)
dummy.coef(aov.res.02) 
grpMeans - grpMeans[1]
options(contrasts=c("contr.sum","contr.poly")) # restore sum-to-zero definition

Which definition of an effect is best? It depends on the circumstances: in some sense it is arbitrary. Fortunately, the \(F\) test in a 1-way ANOVA always evaluates the same null hypothesis (H0: \(\alpha_j = 0\) for all \(j\)) and yields the same \(p\) value regardless of how we define \(\alpha_j\). (However, this is not true in the case of an unbalanced factorial ANOVA, a point that we will return to later in the course).

anova(aov.res.02)

3.3 anova table

The ANOVA table is listed with the anova command. The first line in the ANOVA table lists the change in goodness-of-fit (i.e., SS-residuals) when we force all of the group effects to equal zero. The \(F\) test evaluates the null hypothesis that all of the group effects are zero or, equivalently, that the four conditions are samples from populations with the same mean. Essentially, we are comparing the goodness-of-fit (i.e., SS-residuals) of our full model, which includes and intercept and an effect for each group, and the fit of a reduced model that includes only an intercept. We expect the full model to provide a better fit because it has more parameters. The question is whether the change in goodness-of-fit is unusually large assuming that the true values of the coefficients are zero.

anova(aov.res.01) # standard anova table:

# anova(lm.res.01) # gives same results

In this case, the \(F\) test is significant \(p<0.05\) and therefore we may reject the null hypothesis that the effects are all zero for the alternative hypothesis that they are not all zero.

The data presented in an ANOVA table can be interpreted as a comparison of the goodness-of-fit obtained by a full model, which includes an intercept and effects of condition, and a reduced model that includes only an intercept (i.e., the condition effects are set to zero). To illustrate this point, we will analyze theAOVdata with a reduced linear model that contains only an intercept (i.e., no effect of condition), and compare the goodness-of-fit obtained with the full and reduced models.

# compare directly to reduced model:
aov.reduced.model <- aov(score~1,theAOVdata) # intercept only
anova(aov.reduced.model,aov.res.01) # compare 2 models (compare SS to ANOVA table)

Now calculate the difference between the sum-of-squared residuals for each model. Can you find this difference in the original ANOVA table?

# calculate sum of squared residuals from each model
SS.resid.reduced <- sum( residuals(aov.reduced.model)^2 )
SS.resid.full <- sum( residuals(aov.res.01)^2 )
SS.resid.reduced - SS.resid.full # difference

3.4 association strength & effect size

Association strength and effect size are used to characterize the magnitude of the association between the independent and dependent variables that (in theory) are not dependent on sample size. Association strength refers to the proportion of variance that is accounted for by the grouping variable. One measure, eta-squared (\(\eta^2\)), is simply SS-group divided by SS-total. Another common measure, omega-squared (\(\omega^2\)), is an unbiased estimate of the value of \(\eta^2\) in the population. Omega-squared typically is less than eta-squared, though the difference is small when sample size is large.

# install.packages("effectsize") # use this if you don't have effectsize on your computer
library(effectsize)
eta_squared(aov.res.01) # in effectsize package
omega_squared(aov.res.01) # in effectsize package

Effect size is a standardized measure of differences among means. Cohen’s \(d\), for example, expresses the difference between two group menas relative to the pooled standard deviation. Cohen’s \(f\) is a generalization of Cohen’s \(d\) to cases of 3 or more means.

# install.packages("effectsize") # use this if you don't have effectsize on your computer
# library(effectsize)
cohens_f(aov.res.01)
cohens_f_squared(aov.res.01) # f-squared also is common

One interpretation of Cohen’s \(f\) is that it represents the average standardized deviation of group means from the grand mean:

MS.residuals <- 281 # from ANOVA table
( within.grp.sd <- sqrt(MS.residuals) )
# mean( with(theAOVdata,tapply(score,condition,sd)) ) # same as sqrt(MS.residuals)
( Cohen_d_squared <- ((grpMeans-grandMean)/within.grp.sd)^2 ) # cohen d^2 for each group
sqrt( mean( Cohen_d_squared ) ) # approximate Cohen's f 

You can get measures of association strength and effect size directly from the values in an ANOVA table, which is useful when you are trying to estimate these quantities from previously published data.

( eta2 <- 3226/(3226+21328)) # values from anova table
## [1] 0.1314
sqrt(eta2 / (1-eta2)) # estimate of cohen's f
## [1] 0.3889

3.5 checking assumptions

The analysis of variance assumes that the residuals are distributed normally and that the variance is equal across groups. We can evaluate those assumptions with shapiro.test and bartlett.test.

shapiro.test(residuals(aov.res.01)) # H0: residuals are distributed normally
bartlett.test(score~condition,theAOVdata) # H0: variance is equal across groups

3.6 ANOVA questions

Use the data in acuityLong.df (created in Section 2.1) to answer the following questions. Make sure to use the sum-to-zero definition of effects ( options(contrasts=c("contr.sum","contr.poly")) ).

1. Use the analysis of variance to evaluate the effect of age on acuity and list the ANOVA table.

2. What null hypothesis is being evaluated by the \(F\) test?

3. Calculate omega-squared and Cohen’s \(f\) for condition. What do these values mean?

4. Check the normality and homogeneity of variance assumptions.

5. Show that the SS value for age is the change in the goodness-of-fit when the group effects in the full model are set to zero.

6. Verify that MS-residuals in the full model equals the average within-group variance.