2
$\begingroup$

I have ecological count data across 15 sampling sites. The count data has a Poisson distribution and I have included an offset for "area_searched", which is the m^2 of each sampling site as it varies slightly. I have repeated these measures over 3 time periods (Session_ID). I want to know how I can compare the counts between each of the sessions to each other rather than just the first session. Is it possible to use the emmeans package to do this? If this is not correct, is there another method I can use to compare them between one another? I initially ran this model with the function manyglm() from the mvabund package but it appears that the emmeans package cannot handle this model - is this correct?

My model is:

model <- glm(count ~ Session_ID + offset(log(area_searched)), family = "poisson", data = countdata)

If I do summary(model)

I am only able to compare 3 to 1 and Session 2 to Session 1. I get the following output:

    Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -8.2622     0.1796 -46.002  < 2e-16 ***
Session_ID2   0.2134     0.2304   0.926 0.354383    
Session_ID3  -1.0929     0.3304  -3.308 0.000941 ***

I have used the emmeans package with the following code but I am not sure if it is correct to use with offset() in the model.

A <- emmeans(model, ~Session_ID)
pairs(A)

I get the following output:

 contrast estimate    SE  df z.ratio p.value
 1 - 2      -0.213 0.230 Inf -0.926  0.6237 
 1 - 3       1.093 0.330 Inf  3.308  0.0027 
 2 - 3       1.306 0.313 Inf  4.178  0.0001 

Results are given on the log (not the response) scale. 
P value adjustment: tukey method for comparing a family of 3 estimates

Many thanks!

Improve this question
$\endgroup$
0

2 Answers 2

Reset to default
3
$\begingroup$

The default behavior for emmeans() and its relatives is to use the model predictions, which includes any offset. It seems to be you don't want that in this situation. Fortunately, it is easy to override because there is also an offset argument in the emmeans() function.

B <- emmeans(model, ~Session_ID, offset = 0)
pairs(B)

This will use the same offset in all predictions. An offset of 0 is the same as log(1), i.e. the rate per unit area.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you very much for the response. I tried 'emmeans' with and without 'offset = 0' and it did not change the outcome but I assume that is because my offset values are quite similar to one another for each row of data. Would you be able to explain why I wouldn't want the offset included in this case when doing the Tukey's post-hoc test? $\endgroup$
    – catabolic
    Commented Sep 7, 2020 at 1:16
  • $\begingroup$ If you include the offset, you are comparing predictions. If you don't, you are comparing rates. You can see exactly what offsets are used by listing A@grid; it will have an .offset. column and those will be the average offset used for each prediction. $\endgroup$
    – Russ Lenth
    Commented Sep 7, 2020 at 13:39
0
$\begingroup$

You can do pairwise comparisons between your factor levels in mvabund using the below code:

pairwise_comparisons <- anova.manyglm(mod_name, pairwise.comp = df$factor) pairwise_comparisons

Or in your case:

pairwise_comparisons <- anova.manyglm(model, pairwise.comp = countdata$Session_ID) pairwise_comparisons

Improve this answer
$\endgroup$

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.