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

Thu Jun 9 18:38:51 CEST 2016

You may like the ease with which the lsmeans package can handle this example. There is an 'lstrends' function that estimates slopes for each level of a factor interacting with a covariate; and a 'contrast' function for estimating contarsts thereof. In another situation where more factors are involved, the results may differ somewhat from those in the effects package because I believe the latter weights the estimates by frequency, rather than equally, across levels of uninvolved factors.

> require("lsmeans")

> M2.lst <- lstrends(M2, "int_type", var = "sp_num.center")

> M2.lst
 int_type         sp_num.center.trend         SE df  asymp.LCL  asymp.UCL
 host-parasitoid          -0.67576500 0.14786963 NA -0.9655842 -0.3859459
 plant-disperser           0.10575349 0.35028629 NA -0.5807950  0.7923020
 plant-herbivore          -0.09630775 0.20541014 NA -0.4989042  0.3062887
 plant-pollinator         -0.26938881 0.07695205 NA -0.4202121 -0.1185656

Confidence level used: 0.95 

> contrast(M2.lst, "eff")   ### Sum-to-zero effects
 contrast                   estimate        SE df z.ratio p.value
 host-parasitoid effect  -0.44183798 0.1456253 NA  -3.034  0.0097
 plant-disperser effect   0.33968051 0.2715796 NA   1.251  0.4220
 plant-herbivore effect   0.13761927 0.1769490 NA   0.778  0.5823
 plant-pollinator effect -0.03546179 0.1245606 NA  -0.285  0.7759

P value adjustment: fdr method for 4 tests 

> contrast(M2.lst, "pairwise")   ### Pairwise comparisons
 contrast                             estimate        SE df z.ratio p.value
 host-parasitoid - plant-disperser  -0.7815185 0.3793521 NA  -2.060  0.1663
 host-parasitoid - plant-herbivore  -0.5794572 0.2341922 NA  -2.474  0.0640
 host-parasitoid - plant-pollinator -0.4063762 0.1669175 NA  -2.435  0.0707
 plant-disperser - plant-herbivore   0.2020612 0.4056479 NA   0.498  0.9595
 plant-disperser - plant-pollinator  0.3751423 0.3582628 NA   1.047  0.7216
 plant-herbivore - plant-pollinator  0.1730811 0.2182572 NA   0.793  0.8577

P value adjustment: tukey method for comparing a family of 4 estimates

--
Russell V. Lenth  -  Professor Emeritus
Department of Statistics and Actuarial Science   
The University of Iowa  -  Iowa City, IA 52242  USA   
Voice (319)335-0712 (Dept. office)  -  FAX (319)335-3017

Just because you have numbers, that doesn't necessarily mean you have data.

-----Original Message-----

Message: 1
Date: Tue, 7 Jun 2016 14:31:10 -0300
From: Mariano Devoto <mdevoto at agro.uba.ar>
To: r-sig-mixed-models at r-project.org
Subject: [R-sig-ME] problem extracting main effects from an updated
	glmm
Message-ID:
	<CAJRWjBZd_tL0dUN+5rMgWAeR4Ri2jUzdWVd_2mf_bWnizdJ2oA at mail.gmail.com>
Content-Type: text/plain; charset="UTF-8"

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/