Background: I have count data (beetle count) and I am looking at the effects of a gradient of a treatment on the count data. The gradient is a continuous predictor variable that consists of "7 levels" (i.e., -100% reduction, -80% reduction, -60% reduction, -40% reduction, -20% reduction, 0% reduction and 50% addition). The '0% reduction' means no change, or that is the control. I would like to compare the treatment '-60% reduction' (for example) to '0% reduction' using the GLM output.
How can I use the GLMM output with poisson distribution and log link in R to calculate the % change in count data between '-60% reduction' and '0% reduction'?
This is a sample of the model:
glmmTMB(count_data ~ continuous_predictor + (1|random_effect),
family=poisson(link=log), data=data)
| plot number | treatment | beetle count |
|---|---|---|
| 1 | -60 | 4 |
| 2 | -20 | 13 |
| 3 | 0 | 23 |
| 4 | -100 | 2 |
| 5 | 50 | 10 |
| 6 | -80 | 3 |
| 7 | -40 | 5 |
| 8 | 0 | 14 |
| 9 | -20 | 9 |
| 10 | -60 | 7 |
| 11 | -100 | 1 |
| 12 | -40 | 2 |
Let's make your example reproducible first:
library(glmmTMB)
data <- structure(list(
plot_number = 1:12,
treatment = c(-60L, -20L, 0L, -100L, 50L, -80L,
-40L, 0L, -20L, -60L, -100L, -40L),
beetle_count = c(4L, 13L, 23L, 2L, 10L, 3L, 5L, 14L,
9L, 7L, 1L, 2L)),
class = "data.frame", row.names = c(NA, -12L))
The model you describe with the data you provided looks like this:
model <- glmmTMB(beetle_count ~ treatment + (1|plot_number),
family = poisson(link = log),
data = data)
summary(model)
#> Family: poisson ( log )
#> Formula: beetle_count ~ treatment + (1 | plot_number)
#> Data: data
#>
#> AIC BIC logLik deviance df.resid
#> 68.4 69.8 -31.2 62.4 9
#>
#> Random effects:
#>
#> Conditional model:
#> Groups Name Variance Std.Dev.
#> plot_number (Intercept) 0.1703 0.4127
#> Number of obs: 12, groups: plot_number, 12
#>
#> Conditional model:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 2.366465 0.201081 11.769 < 2e-16 ***
#> treatment 0.015117 0.004148 3.645 0.000268 ***
#> ---
#> Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
What this means is that if we wish to estimate the beetle_count for a given value of treatment, we need to calculate exp(2.366465 + 0.015117 * treatment). Note that when treatment is 0, this simplifies to exp(2.366465) or 10.65964. For a treatment value of -60, the value is exp(2.366465 + 0.015117 * -60) or 4.30357.
So the expected count has dropped from 10.65964 to 4.30357, which means the percentage drop is
100 * ((10.65964 - 4.30357) / 10.65964)
#> [1] 59.62744
Which is almost bang-on 60%
If you want to explore the percentage difference between treatment levels (let's call them treatment_A and treatment_B), the formula simplifies to
100 * (1 - exp(0.015117)^(treatment_A - treatment_B))