\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}