# Log response ratio of means not agreeing with linear mixed model output

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:

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.

• 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. Sep 28, 2018 at 15:01
• 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. Sep 29, 2018 at 1:53
• 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 Sep 29, 2018 at 2:21
• Thanks all. I will plot the estimated means comparisons, using plot(emmobject\$contrasts). It seems too complicated to use the ROM. Sep 29, 2018 at 2:38