[R-sig-ME] query on binomial glm of pollination experiment

Mollie Brooks mo|||eebrook@ @end|ng |rom gm@||@com
Thu Feb 20 21:12:34 CET 2020


Hi Mariano,

I don’t have much time right now, but I was able to improve the model (including convergence) with a few small changes. I think that having nofruits as a separate object may have been causing a problem in the inner-workings of the software. So I moved it to a column of sistrep. I also reordered the treatment levels so that control is the baseline. Then the model appears to run in lme4 with a warning or in glmmTMB without a warning. It looks like some of the coefficients differ between the two models.

sistrep$nofruits <- sistrep$flowers-sistrep$fruits 
sistrep$treatment=factor(sistrep$treatment, levels=c("control", "automan", "cruzman", "tul", "voile"))

library(glmmTMB)
M3 <- glmmTMB(cbind(fruits, nofruits) ~ treatment + (1|plant), data =
sistrep, family = "binomial")
summary(M3)

It could be possible that having better gradient information allows glmmTMB to converge in this case.

BTW, I couldn’t get the getURL command to work, but I downloaded the csv from the web address and read it in without any problems. This may help in case anyone else tries to help.

cheers,
Mollie

> On 20Feb 2020, at 20:38, Mariano Devoto <mdevoto using agro.uba.ar> wrote:
> 
> Dear all,
> 
> I am trying to assess the effect of five pollination treatments on a
> plant's fruit set. The data set is not too good (few samples and unbalanced
> -sadly we lost many field samples to vandalism-, quite a few zeros). I
> counted the n of flowers that received each treatment and the number of
> fruits at the end of the experiment. I am using a binomial model.
> Furthermore, given that several (though not all) plants received all the
> treatments I'd like to include a random factor to account for between-plant
> variation. I tried models with and without a random factor but results
> either look awkward (particularly post hoc comparisons in model M1) or the
> model does not converge (M2).
> 
> I'll be grateful for any advice on how to reliably test the significance of
> the treatments and do multiple pairwise comparisons.
> 
> Thanks in advance for your help!
> 
> Here is the code to load the data and perform the analysis. The data set is
> read from a Google sheet.
> 
> #loading data and curating data
> require(RCurl); require(lme4); require(multcomp)
> my.file <- getURL("
> https://docs.google.com/spreadsheets/d/e/2PACX-1vSu7V5v00cpQ_SFZq5ajS-RUPCgRYf0q249SW3rh9g2q87h7awYOsgBX4xIRDGWQff4Jy6IRnSLuaxH/pub?output=csv
> ")
> sistrep <- read.csv(textConnection(my.file), head=T)
> str(sistrep)
> nofruits <- sistrep$flowers-sistrep$fruits #calculates n of failed
> pollinations
> 
> #plot the data
> prop <- sistrep$fruits/sistrep$flowers
> boxplot(prop ~ treatment, data=sistrep, las=1, ylab="Fruit set (n fruits/n
> flowers)", boxwex=0.6)
> 
> #and here are the models & post hoc comparisons
> M0 <- glm(cbind(fruits, nofruits) ~ 1, data = sistrep, family = "binomial")
> M1 <- glm(cbind(fruits, nofruits) ~ treatment, data = sistrep, family =
> "binomial")
> M2 <- glmer(cbind(fruits, nofruits) ~ treatment + (1|plant), data =
> sistrep, family = "binomial")
> summary(M0)
> summary(M1)
> anova(M0,M1)
> summary(M2)
> 
> summary(glht(M1, mcp(treatment = "Tukey", interaction_average=F)))
> ---------------------------------------------
> 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 5287-0069
> *https://www.researchgate.net/profile/Mariano_Devoto
> <https://www.researchgate.net/profile/Mariano_Devoto>*
> 
> 	[[alternative HTML version deleted]]
> 
> _______________________________________________
> R-sig-mixed-models using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-sig-mixed-models



More information about the R-sig-mixed-models mailing list