[R-sig-ME] problem extracting main effects from an updated glmm

[Fox, John]

Dear Ben and Mariano,

Ben, thanks for figuring out the source of the error before I even had a chance to read Mariano's original message. I think that we can just adopt your fix in the effects package. I'm cc'ing Sandy Weisberg, who wrote mer.to.glm(), in case he sees a problem that I don't.

Best,
 John

> -----Original Message-----
> From: Ben Bolker [mailto:bbolker at gmail.com]
> Sent: Tuesday, June 7, 2016 2:22 PM
> To: r-sig-mixed-models at r-project.org
> Cc: Fox, John <jfox at mcmaster.ca>
> Subject: Re: [R-sig-ME] problem extracting main effects from an updated
> glmm
> 
> 
>   Thanks for presenting a clear reproducible example.
> 
>  Depending on exactly what you need, there are other ways to "extract
> the main effects"  (e.g. fixef(fitted_model),
> coef(summary(fitted_model)), broom::tidy(fitted_model),
> dotwhisker::dwplot(fitted_model)), but I can appreciate the convenience
> of the effects package.
> 
>   The effects package is having trouble because one of its internal
> functions is trying to call glm() with the arguments previously provided
> to glmer, and the format of the 'start' parameter for glmer() is
> inconsistent with the one glm() is expecting (the
> 
>   Here's a crude way to hack the effects package so that it will work
> ...
> 
> ## dump the text of the offending function to a file
> cat(deparse(effects:::mer.to.glm,width.cutoff=200),
>     file="mer.to.glm.R",sep="\n")
> ## NOW EDIT THE FILE IN AN EXTERNAL TEXT EDITOR:
> ## (1) add "mer.to.glm <-" at the top
> ## (2) edit line 32; remove  "start", from argument list to be matched
> ## (3) l. 35, change fixmod -> effects:::fixmod ## read the file back in
> source("mer.to.glm.R")
> ## shove it back into the package namespace
> assignInNamespace("mer.to.glm",mer.to.glm,"effects")
> effects:::mer.to.glm ## print the function to check that it worked
> 
>   Having done that, allEffects(M2) works.
> 
>   This should be pretty easily fixable in the next release of the
> effects package, so you wouldn't have to do the hacking any more.
> 
>   cheers
>     Ben Bolker
> 
> 
> 
> 
> On 16-06-07 01:31 PM, Mariano Devoto wrote:
> > Dear all,
> >
> > we are trying to use a glmm to analyze sampling completeness of
> > interactions in ecological networks based on a literature review of
> > several studies each of which simultaneously sampled a given number of
> networks.
> > Our response variable is the number of interactions observed (Sobs)
> > compared to the number missed (Smiss).
> > Our explanatory variables are:
> > int_type: type of interaction (categorical)
> > n_webs: number of networks studied at the same time in each study
> > (numerical)
> > n_int: number of interactions sampled in each network (numerical)
> > sp_num: number of species in each network We have random effects
> > associated to each study ("datum") which are nested within each
> > research group ("author"). We hope this will account for inherent
> > methodological differences between research groups and uncontrolled
> > ecological background noise for each study.
> > After some exploratory analysis (following protocols in Zuur's books)
> > we also included in the model an interaction term ("int_type:sp_num").
> >
> > After much reading (books, blogs and other online help) this is what
> > we've managed to put together. As this is the first time we are using
> > glmm in our research group, in addition to the particular problem we
> > are having with extracting the model's main effects (see below) we
> > would greatly appreciate any comments/suggestions to improve our
> > general approach. So here's the
> > code:
> >
> > #read data online
> > my.table <- read.csv(file = "
> > http://www.agro.uba.ar//users//mdevoto//mydata.csv")
> >
> > #Initial glmm
> > require(lme4)
> > M0 <- glmer(cbind(Sobs, Smiss) ~ int_type + n_webs + n_int + sp_num +
> > int_type:sp_num + (1 | author/datum), data = my.table, family =
> > binomial)
> >
> > #Rescale explanatory variables after the function suggests to do so
> > my.table$n_webs.center <- scale(my.table$n_webs, center = T, scale =
> > T) my.table$n_int.center <- scale(my.table$n_int, center = T, scale =
> > T) my.table$sp_num.center <- scale(my.table$sp_num, center = T, scale
> > = T)
> >
> > #New model with rescaled variables
> > M1 <- glmer(cbind(Sobs, Smiss) ~ int_type + n_webs.center +
> > n_int.center + sp_num.center + int_type:sp_num.center + (1 |
> > author/datum), data = my.table, family = binomial)
> > summary(M1)
> >
> > # As there seem to be convergence problems, we followed this tutorial
> > to deal with them:
> > https://rstudio-pubs-static.s3.amazonaws.com/33653_57fc7b8e5d484c909b6
> > 15d8633c01d51.html #Restarting the model seems to solve things ss <-
> > getME(M1, c("theta", "fixef"))
> > M2 <- update(M1, start = ss, control = glmerControl(optCtrl =
> > list(maxfun =
> > 20000)))
> > summary(M2)
> >
> > #When we try to extract the main effects is when we run into trouble
> > require(effects)
> > my.effects <- allEffects(M2)
> >
> > We looked extensively online, but can't find a solution to get beyond
> > this point. Any ideas would be most welcome.
> >
> > Best wishes,
> >
> > Mariano
> >
> >
> > *Dr. Mariano Devoto*
> >
> > Profesor Adjunto - Cátedra de Botánica General, Facultad de Agronomía
> > de la UBA Investigador Adjunto del CONICET
> >
> > Av. San Martín 4453 - C1417DSE - C. A. de Buenos Aires - Argentina
> > +5411 4524-8069
> > http://www.agro.uba.ar/users/mdevoto/
> >
> > 	[[alternative HTML version deleted]]
> >
> > _______________________________________________
> > R-sig-mixed-models at r-project.org mailing list
> > https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models
> >