对广义线性模型的学习

1、对广义线性模型（Generalized Linear Model）的学习引言 在学习普通线性模型时就对因变量为离散的情况存有疑问。在统计实验课程研读吴喜之老师的复杂数据一书的第六章时，发现了对离散因变量或者因变量为计数或有序数据时，可采用广义线性模型来处理。因此这燃起了我对于广义线性模型的学习兴趣，通过查阅资料，对此模型有了以下的初步了解。并在对经典方法理论有了一定的了解之后，利用该模型对实际数据进行了处理与分析，同时又用其他方法（包括机器学习等方法）对相同的数据进行了处理，在最后比较了各种方法之间的优缺点。一、数据特点1、横截面数据（Cross-Section Data）：在同一时间，不同统

2、计单位相同统计指标组成的数据列。 Note：与时序数据相比较，其区别在于数据的排列标准不同， 时序数据是按照时间顺序排列的，横截面数据是按照 统计单位排列的。 横截面数据不要求统计对象及其范围相同，但要求统 计的时间相同。#横截面数据即为同一时间截面上的数据2、 横截面数据分析的要点： 异方差问题 由于数据是在某一时期对个体或地域的样本的采集，不同个 体或地域本身就存在差异。 数据的一致性 主要包括变量的样本容量是否一致，样本的取样时期是否一 致，数据的统计标准是否一致。3、 面板数据（Panel Data）：是指在时间序列上取多个截面，对于每一个截面上的数据均为一横截面数据列。 Note：面

3、板数据是一个m*n的数据矩阵，记载的是n个时 间节点上，m个对象的某一数据指标。 其有时间序列和截面两个维度，当这类数据按两个 维度排列时，是排在一个平面上，与只有一个维度 的数据排在一条线上有着明显的不同，整个表格像 是一个面板。 如果从其内在含义上讲，把panel data译为“时间 序列-截面数据”更能揭示这类数据的本质上的特点。4、 广义线性模型主要用于因变量取离散值的情况 当可能值为一切自然数0,1,2,时，多用Poisson分布； 当Y取有限个值（实际是响应可以有有限个状态）时，多项分布 是自然的选择。5、 在很大的程度上可以说，广义线性回归就是针对因变量为有限个值情况的回归分析。

4、但在具体定模型时，需要考虑这有限个状态之间的关系。一种是无序的，即各状态的优劣并无公共的认定。例如外出旅行，有k种交通工具可以选择，其优劣取决于具体情况而并无公认的排序。另一种是有序的，即各状态的优劣次序有公共的认定。如治疗效果、产品质量的分级等。#不同情况建模方法有所不同。2、 广义线性模型的提出广义线性模型的提出源于线性模型在应用上有重要影响的几个缺点：1、 只适用于因变量Y取值为连续的情况。 它特别不适用于分类数据（如Y取0.1为值）。2、 Y的期望E(Y)与自变量X是用线性关系相联系。 选择面太窄，往往与实际情况不符。3、 线性模型的统计推断基本上只适用于误差正态的情形。 在某些Y取值

5、连续的场合，Y的分布是偏态的，如指数分布、 伽马（Gamma）分布等。 广义线性模型的特点正好是对应上面指出的问题：1、 因变量Y可以取连续值或离散值，从常见的应用看，取离散 值的场合更重要。2、 取代，有 函数（其反函数称为联系（或连接）函数（link function） 有较大的选择余地，这样扩大了模型的适用面。3、 Y（q维）有指数型分布 其中，为q维参数向量，是上的有限测度，与 无关（或联系函数使，称自然联系）。 指数型分布是一个适中的选择，一方面它包括了应用上最常 见的一些分布：二项分布、多项分布、Poisson分布，以及 连续型的正态分布、指数分布、伽马分布等。另一方面，这 分布类

6、有很好的分析性质，又便于理论上的研究。3、 广义线性模型设有因变量Y，自变量X，普通线性模型有以下几个特征：1、 （线性：线性指对，而非X）。 Z(X)为X的已知（向量）函数。2、 X，Z(X)，Y都是取值连续的变量，如农作物产量、人的身高 体重之类。3、 Y的分布为正态，或接近正态的分布。广义线性模型从以下几个方面推广： 1、，为一严格单调、充分光滑的函数。 已知，（的反函数）称为联系函数（link function）, 则有。 即不等于，而是的某一函数。 2、X，Z(X)，Y可取连续或离散值，且在应用上更多见的情况为 离散值。如0,1，0,1,2,等。 3、Y的分布属于指数型，正态是其一特

7、例。 4、以下的表格中列出了GLM中常用的几种分布：由上表格中的第二列（Range of y）可以知道，当因变量为对应数据形式时应选择对应的分布来建立模型。5、以下的表格中列出了GLM中常用的几种分布所对应的联系函数：通常称这几种联系函数为标准联系函数，上表中的第三列为偏差。4、 R语言中的模型实现 在R语言中利用stats包中的glm()函数来进行广义线性模型的拟合。和lm函数类似，glm的建模结果可以通过下述的泛型函数进行二次处理，如summary()、coef()、confint()、residuals()、anova()、plot()、predict()。 R提供了一系列广义线性建模工

8、具，从类型上来说包括gaussian，反gaussian，二项式，poisson和gamma模型的响应变量分布以及在响应变量分布没有明确给定时的拟似然(quasi-likelihood)模型。在后者，方差函数(variance function)可以认为是均值的函数，但是在另外一些情况下，该函数可以由响应变量的分布得到。函数glm()的用法：glm(formula, family = gaussian, data, weights, subset, na.action, start = NULL, etastart, mustart, offset, control = list(.), mo

9、del = TRUE, method = "glm.fit", x = FALSE, y = TRUE, contrasts = NULL, .)多数选项与普通线性模型的拟合函数lm()相同，值得注意的是family选项，family即为选择模型的分布，有以下几种选项：binomial(link = "logit")#二项分布gaussian(link = "identity")#正态分布Gamma(link = "inverse")#伽马分布inverse.gaussian(link = "1/mu2&

10、quot;)#反Gaussian分布poisson(link = "log")#泊松分布quasi(link = "identity", variance = "constant")#（quasi-likelihood）#拟家族：响应变量分布没有明确给定时的拟似然模型quasibinomial(link = "logit")#拟二项分布#有过度离散现象时使用：样本观测值变异性过大quasipoisson(link = "log")#拟泊松分布#有过度离散现象时使用：样本观测值变异性过大注：若样本

11、观测值变异性过大，即出现了过度离散现象，此时仍使用二项分布假设就会影响系数检测的显著性。那么补救的方法是使用准二项分布（quasibinomial）。首先要检测样本是否存在过度离散现象，方法是用残差除以残差自由度，若超过1则意味着过度离散。那么将family参数改为quasibinomial。同样，在进行泊松分布也要考虑过度离散现象。其检测方法同样是残差除以其自由度。若确定过度离散存在，则要将family参数设置为准泊松分布(quasipoisson)。在family的分布选项下还有几个常用选型即link和variance，可以用来选择联系函数和方差的形式。Example:glm(y x, f

12、amily = quasi(variance = "mu2", link = "log")5、 建立广义线性模型的实例 1、数据分析：该数据是由美国国家癌症研究所资助的多中心血友病队列研究获得的。该项研究从1978年1月1日到1995年12月31日在16个治疗中心（12个在美国，4个在西欧）跟踪了超过1600个血友病人，该数据一共有2144个观测值及6个变量。下表为变量情况：为了更加直观的分析该数据的特点，截取了原数据中的部分数据行： 变量hiv为分类变量，只有两个选项，1和2；变量factor也为分类变量，有五个选项，1,2,3,4,5；变量year、

13、age和deaths均为整数数据，只有变量py为数量变量。要进行以死亡数即变量deaths作为因变量的回归，由于因变量为整数数据，因此选择广义线性模型来进行拟合。考察因变量中数据的分布情况：发现可将其看作是0,1，或0,1，k的形式，因此我们将采用Poisson对数线性模型（即分布设定为Poisson分布,联系函数设定为对数函数）和多项logit模型（即分布设定为二项分布，联系函数设定为logit函数）两种方法来进行数据的拟合。 2、卡方检验 卡方检验法是在总体X 的分布未知时，根据来自总体的样本，检验关于总体分布的假设的一种检验方法。 由于这个数据的分布信息是未知的，并且我们也不是很容易直观

14、的判断出它的分布信息，因此在这里我们采用卡方检验的方法来判断它的分布信息。 使用卡方检验对总体分布进行检验时，我们先提出原假设: H0：总体X的分布函数为F(x)然后根据样本的经验分布和所假设的理论分布之间的吻合程度来决定是否接受原假设。这种检验通常称作拟合优度检验，它是一种非参数检验。 3、Poisson对数线性模型模型：其中，（i=1,2）代表hiv的两个水平，（j=1,2,5）代表factor的5个水平，代表year（代表year的系数），代表age（代表age的系数），代表py（代表py的系数），代表截距。> ap=glm(deaths.,family='poisson&

15、#39;,data=w)> summary(ap)Call:glm(formula = deaths ., family = "poisson", data = w)Deviance Residuals: Min 1Q Median 3Q Max -2.1139 -0.4316 -0.2209 -0.1026 3.2727 Coefficients: Estimate Std. Error z value Pr(>|z|) (Intercept) -23.135255 1.318652 -17.545 < 2e-16 *hiv2 2.766461 0.20

16、3259 13.611 < 2e-16 *factor2 -0.636420 0.151922 -4.189 2.80e-05 *factor3 -0.403434 0.140538 -2.871 0.0041 * factor4 -0.707524 0.142711 -4.958 7.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.025 5.03e-07 *py 0.033042 0.002845 1

17、1.614 < 2e-16 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1(Dispersion parameter for poisson family taken to be 1) Null deviance: 1892.8 on 2143 degrees of freedomResidual deviance: 1007.6 on 2135 degrees of freedomAIC: 1725.7Number of Fisher Scoring iterations: 6得到的模型拟合结果为：在模型中，定性自变量的各个水平的单独效

18、应是不可估计的，必须加上约束条件，这里的约束条件是每个定性变量第一个水平为0。即效应（hiv1）及（factor1）按照R的默认约束条件都等于0。结果分析：首先，各个变量都很显著，相比较而言factor3和factor5的显著性较差一些。其次，当设定hiv1的效应为0时，hiv2对于死亡数的效应为正，且效应比hiv1的效应大；当设定factor1的效应为0时，factor的其余四个选项对于死亡数的效应均为负，且factor4的效应最大，factor5的效应最小；变量year对于死亡数的影响较大，其余两个变量对其影响较小。由模型拟合结果来分析实际情况，可知hiv为阳性时对血友病有较坏的影响，且影

19、响较大；而在使用凝血因子制剂之后，对于病情均有改善，第二种和第四种制剂对于病情的改善效果较为明显；而变量year对于死亡数的影响明显比age和py的影响大，分析情况可能是因为医疗条件的进步，对于血友病的治疗有了明显的帮助。十折交叉验证对测试集的拟合结果： 4、拟似然（quasi-likelihood）模型 对于所有的族，响应变量的方差依赖于均值并且拥有作为系数(multiplier)的尺度参数。方差对均值的依赖方式是响应分布的一个特性；例如对于poisson分布Var(y)=mu。 对于拟似然估计和推断，我们不是设定精确的响应分布而是设定关联函数和方差函数的形式。因为关联函数和方差函数都依赖于

20、均值。#即拟似然模型为响应变量分布没有明确给定> ap=glm(deaths.,family='quasi',data=w)> summary(ap)Call:glm(formula = deaths ., family = "quasi", data = w)Deviance Residuals: Min 1Q Median 3Q Max -0.8530 -0.2895 -0.0874 0.1447 5.2069 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept)

21、-3.3799953 0.2172690 -15.557 < 2e-16 *hiv 0.3783377 0.0236582 15.992 < 2e-16 *factor -0.0435631 0.0086322 -5.047 4.88e-07 *year 0.0340688 0.0024654 13.819 < 2e-16 *age 0.0169687 0.0032860 5.164 2.64e-07 *py 0.0066769 0.0006726 9.927 < 2e-16 *-Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1

22、1(Dispersion parameter for quasi family taken to be 0.2946817) Null deviance: 795.31 on 2143 degrees of freedomResidual deviance: 630.03 on 2138 degrees of freedomAIC: NANumber of Fisher Scoring iterations: 2由于没有明确的分布，这里并不区分分类变量的各个选项，只给出此变量的效应值，得到的模型拟合结果与Poisson对数线性模型基本一致，具体分析在这里不再赘述。十折交叉验证对测试集的拟合结果

23、： 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 =

24、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,

25、alt.levels = NULL, id.var = NULL, opposite = NULL, drop.index = FALSE, ranked = FALSE, .)参数说明：formula:mlogit提供了条件logit，多项logit，混合logit多种模 型，对于多项logit的估计模型应写为：因变量0|自变量， 如：mode0|income。data：先使用mlogit.data函数使得数据结构符合mlogit函数要求。choice：确定分类变量是什么。shape：如果每一行是一个观测，我们选择wide，如果每一行是表示 一个选择，那么就应该选择long。alt.var：

26、对于shape为long的数据，需要标明所有选择名称。> a=mlogit(deaths0|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)Frequencies of alternatives: 0 1 2 3 4 5 6 0.85494403 0.09888060 0.02891791 0.01305970 0

27、.00279851 0.00093284 0.00046642 nr method21 iterations, 0h:0m:9s g'(-H)-1g = 7.61E-07 gradient close to zero Coefficients : Estimate Std. Error t-value Pr(>|t|) 1:(intercept) -2.6616e+01 2.0816e+00 -12.7863 < 2.2e-16 *2:(intercept) -8.2781e+01 2.1635e+04 -0.0038 0.9969471 3:(intercept) -8.

28、9760e+01 1.8693e+04 -0.0048 0.9961687 4:(intercept) -1.1801e+02 1.9020e+04 -0.0062 0.9950497 5:(intercept) -1.3355e+02 1.2532e+04 -0.0107 0.9914969 6:(intercept) -1.3452e+02 1.4785e+04 -0.0091 0.9927407 1:hiv 2.5000e+00 2.2145e-01 11.2893 < 2.2e-16 *2:hiv 2.3683e+01 1.0818e+04 0.0022 0.9982532 3:

29、hiv 2.3319e+01 9.3463e+03 0.0025 0.9980093 4:hiv 2.1054e+01 9.5101e+03 0.0022 0.9982336 5:hiv 1.7966e+01 6.1673e+03 0.0029 0.9976757 6:hiv 2.2628e+01 7.3922e+03 0.0031 0.9975576 1:factor -4.6267e-03 5.8487e-02 -0.0791 0.9369485 2:factor -1.2710e-01 1.0385e-01 -1.2239 0.2209819 3:factor -4.5252e-01 1

30、.6078e-01 -2.8145 0.0048859 * 4:factor -7.3439e-01 3.7890e-01 -1.9382 0.0525968 . 5:factor -1.5480e+01 2.2125e+03 -0.0070 0.9944175 6:factor 6.6159e-01 1.8911e+00 0.3498 0.7264523 1:year 2.1558e-01 2.1647e-02 9.9590 < 2.2e-16 *2:year 3.5500e-01 4.8365e-02 7.3400 2.136e-13 *3:year 4.3400e-01 7.946

31、9e-02 5.4612 4.728e-08 *4:year 7.7129e-01 2.5454e-01 3.0302 0.0024442 * 5:year 1.1373e+00 7.1766e-01 1.5848 0.1130194 6:year 8.7575e-01 6.3633e-01 1.3762 0.1687448 1:age 1.2110e-01 2.4454e-02 4.9520 7.344e-07 *2:age 1.1098e-01 4.7636e-02 2.3299 0.0198138 * 3:age 1.3162e-01 7.4321e-02 1.7710 0.076561

32、1 . 4:age 1.6933e-01 1.5747e-01 1.0753 0.2822502 5:age 2.5431e-01 3.5224e-01 0.7220 0.4703065 6:age -4.8536e-01 9.4106e-01 -0.5158 0.6060251 1:py 2.6013e-02 5.2107e-03 4.9923 5.966e-07 *2:py 5.9033e-02 9.1636e-03 6.4421 1.179e-10 *3:py 7.8478e-02 1.3123e-02 5.9801 2.230e-09 *4:py 1.0475e-01 3.1720e-

33、02 3.3024 0.0009588 *5:py 1.8260e-01 1.0825e-01 1.6869 0.0916288 . 6:py 1.3230e-01 5.0361e-02 2.6271 0.0086127 * -Signif. codes: 0 * 0.001 * 0.01 * 0.05 . 0.1 1Log-Likelihood: -841.11McFadden R2: 0.28468 Likelihood ratio test : chisq = 669.49 (p.value = < 2.22e-16)从输出结果可以看出，对于不同的死亡数，自变量的系数不同。注：1、

34、多项 Logit模型虽然好用，但从上面的叙述可以看出，多项 Logit 模型最大的限制在于各个类别必须是对等的，因此在可供选择的类别中，不可有主要类别和次要类别混杂在一起的情形。例如在研究旅游交通工具的选择时，可将交通工具的类别粗分为航空、火车、公用汽车、自用汽车四大类，但若将航空类别再依三家航空公司细分出三类而得到总共六个类别，则多项 Logit 模型就不适用，因为航空、火车、公用汽车、自用汽车均属同一等级的主要类别，而航空公司的区别则很明显的是较次要的类别，不应该混杂在一起。2、 多项logit模型的因变量没有所列水平之外的可能，即当分类变量有两个以上的水平且这些水平为仅有的可能时，可以考

35、虑多项logit模型。6、 作为比较：用机器学习的算法模型拟合计数因变量数据 1、随机森林拟合数据的十折交叉验证 2、决策树拟合数据的十折交叉验证 7、 各种方法之间的比较各种方法关于测试集十折交叉验证的NMSE回归方法测试集NMSEPoisson对数线性模型0.7391105拟似然模型0.807315随机森林0.6519371决策树0.730649我们在这里所关注的是测试集的标准化均方误差（NMSE）。显然，对于这个数据，经典的计数模型中Poisson对数线性模型较好；在算法模型中，随机森林较好。但总体来说，按照NMSE从优到劣排序为：随机森林、决策树、Poisson对数线性模型、拟似然模型

36、。附：R语言代码：#1.poisson对数线性模型w=read.csv('hemophilia.csv')ap=glm(deaths.,family='poisson',data=w)summary(ap)AIC(ap)#AIC准则越小越好#十折交叉验证NMSE=rep(0,Z)for(i in 1:Z)m=mmia=glm(deaths.,family="poisson",data=w-m,)y1=predict(a,wm,type="response")NMSEi=mean(wm,D-y1)2)/mean(wm,D-m

37、ean(wm,D)2)(MNMSE=mean(NMSE)#2.拟似然模型w=read.csv('hemophilia.csv')ap=glm(deaths.,family='quasi',data=w)summary(ap)#十折交叉验证NMSE=rep(0,Z)for(i in 1:Z)m=mmia=glm(deaths.,family="quasi",data=w-m,)y1=predict(a,wm,type="response")NMSEi=mean(wm,D-y1)2)/mean(wm,D-mean(wm,D)2)(MNMSE=mean(NMSE)#3.多项logit模型i