22、深入理解泊松回归在SAS中的应用

深入理解泊松回归在SAS中的应用

1. 泊松回归基础与拟合优度评估

在处理计数数据时，泊松回归是一种常用的统计方法。以彩票数据为例，我们可以借助SAS来进行泊松回归分析。以下是拟合泊松回归模型后，用于评估拟合优度的部分指标：
| 评估标准 | 自由度（DF） | 值（Value） | 值/自由度（Value/DF） |
| — | — | — | — |
| 偏差（Deviance） | 16 | 29.4934 | 1.8433 |
| 缩放偏差（Scaled Deviance） | 16 | 29.4934 | 1.8433 |
| 皮尔逊卡方（Pearson Chi - Square） | 16 | 28.2003 | 1.7625 |
| 缩放皮尔逊卡方（Scaled Pearson X2） | 16 | 28.2003 | 1.7625 |
| 对数似然（Log Likelihood） | | 130.3048 | |

泊松回归通过最大化似然函数来估计参数，而非像线性回归那样使用最小二乘法。泊松回归的偏差统计量公式为：
[Deviance = 2\sum_{i}(y_{i}-\lambda_{i}) + y_{i}\log(\frac{\lambda_{i}}{y_{i}})]
其中，(y_{i}) 是观测到的泊松计数。

在泊松分布中，均值和方差通常相等。不过，我们可以通过设置尺度（scale）参数，来拟合方差大于或小于均值的泊松模型。尺度参数由程序员指定，SAS 不会对其进行估计。当指定的尺度参数不等于 1 时，卡方和偏差统计量以及估计参数值的所有标准误差都会被重新缩放。建议的尺度参数值可以从输出中的 “Value/df” 部分获取。

2. obstats 输出解读

在模型语句中指定 “obstats”（观测统计量）选项，能够生成一系列有用的值，包括残差、拟合值等。以下是彩票数据中前四个城镇的部分 “obstats” 输出：
| 观测值（Obs） | 获奖者（winners） | 面积（area） | 人口（pop） | 预测值（Pred） | Xbeta | 标准误差（Std） | Hessian 权重（Hesswgt） | 下限（Lower） | 上限（Upper） | 原始残差（Resraw） | 卡方残差（Reschi） | 偏差残差（Resdev） | 标准化偏差残差（Stresdev） | 标准化卡方残差（Streschi） | 似然残差（Reslik） |
| — | — | — | — | — | — | — | — | — | — | — | — | — | — | — | — |
| 1 | 6 | 6.2 | 17.9 | 6.0282 | 1.7964 | 0.1753 | 6.0282 | 4.27507 | 8.50052 | -0.0282 | -0.01152 | -0.01153 | -0.01278 | -0.01277 | -0.012778 |
| 2 | 11 | 27.9 | 28 | 6.9120 | 1.9332 | 0.0976 | 6.9120 | 5.70806 | 8.36995 | 4.0879 | 1.55490 | 1.43036 | 1.47996 | 1.60882 | 1.488804 |
| 3 | 6 | 33 | 26.2 | 6.4764 | 1.8681 | 0.1238 | 6.4764 | 5.08033 | 8.25625 | -0.4764 | -0.18722 | -0.18959 | -0.19977 | -0.19728 | -0.199532 |
| 4 | 2 | 17.2 | 12.8 | 5.0921 | 1.6277 | 0.1351 | 5.0921 | 3.90706 | 6.63675 | -3.0921 | -1.37029 | -1.56400 | -1.64226 | -1.43885 | -1.624414 |

这些信息可以通过输出语句或输出交付系统（ods）获取。输出语句是捕获 “obstats” 信息的一种方式，而 ods 更加灵活，可用于管理 SAS 过程的输出，包括开启或关闭输出的各个部分，以及保存输出以供后续步骤使用。在相关程序中，ods 语句会捕获这些信息，并创建一个名为 “fitted” 的新数据集，用于程序的后续步骤。

输出中的 “Pred” 表示每个城镇的估计均值 (\hat{\lambda}_{i})，“Xbeta” 是均值对数的估计值，通常被称为线性预测器，其估计标准误差列在 “Std” 列下。通过 “Xbeta” 和 “Std” 可以为每个估计均值创建 95% 的置信区间，计算公式为：
[exp(Xbeta \pm 1.96 \times Std)]

原始残差（Resraw）是观测值与期望值之间的差值：
[Resraw_{i} = y_{i}-\hat{\lambda} {i}]
不过，原始残差虽然直观，但并不是寻找异常值的最佳方法。因为在泊松分布中，均值和方差相等，不同的原始残差会有不同的方差。更好的泊松回归残差是将每个残差标准化后的结果，即卡方残差（Reschi），其定义为：
[Reschi {i} = \frac{y_{i}-\hat{\lambda} {i}}{\sqrt{\hat{\lambda} {i}}}]
这些残差更具可比性，因为它们的方差大致相同。它们被称为卡方残差，是因为每个残差都代表了对熟悉的卡方统计量的个体贡献：
[\chi^{2} = \sum_{i}(Reschi_{i})^{2} = \sum_{i}\frac{(y_{i}-\hat{\lambda} {i})^{2}}{\hat{\lambda} {i}}]

表中给出的其他统计量（Resdev、Stresdev、Streschi、Reslik）也是残差的不同定义，它们的值应该彼此相似，并且与卡方残差相似。Hessian 权重（Hesswgt）是每个观测值的估计方差，等于估计均值。

3. offset 选项的应用

假设我们想将彩票获奖者的平均数量建模为与城镇人口成正比。在其他条件相同的情况下，人口是原来两倍的城镇，我们预期彩票获奖者的数量也会是原来的两倍。在生命率的例子中，我们也希望讨论相对于风险人群的死亡率和发病率。那么在 SAS 中如何实现这一点呢？

如果人口是描述第 (i) 个城镇预期彩票获奖者数量 (\lambda_{i}) 的唯一因素，那么模型可以表示为：
[\lambda_{i} = C \times Population_{i}]
其中，(C) 对于所有城镇具有相同的值，类似于回归系数。

通过对数链接，模型变为：
[\log(\lambda_{i}) = \log(C) + \log(Population_{i})]
这个方程说明了拟合 “成正比” 模型时的问题。我们需要估计 (C)，但 “(\log(Population_{i}))” 项没有需要估计的回归系数，更准确地说，“(\log(Population_{i}))” 前面的回归系数始终等于 1。

以下是一个使用 “offset” 的 SAS 程序示例：

title2 'Fit an offset for mean proportional to population';
data propto;
    /* build a new dataset */
    set lottery;
    /* start with the data of lottery dataset */
    logpop = log(pop);
    /* add the new variable: log population */
run;
proc genmod data=propto;
    model winners = area mill / dist=Poisson obstats offset=logpop;
run;

在这个程序中，数据步骤基于早期的彩票数据创建了一个名为 “propto” 的新数据集。新数据集创建了一个名为 “logpop” 的新变量，它被定义为人口变量的自然对数。

在 “proc genmod” 的模型语句中，列出了因变量（获奖者）、解释变量（面积、磨坊率），指定了泊松分布（dist = Poisson），请求了观测统计量（obstats），并指定人口的对数作为偏移变量。与偏移变量相关的回归系数始终等于 1，因此在这个模型中不需要估计该系数。同样，SAS 也不会对这个变量进行统计显著性检验。

4. 练习部分

以下是一些相关的练习，通过这些练习可以进一步加深对泊松回归的理解：
1. 在均值为 1 的泊松分布中，求事件发生次数为 2 或更少的概率。在 (N = 4) 且 (p = 1/4) 的二项实验中，以及 (N = 8) 且 (p = 1/8) 的二项实验中，分别求事件发生次数为 2 或更少的概率。计算这两个二项实验的均值和方差，观察随着 (N) 增大和 (p) 减小，二项模型的方差是否更接近均值，2 或更少事件发生的概率是否更接近泊松近似。
2. 在航空公司超售航班导致乘客被拒载的例子中，某航班通常会预订多达 194 个座位，但飞机实际只有 185 个物理座位。根据经验，航空公司预计 5% 的乘客不会按时登机。判断泊松分布是否适合这个数据，并解释原因。估计该航班 194 名潜在乘客中 “未出现” 人数的泊松参数 (\lambda)。计算航空公司为该航班可以接受多少预订，并且仍然有 99% 的信心不需要拒载任何乘客。
3. 在之前的示例程序中，我们使用人口的对数作为偏移量，其中人口以千人为单位进行测量。如果以个人为单位测量人口，会发生什么情况？使用两种定义运行程序，并注意输出的差异，解释这些差异产生的原因。
4. 重新审视相关的死亡率数据，将住院死亡人数视为泊松分布的结果。解释为什么这是一个合理的假设，并检查拟合回归模型的残差是否与给定表格中最后一列显示的高低死亡率指示相对应。
5. 分析马萨诸塞州 1854 年各县的精神疾病病例数据。说明病例数量如何随人口或人口对数变化，解释为什么人口密度是这个数据集中一个很好的解释变量。探讨到精神健康中心的距离是解释病例数量还是诊断病例数量，是否存在未诊断或缺失病例的证据。研究倒数距离是否能提供更好的解释或更高的解释价值。使用回归方法单独分析家庭治疗率是否随人口密度或到中心的距离变化，并解释原因。
6. 分析航空公司 “拒载” 乘客的统计数据。思考为什么泊松分布适合建模航班上 “未出现” 乘客的数量，需要考虑哪些假设。不同航空公司的拒载率差异可达 50%，分析估计未来航班 “未出现” 乘客数量需要哪些信息。讨论商务旅客比例高的航班与度假目的地航班的拒载率是否不同，大型飞机飞往大城市的航班与座位较少的小型区域航班的拒载率是否不同，并解释原因。比较自愿被拒载和非自愿被拒载乘客数量的图表，评论使用两种不同尺度的情况，其中一个尺度是另一个的 20 倍。
7. 对于彩票获奖者数据，使用最终模型估计纽黑文市（人口 126,000，面积 21.1 平方英里，磨坊率 37.0）和达勒姆镇（人口 6,000，面积 23.3 平方英里，磨坊率 26.4）预期的彩票获奖者数量。判断这些是否是模型的合适应用，并解释原因。探讨是否有些城镇更幸运，获奖者数量是否与人口成正比，应该如何判断。思考对数尺度下截距的含义，以及其他可能有用的协变量值。判断泊松分布是否适合这个例子，并解释原因。
8. 对加拉帕戈斯群岛的物种数据进行分析。使用泊松分布对每个岛屿上的物种数量进行建模。由于伊莎贝拉岛比其他岛屿大得多，可能会产生较大的影响，可以考虑为伊莎贝拉岛创建一个指示变量。创建新的变量，如变量的乘积（交互项）和比率，并对其进行解释。检查数据中是否存在异常值，分析小岛屿之间是否相似，大岛屿之间的情况如何。判断泊松模型是否适合这个数据，思考岛屿上物种多样性较大是否存在进化优势，这是否能解释数据中的一些异常值，是否还有其他原因导致这些异常值，以及是否应该平等对待所有物种。

通过这些练习，可以更加深入地理解泊松回归在不同场景下的应用和分析方法。在实际应用泊松回归时，需要根据具体问题仔细选择合适的模型和参数，以获得准确可靠的结果。

深入理解泊松回归在SAS中的应用

5. 练习详细分析

5.1 概率计算与分布近似练习

在均值为 1 的泊松分布中，求事件发生次数为 2 或更少的概率，可使用泊松分布概率公式 (P(X = k)=\frac{e^{-\lambda}\lambda^{k}}{k!})，则 (P(X\leq2)=P(X = 0)+P(X = 1)+P(X = 2)=\frac{e^{-1}\times1^{0}}{0!}+\frac{e^{-1}\times1^{1}}{1!}+\frac{e^{-1}\times1^{2}}{2!})。

在 (N = 4) 且 (p = 1/4) 的二项实验中，使用二项分布概率公式 (P(X = k)=C_{N}^{k}p^{k}(1 - p)^{N - k})，则 (P(X\leq2)=P(X = 0)+P(X = 1)+P(X = 2)=C_{4}^{0}(\frac{1}{4})^{0}(1-\frac{1}{4})^{4}+C_{4}^{1}(\frac{1}{4})^{1}(1-\frac{1}{4})^{3}+C_{4}^{2}(\frac{1}{4})^{2}(1-\frac{1}{4})^{2})。

同理，在 (N = 8) 且 (p = 1/8) 的二项实验中，(P(X\leq2)=P(X = 0)+P(X = 1)+P(X = 2)=C_{8}^{0}(\frac{1}{8})^{0}(1-\frac{1}{8})^{8}+C_{8}^{1}(\frac{1}{8})^{1}(1-\frac{1}{8})^{7}+C_{8}^{2}(\frac{1}{8})^{2}(1-\frac{1}{8})^{6})。

二项分布的均值为 (Np)，方差为 (Np(1 - p))。随着 (N) 增大和 (p) 减小，二项模型的方差会更接近均值，2 或更少事件发生的概率也会更接近泊松近似，因为当 (N) 很大，(p) 很小时，二项分布近似于泊松分布。

5.2 航空公司超售问题

对于航空公司超售航班的问题，泊松分布适合建模 “未出现” 乘客的数量，因为乘客是否出现可看作独立事件，且在一定时间和空间内，“未出现” 乘客的平均数量相对稳定。泊松参数 (\lambda = 194\times0.05 = 9.7)。

设航空公司可以接受 (n) 个预订，要使有 99% 的信心不需要拒载任何乘客，即 (P(X\leq185 - n)\geq0.99)，可通过泊松分布表或统计软件求解 (n)。

5.3 人口测量单位对偏移量的影响

当以个人为单位测量人口时，运行之前的程序，输出会有所不同。因为偏移量的值发生了变化，导致模型中的线性预测器和估计参数也会改变。但从本质上来说，模型的结构和原理是不变的，只是参数的数值会根据人口测量单位的不同而调整。

5.4 死亡率数据的泊松回归

将住院死亡人数视为泊松分布的结果是合理的，因为死亡事件的发生相对独立，且在一定人群和时间内，死亡的平均发生率相对稳定。检查拟合回归模型的残差，可通过计算残差（如原始残差、卡方残差等），观察残差的分布是否与给定表格中最后一列显示的高低死亡率指示相对应，若残差分布合理，则说明模型能够较好地拟合数据。

5.5 精神疾病病例数据分析

病例数量通常会随着人口或人口对数的增加而增加，因为人口越多，潜在的病例数量可能就越多。人口密度是一个很好的解释变量，因为它反映了人口的集中程度，人口密度高的地区可能更容易传播疾病，也更容易获得医疗资源，从而影响病例数量。

到精神健康中心的距离可能会影响病例的诊断和报告，距离远的地区可能存在未诊断或缺失病例的情况。倒数距离可能会提供更好的解释，因为它可以突出近距离的影响。使用回归方法分析家庭治疗率与人口密度或到中心的距离的关系，可建立回归模型，观察回归系数的显著性和正负性，以确定它们之间的关系。

5.6 航空公司拒载统计分析

泊松分布适合建模航班上 “未出现” 乘客的数量，需要考虑的假设包括乘客是否出现相互独立，以及在一定时间段内 “未出现” 乘客的平均数量相对稳定。

估计未来航班 “未出现” 乘客数量需要的信息包括历史航班的 “未出现” 率、航班的类型（商务航班或度假航班）、乘客的来源地等。商务旅客比例高的航班拒载率可能较低，因为商务旅客行程相对固定；大型飞机飞往大城市的航班与小型区域航班的拒载率可能不同，因为航班的乘客构成、航线特点等因素不同。

比较自愿被拒载和非自愿被拒载乘客数量的图表，使用两种不同尺度会使数据的变化趋势看起来不同。一个尺度是另一个的 20 倍，可能会夸大或缩小某一数据的变化，在分析时需要注意这种尺度差异对结果的影响。

5.7 彩票获奖者数据建模

使用最终模型估计纽黑文市和达勒姆镇预期的彩票获奖者数量，将两市的人口、面积和磨坊率等数据代入模型即可。判断这些是否是模型的合适应用，需要考虑模型的适用范围和数据的代表性。如果模型是基于与这两个城市情况相似的数据建立的，那么估计结果可能是合理的；否则，可能需要进一步调整模型。

探讨是否有些城镇更幸运，可通过观察残差是否存在异常值来判断。获奖者数量是否与人口成正比，可通过分析模型中人口相关参数的显著性和大小来确定。对数尺度下截距的含义是当所有解释变量都为 0 时，响应变量的对数均值。其他可能有用的协变量值包括城镇的经济发展水平、教育程度等。判断泊松分布是否适合这个例子，需要检查数据是否满足泊松分布的假设，如事件发生的独立性和均值与方差相等。

5.8 加拉帕戈斯群岛物种数据分析

使用泊松分布对每个岛屿上的物种数量进行建模，由于伊莎贝拉岛比其他岛屿大得多，可能会对模型产生较大影响，为其创建指示变量可以控制这种影响。创建新的变量，如面积与距离的乘积、物种数量与面积的比率等，可以挖掘变量之间的潜在关系。

检查数据中是否存在异常值，可通过观察残差或使用统计方法（如箱线图）来判断。小岛屿之间可能由于生态环境相似而具有相似的物种数量，大岛屿之间可能由于面积、环境等因素差异较大而不同。

判断泊松模型是否适合这个数据，需要检查数据的均值和方差是否接近，以及事件发生的独立性。岛屿上物种多样性较大可能存在进化优势，因为更多的物种可以提供更多的生态位和资源利用方式，这可能解释数据中的一些异常值。其他导致异常值的原因可能包括人类活动、自然灾害等。是否应该平等对待所有物种，需要根据研究目的和数据特点来决定，如果不同物种对生态系统的影响不同，可能需要进行加权处理。

6. 总结与流程图

通过以上对泊松回归在 SAS 中的应用及相关练习的分析，我们可以总结出以下步骤：
1. 数据准备：收集和整理相关数据，确保数据的质量和完整性。
2. 模型建立：选择合适的解释变量，指定泊松分布，使用 SAS 的 proc genmod 过程建立模型。
3. 模型评估：通过拟合优度指标（如偏差、皮尔逊卡方等）评估模型的拟合效果。
4. 残差分析：计算不同类型的残差（如原始残差、卡方残差等），检查残差的分布，寻找异常值。
5. 模型调整：根据残差分析和实际情况，调整模型的参数或添加新的解释变量。
6. 结果应用：使用拟合好的模型进行预测和推断，解释模型的参数和结果。

以下是一个 mermaid 格式的流程图，展示了泊松回归分析的主要流程：

graph LR
    A[数据准备] --> B[模型建立]
    B --> C[模型评估]
    C --> D{模型是否合适?}
    D -- 是 --> E[残差分析]
    D -- 否 --> B
    E --> F{是否存在异常值?}
    F -- 是 --> G[模型调整]
    F -- 否 --> H[结果应用]
    G --> B

通过这个流程图，我们可以清晰地看到泊松回归分析的整个过程，从数据准备到结果应用，每个步骤都相互关联，需要不断地评估和调整，以确保模型的准确性和可靠性。在实际应用中，我们可以根据具体问题灵活运用这些步骤，更好地理解和分析数据。