泊松回归分析(Poisson Regression Analysis)——Python软件实现

关键词：Python; 泊松回归; Poisson回归; 等离散

一、案例介绍

某临床医师对39名有胸闷症状的非器质性心脏病男性患者的24小时早搏数(beat)进行了临床研究记录，每个患者的研究影响因素包括是否喝浓茶(tea，0为不喝浓茶，1为喝浓茶)、是否吸烟(smoke，0为不吸烟，1为吸烟)。请利用该资料对24小时早搏数的影响因素进行分析，部分数据见图1。本案例数据可从“附件下载”处下载。

二、问题分析

本案例的分析目的是了解有胸闷症状的非器质性心脏病男性患者24小时早搏数的影响因素。了解单位时间、单位面积或单位空间内某事件发生数的影响因素，可以考虑使用Poisson回归分析。但需要满足以下6个条件：

条件1：观察变量为计数变量。本研究中心脏病男性患者的24小时早搏数为计数变量，该条件满足。

条件2：观察变量的发生相互独立。本研究中各研究对象的每次早搏事件发生都是独立的，不存在互相干扰的情况，该条件满足。

条件3：至少有1个自变量，可以是分类变量，也可以是连续变量。本研究中有两个分类自变量，分别为是否喝浓茶和是否吸烟，该条件满足。

条件4：观察变量不存在显著的异常值，该条件需要通过软件分析后判断。

条件5：观察变量服从Poisson分布，即满足等离散性，表现为计数值的平均值(近似)等于方差。该条件可以通过数据特征进行初步判断，本研究的观察变量为心脏病男性患者24小时早搏数(计数资料)，从专业知识可知，早搏发生频数较低，各单位时间内的发生情况相互独立，基本满足Poisson分布的条件。同时还也可以结合软件分析进行判断。

条件6：自变量之间无多重共线性。该条件需要通过软件分析后判断。

三、软件操作及结果解读

(一) 导入数据

import pandas as pd #导入pandas包
df = pd.read_csv(r"泊松回归.csv")
pd.set_option('display.max_rows', 10)
df

在数据栏目中可以查看全部数据情况，数据集中共有4个变量和39个观察数据，4个变量分别为个体编号(ID)、是否喝浓茶(tea)、是否吸烟(smoke)及24小时早搏数(beat)。

如果数据集较大也可使用如下命令查看数据框结构(图2)：

df.info()   #查看数据结构

(二) 适用条件判断

1.条件4判断(异常值检测)

(1) 软件操作

import matplotlib.pyplot as plt#绘制箱线图
plt.boxplot(df.beat)
plt.show()

(2) 结果解读

异常值通过箱线图和专业知识进行判断，箱线图(图3)提示存在一个异常值。查看数据表可以发现在早搏数的第14位数值为17，依据专业可判定该值可以保留。综上，本案例未发现需要处理的异常值，满足条件4。

2.条件5判断(等离散性检验)

(1) 软件操作

print(df.beat.describe())
import numpy as np
print(np.var(df.beat))

(2) 结果解读

描述性分析结果(图4)显示早搏数的均数为7.308，方差为10.418，均数近似等于方差，提示数据满足等离散性。数据满足条件5。

3.条件6判断(多重共线性诊断)

(1) 软件操作

变量之间的共线性诊断可以采用容忍度法(Tolerance)、方差膨胀因子法(VIF)、相关系数法和条件数法。本案例采用方差膨胀因子法来查看变量之间的多重共线性，可在线性回归后计算方差膨胀因子。

from statsmodels.stats.outliers_influence import variance_inflation_factor
df.loc[:,'intercept'] = 1
vif = [variance_inflation_factor(df.drop(['ID','beat'],axis = 1), i) for i in range(df.drop(['ID','beat'],axis = 1).shape[1])]
pd.DataFrame(data = zip(df.drop(['ID','beat'],axis = 1),vif),columns = ['变量','vif'])

(2) 结果解读

方差膨胀因子(VIF)结果(图5)显示，2个自变量的VIF均远小于10，提示变量间不存在严重的多重共线性，可知满足条件6。

(三) 模型拟合

1. 软件操作

from statsmodels.formula.api import poisson
model_pos = poisson(formula='''beat~tea+smoke''',data=df)
results_pos = model_pos.fit()
results_pos.summary()#建立泊松回归模型

#计算变量tea效应量的95%CI#

import math
IRR1=math.exp(0.4231)
IRR1CIlow=math.exp(0.187)
IRR1CIhigh=math.exp(0.659)
print(IRR1)
print(IRR1CIlow)
print(IRR1CIhigh)

#计算变量smoke效应量的95%CI#

IRR2=math.exp(0.6084)
IRR2CIlow=math.exp(0.335)
IRR2CIhigh=math.exp(0.882)
print(IRR2)
print(IRR2CIlow)
print(IRR2CIhigh)

2. 结果解读

Poisson回归分析的结果(图6)中列出了截距和自变量的“coef (非标准化系数)”、“std err (标准误)”、模型变量系数的95%CI、统计量z值及P值。并计算了模型变量系数的IRR值及其95%CI (图7~图8)。

回归分析结果显示(图6)，变量“是否吸烟”的P值＜0.001，说明吸烟患者和不吸烟患者24小时早搏数的差异存在统计学意义；变量“是否喝浓茶”的P值＜0.001，则表明喝浓茶和不喝浓茶的患者，24小时早搏数的差异存在统计学意义；从系数的值来看，吸烟和喝浓茶对患者24小时早搏数均有正向影响。

变量tea效应量计算结果(图7)显示，喝浓茶患者早搏的概率是不喝浓茶患者的1.527倍(95%CI：1.206~1.933)，变量smoke效应量计算结果(图8)显示，吸烟患者早搏的概率是不吸烟患者的1.837倍(95%CI：1.397~2.415)。

四、结论

本研究采用Poisson回归分析患者胸闷症状的非器质性心脏病男性患者24小时早搏数的影响因素。观察变量早搏数为计数变量，且单位时间内其发生相互独立，经过软件判断服从Poisson分布，不存在需要处理的异常值，且自变量之间不存在严重共线性，数据满足Poisson回归分析的条件。

结果显示，吸烟患者和不吸烟患者24小时早搏数存在差异(P＜0.001)，吸烟患者早搏的概率是不吸烟患者的1.837倍(95%CI：1.397~2.415)，喝浓茶患者早搏的概率是不喝浓茶患者的1.527倍(95%CI：1.206~1.933)。

五、知识小贴士

(一) Poisson分布

• 人类稀有疾病或一些卫生事件，如恶性肿瘤、某地在一个月内因交通事故死亡人数、1 mL水中大肠杆菌数等计数资料，具有发病率低或者不像二项分布资料有分母能计算比例等特点。因此，这些事件数的多少除了取决于事件的实际发生数，还取决于计数时研究者所观察的范围，即观察多长时间、多大人群、多大面积等。使用发病密度等密度指标描述这些事件的群体特征比较合适。对于此类罕见事件的发生，如果事件之间彼此相互独立，观察样本含量较大时，则具有平均计数等于方差的特点。这类事件的发生次数往往服从Poisson分布。

(二) Poisson回归

• Poisson回归主要用于单位时间、单位面积、单位空间内某事件发生数的影响因素分析。在进行稀有事件等计数资料的影响因素分析时，应首先对资料的过离散情况进行判断或检验分析，然后选择正确的回归模型分析，才能得到正确的结果。一般先从专业方面判断，然后用统计学方法检验资料是否存在过离散现象。如果资料存在过离散情况，选用负二项回归模型；如果资料无过离散情况，选用Poisson回归模型分析。Poisson回归和多重线性回归、Logistics回归、Cox比例风险模型、负二项回归等都是医学领域中应用最多的广义线性模型之一。