\documentclass{article}
\usepackage{natbib}
\usepackage{sober}
\usepackage{color}
\usepackage{hyperref}
\usepackage{url}
\usepackage{listings}
\usepackage[utf8]{inputenc}
\newcommand{\code}[1]{{\tt #1}}
\title{GLMM on symbiont effects on coral predation}
\author{Ben Bolker}
\date{\today}
\begin{document}
\bibliographystyle{ESA1009}
\newcommand{\todo}[1]{{\color{red} \emph{#1}}}
\maketitle
\SweaveOpts{height=5,width=10}
\setkeys{Gin}{width=\textwidth}
<>=
options(SweaveHooks=list( fig=function() par(mfrow=c(1,2)) ))
@
<>=
options(SweaveHooks=list( fig=function() par(mar=c(5,11,4,9)+0.1)) )
@
\section{Preliminaries}
The purpose of this document is to explore/explain some
of the nitty-gritty details of fitting, and making inferences
from, GLMMs in R. It may turn into a ``real paper'' at some
point; we'll see.
The data are from Adrian Stier and Sea McKeon, and represent
trials of \emph{Culcita} sea stars attacking coral
that harbor differing combinations of protective symbionts (crabs
and shrimp). The design is a randomized complete block design
with some replication (2 replications per treatment per block,
4 treatments (no symbionts, crabs alone, shrimp alone,
both crabs and shrimp), each of these units of 8 repeated in
10 blocks). The data represent a reasonably typical small-scale
ecological field experiment.
Some issues that do \emph{not} arise in this example:
\begin{itemize}
\item extremely large data sets (requiring concerns for
computational efficiency and stability)
\item extremely small (5--6 or less) numbers of random
effect units that have to be dealt with in some special
way (treatment as fixed effects or supply of a Bayesian
prior \citep{gelman_prior_2006})
\item crossed random effects
\item overdispersion
\item spatial or temporal correlations that decrease
with separation in space or time
\item other interesting correlation structures (phylogenetic?)
\item lack of balance
\end{itemize}
<>=
source("glmmfuns.R")
@
The basic data can be reduced, for the purposes
of this exercise, to a single treatment (\code{ttt})
[which consists of combinations of different
symbionts: crab, shrimp, both or neither, but never
mind]; a binary response (\code{predation});
and a blocking factor (\code{block}).
Data entry/summary:
<>=
x=read.csv("culcitalogreg.csv",
colClasses=c(rep("factor",2),
"numeric",
rep("factor",6)))
## abbreviate slightly
levels(x$ttt.1) <- c("none","crabs","shrimp","both")
@
Adjust contrasts for the treatment,
to compare (1) no-symbiont vs symbiont cases,
(2) crabs vs shrimp, (3) effects of a single pair/type
of symbionts vs effects of having both:
<>=
contrasts(x$ttt)=
matrix(c(3,-1,-1,-1,
0,1,-1,0,
0,1,1,-2),
nrow=4,
dimnames=list(c("none","C","S","CS"),
c("symb","C.vs.S","twosymb")))
@
\section{Plots}
Two ways to plot the data:
1. A basic barplot that loses information on
blocks but gives the distribution of numbers
of predation events per block:
<>=
library(reshape)
m <- melt(x[,c(1,4,3)],id.vars=1:2)
m3 <- recast(m,ttt.1+block~variable,fun.aggregate=sum)
p <- with(m3,table(predation,ttt.1))
op <- par(las=1,cex=1.5,bty="l")
bb <- barplot(p,ylab="Number of blocks",xlab="Symbionts",
main="Number of predation events",
col=gray(c(1,0.8,0.5)))
## the rest of this is fanciness to plot the values
## within the bar segments (rather than just
## a plain old legend)
bloc <- apply(p,2,cumsum)
midb <- rbind(bloc[1,]/2,
(bloc[1,]+bloc[2,])/2,
(bloc[2,]+bloc[3,])/2)
text(bb[col(p)][p>0],midb[p>0],c(0,1,2)[row(p)][p>0])
par(op)
@
2. A (too?) clever strip plot that shows
the pattern for each block:
<>=
library(lattice)
ctab <- with(x,as.data.frame.table(tapply(predation,list(block,ttt.1),sum)))
names(ctab) <- c("block","ttt","pred")
ctab$rblock <- with(ctab,
reorder(block,pred))
ctab$jpred <- ctab$pred+seq(-0.1,0.1,length.out=10)[ctab$rblock]
print(stripplot(jpred~ttt,groups=block,data=ctab,
type=c("p","l"),ylab="Predation events by block",
scales=list(y=list(tick.number=3))))
@
\todo{A simpler plot, or one based on model fits?}
\section{Basic fits}
\todo{should I try to test the \code{ttt:block} interaction?
Probably will be a mess since within-block/treatment replication
is so limited (two binary samples), but worth it for completeness?}
\subsection{glm (base R)}
Fit without blocking (i.e., no random effects):
<>=
mod1=glm(predation~ttt,binomial,data=x)
@
\subsection{glmer (lme4)}
Fitting with the default Laplace approximation and
with adaptive Gauss-Hermite quadrature (abbreviated
sometimes as AGQ and sometimes as (A)GHQ):
<>=
library(lme4)
## Laplace
mod2 <- glmer(predation~ttt+(1|block),family=binomial,data=x)
## AGQ
mod2B <- glmer(predation~ttt+(1|block),family=binomial,
nAGQ=8,data=x)
## mod2C <- glmer(predation~ttt+(ttt|block),family=binomial,
## nAGQ=8,data=x)
@
\subsection{glmmML (glmmML)}
The \code{glmmML} function is in a separate package
(the \code{glmmML} package, surprisingly enough).
It handles only random intercept models with a single
blocking factor
(I chose \code{nAGQ=8} above to match the default
number of AGQ points (8) for \code{glmmML} with
\code{method="ghq"}).
<>=
library(glmmML)
## Laplace (default)
mod3=glmmML(predation~ttt,family=binomial,data=x,cluster=block)
## AGQ (=GHQ) with 8 points
mod3B=glmmML(predation~ttt,family=binomial,data=x,cluster=block,
method="ghq")
@
Fitting with Laplace gives a warning:
\code{Hessian non-positive definite. No variance!}
This means it can fit the model but fails to get
(approximate) standard errors for the parameters.
We'll decide later whether we want to worry about this.
\subsection{glmm (repeated)}
Jim Lindsey's \code{repeated} package (available
not from CRAN, but from \url{http://www.commanster.eu/rcode.html}) contains the \code{glmm} function:
<>=
library(repeated)
mod7=glmm(cbind(predation,1-predation)~ttt,family=binomial,
data=x,nest=block,points=8)
mod7.ci = confint(mod7) ## not reliable!
@
Have to specify predation as $(k,N-k)$, \code{glmm}
apparently doesn't recognize binary data as such.
\code{glmm} has a \code{confint} method, which is supposed
to compute profile confidence intervals, but it gives
ridiculously optimistic results, so we'll fall back
on the estimated standard errors for CIs below.
\subsection{glmmADMB}
<<>>=
library(glmmADMB)
@
<>=
mod4=glmm.admb(predation~ttt,random=~1,group="block",
data=x,family="binomial",link="logit")
@
<>=
mod4B=glmm.admb(predation~ttt,random=~1,group="block",
data=x,family="binomial",link="logit",impSamp=2)
@
Can supposedly improve the Laplace fit with importance
sampling, but I didn't really understand the meaning of
the {\tt impSamp} parameter from the
documentation, and got errors when
I tried setting it $>0$
(after step 2: \code{The function maximizer failed}).
(This is mentioned, but not discussed much in, the \code{glmm.admb}
help page: probably need to see ref therein:
\cite{skaug_automatic_2004}, which
has a line about importance sampling, further referencing
\cite{skaug_automatic_2002}.
\subsection{MCMCglmm}
Would have tried \code{MCMCglmm}, but it is
impossible to run \code{MCMCglmm} models \emph{without}
individual-level random effects --- which in this case
(with binary data) causes convergence problems.
<>=
library(MCMCglmm)
mod5 <- MCMCglmm(predation~ttt,random=~block,rcov=~units,
family="categorical",data=xaug,verbose=FALSE)
## coef
colMeans(mod5$Sol)
apply(mod5$Sol,2,posterior.mode)
colMeans(mod5$VCV)
## somewhat odd answers.
## check convergence etc etc etc
## cannot (apparently) fit WITHOUT residual variance term??? (rcov=0)
@
\subsection{glmmBUGS}
(Have to do this next bit carefully:
both of the next two approaches load the
\code{nlme} package (used by \code{glmmPQL}).
Since \code{nlme} and \code{lme4} conflict,
we have to unload/reload the packages.)
<>=
detach("package:lme4")
library(MASS)
library(nlme)
@
The \code{glmmBUGS} package is supposed to translate
a mixed-model statement into BUGS code which can then
be run via \code{R2WinBUGS}.
This doesn't quite work --- first because the binomial-family
specification doesn't get translated correctly (a
\code{dbin(x)} specification should read
\code{dbin(x,1)}), and then from an unknown BUGS
problem \ldots
<>=
library(glmmBUGS)
gg1 <- glmmBUGS(predation~ttt, data=x, effects="block",
family="binomial")
startingValues <- gg1$startingValues
source("getInits.R")
library(R2WinBUGS)
gg2 <- bugs(gg1$ragged,getInits,
parameters.to.save = names(getInits()),
working.directory=getwd(),
model.file="model.bug")
@
\subsection{glmmPQL}
Also worth trying with
{\tt glmmPQL} (which uses a less reliable algorithm,
but is more flexible etc.):
<>=
mod5=glmmPQL(predation~ttt,random=~1|block,family=binomial,data=x)
pqlfix <- fixef(mod5)
pqlfix.se <- summary(mod5)$tTable[,"Std.Error"]
pqlran <- as.numeric(VarCorr(mod5)["(Intercept)","StdDev"])
@
<>=
detach("package:nlme")
library(lme4)
@
\subsection{BUGS}
Or, we can simply write our own \code{BUGS} file:
\lstset{numbers=left,
numberstyle=\small,
basicstyle=\ttfamily}
\lstinputlisting{culcita_bugs.txt}
<<>>=
library(R2WinBUGS)
@
<>=
mm <- model.matrix(predation~ttt,data=x)
bugsdata <- list(mm=mm,N=nrow(x),
nblock=length(levels(x$block)),
block=as.numeric(x$block),
predation=x$predation)
bugsinit <- list(list(b=c(0.5,rep(0,3)),
eps = rep(0,10),
tau.block=1),
list(b=c(0.5,rep(-0.5,3)),
eps = rep(0,10),
tau.block=1),
list(b=c(0.5,rep(0.5,3)),
eps = rep(0,10),
tau.block=1))
b1 <- bugs(data=bugsdata,inits=bugsinit,
parameters=c("b","eps","tau.block","sd.block"),
model.file="culcita_bugs.txt",
n.chains=3,
n.iter=4000,
working.directory=getwd())
@
We can get 95\% credible intervals \ldots
<>=
library(emdbook)
b2 <- as.mcmc.list(b1)
b3 <- lump.mcmc.list(b2)
head(coda::HPDinterval(b3),4)
@
\ldots or examine the posterior densities (of
block effect expressed as precision ($\tau$) or
standard deviation ($\sigma$) \ldots
<>=
b4 <- lapply(b2,function(x) {x[,c("sd.block","tau.block")]})
class(b4) <- "mcmc.list"
print(densityplot(b4,from=0))
## z <- data.frame(do.call(rbind,b4),chain=rep(1:3,each=334))
## densityplot(~tau.block,groups=chain,data=z,
## panel=function(...) {
## panel.densityplot(...)
## panel.curve(dgamma(x,0.01,0.01),col="gray")
## },
## from=0)
@
\subsection{glmmAK}
The \code{glmmAK} package is a lesser-known alternative.
Its range of allowed models is small (logistic and
Poisson regression, plus multinomial/ordinal logistic
regression), and its interface is a bit quirky \ldots
I think the specification below gives a model that
is equivalent to the BUGS model above \ldots
<>=
library(glmmAK)
mod6 = cumlogitRE(y=x$predation,x=mm[,-1],intcpt.random=TRUE,
cluster=x$block,prior.fixed=list(mean=0,var=100),
prior.random=list(Ddistrib="gamma",
Dshape=0.01,DinvScale=100),
nsimul=list(niter=6000,nthin=40,nburn=2000))
@
<<>>=
mod6fixchain <- read.table("betaF.sim",header=TRUE)
mod6rechain <- read.table("varR.sim",header=TRUE)
b6 = as.mcmc(cbind(mod6fixchain,sqrt(mod6rechain[,1])))
@
\todo{SAS, ASREML/Genstat, others?}
\subsection{Comparisons}
(Ugly code suppressed here --- look at the original
Sweave file if you want the gory details on extracting
coefficients, standard errors, etc. from the model fits.)
<>=
fixefmat <- rbind(glm=coef(mod1),
glmmPQL=pqlfix,
glmer.Laplace=fixef(mod2),
glmmML.Laplace=coef(mod3),
glmmADMB.Laplace=mod4$b,
glmer.AGQ=fixef(mod2B),
glmmML.AGQ=coef(mod3B),
glmm_repeated.AGQ=coef(mod7)[1:4],
BUGS.MCMC=b1$mean$b,
glmmAK.MCMC=colMeans(b6[,1:4]))
sdvec <- c(0,pqlran,
attr(VarCorr(mod2)[[1]],"stddev"),
mod3$sigma,
sqrt(mod4$S),
attr(VarCorr(mod2B)[[1]],"stddev"),
mod3B$sigma,
coef(mod7)[5],
b1$mean$sd.block,
mean(b6[,5]))
(compmat <- round(cbind(fixefmat,sdran=sdvec),3))
@
Conclusions: ignoring the random effect
underestimates all of the fixed effect sizes
(this is more or less as expected); PQL underestimates
the standard deviation of the random effect
by about 50\% relative to Laplace/AGQ. Laplace
and AGQ give very similar answers:
the Laplace results are identical across all
3 packages tried.
AGQ makes a small difference,
but even the direction of change depends
on which implementation you use.
(I wouldn't say that any of these
differences are large enough to concern us in practice
in this case). The two MCMC-based methods, \code{BUGS}
and \code{glmmAK}, give substantially higher estimates
of the fixed and random effects (this is not due to
asymmetry of the posterior distributions/differences
between the mode and the mean of the posterior).
By analogy with randomized complete block designs
for LMMs, the ``denominator df'' for this problem
should be 27 = $(10-1)\times(4-1)$: see for example \url{http://www.tfrec.wsu.edu/ANOVA/RCBsub.html}.
Add confidence intervals, crudely based on Wald
intervals with 27 df (except for BUGS/\code{glmmAK} results,
which are credible intervals):
(\code{glmmADMB} gives a standard deviation of the
variance estimate, but $\hat \sigma^2=\Sexpr{round(mod4$S,1)}$
and its standard deviation is supposed to be
\Sexpr{round(mod4$sd_S,1)}, so I can't really use
it to construct sensible confidence intervals.)
<>=
fac <- qt(0.975,27)
fixef.low <- rbind(glm=coef(mod1)-fac*coef(summary(mod1))[,"Std. Error"],
glmmPQL=pqlfix-fac*pqlfix.se,
glmer.Laplace=fixef(mod2)-fac*
summary(mod2)@coefs[,"Std. Error"],
glmmML.Laplace=rep(NA,4),
glmmADMB.Laplace=mod4$b-fac*mod4$stdbeta,
glmer.AGQ=fixef(mod2B)-fac*
summary(mod2B)@coefs[,"Std. Error"],
glmmML.AGQ=coef(mod3B)-fac*mod3B$coef.sd,
glmm.AGQ=coef(mod7)[1:4]-fac*
coef(summary(mod7))[1:4,"Std. Error"],
BUGS.MCMC=coda::HPDinterval(b3)[1:4,"lower"],
glmmAK.MCMC=coda::HPDinterval(b6)[1:4,"lower"])
fixef.hi <- rbind(glm=coef(mod1)+fac*coef(summary(mod1))[,"Std. Error"],
glmmPQL=pqlfix+fac*pqlfix.se,
glmer.Laplace=fixef(mod2)+fac*
summary(mod2)@coefs[,"Std. Error"],
glmmML.Laplace=rep(NA,4),
glmmADMB.Laplace=mod4$b+fac*mod4$stdbeta,
glmer.AGQ=fixef(mod2B)+fac*
summary(mod2B)@coefs[,"Std. Error"],
glmmML.AGQ=coef(mod3B)+fac*mod3B$coef.sd,
glmm_repeated.AGQ=coef(mod7)[1:4]+fac*
coef(summary(mod7))[1:4,"Std. Error"],
BUGS.MCMC=coda::HPDinterval(b3)[1:4,"upper"],
glmmAK.MCMC=coda::HPDinterval(b6)[1:4,"upper"])
sdvec.low <- c(0,
NA,
NA,NA,
NA,
mod3B$sigma-fac*mod3B$sigma.sd,
coef(mod7)[5]-fac*coef(summary(mod7))[5,"Std. Error"],
coda::HPDinterval(b3)["sd.block","lower"],
coda::HPDinterval(b6)[5,"lower"])
sdvec.hi <- c(0,NA,
NA,NA,
NA,
NA,
mod3B$sigma+fac*mod3B$sigma.sd,
coef(mod7)[5]+fac*coef(summary(mod7))[5,"Std. Error"],
coda::HPDinterval(b3)["sd.block","upper"],
coda::HPDinterval(b6)[5,"upper"])
compmat.low <- cbind(fixef.low,sdran=sdvec.low)
compmat.hi <- cbind(fixef.hi,sdran=sdvec.hi)
@
<>=
m <- melt(compmat)
low <- melt(compmat.low)
hi <- melt(compmat.hi)
m$X1 <- factor(m$X1,levels=rownames(compmat))
## levels(m$X1) <- c("glm","PQL","Laplace (1)","Laplace (2)",
## "Laplace (3)","AGQ (1)","AGQ (2)","BUGS")
m$X2 <- factor(m$X2,levels=colnames(compmat))
m$low <- low$value
m$hi <- hi$value
## trellis.par.set(add.text=list(cex=0.8))
## print(xyplot(X1~value|X2,
## data=m,pch=16,
## ## scales=list(x=list(relation="free")),
## scales=list(y=list(alternating=FALSE)),
## xlab="",ylab="",
## panel=function(...) {
## panel.xyplot(...)
## panel.abline(v=0)
## }))
## trellis.par.set(theme=orig) ## restore settings
<>
library(ggplot2)
print(ggplot(m,aes(x=value,y=X1))+geom_point()+
geom_errorbarh(aes(xmin=low,xmax=hi))+
facet_wrap(~X2,scales="free_x")+geom_vline(x=0,col="red")+
scale_x_continuous("effect size",expand=c(0.1,0)))
<>
@
\section{Inference}
\subsection{Wald $Z$ tests}
These are given (for fixed effects) by both packages
as \verb+Pr(>|z|)+ in
the standard model print-out (they're missing for
glmmML/Laplace because there was a computational
problem) and suggest that there is an effect
of having any symbionts at all (\code{ttt1}),
but that crabs are not different from shrimp
(\code{ttt2}) and that having two symbionts
isn't any better than having one (\code{ttt3}).
Replicating the results for the
glmmML code here (this appears
in \code{summary(mod3B)}, but it
doesn't look like this function call
returns anything useful --- it just
prints the object in the usual way).
<>=
csd <- mod3B$coef.sd
z.table <- cbind(coef(mod3B),csd,
coef(mod3B)/csd,
pnorm(abs(coef(mod3B))/csd,
lower.tail=FALSE)*2)
colnames(z.table) <- c("coef","se(coef)","Z","Pr(>|Z|)")
round(z.table,4)
@
You can get equivalent fits from the
\code{glmer} fit (\verb+summary(mod2B)@coefs+)
or the \code{glmm} fit (\code{coef(summary(mod7))}).
The main problems
with this analysis are (1) it doesn't
take account of the residual degrees of
freedom/uncertainty in standard error estimates
(asymptotic resullts) and (2) it's a fairly crude,
Wald-based estimate --- not clear whether the null distribution
of the parameter estimates divided by their standard
errors is really normal $Z$ (or $t$) in this case
(3) it tests parameters (contrasts) one at a time.
(I've written a \code{waldF()} function to compute
Wald $F$ tests --- it can be found in \code{glmmfuns.R}).
\subsection{Wald $t$ tests}
One of the harder parts of GLMMs
(in the traditional Wald- or $F$-testing approach)
is figuring out the appropriate degrees of freedom.
As mentioned above,
I will use $(10-1)\times(4-1)=27$ df as the
denominator df (may check later that this
actually fits the null Wald distribution \ldots).
(For the same model,
\code{nlme}'s ``inner-outer'' approach
\citep{pinheiro_mixed-effects_2000} gives
37 df instead --- but this is only the difference
between a critical value of
\Sexpr{round(qt(0.975,27),3)} for 27 df
and \Sexpr{round(qt(0.975,37),3)} for 37 df,
vs. the magic 1.96 for a $Z$ test.)
None of the tests we did cross
the magic $p<0.05$ line in either
direction under these circumstances
(results are similar for \code{glmer}):
<>=
round(Z.to.t(z.table,27),3)
@
\code{Z.to.t()} is a small function to hack
$Z$ tables (provided e.g. by \code{glmer} and \code{glmm}
into $t$ tables, \emph{provided that you can specify
the appropriate ``denominator'' degrees of freedom}.
This is of course not using any of the fancy
technology (which Doug Bates distrusts) for
computing corrected degrees of freedom
(Satterthwaite or Kenward-Roger) --- but which
may not be necessary here in any case.
\subsection{LRT}
It is possible, but strongly advised against, to
use Likelihood Ratio tests to make inferences about
fixed effects \citep{pinheiro_mixed-effects_2000} --- %
it might be the recommended way to
make inferences about random effects.
(Fabian Scheipl's \code{RLRSim} package is probably
the best way to make inferences on significance
of random effects in simple LMMs, but it
handles neither GLMMs nor more complex (crossed,
multi-random-factor, etc.) LMMs.)
One problem with testing the random effects here
is that \code{glmer} can't fit models with no random effects
at all. In principle we can compare the log-likelihood
from glmer (with block effects) to that from glm
(without block effects). In practice we have to
be \emph{very careful} when comparing log-likelihoods
calculated by different functions, because they can
incorporate (or drop) different additive constants
(which makes no difference when we compare results
computed by the same function, but can screw us up
when we compare results computed by different functions).
Ways to check/test this: (1) look at the code for
both functions (ugh); (2) simulate data with a very
small or zero block effect, fit it with glmer and glm
(presumably getting a very small/zero estimate of the
random effect and hence similar log-likelihoods),
and compare; (3) set the variance
parameter in the glmer fit to zero and re-evaluate
(this is what I do below).
<>=
##source("glmmprof.R")
## source("glmersim-funs.R") ## needed for simulate.lmer
## get deviance of model at zero random effects
## watch out, this function may (because of an lme4 bug)
## actually **modify the original fitted object** (!)
zerodev <- function(mm) {
class(mm) <- class(mm) ## HACK: "touch" model object to force copy
varprof(mm,0,0,1)[["ML"]]
}
fixefsd <- function(mm) {
sqrt(diag(vcov(mm)))
}
@
Using the \code{varprof} function (too ugly
to show here --- defined in \code{glmmfuns.R})
which computes a deviance profile
by varying the value of the random-effects standard deviation,
we can show the profile.
<>=
profres <- varprof(mod2B)
@
<>=
sdest <- attr(VarCorr(mod2B)[[1]],"stddev")
sp1 <- with(subset(as.data.frame(profres),sdsdest),
splines::interpSpline(sd, ML, na.action=na.omit))
bsp2 <- splines::backSpline(sp2)
upper <- predict(bsp2,-2*logLik(mod2B)+qchisq(0.9,1))$y
@
\todo{since I'm not that interested in inference
on random effects in this case, this section
is a little dusty/could use some work}
In the figure below,
the solid line is the critical level for the LRT based on
a $\chi^2$ distribution with 1 df, the dashed line is for
$\chi^2_{1/2}$ (see ``boundary effects'' in the GLMM paper
for why $\chi^2_{1/2}$ may be more appropriate).
The $p$ value based on this distribution is
% would be: SEXPR
% omitted because of lme4 model-munging problems? fixed now?
\Sexpr{round(c(pchisq(zerodev(mod2B)-(-2*logLik(mod2B)),1,lower.tail=FALSE)/2),2)}.
<>=
plot(profres[,"sd"],profres[,"ML"],
xlab="random-effects SD",ylab="deviance")
crit1 <- -2*logLik(mod2B)+qchisq(0.95,1)
crithalf <- -2*logLik(mod2B)+qchisq(0.9,1)
abline(h=c(crit1,crithalf),lty=1:2)
text(0,crit1,expression(chi[1]^{2}),pos=3)
text(0,crithalf,expression(chi[1/2]^{2}),pos=1)
segments(lower,crithalf,upper,crithalf,lwd=2)
points(0,-2*logLik(mod1),pch=16,col=2)
@
The deviance calculated by glmer for $\sigma=0$
is \Sexpr{round(profres[1,"ML"],2)}; the
deviance calculated by glm is
\Sexpr{round(-2*logLik(mod1),2)}
[shown in red in the figure], so the
computed log-likelihoods do \emph{not} match ---
we were right to be careful. (That, or I made
some other mistake here.)
The 95\% profile confidence limits on the
random-effects sd, based on
$\chi^2_{1/2}$ (see below for why I think
this is the best answer) are
\{\Sexpr{round(lower,2)},\Sexpr{round(upper,2)}\}
(see ugly invisible code above,
taken out of one of the \code{mle} profiling functions,
for how to get these numbers).
(This value takes boundary effects into account,
but not finite sample size --- from that point
of view, it may be a little bit conservative.)
We can also get likelihood profiles
for the fixed effects. The basic recipe is
to leave one fixed effect parameter (the one
being profiled) out of the model;
use the \code{offset} argument to \code{glmer}
to set it to a particular (fixed) value; and
re-fit the model. For factors this is a little
bit trickier because it's hard to leave one
parameter out of the model. The way I ended
up doing this was a bit of a hack --- treating
the model as a regression analysis on the individual
columns of the model matrix (it would be easier
if one could pass a model matrix directly).
Right now the code is fairly ugly and fairly specific
to this problem, it could certainly be
cleaned up and generalized \ldots
Plots of the deviance profiles for each fixed-effect
parameter. Red lines are quadratic approximations
based on the local curvature, as used in the Wald $Z$
tests; horizontal lines are the LRT cutoff --- for
quadratic approximations, equivalent to the Wald $Z$
95\% critical value.
<>=
vals <- fixprof(mod2B)
@
<>=
op <- par(mfrow=c(2,2),mgp=c(2,1,0),mar=c(3,3,1,1))
mindev <- (-2*logLik(mod2B))
for (i in 1:4) {
with(vals[[i]],
{
plot(parval,ML,xlab=names(vals)[i],ylab="deviance",
ylim=c(65,78))
curve(((x-fixef(mod2B)[i])/fixefsd(mod2B)[i])^2+mindev,
add=TRUE,col=2)
abline(h=mindev+qchisq(0.95,1))
})
}
par(op)
@
The profile for the intercept ran into some
numerical difficulties beyond about 3.5 --- these
are simply not shown on the plot.
\todo{compare likelihood profile
confidence intervals
(based on asymptotic $\chi^2_1$, or on
a scaled $F_{1,27}$ or $F_{1,??}$ ?)
with Wald $Z$/$t$}
\subsection{AIC}
\todo{Finish this!}
<>=
## define logLik and AIC for glmmML objects
logLik.glmmML <- function(x) {
loglik <- (-x$deviance)/2
attr(loglik,"df") <- length(coef(x))
loglik
}
AIC.glmmML <- function(x) x$aic
@
Both \code{glmmML} and \code{glmer} fits will give AIC
values, although it is a bit tricky
to get them --- I had to write the extractor
methods for \code{glmmML} myself.
<>=
library(stats4)
sapply(list(mod2,mod3,mod2B,mod3B),AIC)
sapply(list(mod2,mod3,mod2B,mod3B),logLik)
@
Same issues about comparing log-likelihoods
from different functions apply here, too, but
with additional complication -- parameter counting
(as well as likelihood constants) can differ
among functions, e.g. whether one includes an
implicitly estimated variance term or not.
So you have to be careful. Once you have the
log-likelihoods, though, it is fairly easy to
figure out how a given function is counting --- easier
than working out the log-likelihood comparisons
(see above).
\section{Randomization/simulation approaches}
The idea here is to simulate from the null model
(zero block effect, or zero for some treatment effect
or combination of effects), refit the full
model, get some summary statistics -- the log-likelihood
or deviance (or Z or t statistic) -- and repeat.
For many (1000 or so) simulations these statistics
will give the null distribution, which you can then
compare against the observed values.
In general, this is fairly simple (not necessarily
even worth packing into/hiding behind a function):
fit the reduced model; use
\code{simulate} to generate new realizations from
the reduced model; for each realization,
use \code{refit} to refit
the full model and (possibly) the
reduced model to the data simulated from the
reduced model; and compute some sort of summary
statistics (deviance/log-likelihood difference,
$Z/t$ statistic, ?).
\subsection{Random effects}
\todo{Also a little dusty --- redo along the
lines of the fixed effect example below?}
For dropping the only random
effect from the model we can do almost the
same thing, but we have a little problem ---
we need to simulate from a
\code{glm} fit instead of a reduced \code{glmer} fit,
because \code{glmer} can't fit models with
no random effects.
(In this case I had to write \code{simulate.glm}
since it didn't exist: it's in \code{glmmprof.R}.)
As promised, the code is actually pretty
simple.
The extra complexities are:
\begin{itemize}
\item a little bit of code to load previously
computed results from a file/save the results
to a file (since this is computationally intensive);
\item because \code{glm} computes likelihood differently
from \code{glmer}, we use the \code{zerodev} function (defined above)
to return the value of the deviance when a
model is refitted with the random-effects standard
deviation forced equal to zero --- rather than
what we would usually do, refitting a \code{glm}
model to the new realization and computing the
deviance that way
\end{itemize}
<>=
if (file.exists("glmersim1.RData")) {
load("glmersim1.RData")
} else {
svals <- simulate(mod1,nsim=1000)
t0 <- system.time(sim1 <- apply(svals,2,
function(x) {
r <- refit(mod2,x)
zerodev(r)- deviance(r)
}))
save("svals","t0","sim1",file="glmersim1.RData")
## restore random numbers??
}
@
(this took \Sexpr{round(t0[1],1)} seconds on my laptop).
The following figure compares the $\chi^2_1$ distribution
with the actual distribution of deviance differences:
<>=
library(lattice)
q1 <- function(x) qchisq(x,1)
print(qqmath(~sim1,distribution=q1,
prepanel = prepanel.qqmathline,
panel = function(x, distribution, ...) {
panel.qqmathline(x,distribution=distribution, ...)
panel.qqmath(x,distribution=distribution,...)
}))
@
The next figure compares the
observed null distribution to the $\chi^2_{1/2}$ distribution.
(The quantiles of this distribution are equivalent
to $\chi^2_1(1-2\alpha)$, e.g. the 95\% quantile
is equal to the 90\% quantile for $\chi^2_1$.
The \code{qchibarsq} function from the \code{emdbook}
package is a convenience function that implements
this distribution.)
The solid line, which looks really
bad, intersects the (0.25,0.75) quantiles ---
this distribution is extremely skewed, so most
of the points actually fall in the lower left-hand
corner. The dashed line (admittedly cherry-picked)
intersects the (0.25,0.95) quantiles and seems OK
(\emph{this varies by realization --- it looked better
for a previous set of simulation results \ldots})
<>=
library(emdbook)
qhalf <- function(x) qchibarsq(x,1)
print(qqmath(~sim1,distribution=qhalf,
prepanel = prepanel.qqmathline,
panel = function(x, distribution, ...) {
panel.qqmathline(x,distribution=distribution, ...)
panel.qqmathline(x,distribution=distribution, p=c(0.25,0.95), lty=2,...)
panel.qqmath(x,distribution=distribution,...)
}))
@
The critical values/quantiles from the simulated and
the $\chi^2_{1/2}$ distribution:
<>=
cbind(sim=quantile(sim1,c(0.9,0.95,0.99)),
qhalf=qchibarsq(c(0.9,0.95,0.99)))
obsdev <- c(zerodev(mod2B)-(-2*logLik(mod2B)))
@
Estimated $p$ values for the random effect
from simulated and $\chi^2_{1/2}$
(both very very small):
<>=
pchibarsq(obsdev,lower.tail=FALSE)
mean(c(sim1,obsdev)>=obsdev)
@
This is essentially 1/1001 --- it is
general practice (ref???) to include
the observed value in the set of null-hypothesis
values, which means that the minimum
$p$-value achieved in this way is $1/(N+1)$.
\subsection{Fixed effects}
\section*{Permutation}
A basic function (maybe there is a better way?)
to permute predation events within blocks:
<>=
permfun <- function(x) {
r <- do.call("rbind",
lapply(split(x,x$block),
function(z) {
z$predation <- sample(z$predation)
z
}))
## to use refit(), we have to make sure that
## we reassemble into the proper order matching the original
## data -- better/cleaner way to do this?
r[order(r$ttt,as.numeric(as.character(r$block))),]
}
@
Check that it preserves the structures we want:
<>=
## experimental design --- ttt * block
with(permfun(x),table(block,ttt))
## total predation events
c(sum(x$predation),sum(permfun(x)$predation))
## predation events per block
rbind(tapply(x$predation,list(x$block),sum),
tapply(permfun(x)$predation,list(x$block),sum))
@
A function to fit the full and reduced model to a
new predation vector and return the difference in deviance and the Wald statistics:
<>=
statfun <- function(x,mod0,mod1) {
mod1 = try(refit2(mod1,x))
mod0 = try(refit2(mod0,x))
if (inherits(mod0,"try-error")) {
if (inherits(mod1,"try-error")) {
## both bad
return(rep(NA,5))
} else return(c(NA,summary(mod1)@coefs[,"z value"]))
}
likratio <- c(-2*(logLik(mod1)-logLik(mod0)))
wald <- summary(mod1)@coefs[,"z value"]
c(likratio,wald)
}
@
A function to do Monte Carlo simulation of
the predation events from the null model:
<>=
mcfun <- function(x) {
my.mer.sim(mod0)
}
@
Run permutations:
<>=
source("~/projects/glmm/glmmfuns.R")
permvals <- as.data.frame(t(replicate(2000,
statfun(permfun(x)$predation,
mod0,mod1))))
@
<>=
## junk
myrep <- function (n, expr, simplify = TRUE) {
pb <- txtProgressBar()
sapply(seq(n),
function(i,...) {
setTxtProgressBar(pb,i/n)
eval.parent(substitute(function(...) expr))
},
simplify = simplify)
}
@
Run MC simulations:
<>=
mcvals <- as.data.frame(t(replicate(4000,
statfun(mcfun(x,mod0,mod1)))))
@
<>=
load("culcita_mcvals_1.RData")
load("culcita_permvals_1.RData")
@
Even after we throw out fits with
warning messages there are 2 or 3 really
awful outliers, but we will just ignore
these (doesn't make too much difference
to $p$-values anyway):
<>=
ss <- sort(abs(mcvals$ML))
sn <- sign(mcvals$ML)[order(abs(mcvals$ML))]
plot(sort(abs(mcvals$ML)),col=sn+2)
@
Analyze MC results:
<>=
## see note at bottom of ?hist examples about
## comparing distributions with histograms ...
mod1 = glmer(predation~ttt2+(1|block),data=x,family=binomial)
mod0 = update(mod1,.~.-ttt2)
mcdev <- -mcvals$ML[abs(mcvals$ML)<50]
permdev <- na.omit(-permvals$ML)
obsdev <- c(2*(logLik(mod1)-logLik(mod0)))
@
<>=
op <- par(mar=c(5,4,2,5)+0.1)
hist(mcdev,breaks=200,col="gray",freq=FALSE,las=1,
xlab="Deviance",main="")
legend("top",
c("MC distrib",
"perm. distrib",
expression(chi[3]^2)),
fill=c("gray","red",NA),
lty=c(NA,NA,1),
col=c(NA,NA,"cyan"))
abline(v=obsdev,col=2,lwd=2)
curve(dchisq(x,3),col=5,lwd=2,add=TRUE,n=301)
par(new=TRUE)
hist(permdev,breaks=200,col="red",freq=FALSE,las=1,
axes=FALSE,xlim=c(0,20),xlab="",ylab="",main="")
axis(side=4,las=1)
@
% mean of F distribution -- n_2/(n_2-2) for n_2>2
% mean of chi-sq = n
% how do we match up chi-sq(num.df) with F(num.df,den.df) ?
% dist(chi-sq/n) = dchisq(x*n,n)*n
%
<>=
## playing around with nchisq/F equivalence
## distribution of (chi-sq/n)
dnchisq <- function(x,df,...) {
dchisq(x*df,df,...)*df
}
## distribution of (F*n_1)
dnf <- function(x,df1,df2,...) {
df(x/df1,df1,df2,...)/df1
}
pnf <- function(x,df1,df2,...) {
pf(x/df1,df1,df2,...)
}
@
<>=
curve(dnchisq(x,3))
curve(df(x,3,1000),add=TRUE,col=2)
curve(df(x,3,100),add=TRUE,col=3)
curve(df(x,3,20),add=TRUE,col=4)
curve(df(x,3,10),add=TRUE,col=5)
curve(df(x,3,5),add=TRUE,col=6)
curve(dchisq(x,3),from=0,to=8)
curve(dnf(x,3,100),add=TRUE,col=2)
curve(dnf(x,3,5),add=TRUE,col=6)
curve(pchisq(x,3),from=0,to=8)
curve(pnf(x,3,100),add=TRUE,col=2)
curve(pnf(x,3,5),add=TRUE,col=6)
@
The gray bars are the observed Monte Carlo
distribution for simulations of the null hypothesis ---
<>=
qqplot(mcdev, qchisq(ppoints(mcdev), df = 3),
type="l")
abline(0,1, col = 2)
@
Permutation distribution only has
\Sexpr{length(unique(round(permdev,7)))} unique
values (once we account for numeric errors
on the order of 1e-9 from fitting models with
data in different orders \ldots)
<<>>=
mean(obsdev<=pchisq(obsdev,3,lower=FALSE)) ## LRT: approx 8.3e-5
@
It turns out that the LRT $p$-value is
anticonservative in this finite-sample case, as previously reported
for typical small-sample LMMs by \cite{pinheiro_mixed-effects_2000}.
That is, the reported or nominal
$p$ value ($y$ axis) is generally below the Monte
Carlo $p$ value ($x$ axis) (or in other words, the
(black) $\chi^2_3$ line is always below the (gray) 1:1 line
in the figure below.
<>=
ndev <- length(mcdev)
matplot(1:ndev/(ndev+1),
cbind(pchisq(rev(sort(mcdev)),3,lower.tail=FALSE),
pnf(rev(sort(mcdev)),3,60,lower.tail=FALSE),
pnf(rev(sort(mcdev)),3,40,lower.tail=FALSE),
pnf(rev(sort(mcdev)),3,27,lower.tail=FALSE)),
type="s",
col=c(1,2,4,5),lty=1,
xlim=c(0.001,1),ylim=c(0.001,1),
log="xy",las=1,bty="l",
xlab="True p-value",ylab="")
mtext("Nominal p-value",side=2,line=3.5,xpd=NA)
abline(a=0,b=1,col="gray")
brk <- c(0.001,0.005,0.01,0.05,0.1)
abline(h=brk,v=brk,col="gray",lty=3)
legend("bottomright",c(expression(chi[3]^2),
"scaled F(3,60)",
"scaled F(3,40)",
"scaled F(3,27)"),
lty=1,col=c(1,2,4,5),y.intersp=1.25)
u <- par("usr")
@
Also drawn in this figure are the $F$ distributions
(scaled to match the $\chi^2_3$ distribution in the
large-denominator-df limit) for varying denominator
df. The line based on $df_2=27$ (expectation from
classic ANOVA tables). The lines for $df_2=40$ and
$df_2=60$ are closer, but none of them seem to match
precisely --- as Bates would say, we don't necessarily
have a reason to believe that the null distribution
of deviance differences is $F$ distributed in any case.
\todo{do the same kinds of checks for Wald statistics
and for likelihood ratio tests for dropping individual
contrasts}
\section{To do}
\begin{itemize}
\item more AIC stuff, model averaging?
\item incorporate profile confidence limits (backspline) in varprof
code? profiling code could be tested and polished quite a bit:
(1) better ways to manipulate formulas and model matrices;
(2) test on other models, more robust, work out what to do about
zero-intercept models, binomial responses with multiple columns,
etc.
(3) make profile objects match structure of glm/mle profile
objects, so we can use the code for plotting, confidence intervals
etc.
\item structured bootstrapping (by block) example?
(Can easily bootstrap on blocks: however, this will be
a little slower than resampling just the response vector;
we will need to refit every time. Sampling design will
be such that we will have unbalanced data --- fewer blocks,
some blocks overrepresented.) Start by splitting data by
block, sample blocks with replacement, replicate as necessary,
squash back together, fit \ldots
\item ASReML, SAS results ??? (MLWin, Stata, ???)
\item make glmmBUGS work?
\item can trace of hat matrix be resurrected??
\item work with simulated data of a similar size?
\end{itemize}
\bibliography{glmm}
\end{document}