对广义线性模型的学习.docx

对广义线性模型的学习.docx

《对广义线性模型的学习.docx》由会员分享，可在线阅读，更多相关《对广义线性模型的学习.docx（18页珍藏版）》请在冰点文库上搜索。

对广义线性模型的学习

对广义线性模型（GeneralizedLinearModel）的学习

引言在学习普通线性模型时就对因变量为离散的情况存有疑问。

在统计实验课程研读吴喜之老师的《复杂数据》一书的第六章时，发现了对离散因变量或者因变量为计数或有序数据时，可采用广义线性模型来处理。

因此这燃起了我对于广义线性模型的学习兴趣，通过查阅资料，对此模型有了以下的初步了解。

并在对经典方法理论有了一定的了解之后，利用该模型对实际数据进行了处理与分析，同时又用其他方法（包括机器学习等方法）对相同的数据进行了处理，在最后比较了各种方法之间的优缺点。

一、数据特点

1、横截面数据（Cross-SectionDat）a：

在同一时间，不同统计单位相同统计指标组成的数据列。

Note：

①与时序数据相比较，其区别在于数据的排列标准不同，时序数据是按照时间顺序排列的，横截面数据是按照统计单位排列的。

②横截面数据不要求统计对象及其范围相同，但要求统计的时间相同。

#横截面数据即为同一时间截面上的数据

2、横截面数据分析的要点：

①异方差问题由于数据是在某一时期对个体或地域的样本的采集，不同个

体或地域本身就存在差异。

②数据的一致性

主要包括变量的样本容量是否一致，样本的取样时期是否一致，数据的统计标准是否一致。

3、面板数据（PanelData）：

是指在时间序列上取多个截面，对于每一个截面上的数据均为一横截面数据列。

Note：

①面板数据是一个m*n的数据矩阵，记载的是n个时间节点上，m个对象的某一数据指标。

②其有时间序列和截面两个维度，当这类数据按两个维度排列时，是排在一个平面上，与只有一个维度的数据排在一条线上有着明显的不同，整个表格像是一个面板。

③如果从其内在含义上讲，把paneldata译为“时间序列-截面数据”更能揭示这类数据的本质上的特点。

4、广义线性模型主要用于因变量取离散值的情况当可能值为一切自然数0,1,2,,,时，多用Poisson分布；当Y取有限个值（实际是响应可以有有限个状态）时，多项分布是自然的选择。

5、在很大的程度上可以说，广义线性回归就是针对因变量为有限个

值情况的回归分析。

但在具体定模型时，需要考虑这有限个状态之间的关系

一种是无序的，即各状态的优劣并无公共的认定。

例如外出旅行，有k种交通工具可以选择，其优劣取决于具体情况而并无公认的排序。

另一种是有序的，即各状态的优劣次序有公共的认定。

如治疗效果、产品质量的分级等。

#不同情况建模方法有所不同。

二、广义线性模型的提出

广义线性模型的提出源于线性模型在应用上有重要影响的几个缺点：

1、只适用于因变量Y取值为连续的情况。

它特别不适用于分类数据（如Y取0.1为值）。

2、Y的期望E（Y）与自变量X是用线性关系E（Y）Z（X）相联系。

选择面太窄，往往与实际情况不符。

3、线性模型的统计推断基本上只适用于误差正态的情形。

在某些Y取值连续的场合，Y的分布是偏态的，如指数分布、伽马（Gamma）分布等。

广义线性模型的特点正好是对应上面指出的问题：

1、因变量Y可以取连续值或离散值，从常见的应用看，取离散值的场合更重要。

2、取代E（Y）ZT（X），有E（Y）h（ZT（X））函数h（其反函数g称为联系（或连接）函数（linkfunction））有较大的选择余地，这样扩大了模型的适用面。

3、Y（q维）有指数型分布Y~exp（Yb（））d（Y）其中，θ为q维参数向量，μ是R上的σ有限测度，μ与θ无关（或联系函数gh1使ZT（X），称自然联系）。

指数型分布是一个适中的选择，一方面它包括了应用上最常见的一些分布：

二项分布、多项分布、Poisson分布，以及连续型的正态分布、指数分布、伽马分布等。

另一方面，这分布类有很好的分析性质，又便于理论上的研究。

三、广义线性模型

设有因变量Y，自变量X，普通线性模型有以下几个特征：

1、E（Y）Z（X）（线性：

线性指对β，而非X）。

Z（X）为X的已知（向量）函数。

2、X，Z（X），Y都是取值连续的变量，如农作物产量、人的身高体重之类。

3、Y的分布为正态，或接近正态的分布。

广义线性模型从以下几个方面推广：

1、E（Y）h（ZT（X）），h为一严格单调、充分光滑的函数。

h已知，gh1（h的反函数）称为联系函数（linkfunction）,则有g（）ZT。

即E（Y）不等于ZT（X），而是ZT（X）的某一函数。

2、X，Z（X），Y可取连续或离散值，且在应用上更多见的情况为离散值。

如｛0,1｝，｛0,1,2,,,｝等。

3、Y的分布属于指数型，正态是其一特例。

4、以下的表格中列出了GLM中常用的几种分布：

由上表格中的第二列（Rangeofy）可以知道，当因变量为对应数据形式时应选择对应的分布来建立模型。

5、以下的表格中列出了GLM中常用的几种分布所对应的联系函数：

通常称这几种联系函数为标准联系函数，上表中的第三列为偏差。

四、R语言中的模型实现

在R语言中利用stats包中的glm（）函数来进行广义线性模型的拟合。

和lm函数类似，glm的建模结果可以通过下述的泛型函数进行二次

处理，如summary（）、coef（）、confint（）、residuals（）、anova（）、plot（）、predict（）。

R提供了一系列广义线性建模工具，从类型上来说包括gaussian，反gaussian，二项式，poisson和gamma模型的响应变量分布以及在响应变量分布没有明确给定时的拟似然（quasi-likelihood）模型。

在后者，方差函数（variancefunction）可以认为是均值的函数，但是在另外一些情况下，该函数可以由响应变量的分布得到。

函数glm（）的用法：

glm（formula,family=gaussian,data,weights,subset,na.action,start=NULL,etastart,mustart,offset,

control=list（...）,model=TRUE,method="glm.fit",x=FALSE,y=TRUE,contrasts=NULL,...）

多数选项与普通线性模型的拟合函数lm（）相同，值得注意的是family选项，family即为选择模型的分布，有以下几种选项：

#拟家族：

响应变量分布没有明确给定时的拟似然模型

#有过度离散现象时使用：

样本观测值变异性过大

quasipoisson（link="log"）#拟泊松分布

#有过度离散现象时使用：

样本观测值变异性过大注：

若样本观测值变异性过大，即出现了过度离散现象，此时仍使用二项分布假设就会影响系数检测的显著性。

那么补救的方法是使用准二项分布（quasibinomial）。

首先要检测样本是否存在过度离散现象，方法是用残差除以残差自由度，若超过1则意味着过度离散。

那么将family参数改为quasibinomial。

同样，在进行泊松分布也要考虑过度离散现象。

其检测方法同样是残差除以其自由度。

若确定过度离散存在，则要将family参数设置为准泊松分布（quasipoisson）。

在family的分布选项下还有几个常用选型即link和variance，可以用来选择联系函数和方差的形式。

Example:

glm（y~x,family=quasi（variance="mu^2",link="log"））

五、建立广义线性模型的实例

1、数据分析：

该数据是由美国国家癌症研究所资助的多中心血友病队列研究获得

的。

该项研究从1978年1月1日到1995年12月31日在16个治疗

中心（12个在美国，4个在西欧）跟踪了超过1600个血友病人，该

数据一共有2144个观测值及6个变量。

下表为变量情况：

为了更加直观的分析该数据的特点，截取了原数据中的部分数据行：

变量hiv为分类变量，只有两个选项，1和2；变量factor也为分类变量，有五个选项，1,2,3,4,5；变量year、age和deaths均为整数数据，只有变量py为数量变量。

要进行以死亡数即变量deaths作为因变量的回归，由于因变量为整数数据，因此选择广义线性模型来进行拟合。

考察因变量中数据的分布情况：

k}的形式，因此我们

发现可将其看作是{0,1，,,}或{0,1，,,，将采用Poisson对数线性模型（即分布设定为Poisson分布,联系函数设定为对数函数）和多项logit模型（即分布设定为二项分布，联系函数设定为logit函数）两种方法来进行数据的拟合。

2、卡方检验

卡方检验法是在总体X的分布未知时，根据来自总体的样本，检验关于总体分布的假设的一种检验方法。

由于这个数据的分布信息是未知的，并且我们也不是很容易直观的判断出它的分布信息，因此在这里我们采用卡方检验的方法来判断它的分布信息。

使用卡方检验对总体分布进行检验时，我们先提出原假设:

H0：

总体X的分布函数为F（x）然后根据样本的经验分布和所假设的理论分布之间的吻合程度来决定是否接受原假设。

这种检验通常称作拟合优度检验，它是一种非参数检验。

3、Poisson对数线性模型

模型：

其中，i（i=1,2）代表hiv的两个水平，j（j=1,2,,,,5）代表factor的5个水平，x1代表year（1代表year的系数），x2代表age（2代表age的系数），x3代表py（3代表py的系数），0代表截距。

>ap=glm（deaths~.,family='poisson',datw）

a=

>summary（ap）

Call:

glm（formula=deaths~.,family="poisson",data=w）

DevianceResiduals:

Min1QMedian3QMax

-2.1139-0.4316-0.2209-0.10263.2727

Coefficients:

EstimateStd.ErrorzvaluePr（>|z|）

（Intercept）

-23.135255

1.318652-

17.545

<2e-16***

hiv2

2.766461

0.203259

13.611

<2e-16***

factor2

-0.636420

0.151922

-4.1892.80e-05***

factor3

-0.403434

0.140538

-2.871

0.0041**

factor4

-0.707524

0.142711

-4.9587.13e-07***

factor5

-0.371585

0.146238

-2.541

0.0111*

year

0.211047

0.014090

14.979

<2e-16***

age

0.077867

0.015495

5.0255.03e-07***

py

0.033042

0.002845

11.614

<2e-16***

Signif.codes:

0‘***

'0.001

0.01‘*'0.05‘.'0.1‘'1

（Dispersionparameterforpoissonfamilytakentobe1）

Nulldeviance:

1892.8on2143degreesoffreedomResidualdeviance:

1007.6on2135degreesoffreedomAIC:

1725.7

NumberofFisherScoringiterations:

6

得到的模型拟合结果为：

ln（）23.1402.7700.640.400.710.370.21x10.08x20.03x3在模型中，定性自变量的各个水平的单独效应是不可估计的，必须加上约束条件，这里的约束条件是每个定性变量第一个水平为0。

即效

应1（hiv1）及1（factor1）按照R的默认约束条件都等于0。

结果分析：

首先，各个变量都很显著，相比较而言factor3和factor5的显著性

较差一些其次，当设定hiv1的效应为0时，hiv2对于死亡数的效应为正，且效应比hiv1的效应大；当设定factor1的效应为0时，factor的其余四个选项对于死亡数的效应均为负，且factor4的效应最大，factor5的效应最小；变量year对于死亡数的影响较大，其余两个变量对其影响较小。

由模型拟合结果来分析实际情况，可知hiv为阳性时对血友病有较坏的影响，且影响较大；而在使用凝血因子制剂之后，对于病情均有改善，第二种和第四种制剂对于病情的改善效果较为明显；而变量year对于死亡数的影响明显比age和py的影响大，分析情况可能是因为医疗条件的进步，对于血友病的治疗有了明显的帮助。

十折交叉验证对测试集的拟合结果：

4、拟似然（quasi-likelihood）模型对于所有的族，响应变量的方差依赖于均值并且拥有作为系数

（multiplier）的尺度参数。

方差对均值的依赖方式是响应分布的一个特性；例如对于poisson分布Var（y）=mu。

对于拟似然估计和推断，我们不是设定精确的响应分布而是设定关联函数和方差函数的形式。

因为关联函数和方差函数都依赖于均

#即拟似然模型为响应变量分布没有明确给定

>ap=glm（deaths~.,family='quasi',data=w）

>summary（ap）

Call:

glm（formula=deaths~.,family="quasi",data=w）

DevianceResiduals:

Min1QMedian3QMax

-0.8530-0.2895-0.08740.14475.2069

Coefficients:

EstimateStd.ErrortvaluePr（>|t|）

（Intercept）

-3.37999530.2172690-15.557<2e-16***

hiv

0.3783377

0.0236582

15.992<2e-16***

factor

-0.0435631

0.0086322

-5.0474.88e-07***

year

0.0340688

0.0024654

13.819<2e-16***

age

0.0169687

0.0032860

5.1642.64e-07***

py

0.0066769

0.0006726

9.927<2e-16***

（Dispersionparameterforquasifamilytakentobe0.2946817）

Nulldeviance:

795.31on2143degreesoffreedomResidualdeviance:

630.03on2138degreesoffreedomAIC:

NA

NumberofFisherScoringiterations:

2

由于没有明确的分布，这里并不区分分类变量的各个选项，只给出此变量的效应值，得到的模型拟合结果与Poisson对数线性模型基本一致，具体分析在这里不再赘述。

十折交叉验证对测试集的拟合结果：

5、多项logit模型

模型：

多项logit模型在类别上仍可归为广义线性模型，是二项分布的logistic回归向多项分布的推广，但是在R语言的glm（）函数中只能进行二项分布的回归，而无法进行多项分布的回归。

因此我们利用R语言mlogit包中的mlogit（）函数来进行模型的拟合。

函数mlogit的用法：

mlogit（formula,data,subset,weights,na.action,start=NULL,

alt.subset=NULL,reflevel=NULL,

nests=NULL,un.nest.el=FALSE,unscaled=FALSE,

heterosc=FALSE,rpar=NULL,probit=FALSE,

R=40,correlation=FALSE,halton=NULL,

random.nb=NULL,panel=FALSE,estimate=TRUE,

seed=10,...）

mlogit.data（data,choice,shape=c（"wide","long"）,varying=NULL,

sep=".",alt.var=NULL,chid.var=NULL,alt.levels=NULL,

id.var=NULL,opposite=NULL,drop.index=FALSE,

ranked=FALSE,...）

参数说明：

formula:

mlogit提供了条件logit，多项logit，混合logit多种模型，对于多项logit的估计模型应写为：

因变量~0|自变量，如：

mode~0|income。

data：

先使用mlogit.data函数使得数据结构符合mlogit函数要求。

choice：

确定分类变量是什么。

shape：

如果每一行是一个观测，我们选择wide，如果每一行是表示一个选择，那么就应该选择long。

alt.var：

对于shape为long的数据，需要标明所有选择名称。

>a=mlogit（deaths~0|hiv+factor+year+age+py,data=w1）

>summary（a）

Call:

mlogit（formula=deaths~0|hiv+factor+year+age+py,data=w1,method="nr",print.level=0）

Frequenciesofalternatives:

0123456

0.854944030.098880600.028917910.013059700.002798510.000932840.00046642nrmethod

21iterations,0h:

0m:

9sg'（-H）^-1g=7.61E-07gradientclosetozero

Coefficients:

EstimateStd.Errort-valuePr（>|t|）

1:

（intercept）-2.6616e+012.0816e+00-12.7863<2.2e-16***

2:

（intercept）-8.2781e+012.1635e+04-0.00380.9969471

3:

（intercept）-8.9760e+011.8693e+04-0.00480.9961687

4:

（intercept）-1.1801e+021.9020e+04-0.00620.9950497

5:

（intercept）-1.3355e+021.2532e+04-0.01070.9914969

6:

（intercept）-1.3452e+021.4785e+04-0.00910.99274071:

hiv2.5000e+002.2145e-0111.2893<2.2e-16***2:

hiv2.3683e+011.0818e+040.00220.9982532

3:

hiv

2.3319e+01

9.3463e+03

0.00250.9980093

4:

hiv

2.1054e+01

9.5101e+03

0.00220.9982336

5:

hiv

1.7966e+01

6.1673e+03

0.00290.9976757

6:

hiv

2.2628e+01

7.3922e+03

0.00310.9975576

1:

factor

-4.6267e-03

5.8487e-02

-0.07910.9369485

2:

factor

-1.2710e-01

1.0385e-01

-1.22390.2209819

3:

factor

-4.5252e-01

1.6078e-01

-2.81450.0048859**

4:

factor

-7.3439e-01

3.7890e-01

-1.93820.0525968.

5:

factor

-1.5480e+01

2.2125e+03-0.00700.9944175

6:

factor

6.6159e-01

1.8911e+00

0.34980