泊松回归[一般对数线性分析过程]

转自个人微信公众号【Memo_Cleon】的统计学习笔记：泊松回归：一般对数线性分析过程。

【1】计数资料[count data]

“哮喘患者一个月内急性发作次数”是什么资料？

表中[是否发作]和[发作次数]的数据类型一样吗？显然不一样。[是否发作]属于分类资料（categorical data），其值是“是”或者“否”，没有单位，同表中的[组别]、[性别]性质是一样的，而[发作次数]是靠清点计数得到的，虽然它有自己的单位（本例是次），但与连续资料又有不同，它只能是非负的整数，在统计上这类资料属于计数资料[count data]。我们也可以将[发作次数]转换成[等级]，这种资料就变成了分类资料中的等级资料（有序多分类资料）。

分类资料最常用的单因素分析方法是卡方分析，多因素分析可以考虑logistic回归、对数线性回归，计数资料最常用的就是Poisson回归，即泊松回归，全称应该是泊松对数线性回归模型。泊松回归常用来研究单位时间/单位面积/单位空间内某事件的发生数的影响因素。

【2】Poisson分布

像线性回归要求满足正态分布一样，泊松回归要求满足泊松分布，泊松分布的概率密度函数为：

表示单位时间/单位面积/单位空间内某事件发生k次的概率。λ是这个分布唯一的参数，λ是该分布的均值，也是该分布的方差，即均值等于方差，λ越大Poisson分布越逼近正态分布。

【3】Poisson回归

如果某一事件的的观测频数服从一个均值为λ的Poisson分布，则Poisson回归模型如下：

由上述公式很容易得出偏回归系数β的实际意义：β₀表示各个“自变量取值为0”时观测频数的自然对数值，注意这里的取值为0的水平指的是被设为参照的那个水平，在SPSS中默认高水平赋值的那个水平为参照水平，跟自变量实际的value赋值不一致。如果模型是非饱和模型，结果会与实际计数计算的有差别，但如果是饱和模型则完全一致。β_i则表示其他自变量取值不变时，自变量x_i每改变一个单位，所引起的观测频数λ的自然对数值的改变量，β_i为正表示自变量每增加一个单位观测频数λ会增加，βi为负值表示自变量每增加一个单位观测频数λ会减少。

当不同观察单位的发生事件观测基数不一样时，需要将发生数除以相应单位的总观测数，用于抵消观察单位数的不同对结果造成的的影响。

其中lnn称为偏移量（offset）。

此时β₀表示参照水平的事件发生率的自然对数，exp(β₀)表示参照水平的事件发生率。β_i则表示其他自变量取值不变时，自变量xi每改变一个单位，所引起的事件发生率自然对数值的改变量。

我们把上面的公式换种像我这样的数学小白更容易理解的表达方式：

但其他变量保持不变，自变量xi改变1个单位，取值由xi变为xi+1，发生率比值（incidence-rate ratios，IRR，意同RR）:

这样系数βi的解释可以转换成更具有实际意义的IRR，即exp(βi)，意为自变量x_i每改变一个单位，事件发生率是参照水平的exp(βi)倍，或者说相比参照水平，事件发生率会改变[exp(βi)-1]倍。

个人以为SPSS中的Poisson回归操作界面并不友好，对话框中去掉了因变量的这个选框，大概理由是结局变量就是单元格的计数，所以不需要选择，但却需要专门的数据格式来匹配这个模型，如果你习惯了在对话框中确定结局变量和解释变量，对此可能会有些不适应。

【4】示例1：英国男性医生冠心病死亡与吸烟、年龄的关系。引自自stata16帮助文件。

数据录入：smoke（0=不吸烟，1=吸烟），agecat（0=35-44岁，1=45-54岁，2=5-64岁，3=65-74岁，5=75-84岁），deaths表示存活人数变量，pyears表示人年数，相当于总观察数。要注意这种输入格式跟我们以前常用的交叉表格式是不一样的。

我们熟悉的交叉表数据录入格式是下面这种形式，但下面这种数据格式适用于广义线性模型，并不适合一般对数线性分析中的Poisson模型。

数据加权：Data>>Weight Cases…，激活加权变量，选入变量死亡数[deaths]

泊松回归：Analyze>>Loglinear>>General…

[因子]：选入变量smoke、agecat；

[单元格结构]：选入观察基数变量pyears。单元格结构选入单元格的观察单位数，用于计算偏移量，去除基数的不同对模型造成的影响；

【模型】自定义模型，选入smoke和agecat的主效应；继续；

【选项】选中参数估计；δ取值0-1，默认为0.5，模型会将所有单元格的频数或者计数加上该值进行参数估计，以避免单元格频数或计数为0时引起的计算问题。数据量较少时会略微影响参数的估计值，但不会影响统计检验结果。如果数据不存在空单元格可以将其改为0。继续；

OK

结果：

（1）数据信息及模型收敛情况

（2）拟合优度检验：两种方法的P<0.05，表明该模型与饱和模型存在差异，也就是该模型还需要进一步纳入别的选项。由于两个等级变量（smoke、age）的饱和模型包含3项，即smoke、agecat主效应以及两者的交互项smoke*agecat，只含有smoke和agecat主效应的简约模型与饱和模型相比只缺少交互项，拟合优度检验具有统计学意义表明交互项是存在的，应该纳入模型。

（3）单元格计数与残差：显示实际计数、理论计数、残差、标准化残差、校正残差即偏差等信息。

（4）参数估计：这是回归模型中最为关键的结果。

校正年龄因素后，与抽烟的男医生相比，不抽烟的男医生的冠心病死亡风险更低，不抽烟是抽烟死亡风险的0.70倍【IRR=exp(-0.355)=0.70】，或者说抽烟是不抽烟死亡风险的1.43倍【IRR=exp(0.355)=1.43，或IRR=1/exp(0.355)】。

校正吸烟因素后，随着年龄的增加，死亡风险也在增加。75-84岁的男医生的冠心病死亡风险分别是35-44岁、45-54岁、55-64岁、65-74岁的男医生的40.45倍、9.17倍、2.92倍和1.42倍。

其他任意两水平的IRR可用exp(两者系数相减值）计算，统计学检验可通过95%CI是否有重叠进行大体评估。比如45-54岁相比35-44岁的冠心病死亡IRR=exp(-2.216+3.700)=4.41，两者95%CI没有重叠部分，表明具有统计学意义。

（5）相关参数估计与协方差参数估计：显示相关系数与协方差矩阵。

（6）默认的几个散点图：实际计数、理论计数与校正残差的散点图矩阵；校正残差的正态Q-Q图；校正残差的去趋势Q-Q图。从散点图矩阵和正态Q-Q图来看，扣除期望值后的残差基本满足正态分布，模型拟合的还算可以，但去趋势Q-Q图显示，还是存在一定的趋势，这就意味着残差还有未能被模型中的变量解释的部分。

注：nondetrended normal plots:The observed values are plotted against the expected values.detrended plots: The observed values are plotted against the deviations from the expected values.QQ图上点偏离斜线的距离肉眼观察可能会误导我们，感觉像是在斜线周边。去势QQ图可以更清楚得看出实际偏离有多大。

拟合优度检验结果表明，当前简约模型与饱和模型存在差异，我们可以在【模型】中构建新的模型进行分析。当前简约模型与饱和模型只差一个交互项，我们也可以在原来的简约模型中添加缺少的交互项。我们也选择的饱和模型，饱和模型是默认的分析模型。为更好的理解各系数的意义，也可以在【选项】中选中[设计矩阵(Design matrix)]。

饱和模型中各单元格的理论计数与实际计数完全一致，因此各种残差均为0，另外，All residuals are zero under the model. Therefore the requested charts will not be created，各种残差图也不再创建。

我们看一下饱和模型的参数估计结果：交互作用的确有统计学意义。本例吸烟与年龄的交互作用具有统计学意义，偏回归系数值表明吸烟对冠心病死亡的影响随着年龄的增大而增大。交互作用具有统计学意义，无论吸烟与年龄的主效应有无统计学意义都应该纳入模型。

由于交互作用的存在，单独分析smoke与age的主效应已无多大实际意义，我们需要进一步分析单独效应（按某个因素不同水平下另外一个因素的效应）。

我们在前面讲过，常数项代表的是参照水平的事件发生率的对数，或者说exp(常数项）=参照水平的事件发生率，饱和模型时与实际计数计算结果完全一致。本例exp（-3.949）=1.93%，参照水平【吸烟（smoke=1）、年龄在75-84岁（agecat=5）】的死亡人数102例，观察数5317，死亡发生率=102/5317=1.92%。你可能会问，为什么不完全一样呢？！你可能已经知道了，在Poisson对话框的【选项】中的δ我们使用了默认的在原单元格的计数值上增加0.5以避免0值单元格带来的计算问题，实际上本例各单元格并没有0值，如果将δ修改为0结果见下面第二张表，常数项估计值为-3.954，exp(-3.954)=1.92%。

单独效应可以大体计算如下：

• 不同年龄的男医生

对于35-44岁的男医生，不吸烟发生冠心病死亡的发生率是exp(-3.954+0.100--3.447-1.847)，吸烟发生冠心病死亡的发生率是exp(-3.954+0-3.447+0)，吸烟IRR=exp(-3.954+0-3.447+0)/exp(-3.954+0.100-3.447-1.847)，即吸烟发生死亡的风险是不吸烟的[1/exp(0.100-1.847)]=5.74倍；

同样的，对于45-54岁的男医生，吸烟发生死亡的风险是不吸烟的[1/exp(0.100-0.860)]=2.14倍；

对于55-64岁的男医生，吸烟发生死亡的风险是不吸烟的[1/exp(0.100-0.484)]=1.47倍；

对于65-74岁的男医生，吸烟发生死亡的风险是不吸烟的[1/exp(0.100-0.405)]=1.36倍；

对于75-84岁的男医生，吸烟发生死亡的风险是不吸烟的[1/exp(0.100-0)]=0.90倍。

• 吸烟和不吸烟的男医生

对于不吸烟的男医生，75-84岁的男医生发生冠心病死亡的风险是35-44岁的男医生[1/exp(-3.447-1.847)]=199.14倍，同样可以获得不吸烟男医生各年龄段的IRR，感兴趣的可以自己去计算一下。

对于吸烟的男医生，75-84岁的男医生发生冠心病死亡的风险是35-44岁的男医生[1/exp(-3.447+0)]=31.41倍，同样可以获得吸烟男医生各年龄段的IRR，也不再累述。

通过参数估计表可以对简单效应做上述的大体判断，但我们无法判断各个效应是否有统计学意义，而且计算起来啰里啰嗦的，如果想明确是否有统计学意义，可以通过分割文件分组比较来完成。

Data>>Spilit Files…，操作如下：

结果如下表。具体解读就不再赘述了，通过吸烟与否进行文件拆分比较各年龄段的IRR感兴趣的也可以自己去操作一下看看。

转自个人微信公众号【Memo_Cleon】的统计学习笔记：泊松回归：一般对数线性分析过程。

… E N D …