You should always examine the data. If you are worried about data snooping, then make sure that you decide beforehand (and if necessary write down) which analyses you plan do to, and why. Proceeding with analysis without doing at least cursory textual and graphical summaries of the data is insane.
From this we can see that the data set is extremely simple (a single numerical response and a single categorical predictor). The number of operators is so small that we can immediately see the experimental design (5 measurements for each of 4 operators), but if necessary we could use
to get a full enumeration of the number of observations per operator (in this case, we could also increase the value of maxsum (default 7), which specifies how many levels of factors to report, e.g. summary(pulp,maxsum=10): see ?summary.data.frame). For more complicated data sets, table is really useful for figuring out the experimental design: how treatments are distributed among blocks, whether the design is balanced, etc..
We could go slightly further in getting summaries:
The command stem(pulp$bright) gives a text-based approximation of a histogram of the data (something of an anachronism at this point, but cute).
At this point it really makes much more sense to turn to graphical summaries.
Try
to see what you get. Were you surprised? Why is this form of plotting suboptimal for a data set with only 5 observations per group, even though it would make sense for a larger data set …? You could add the individual points to the plot:
By now you should have noticed that some points have duplicated numerical values. The quick and dirty way to deal with this is to jitter the data, e.g.:
Here is a ggplot version. Scaling the points according to the number of overplotted points is easy: just use stat_sum(). The complexity below (aes(size=..n..) and the scale_size_continuous stuff) is there to make the plot a bit prettier.
As a side note, the fact that there are duplicated data points means that the data set cannot technically be normal: this would be unlikely to be a big deal unless we were trying to do something very sophisticated with the data set (and probably nothing we should be trying to do with only 20 observations in 4 groups).
In such a simple case (and with a small number of treatment levels), there is going to be very little practical difference between treating operator as fixed vs. random, but the minor differences are somewhat educational.
Fitting the model as a fixed effect using lm, and using two ways of testing (aov and anova):
We can also use anova to explicitly test one model against its reduced version (again with the same results):
The aov() function allows a simple form of random effects estimation, via the Error specification:
(We could also have expressed the formula as an “intercept-only” model via bright~1+Error(operator).) There are no fixed effects at all in this model, so the summary shows only the error SSQ and MSQ.
We can (not particularly conveniently) test the hypothesis that the among-operator variance is zero by using pf(...,lower.tail=FALSE) to compute the appropriate upper tail area:
In this case the p-values are identical whether we treat operator as a random or a fixed effect …
I did the test above by looking at the results from summary and plugging the appropriate numbers into pf. As a general practice, it is better to specify this test based on digging the values out from inside the s1 object, but this is actually a tremendous nuisance because of the way the results are organized internally: it takes quite a bit of digging around with str to figure out what to extract.
The bottom line is that, while you can do ANOVA on classical balanced designs in R (see MASS chapter 10 for more details), it is not very convenient to do so.
In particular, the classical machinery provided by R is not well geared to doing F tests of random effects; most of the testing framework assumes that you want to test particular fixed effects with appropriate error terms. (My guess about this omission is that, in the classical experimental design framework, the random effect terms (e.g. experimental block) are all nuisance variables that are neither of interest in themselves nor dispensable if they turn out to be non-significant. In contrast, in genetics — the other historical source of random and mixed effects ANOVA — the random effects are of great interest …)
Rule of thumb: if you have a design that is more than slightly complicated (e.g. unbalanced, non-nested, etc.), then you should use lme (nlme package) or lmer (lme4 package).
The summary extracts summary statistics about the fit (AIC, BIC, and log-likelihood); standard deviations (and correlations in more complex cases) for each random effect fitted and for the residual standard error; a standard table of fixed-effect parameters with the value, standard error, and Wald statistics (and p values); a summary of residuals (good for a preliminary assessment of whether the residuals are badly skewed or otherwise wonky); and basic statistics on the number of levels per group (also good for sanity-checking your results).
The standard accessors (fitted, residuals, logLik, AIC, predict etc.) work for lme models, with some wrinkles (for each method, see the appropriate lme-specific help page, e.g. ?fitted.lme for fitted, to see specifically what it does when applied to lme objects). There are a few accessors which are new or different.
Exercise: Try out these accessors on the model fit above. Convince yourself that (1) the residual variance is the same as estimated above by the classical approach; (2) that the estimated variance among operators is not the same; (3) that the classical estimate of the among-operator variance is within the 95% confidence intervals estimated by lme.
Now let’s use ranef to demonstrate the shrinkage of the estimates (predictions) of operator effects:
The slope of the
regression is 0.76.
If we preferred, we could do this in terms of the predicted deviations from the mean by operator rather than the predicted values:
It turns out to be a nuisance to test the significance of the only random effect in a model fitted via lme. Normally you can use anova to test the differences between a full and a reduced model, but we can’t use lme to fit the model without the random effect of observer, because lme insists on at least one random effect.
An alternative is to use the RLRsim package, which uses very efficient simulation to test the effect of dropping a single random effect:
If you run this a few times, you’ll get slightly different answers each time, but all should agree reasonably closely with the classical answer (p ≈ 0.02) above.
We can also compute the intra-class correlation (VarCorr produces a slightly complicated object, so we collapse it to a numeric vector before trying to manipulate it):
In lme4, random effects are expressed in a single formula along with fixed effects, as terms of the form (effect|group) (i.e. like the terms in the random= argument to lme, but surrounded by parentheses).
This gives very similar results to lme. One of the main differences is that in the summary, lmer does not give p values …
And now for a brief exploration of a (much) more complicated data set …
In this data set, the response variable is the number of fruits (i.e. seed capsules) per plant. The number of fruits produced by an individual plant (the experimental unit) was hypothesized to be a function of fixed effects, including nutrient levels (low vs. high), simulated herbivory (none vs. apical meristem damage), region (Sweden, Netherlands, Spain), and interactions among these. Fruit number was also a function of random effects including both the population and individual genotype. Because Arabidopsis is highly selfing, seeds of a single individual served as replicates of that individual. There were also nuisance variables, including the placement of the plant in the greenhouse, and the method used to germinate seeds. These were estimated as fixed effects but interactions were excluded.
Running
(try it for yourself!) reveals that we have only 2–4 populations per region and 2–3 genotypes per population. However, we also have 2–13 replicates per genotype per treatment combination (four unique treatment combinations: 2 levels of nutrients by 2 levels of simulated herbivory). Thus, even though this was a reasonably large experiment (625 plants), there were a very small number of replicates with which to estimate variance components, and many more potential interactions than our data can support. Therefore, judicious selection of model terms, based on both biology and the data, is warranted. We note that we don’t really have enough levels per random effect, nor enough replication per unique treatment combination. Therefore, we decide to ignore the fixed effect of “region”, but now recognize that populations in different regions are widely geographically separated.
A quick glance at the mean-variance relationship:
Note that there actually three points at (0,0) (the three genotype by nutrient by clipping treatments with no fruits at all).
The linear regression indicates a relationship of var ≈ V 0 ⋅ mean1.7 (depending on details like whether we fit log(y) ∝ log(x) or log(1 + y) ∝ log(1 + x)). We will hope that adding an observation-level random effect will take care of this overdispersion.
Here is an ambitious model fit that includes all the fixed effects (nutrient × clipping interaction, experimental design details, as well as random effects for the interactions of the treatments at both the among-population and the within-population, among-genotype level.
Pay special attention to the (amd*nutrient|popu) and (amd*nutrient|gen) terms; these allow the effects of clipping (amd), nutrient addition, and their interaction to vary across populations and across genotypes. Genotypes are implicitly nested within populations; that is, they have unique labels, but each label occurs only within a single population. If they did not have not have unique labels — if we had labelled genotypes 1, 2, 3 …within each population — then we would have to specify the genotype effect as (amd*nutrient|popu:gen). We might be able to code the two random effects together as (amd*nutrient|popu/gen), but I didn’t try it.
You can run this model yourself (it took about a minute to run on my system); I ran this model and a variety of others as a batch, so you can load the results and access the value as mf_fits$full:
We have “successfully” fitted the full model — that is, we didn’t get an errors, and none of the random-effects variances were reduced to zero. However …when we look at the random-effects section, we see that the correlations between the random effects components are all -1 or 1, indicating that we have overfitted the model.
There is not very much guidance available in the literature for deciding how to cut down such complex models to a manageable size. In all, I fitted 14 models to the data: the full model, plus all combinations of (amd only, nutrient only, intercept only) varying at the population or genotype level (9 models), plus a model with only an intercept effect of genotype or population (no random effects at the other level), or an observation-level random effect only (no effects of population or genotype).
In previous attempts, we had used AIC as a quick-and-dirty metric of model fit. The first look at the AIC values shows that we definitely need something about population and/or genotype included in the model (lower AIC is better, and a rule of thumb is that models > 10 AIC units higher than the best model are not worth further consideration)
And here is a close-up on all models that are not more than 10 AIC units higher than the best model:
The models with variation in intercept only (at both genotype and population levels) or with an intercept model for population and variation in nutrient effect across genotypes look reasonable.
Unfortunately, even though the model is rated as being the best, we still have perfect negative correlation between the variation in intercept (growth at low nutrient) and the effect of increased nutrient acorss genotypes:
For what it’s worth, including or excluding the random effect of nutrient at the genotype level has very little effect on the estimates of the fixed effects:
In either case, we can conclude that nutrient and clipping have opposite effects (positive and negative, respectively) on the number of fruits, with nutrients being stronger than clipping; that the experimental design variables (status and rack) do make a difference: and that nutrients ameliorate the effect of clipping, although not quite significantly.
The parametric bootstrap takes about 12 seconds per replicate, so you shouldn’t try to run the following code — the results are read in from the batch file.
(Don’t actually do this unless you want to spend a long time!)
Very close
…a slightly more sophisticated way to way view this, somewhat similar
to a q-q plot, is to plot the p-value for the parametric boostrap on the
x axis against the nominal value (from the χ12 distribution) on the y
axis.
That the
black (observed) line is below the (dark gray) 1:1 line indicates that the
standard likelihood ratio test is actually conservative in this particular
case (the nominal p-value is larger than the parametric bootstrap value).
The red dot shows the observed p-value for the data (nominal=0.063,
PB=0.073).
However, it’s a bit disappointing (if we are interested in the evolutionary consequences of variation across genotypes) that we don’t seem to have quite enough power to estimate covariation of traits across genotypes:
Exercise: Analyze one of the relevant data sets in
or a (small- to medium-size!) data set of your choice. The data must contain at least one predictor that can plausibly be treated as a random effect.