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

[Ben Bolker]

  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_57fc7b8e5d484c909b615d8633c01d51.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
>