>=
pb1 <- PBrefdist(mp4,mp5)
@
Once we have the null distribution, doing inference
with it is quick.
<<>>=
PBmodcomp(mp4,mp5,ref=pb1)
@
This compares the likelihood ratio and
parametric bootstrap tests: the parametric
bootstrap test is slightly more conservative,
but both tests have $0.05>=
overdisp_fun(mp5)
@
If we did \code{locscaleplot(mp5)} we would see
very much the same pattern as in Figure~\ref{fig:locscale}.
In addition, we can check the normality of the random effects and find they are reasonable (Fig. \ref{fig:reN}).
We use \code{ldply} from the \code{reshape} package to collapse the list of random effect values into a data
frame (we might have to do something different if there were more than one random effect within each level,
e.g. a model including \code{(nutrient|gen)}). The fancy \code{panel} code in the figure adds a reference
line to the Q-Q plot: see \code{?qqmath}.
<>=
reStack <- ldply(ranef(mp5))
print( qqmath( ~`(Intercept)`|.id, data=reStack, scales=list(relation="free"),
prepanel = prepanel.qqmathline,
panel = function(x, ...) {
panel.qqmathline(x, ...)
panel.qqmath(x, ...)
},
layout=c(3,1)))
@
\begin{figure}[ht]
\centering
\includegraphics[width=.6\linewidth]{reN}
\caption{Quantile-quantile plot of the genotype random effects of the final model.}
\label{fig:reN}
\end{figure}
\section{Conclusions}
Our final model includes fixed effects of nutrients
and clipping, as well as the nuisance variables \code{rack} and \code{status};
observation-level random effects to account for overdispersion;
and variation in overall fruit set at the population and genotype levels.
However, we don't (apparently) have quite enough information to estimate the
variation in clipping and nutrient effects, or their interaction,
at the genotype or population levels. There is a strong overall positive effect
of nutrients and a slightly weaker negative effect of clipping.
The interaction between
clipping and nutrients is only weakly supported (i.e. the $p$-value is not
very small), but it is positive and about the same magnitude as the
clipping effect, which is consistent
with the statement that ``nutrients cancel out the effect of herbivory''.
To test this equivalence formally, we could set the analysis up as a one-way
layout (with levels [\code{Nutrient1:Unclipped},
\code{Nutrient1:Clipped},
\code{Nutrient8:Clipped},
\code{Nutrient8:Unclipped}]); the contrast between
\code{Nutrient1:Unclipped} and \code{Nutrient8:Clipped} would then
test the hypothesis (see
\url{http://glmm.wikidot.com/local--files/examples/culcita_glmm.pdf} for a
similar example).
This concludes part~1. To see the comparisons of the same analysis with
a variety of different estimation approaches using different R packages, see part~2.
Save results to use in part~2:
<<>>=
save("dat.tf","mp4",file="Banta_part1.RData")
@
\bibliographystyle{ESA1009}
\bibliography{glmm}
\end{document}