1

背景：我有计数数据（甲虫计数），我正在研究处理梯度对计数数据的影响。梯度是一个连续的预测变量，由“7个级别”组成（即-100%减少、-80%减少、-60%减少、-40%减少、-20%减少、0%减少和50%增加） . '0% 减少'意味着没有变化，或者这就是控制。我想使用 GLM 输出将处理“-60% 减少”（例如）与“0% 减少”进行比较。

如何在 R 中使用具有泊松分布和日志链接的 GLMM 输出来计算“减少 -60%”和“减少 0%”之间计数数据的百分比变化？

这是模型的示例：

glmmTMB(count_data ~ continuous_predictor + (1|random_effect),
        family=poisson(link=log), data=data)
地块编号 治疗 甲虫数
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
4

1 回答 1

1

让我们首先让您的示例可重现：

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))

您使用提供的数据描述的模型如下所示：

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

这意味着如果我们希望估计beetle_count给定值的treatment，我们需要计算exp(2.366465 + 0.015117 * treatment)。请注意，当处理为 0 时，这将简化为exp(2.366465)10.65964。对于 -60 的处理值，该值为exp(2.366465 + 0.015117 * -60)4.30357。

所以预期计数从 10.65964 下降到 4.30357，这意味着下降百分比是

100 * ((10.65964 - 4.30357) / 10.65964)
#> [1] 59.62744

这几乎是 60%

如果您想探索治疗水平之间的百分比差异（我们称它们为treatment_A和treatment_B），公式简化为

100 * (1 - exp(0.015117)^(treatment_A - treatment_B))
于 2022-02-01T10:32:42.237 回答