2
$\begingroup$

I have run the following glmm:

mod<-glmer(data=newdata, total_flr_vis ~ treatment + flr_num + (1|individual), family=poisson)

and get this output from summary()

Generalized linear mixed model fit by maximum likelihood (Laplace Approximation) ['glmerMod']
Family: poisson  ( log )
Formula: total_flr_vis ~ treatment + flr_num + (1 | individual)
Data: newdata

 AIC      BIC   logLik deviance df.resid 
706.8    718.9   -349.4    698.8      148 

Scaled residuals: 
Min      1Q  Median      3Q     Max 
-3.1461 -0.8155 -0.3522 -0.2022 14.1669 

Random effects:
Groups     Name        Variance Std.Dev.
individual (Intercept) 4.066    2.016   
Number of obs: 152, groups:  individual, 42

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -1.610176   0.432666  -3.722 0.000198 ***
treatmentR  -0.457492   0.121054  -3.779 0.000157 ***
flr_num      0.037063   0.005064   7.319  2.5e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Correlation of Fixed Effects:
             (Intr) trtmnR
  treatmentR -0.117       
  flr_num    -0.416  0.040

But I get the following plot for treatment using plot(allEffects(mod))

Effects Plots

I don't understand why the effects plot shows overlapping error bars while the summary() output tells me that the effect of treatment is significant. Is there a problem with the model, or is it the plot? How can I troubleshoot this?

Here is the residual plot (which I got using plot(mod))

Residual plot

I'm not totally sure how to interpret this plot, but it does not look random to me, thus I suspect that there is something wrong with the model.

I am happy to post data if someone can tell me where to do that.

Note: I posted this on Stack Overflow and someone suggested that I am misinterpreting the first figure. This is very likely, but, unfortunately, the user did not suggest an alternative way to interpret the figure. If anyone here can do so, I would be very grateful.

Any help would be very welcome.

$\endgroup$
2
$\begingroup$

Your model has a random effect for individual which effects the marginal distribution of your observations when split into treatment groups, which is what you see in the first plot (I think, because you don't describe the data.) In the reproducible example below, there is a very obvious effect of the factor 'x', but it there is a lot of overlap in the 'allEffects' plot because this plot seems to average over the grouping factor.

The short answer is that effective blocking reduces error variance.

effect <- 3
ngroups <- 4
nper <- 5
total <- ngroups*nper
groups <- rnorm(ngroups, 0, 3)
offset <- 2
library(dplyr)
df <- data.frame(x = as.factor(rep(c(0, 1), each=ngroups/2*nper)),
             groupid = as.factor(rep(1:ngroups, times=nper))) %>%
  mutate(y = rpois(total, exp(offset + effect*(x==1) + groups[groupid] + rnorm(total))))
library(lme4)
library(effects)
fit <- glmer(y ~ x + (1|groupid), family=poisson("log"), data=df)
summary(fit)
plot(allEffects(fit))

library(ggplot2)
ggplot(df, aes(x=groupid, y=y, fill=x)) + geom_boxplot(position="dodge") + 
scale_y_continuous(trans="log")

enter image description here enter image description here

$\endgroup$
  • $\begingroup$ Thank very much @Eric Mittman. I used individual as a random effect to account for the fact that this was a paired design. Each individual received both treatments simultaneously on different stalks. Should I specify the model differently, or simply plot as you suggest? $\endgroup$ – JKO Apr 13 '17 at 15:10
  • $\begingroup$ I don't have reason to think there is a problem with the plot, provided that it is interpreted correctly. I'm not sure what exactly you want to convey; if it is just that the standard error for the difference is smaller than the estimate, I don't think a plot is necessary. Something like the second plot of mine could be persuasive. I can't speak to the appropriateness of your model; the one that I fit to the dummy data here matches exactly the true data generating process without the overdispersion ("rnorm(total)", which is easy to deal with; see: stats.stackexchange.com/a/9670/87365). $\endgroup$ – HStamper Apr 13 '17 at 16:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.