2
$\begingroup$

I have a linear mixed model structured like so:

Richness~Time*Treatment+(1|Site)

Time has two levels (Pre and Post) as does Treatment (Treatment and Control). Now, richness is count data and so I'm using a poisson distribution with a log link (residual plot looks fine).

I am predicting there will no difference between Pre Treatment and control samples (before a fire, in this case), but there will be a difference in the Post samples.

The model gives the following output:

FIXED EFFECTS:
                               Est.    S.E.   z val.       p   
(Intercept)                 0.40905 0.21488  1.90365 0.05696  .
TreatmentTreatment         -0.50263 0.18042 -2.78584 0.00534 **
TimePre                     0.29058 0.15135  1.91993 0.05487  .
TreatmentTreatment:TimePre  0.54499 0.22800  2.39028 0.01684  *

and doing a pairwise comparison with emmeans:

contrast                          estimate        SE  df z.ratio p.value
Control,Post - Treatment,Post   0.50262890 0.1804227 Inf   2.786  0.0274
Control,Post - Control,Pre     -0.29058074 0.1513496 Inf  -1.920  0.2196
Control,Post - Treatment,Pre   -0.33293869 0.1451375 Inf  -2.294  0.0994
Treatment,Post - Control,Pre   -0.79320964 0.1758399 Inf  -4.511  <.0001
Treatment,Post - Treatment,Pre -0.83556759 0.1705221 Inf  -4.900  <.0001
Control,Pre - Treatment,Pre    -0.04235795 0.1393992 Inf  -0.304  0.9903

Or to do the contrasts differently:

Time = Post:
contrast               estimate        SE  df z.ratio p.value
Control - Treatment  0.50262890 0.1804227 Inf   2.786  0.0053

Time = Pre:
contrast               estimate        SE  df z.ratio p.value
Control - Treatment -0.04235795 0.1393992 Inf  -0.304  0.7612

So far so good. There is clearly no difference before the fire but there is a difference after the fire.

Now, the problem. It's always better to display actual data rather than modelled data (at least in my opinion), and to do that I plotted the log response ratio (using escalc(measure="ROM") from the metafor R package). I averaged richness across sites, calculated SD and used the number of sites as n. I then calculated the difference between Treatment and Control sites, for Pre and Post separately, with 95% CIs: enter image description here

This does not agree with the statistical results above - it shows no difference in pre or post, or at best, a difference in the Pre samples.

I suspect the problem could lie in at least two places. The first is perhaps the most important.

1) Differences between sites are not accounted for in the ROM. Is there a way to weight the sites' richness prior to calculating the overall mean and SD for use in the ROM?

2) I only have a small number of sites (n=5), which can bias the ROM, and there is a fix for that, but I cannot find a way to implement it in R: https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/14-2402.1

Any assistance or insights appreciated, or if this is a duplicate.

$\endgroup$
4
  • $\begingroup$ Why not try to use the effects package to get the effects you are interested in? See summary(allEffects(model)) for example. The contrasts you tested are comparing log expected means via their differences, so it would make sense to compute the expected means themselves and visualize them to get a correspondence between the contrasts and the expected means. $\endgroup$ Sep 28, 2018 at 15:01
  • $\begingroup$ You can use plot(emmeans(...)) to display confidence intervals. If you want to plot the data, just use a scatter plot. You can add the CIs to that plot by using as.data.frame(emmeans(...)) as a data source. $\endgroup$
    – Russ Lenth
    Sep 29, 2018 at 1:53
  • $\begingroup$ Better yet, see the example here: cran.r-project.org/web/packages/emmeans/vignettes/… for a way to add the data to the interaction-plot-style display produced by emmeans $\endgroup$
    – Russ Lenth
    Sep 29, 2018 at 2:21
  • $\begingroup$ Thanks all. I will plot the estimated means comparisons, using plot(emmobject$contrasts). It seems too complicated to use the ROM. $\endgroup$ Sep 29, 2018 at 2:38

1 Answer 1

3
$\begingroup$

There are two extra points:

  1. Because you’re using a mixed effects Poisson regression with a log (nonlinear) link function, the coefficients you obtain have an intepretation conditional on the random effect. I.e., an interpretation conditional on site. This seems to be different than the plot you produce where you average over sites. Hence, it is a bit as comparing apples with oranges.
  2. Another potential reason for the differences is missing data. Namely, the mixed model will provide you correct inferences when the missing data mechanism is missing at random. Whereas, typically, plots on the observed data are only valid under missing completely at random.
$\endgroup$
1
  • $\begingroup$ Yes, thinking about it, my question is essentially asking how to model mean values for the ROM, which is kind-of pointless - I might as well stick to the estimated means. I wonder though, is there a way to modify the calcuation of confidence intervals to reflect a poisson log-link distribution? That is possibly the main problem - the confidence intervals are assuming a normal distribution, since it is calculated from ordinary standard deviation. Richness never follows a normal distribution, so even if the control and treatment have equal variance, the CIs rely on incorrect assuptions. $\endgroup$ Sep 29, 2018 at 19:38

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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