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

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

【1】计数资料[count data]

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

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

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

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

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

【3】Poisson回归

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

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

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

其中lnn称为偏移量(offset)。

此时β0表示参照水平的事件发生率的自然对数,exp(β0)表示参照水平的事件发生率βi则表示其他自变量取值不变时,自变量xi每改变一个单位,所引起的事件发生率自然对数值的改变量。
我们把上面的公式换种像我这样的数学小白更容易理解的表达方式:

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

这样系数βi的解释可以转换成更具有实际意义的IRR,即exp(βi)意为自变量xi每改变一个单位,事件发生率是参照水平的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 …

(0)

相关推荐