2022年对广义线性模型的学习 .pdf

《2022年对广义线性模型的学习 .pdf》由会员分享，可在线阅读，更多相关《2022年对广义线性模型的学习 .pdf（20页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

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

2、间，不同统计单位相同统计指标组成的数据列。Note： 与时序数据相比较， 其区别在于数据的排列标准不同，时序数据是按照时间顺序排列的，横截面数据是按照统计单位排列的。横截面数据不要求统计对象及其范围相同，但要求统计的时间相同。#横截面数据即为同一时间截面上的数据2、横截面数据分析的要点：异方差问题名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 1 页，共 20 页 - - - - - - - - - 由于数据是在某一时期对个体或地域的样本的采集，不同个体或地域本身就存在差异。数据的一

3、致性主要包括变量的样本容量是否一致，样本的取样时期是否一致，数据的统计标准是否一致。3、面板数据（ Panel Data ）：是指在时间序列上取多个截面，对于每一个截面上的数据均为一横截面数据列。Note：面板数据是一个m*n 的数据矩阵，记载的是n 个时间节点上， m 个对象的某一数据指标。其有时间序列和截面两个维度，当这类数据按两个维度排列时，是排在一个平面上，与只有一个维度的数据排在一条线上有着明显的不同，整个表格像是一个面板。如果从其内在含义上讲，把panel data译为“时间序列-截面数据”更能揭示这类数据的本质上的特点。4、广义线性模型主要用于因变量取离散值的情况当可能值为一切自

4、然数0,1,2, 时，多用Poisson分布；当 Y取有限个值（实际是响应可以有有限个状态）时，多项分布是自然的选择。5、在很大的程度上可以说，广义线性回归就是针对因变量为有限个值情况的回归分析。但在具体定模型时，需要考虑这有限个状态之间的关系。名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 2 页，共 20 页 - - - - - - - - - 一种是无序的，即各状态的优劣并无公共的认定。 例如外出旅行，有 k 种交通工具可以选择，其优劣取决于具体情况而并无公认的排序。另一种是有

5、序的， 即各状态的优劣次序有公共的认定。如治疗效果、产品质量的分级等。#不同情况建模方法有所不同。二、广义线性模型的提出广义线性模型的提出源于线性模型在应用上有重要影响的几个缺点：1、只适用于因变量Y取值为连续的情况。它特别不适用于分类数据（如Y取 0.1 为值）。2、Y的期望 E(Y) 与自变量 X是用线性关系相联系。选择面太窄，往往与实际情况不符。3、线性模型的统计推断基本上只适用于误差正态的情形。在某些 Y取值连续的场合， Y的分布是偏态的，如指数分布、伽马（Gamma）分布等。广义线性模型的特点 正好是对应上面指出的问题：1、因变量 Y可以取连续值或离散值，从常见的应用看，取离散值的场

6、合更重要。2、取代，有)()(XZhYET函数h（其反函数g称为联系（或连接）函数 （link function） ）有较大的选择余地，这样扩大了模型的适用面。)()(XZYET)()(XZYET名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 3 页，共 20 页 - - - - - - - - - 3、Y（q 维）有指数型分布)()(exp(YdbYY其中，为 q 维参数向量， 是qR上的有限测度， 与 无关（或联系函数1hg使)(XZT，称自然联系）。指数型分布是一个适中的选择，

7、一方面它包括了应用上最常见的一些分布：二项分布、多项分布、Poisson 分布，以及连续型的正态分布、指数分布、伽马分布等。另一方面，这分布类有很好的分析性质，又便于理论上的研究。三、广义线性模型设有因变量 Y，自变量 X，普通线性模型有以下几个特征：1、（线性：线性指对 ，而非 X）。Z(X) 为 X的已知（向量）函数。2、X，Z(X)，Y都是取值连续的变量，如农作物产量、人的身高体重之类。3、Y的分布为正态，或接近正态的分布。广义线性模型从以下几个方面推广：1、)()(XZhYET，h为一严格单调、充分光滑的函数。h已知，1hg（h的反函数）称为联系函数（link function ）,

8、则有TZg)(。即)(YE不等于)(XZT，而是)(XZT的某一函数。2、X，Z(X) ，Y可取连续或离散值，且在应用上更多见的情况为离散值。如 0,1，0,1,2,等。3、Y的分布属于指数型，正态是其一特例。)()(XZYET名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 4 页，共 20 页 - - - - - - - - - 4、以下的表格中列出了GLM中常用的几种分布：由上表格中的第二列（ Range of y）可以知道，当因变量为对应数据形式时应选择对应的分布来建立模型。5

9、、以下的表格中列出了GLM中常用的几种分布所对应的联系函数：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 5 页，共 20 页 - - - - - - - - - 通常称这几种联系函数为标准联系函数，上表中的第三列为偏差。四、 R语言中的模型实现在R语言中利用 stats 包中的 glm() 函数来进行广义线性模型的拟合。和 lm 函数类似， glm 的建模结果可以通过下述的泛型函数进行二次名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - -

10、 - - - - - 名师精心整理 - - - - - - - 第 6 页，共 20 页 - - - - - - - - - 处理，如 summary()、coef() 、confint()、residuals()、anova() 、plot()、predict()。 R提供了一系列广义线性建模工具， 从类型上来说包括gaussian ，反 gaussian ，二项式，poisson 和 gamma 模型的响应变量分布以及在响应变量分布没有明确给定时的拟似然(quasi-likelihood)模型。在后者，方差函数 (variance function)可以认为是均值的函数，但是在另外一些情况

11、下，该函数可以由响应变量的分布得到。函数 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即为选择模型的分布，有以下几种选项：binom

12、ial(link = logit)#二项分布gaussian(link = identity)#正态分布Gamma(link = inverse)#伽马分布inverse.gaussian(link = 1/mu2)#反 Gaussian 分布poisson(link = log)#泊松分布quasi(link = identity, variance = constant)#（quasi-likelihood）#拟家族：响应变量分布没有明确给定时的拟似然模型quasibinomial(link = logit)#拟二项分布名师资料总结 - - -精品资料欢迎下载 - - - - - - -

13、- - - - - - - - - - - 名师精心整理 - - - - - - - 第 7 页，共 20 页 - - - - - - - - - #有过度离散现象时使用：样本观测值变异性过大quasipoisson(link = log)#拟泊松分布#有过度离散现象时使用：样本观测值变异性过大注：若样本观测值变异性过大，即出现了过度离散现象，此时仍使用二项分布假设就会影响系数检测的显著性。那么补救的方法是使用准二项分布（quasibinomial）。首先要检测样本是否存在过度离散现象，方法是用残差除以残差自由度，若超过1 则意味着过度离散。那么将family参数改为quasibinomial

14、。同样，在进行泊松分布也要考虑过度离散现象。其检测方法同样是残差除以其自由度。若确定过度离散存在，则要将family参数设置为准泊松分布(quasipoisson)。在 family的分布选项下还有几个常用选型即link和 variance ，可以用来选择联系函数和方差的形式。Example:glm(y x, family = quasi(variance = mu2, link = log)五、建立广义线性模型的实例1、数据分析：该数据是由美国国家癌症研究所资助的多中心血友病队列研究获得的。该项研究从 1978 年 1 月 1 日到 1995 年 12 月 31 日在 16 个治疗中心（ 1

15、2 个在美国， 4 个在西欧）跟踪了超过1600 个血友病人，该数据一共有 2144 个观测值及 6 个变量。下表为变量情况：名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 8 页，共 20 页 - - - - - - - - - 为了更加直观的分析该数据的特点，截取了原数据中的部分数据行：变量 hiv 为分类变量，只有两个选项， 1 和 2；变量 factor也为分类变量，有五个选项， 1,2,3,4,5；变量 year、age 和 deaths 均为整数数据，只有变量 py 为数

16、量变量。要进行以死亡数即变量deaths 作为因变量的回归，由于因变量为整数数据，因此选择广义线性模型来进行拟合。考察因变量中数据的分布情况：发现可将其看作是 0,1 ，, 或0,1 ，, ，k 的形式，因此我们名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 9 页，共 20 页 - - - - - - - - - 将采用 Poisson 对数线性模型（即分布设定为Poisson 分布, 联系函数设定为对数函数）和多项logit模型（即分布设定为二项分布，联系函数设定为 logit函

17、数）两种方法来进行数据的拟合。2、卡方检验卡方检验法是在总体X 的分布未知时，根据来自总体的样本，检验关于总体分布的假设的一种检验方法。由于这个数据的分布信息是未知的，并且我们也不是很容易直观的判断出它的分布信息， 因此在这里我们采用卡方检验的方法来判断它的分布信息。使用卡方检验对总体分布进行检验时，我们先提出原假设: H0：总体 X的分布函数为 F(x) 然后根据样本的经验分布和所假设的理论分布之间的吻合程度来决定是否接受原假设。这种检验通常称作拟合优度检验，它是一种非参数检验。3、Poisson 对数线性模型模型：其中，i（i=1,2 ）代表 hiv 的两个水平，j（j=1,2, ,5 ）

18、代表factor的 5 个水平，1x代表 year（1代表 year 的系数），2x代表 age（2代表 age 的系数），3x代表 py（3代表 py 的系数），0代表截距。 ap= glm(deaths.,family=poisson,data=w) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 10 页，共 20 页 - - - - - - - - - summary(ap) Call: glm(formula = deaths ., family = poisson, da

19、ta = 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.203259 13.611 2e-16 * factor2 -0.636420 0.151922 -4.189 2.80e-05 * factor3 -0.403434 0.140538 -2.8

20、71 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 11.614 ap=glm(deaths.,family=quasi,data=w) summary(ap) Call: glm(formula = deaths ., family = quasi, d

21、ata = 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) -3.3799953 0.2172690 -15.557 2e-16 * hiv 0.37833770.0236582 15.992 2e-16 * factor -0.04356310.0086322 -5.047 4.88e-07 * year 0.03406880.0024654 13.819

22、 2e-16 * age 0.01696870.0032860 5.164 2.64e-07 * py 0.00667690.0006726 9.927 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.

23、09888060 0.02891791 0.01305970 0.00279851 0.00093284 0.00046642 nr method 21 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.99694

24、71 3:(intercept) -8.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.0

25、022 0.9982532 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 15 页，共 20 页 - - - - - - - - - 3: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.9975

26、576 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.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.72

27、64523 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.9469e-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.2110

28、e-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.0765611 . 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

29、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-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 1 Log-Likel

30、ihood: -841.11 McFadden R2: 0.28468 Likelihood ratio test : chisq = 669.49 (p.value = 2.22e-16) 从输出结果可以看出，对于不同的死亡数，自变量的系数不同。注： 1、多项 Logit模型虽然好用，但从上面的叙述可以看出， 多项 Logit 模型最大的限制在于各个类别必须是对等的，因此在可供选择的类别中，不可有主要类别和次要类别混杂在一起的情形。例如在研究旅游交通工具的选择时，可将交通工具的类别粗分为航空、火车、公用汽车、名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - -

31、 - - - - - - - 名师精心整理 - - - - - - - 第 16 页，共 20 页 - - - - - - - - - 自用汽车四大类， 但若将航空类别再依三家航空公司细分出三类而得到总共六个类别，则多项 Logit 模型就不适用，因为航空、火车、公用汽车、自用汽车均属同一等级的主要类别，而航空公司的区别则很明显的是较次要的类别，不应该混杂在一起。2、多项 logit模型的因变量没有所列水平之外的可能，即当分类变量有两个以上的水平且这些水平为仅有的可能时，可以考虑多项 logit模型。六、作为比较：用机器学习的算法模型拟合计数因变量数据1、随机森林拟合数据的十折交叉验证2、决策

32、树拟合数据的十折交叉验证七、各种方法之间的比较各种方法关于测试集十折交叉验证的NMSE 回归方法测试集 NMSE Poisson 对数线性模型0.7391105 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 17 页，共 20 页 - - - - - - - - - 拟似然模型0.807315 随机森林0.6519371 决策树0.730649 我们在这里所关注的是测试集的标准化均方误差（NMSE ）。显然，对于这个数据， 经典的计数模型中Poisson 对数线性模型较好； 在算法

33、模型中，随机森林较好。但总体来说，按照NMSE 从优到劣排序为：随机森林、决策树、 Poisson 对数线性模型、拟似然模型。附：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=mmi a=glm(deaths.,family=poisson,data=w-m,) y1=predict(a,wm,type=response) NMSEi=

34、mean(wm,D-y1)2)/mean(wm,D-mean(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=mmi a=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) 名师资料

35、总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 18 页，共 20 页 - - - - - - - - - (MNMSE=mean(NMSE) #3. 多项 logit模型install.packages(mlogit) library(mlogit) w=read.csv(hemophilia.csv) w1=mlogit.data(w,choice=deaths,shape=wide) a=mlogit(deaths0|hiv+factor+year+age+py,data=w1) su

36、mmary(a) #十折交叉验证NMSE=rep(0,Z) for(i in 1:Z) m=mmi a=mlogit(deaths0|hiv+factor+year+age+py,data=w1-m,na.action=F) y1=predict(a,w1m,type=response) NMSEi=mean(w1m,D-y1)2)/mean(w1m,D-mean(w1m,D)2) (MNMSE=mean(NMSE) #4. 随机森林拟合数据的十折交叉验证install.packages(randomForest) library(randomForest) Fold=function(Z=1

37、0,w,D,seed=7777) n=nrow(w) d=1:n;dd=list() e=levels(w,D) T=length(e);set.seed(seed) for(i in 1:T) d0=dw,D=ei;j=length(d0) ZT=rep(1:Z,ceiling(j/Z)1:j id=cbind(sample(ZT,length(ZT),d0);ddi=id mm=list();for(i in 1:Z)u=NULL; for(j in 1:T)u=c(u,ddjddj,1=i,2) mmi=u return(mm) for(i in 1:2)w,i=factor(w,i)#

38、将定性变量用factor()标出H=2;Z=10;n=nrow(w)#Z为折数 H=2为选择序号为2 的分类变量mm=Fold(Z,w,H,8888) D=6#因变量位置set.seed(1010) NMSE=rep(0,Z) for(i in 1:Z) m=mmi A=randomForest(deaths.,data=w-m, importance=TRUE,proximity=TRUE) y1=predict(A,wm,) NMSEi=mean(wm,D-y1)2)/mean(wm,D-mean(wm,D)2) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - -

39、 - - - - - - - - - 名师精心整理 - - - - - - - 第 19 页，共 20 页 - - - - - - - - - (MNMSE=mean(NMSE) #5. 决策树拟合数据的十折交叉验证library(rpart) NMSE=rep(0,Z) for(i in 1:Z) m=mmi a=rpart(deaths.,data=w-m,) y1=predict(a,wm,) NMSEi=mean(wm,D-y1)2)/mean(wm,D-mean(wm,D)2) (MNMSE=mean(NMSE) 名师资料总结 - - -精品资料欢迎下载 - - - - - - - - - - - - - - - - - - 名师精心整理 - - - - - - - 第 20 页，共 20 页 - - - - - - - - -