Skip to content

problem extracting main effects from an updated glmm

3 messages · Mariano Devoto, Ben Bolker, John Fox

Original
Mariano Devoto
# 
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/
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:
John Fox
# 
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