[R-sig-ME] problem extracting main effects from an updated glmm
Lenth, Russell V
russell-lenth at uiowa.edu
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/
More information about the R-sig-mixed-models
mailing list