[R-sig-ME] Plotting residuals for GLMER model and zero counts

[Ben Bolker]

  At present this is only in the development (1.1-1) version, which is
on Github.  If you have development tools (compiler etc.) installed, then

library(devtools)
install_packages("lme4",user="lme4")

  If not, there may be binaries at
http://lme4.r-forge.r-project.org/repos (I don't remember ...)

  cheers
    Ben


On 13-11-14 02:30 PM, Malcolm Fairbrother wrote:
> Dear Ben,
> 
> I was looking at the nifty code you sent around on 6 Nov (see below),
> but on my system I'm getting stuck at:
> 
>> s <- simulate(form[-2],  ## use one-sided formula
> +     newdata=d,
> +     newparams=params,
> +     family=binomial)
> Error in UseMethod("simulate") : 
>   no applicable method for 'simulate' applied to an object of class
> "formula"
> 
> I just checked, and I think I've got the latest version of lme4
> (lme4_1.0-5, right?). I presume I'm missing something simple...?
> 
> Thanks,
> Malcolm
> 
> 
> 
> 
> 
>     Date: Wed, 6 Nov 2013 21:33:20 +0000 (UTC)
>     From: Ben Bolker <bbolker at gmail.com <mailto:bbolker at gmail.com>>
>     To: r-sig-mixed-models at r-project.org
>     <mailto:r-sig-mixed-models at r-project.org>
>     Subject: Re: [R-sig-ME] Plotting residuals for GLMER model and zero
>             counts
>     Message-ID: <loom.20131106T221303-818 at post.gmane.org
>     <mailto:loom.20131106T221303-818 at post.gmane.org>>
>     Content-Type: text/plain; charset=us-ascii
> 
>     Francesco <fbromano at ...> writes:
> 
>     > Dear R-ers,
>     >
>     > although a number of options exist out there to plot logit models, I
>     > can't seem
>     > to find one that works for glmer. What I would like to do is, in
>     light of
>     > having found no interaction between the two fixed effects (which
>     answers my
>     > research question), look at a plot to tell how well the model fits
>     the data.
>     >
>     > The package LMERconveniencefunctions ver 2.0
>     > has a nice plot for lmer but this won't work on my model because R
>     > automatically
>     > asks me to use glmer. It has a factor DV (Correct), two fixed
>     factor IVs
>     > (Group and Syntax),
>     > and two random effect (ID and item).
>     >
>     > Here's the output of calling:
>     >
>     > model<- glmer(Correct ~ Group * Syntax + (Syntax + 1 | ID) +
>     (Group + 1
>     > | item), data=...., family=binomial)
> 
>     I was going to ask you to please generate a (small) reproducible example
>     (e.g. see http://tinyurl.com/reproducible-000 ), but I couldn't resist
>     the urge to show off a new feature of the development branch of lme4,
>     which is that you can use simulate() with a regular (g)lmer formula
>     to simulate data corresponding to a given mixed model ...
> 
>     library(lme4)
>     set.seed(101)
>     ## generate factorial combinations
>     d <- expand.grid(Group=c("ns","nns"),Syntax=c("'s","of"),
>           ID =factor(1:38),item=factor(1:16))
>     ## subsample randomly down to actual data size
>     d <- d[sample(nrow(d),size=451),]
>     ## check ...
>     with(d,table(interaction(Group,Syntax)))
>     ## parameters (approximated from your output)
>     params <- list(fixef=c(-9,-2,5,-1.5),theta=c(7,-1,7,3,-0.5,6))
>     form <- Correct ~ Group * Syntax + (Syntax + 1 | ID) + (Group + 1| item)
>     s <- simulate(form[-2],  ## use one-sided formula
>         newdata=d,
>         newparams=params,
>         family=binomial)
>     ## turn response into a factor
>     d$Correct <- factor(s[[1]],labels=c("InC","Cor"))
>     fit1 <- glmer(Correct ~ Group * Syntax +
>         (Syntax + 1 | ID) + (Group + 1| item),
>         data=d,
>         family=binomial)
>     dAug <- data.frame(d,res=residuals(fit1,"pearson"))
>     library(lattice)
>     plot(fit1)  ## the default plot is not too useful
>     ## or
>     bwplot(res~Group:Syntax,data=dAug)
>     plot(fit1,Group:Syntax~resid(.))
>     ## how bad does the fit get if we leave out the interaction?
>     fit2 <- update(fit1,.~.-Group:Syntax)
>     plot(fit2,Group:Syntax~resid(.))
>     anova(fit1,fit2)
>     ## you can also separate by item etc. pretty easily ...
>     plot(fit2,Group:Syntax~resid(.)|item)
>