# Post Hoc Test of interaction factor in binomial glmm with proportions

I've performed a binomial glmm, because my data are proportions of a species in a sample of +/- 100 Individuals. I test the interaction of two factors and use car::Anova to get the p-values. My random factor is the ID of the field (subject) that was sampled. I sampled 12 fields from 4 different classes (factor1), therefore i think the levels of the random factor is 12. The reason I use the random factor is that i want to correct for repeated measurements. Each field was measured only two times. These two timepoints are the levels of factor 2.

my model:

binomial.glmm <- glmer(cbind(species1, not species1) ~ factor1*factor2 + (1|field), family=binomial(link="logit"), data)


On the one side, I'm using the glht() function from the multcomp package to perform a post-Hoc Tukey test with bonferroni adjustment (all pairwise comparisons). And on the other side, I addtionally plot fitted values with confidence intervals.

My problem ist, that A: the a,b,c,... letters from the post hoc test do not make sense in my opinion. And B: The confidence intervals are incredibly large.

Now I'm thinking that the umber of real replicates (3) in each class is too small? Could that be the reason? Am i siply unable to test for the interaction effect?

I get the confidence intervals with the following R code:

testdata = expand.grid(factor1=unique(data$factor1), factor2 = unique(data$factor2))

X <- model.matrix(~ factor1*factor2, data = testdata)
testdata$fit <- X %*% fixef(binomial.glmm) testdata$SE <- sqrt(  diag(X %*%vcov(binomial.glmm) %*% t(X))  )
testdata$upr=testdata$fit+1.96*testdata$SE testdata$lwr=testdata$fit-1.96*testdata$SE


Then i plot the fit, upr and lwr, after calculating the exp(), using ggplot:

ggplot(testdata, aes(x = factor1, y = exp(fit))) +
#geom_bar(stat="identity",position = position_dodge(1), col="454545", size=0.15, fill="grey") +
geom_point(aes(x=as.numeric(factor1)+0.3),pch=23, bg="aquamarine2") +
geom_errorbar(aes(x=as.numeric(factor1)+0.3, ymin = exp(lwr), ymax = exp(upr)),position = position_dodge(1),col="black",width=0.15, size=0.15) +
geom_boxplot(aes(y=response), data=data) +
facet_grid(.~factor2) +
geom_hline(xintercept = 1, size=0.15) +
ylab("Species1?") +
xlab("Factor1") +
scale_x_discrete(labels=c("A", "B", "C", "D"))


Here are the plots:  Results from the post Hoc tukey Test, unfortuneately i wasn't able to make a pretty plot. The letters can't be true!!??

Predicitions plots. The blue points represent the fit with confidence interval, The boxplots are obtained from the real data

Interestingly, the prediction plots get much better, when i use a lmm with arcsine transformed fractions of Species1 as response. However everybody argues that this would not be gould practice,...

The residual plots to validate the models are a bit misbehaved, although i included the random term (1|ID) with ID = 1:nrow(data). The problem is, that if this does not improve the model, nothing could improve the model and at least there is only one single package (nparLD) that is able to perform non-parametric test with interactions with data that have repeated measurements. And i I'm afraid that procedure would become too difficult to explain in a matherial and methods section.

• You have a logit model, so the correct back-transformation is not the exponential function! ("logit" is not the same as "log") – Russ Lenth Nov 10 '15 at 17:34

You might try the lsmeans package, as it makes some of this stuff easier and clearer. To get the results on the back-transformed scale, do

require(lsmeans)
lsm <- lsmeans(binomial.glmm, ~ factor1 * factor2)
summary(lsm, type = "response")


To see the results graphically, do

plot(lsm, by = "factor2", intervals = TRUE, type = "response")


or, for an interaction-plot style,

lsmip(lsm, factor1 ~ factor2, type = "response")


In a mixed model like this, it is often the case that the SEs of these least-squares means (AKA predictions) will be much larger than those of some or all of the pairwise differences, because the between-subjects variations cancel out in those comparisons. So the displayed CIs can be very misleading for comparing the predictions (and you shouldn't use CIs to do comparisons in any case).

To get the Tukey-adjusted comparisons, do

summary(pairs(lsm), type = "response")


(This actually computes the differences on the logit scale, then back-transforms, so that the results are odds ratios. If you want differences of proportions instead, do pairs(regrid(lsm)) instead.)

• Thank you very much for the Answer. I could cry for joy, since month I try to figure out how to analyse my Data. Again i make the experience that i gain the most from Answers like yours. I'm no mathematician, so i would like to ask you if you could recommend some books, homepages, blogs, tutorials that can help to learn that stuff autodidactically? Since I searched the furthest corners of the internet to get the code I have used, but haven't heard about lsmeans before. If you got time, could you take look at my questions, i posted in the Answer section beneath? – Pharcyde Nov 11 '15 at 8:25

I use the "Answer" option, since the comments leave me with too less space.

I'm afraid, i need to have some points cleared:

1. Then I believe the low number of replicates is not really a problem in GLMM. And the low number of two repeated measurements do not influence the reliability of the random factor, if enough subjects were measured.

2. Could you please explain what you mean by "because the between subject variation cancel" has it sth to do with the "between" and "withnín" vocabulary known from repeated measures ANOVA?

3. If i should not use CIs, then this is also true for the plot(lsm) graphics? Hence i do not compare the factors with regard to significant differences. I always thought that if two CIs do not overlap, this would be the best indictor of significance. And I noticed, that the values in the plot(lsm) graphics, are similar to the ranges of the proportions of my data.

4. Can i use your code with poisson models, too?

5. I used the lmerTest package with your code for "lmer()"/gaussian lmms and I think the plots tell me which coefficients are significantly different from zero. I searched the internet and found many posts that argue that in the gaussian lmm case the multcomp::glht() function is appropriate for multicomparisons.

5.1 Comment: I used your code and found that lsmeans does not work for "glmer()" - models if the lmerTest package is loaded!

Thank you so much for your worthful time!!

• 1. Low number of reps is ALWAYS a problem. 3. Don't compare CIs ever. Use results of actual tests of pairwise comparisons. 4. Yes. 5. True. If it's loaded, use lsmeans::lsmeans(...) – Russ Lenth Nov 11 '15 at 16:26
• 2. I don't know what's true of your data, but when there is substantial variation among subjects, then you get better resolution comparing treatments given to the same subjects than you get comparing treatments given to different subjects. – Russ Lenth Nov 11 '15 at 16:29
• I've got longitudinal data from an ecological field stuy, my subjects are fields of different age. Since I can't manipulate the fields, i can't have the same "treatment" (=age of the field) on every field. Since I also use time (repeated measurements) as fixed factor, I always wondered if it is necessary to use the field(subject)-ID as random factor, at all. – Pharcyde Nov 12 '15 at 7:44
• Yes, I think you probably keep field in there as a random factor. Age is between-fields, and time is within-fields. – Russ Lenth Nov 12 '15 at 13:51