R - 泊松回归

R 语言中的现代泊松回归

泊松回归（Poisson Regression）是一种广义线性模型（GLM），用于响应变量表示计数数据的情况。这包括每天的客户投诉数量、某个路口的交通事故数量或某个栖息地中发现的物种数量等。计数必须是非负整数（0, 1, 2, …）。

核心假设是响应变量服从泊松分布（Poisson distribution）。该模型预测预期计数的对数是预测变量的线性组合：

log(E[y]) = β₀ + β₁x₁ + β₂x₂ + … + βₙxₙ

y 是计数响应变量。
E[y] 是预期（平均）计数。
β 值是模型系数。
x 值是预测变量。

`glm()` 函数

在 R 语言中，拟合泊松回归的主要函数是 glm()。

语法

glm(formula, data, family = poisson)

formula：定义了关系，例如 counts ~ predictor1 + predictor2。
data：包含变量的数据框。
family：指定误差分布和连接函数。对于泊松回归，我们使用 family = poisson，这意味着使用对数连接函数。

示例：分析织机中的经纱断裂

我们将使用内置的 warpbreaks 数据集。它包含基于 wool（羊毛）类型（A 或 B）和 tension（张力）水平（L、M、H）的织造过程中纱线断裂数量的数据。我们的目标是建模羊毛类型和张力如何影响断裂数量。

# 加载用于数据处理和模型解释的现代库
# install.packages(c("tidyverse", "broom"))
library(tidyverse)
library(broom)

# 检查数据
data(warpbreaks)
glimpse(warpbreaks)

glimpse() 输出显示了我们的变量：

Rows: 54
Columns: 3
$ breaks  <dbl> 26, 30, 54, 25, 70, 52, 51, 26, 67, 18, ...
$ wool    <fct> A, A, A, A, A, A, A, A, A, B, ...
$ tension <fct> L, L, L, L, L, L, M, M, M, L, ...

拟合和解释回归模型

# 拟合泊松广义线性模型。该公式将 'breaks' 建模为 'wool' 和 'tension' 的函数。
poisson_model <- glm(breaks ~ wool + tension, data = warpbreaks, family = poisson)

# 获取模型的综合摘要
summary(poisson_model)

摘要输出是理解模型的关键：

Call:
glm(formula = breaks ~ wool + tension, family = poisson, data = warpbreaks)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-3.6871  -1.6503  -0.4269   1.1902   4.2616

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  3.69196    0.04541  81.302  < 2e-16 ***
woolB       -0.20599    0.05157  -3.994 6.49e-05 ***
tensionM    -0.32132    0.06027  -5.332 9.73e-08 ***
tensionH    -0.51849    0.06396  -8.107 5.21e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for poisson family taken to be 1)

    Null deviance: 297.37  on 53  degrees of freedom
Residual deviance: 210.39  on 50  degrees of freedom
AIC: 493.06

Number of Fisher Scoring iterations: 4

解释系数

在泊松模型中，Estimate 列显示了预期计数的对数的变化。为了使其更直观，我们必须对系数进行指数化：exp(Estimate)。

基线：(Intercept) 对应于基线组：wool ‘A’ 和 tension ‘L’（每个因子的第一个水平）。该组的断裂对数计数为 3.69。
woolB：系数为 -0.206。这意味着在张力不变的情况下，羊毛 B 的断裂对数计数比羊毛 A 低 0.206。乘法效应是 exp(-0.206) = 0.814。因此，羊毛 B 的断裂数量比羊毛 A 减少了约 18.6%。
tensionM 和 tensionH：类似地，中等张力和高张力与低张力相比，断裂数量更少。例如，高张力（tensionH）与 exp(-0.518) = 0.596 相关联，即比低张力减少了约 40.4% 的断裂。
显著性：所有预测变量的 p 值（Pr(>|z|)）都很低，这表明羊毛类型和张力水平都对断裂数量具有统计学上显著的影响。

使用 `broom` 包进行现代解释

使用 broom::tidy() 简化了提取和使用这些系数的过程。

# 整理模型输出
tidy_model <- tidy(poisson_model)

# 添加指数化系数以便于解释
tidy_model %>%
  mutate(rate_ratio = exp(estimate))

# 输出：
# A tibble: 4 × 6
#  term        estimate std.error statistic      p.value rate_ratio
#  <chr>          <dbl>     <dbl>     <dbl>        <dbl>      <dbl>
#1 (Intercept)    3.69     0.0454     81.3  1.33e-164        40.1
#2 woolB         -0.206    0.0516     -3.99 6.49e-  5         0.814
#3 tensionM      -0.321    0.0603     -5.33 9.73e-  8         0.725
#4 tensionH      -0.518    0.0640     -8.11 5.21e- 16         0.595

实际考量与最佳实践

过度分散（Overdispersion）：泊松模型的一个关键假设是计数数据的均值等于其方差。然而，在实际情况中，方差常常大于均值（这种情况称为过度分散）。您可以通过检查残差离差（Residual Deviance）是否远大于自由度来发现这一点（在我们的例子中，210.39 >> 50，这表明存在过度分散）。
替代模型：如果存在过度分散，**负二项回归（Negative Binomial Regression）**是更好的选择。您可以使用 MASS 包中的 family = "negbin" 来拟合它（glm.nb() 是一个方便的封装函数）。
暴露量（Exposure）/偏移量（Offset）：如果计数是在不同时期收集的（例如，每月事故数 vs. 每年事故数），您应该在模型中包含一个暴露量或偏移量项来规范化比率。例如：glm(counts ~ predictor + offset(log(time)), family = poisson)。