对广义线性模型的学习

对广义线性模型(GeneralizedLinearModeI)的学习引言在学习普通线性模型时就对因变量为离散的情况存有疑问。在统计实验课程研读吴喜之老师的《复杂数据》一 书 关于书的成语关于读书的排比句社区图书漂流公约怎么写关于读书的小报汉书pdf 的第六章时，发现了对离散因变量或者因变量为计数或有序数据时，可采用广义线性模型来处理。因此这燃起了我对于广义线性模型的学习兴趣，通过查阅资料，对此模型有了以下的初步了解。并在对经典 方法 快递客服问题件处理详细方法山木方法pdf计算方法pdf华与华方法下载八字理论方法下载 理论有了一定的了解之后，利用该模型对实际数据进行了处理与 分析 定性数据统计分析pdf销售业绩分析模板建筑结构震害分析销售进度分析表京东商城竞争战略分析 ，同时又用其他方法(包括机器学习等方法)对相同的数据进行了处理，在最后比较了各种方法之间的优缺点。—、数据特点k横截面数据(Cross-SectionDat)a:在同一时间，不同统计单位相同统计指标组成的数据列。Note:①与时序数据相比较，其区别在于数据的排列 标准 excel标准偏差excel标准偏差函数exl标准差函数国标检验抽样标准表免费下载红头文件格式标准下载 不同，时序数据是按照时间顺序排列的，横截面数据是按照统计单位排列的。②横截面数据不要求统计对象及其范围相同，但要求统计的时间相同。#横截面数据即为同一时间截面上的数据2、横截面数据分析的要点:①异方差问 题 快递公司问题件快递公司问题件货款处理关于圆的周长面积重点题型关于解方程组的题及答案关于南海问题 由于数据是在某一时期对个体或地域的样本的采集，不同体或地域本身就存在差异。②数据的一致性主要包括变量的样本容量是否一致，样本的取样时期是否一致，数据的统计标准是否一致。3、面板数据(PanelData):是指在时间序列上取多个截面，对于每一个截面上的数据均为一横截面数据列。Note:①面板数据是一个m*n的数据矩阵，记载的是n个时间节点上，m个对象的某一数据指标。②其有时间序列和截面两个维度，当这类数据按两个维度排列时，是排在一个平面上，与只有一个维度的数据排在一条线上有着明显的不同，整个表格像是一个面板。③如果从其内在含义上讲，把paneldata译为“时间序列-截面数据”更能揭示这类数据的本质上的特点O4、广义线性模型主要用于因变量取离散值的情况当可能值为一切自然数0,1,2,,,时，多用Poisson分布；当丫取有限个值(实际是响应可以有有限个状态)时，多项分布是自然的选择。5、在很大的程度上可以说，广义线性回归就是针对因变量为有限个值情况的回归分析。但在具体定模型时，需要考虑这有限个状态之间的关系一种是无序的，即各状态的优劣并无公共的认定。例如外出旅行，有k种交通工具可以选择，其优劣取决于具体情况而并无公认的排序。另一种是有序的，即各状态的优劣次序有公共的认定。如治疗效果'产品质量的分级等。#不同情况建模方法有所不同。二、广义线性模型的提出广义线性模型的提出源于线性模型在应用上有重要影响的几个缺占?/KW■1、只适用于因变量丫取值为连续的情况。它特别不适用于分类数据(如丫取0.1为值)。2、丫的期望E(Y)与自变量X是用线性关系？•••相联系。选择面太窄，往往与实际情况不符。3、线性模型的统计推断基本上只适用于误差正态的情形。在某些丫取值连续的场合，丫的分布是偏态的，如指数分布、伽马(Gamma分布等。广义线性模型的特点正好是对应上面指出的问题：因变量丫可以取连续值或离散值，从常见的应用看，取离散值的场合更重要。2、取代?小，有E(Y)h(ZT(X))函数h(其反函数g称为联系(或连接)函数(linkfunction))有较大的选择余地，这样扩大了模型的适用面。3、Y(q维)有指数型分布Y~exp(Yb())d(Y)其中，0为q维参数向量，U是R上的。有限测度，U与e无关(或联系函数1Tgh1使ZT(X),称自然联系)。指数型分布是一个适中的选择，一方面它包括了应用上最常见的一些分布：二项分布多项分布、Poisson分布，以及连续型的正态分布'指数分布、伽马分布等。另一方面，这分布类有很好的分析性质，又便于理论上的研究。三、广义线性模型设有因变量丫，自变量X,普通线性模型有以下几个特征：1、—(线性：线性指对P,而非X)。Z(X)为X的已知(向量)函数。2、X,Z(X),丫都是取值连续的变量，如农作物产量、人的身高体重之类。3、丫的分布为正态，或接近正态的分布。广义线性模型从以下几个方面推广：1、E(Y)h(ZT(X)),h为一严格单调、充分光滑的函数。h已1知，gh(h的反函数)称为联系函数(linkfunction),则有g()ZTo即E(Y)不等于ZT(X),而是ZT(X)的某一函数。2、X,Z(X),Y可取连续或离散值，且在应用上更多见的情况为离散值。如{0,1},{0,12,,}等。3、丫的分布属于指数型，正态是其一特例。4、以下的表格中列出了GLM中常用的几种分布由上表格中的第二列(Rangeofy)可以知道，当因变量为对应数据形式时应选择对应的分布来建立模型。5、以下的表格中列出了GLM中常用的几种分布所对应的联系函数通常称这几种联系函数为标准联系函数，上表中的第三列为偏差四、R语言中的模型实现在R语言中利用stats包中的glmO函数来进行广义线性模型的拟合。和Im函数类似，glm的建模结果可以通过下述的泛型函数进行二次处理，女口summary()xcoef()、confint()、residuaIs()、anova()、plot()、predict()。R提供了一系列广义线性建模工从类型上来说包括gaussian，反gaussian，二项式，poisson禾口gamma模型的响应变量分布以及在响应变量分布没有明确给定时的拟似然(quasi-IikeIihood)模型。在后者，方差函数(variancefunction)可以认为是均值的函数，但是在另外一些情况下，该函数可以由响应变量的分布得到。函数glm()的用法：gIm(formula,famiIy二gaussian,data,weights.subset,na.action,start二NULL,etastart,mustart,offset,control=Iist(...),modeI=TRUE,method="glm.fit",x=FALSE,y=TRUE,contrasts=NULL,???)多数选项与普通线性模型的拟合函数lm()相同，值得注意的是family选项，family即为选择模型的分布，有以下几种选项：binomiaI(Iink="logit")#二项分布gaussian(link二Nidentity)”#正态分布Gamma(link二"inverse11)#伽马分布inverse.gaussian(link="1/mu'2")#反Gaussian分布poisson(link=="log")#泊松分布quasi(link="identity",variance="constant")#(CjUdSi—1ik?1!hood#拟二项分布quasibinomiaI(link二"logit”)#拟家族：响应变量分布没有明确给定时的拟似然模型quasibinomial(link二"logit”)#拟二项分布#有过度离散现象时使用：样本观测值变异性过大quasipoisson(link二"log")#拟泊松分布#有过度离散现象时使用：样本观测值变异性过大注：若样本观测值变异性过大，即出现了过度离散现象，此时仍使用二项分布假设就会影响系数检测的显著性。那么补救的方法是使用准二项分布(quasibinomial)。首先要检测样本是否存在过度离散现象，方法是用残差除以残差自由度，若超过1则意味着过度离散。那么将family参数改为quasibinomiaI°同样，在进行泊松分布也要考虑过度离散现象。其检测方法同样是残差除以其自由度。若确定过度离散存在，则要将family参数设置为准泊松分布(quasipoisson)o在family的分布选项下还有几个常用选型即link和variance,可以用来选择联系函数和方差的形式。ExampIe:gIm(y〜x,family=quasi(varianee="mu'2",link="Iog"))五、建立广义线性模型的实例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作为因变量的回归，由于因变量为整数数据，因此选择广义线性模型来进行拟合。考察因变量中数据的分布情况：发现可将其看作是｛0,1,,,｝或k｝的形式，因此我们｛0,15,,,将采用Poisson对数线性模型（即分布设定为Poisson分布,联系函数设定为对数函数）和多项logit模型（即分布设定为二项分布，联系函数设定为logit函数）两种方法来进行数据的拟合。2、卡方检验卡方检验法是在总体X的分布未知时，根据来自总体的样本，检验关于总体分布的假设的一种检验方法。由于这个数据的分布信息是未知的，并且我们也不是很容易直观的判断出它的分布信息，因此在这里我们采用卡方检验的方法来判断它的分布信息。使用卡方检验对总体分布进行检验时，我们先提出原假设：H0:总体X的分布函数为F（x）然后根据样本的经验分布和所假设的理论分布之间的吻合程度来决定是否接受原假设。这种检验通常称作拟合优度检验，它是一种非参数检验3、Poisson对数线性模型模型:其中，i（i=1,2）代表hiv的两个水平，」（j=1,2,,,,5）代表factor的5个水平，xi代表year（丨代表year的系数）X2代表age（2代表age的系数）,X3代表py（3代表py的系数)，。代表截距。>ap=glm(deaths;family=fpoisson'fdatw)a>summary(ap)Call:IIglm(formula—二deaths〜family—二poissorT.data—二w)DevianceResiduals:Min1QMedian3QMax-2.1139-0.4316-0.2209-0.10263.2727Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)-23.1352551.31865217.545<2e-16***hiv22.7664610.20325913.611<2e-16***factor2-0.6364200.151922-4.1892.80e-05***factor3-0.4034340.140538-2.8710.0041**factor4-0.7075240.142711-4.9587.13e-07***factor5-0.3715850.146238-2.5410.0111*year0.2110470.01409014.979<2e-16***age0.0778670.0154955.0255.03e-07***py0.0330420.00284511.614<2e-16***Signif.codes:06***10.0010.010.(Dispersionparameterforpoissonfamilytakentobe1)NulIdevianee:1892.8on2143degreesoffreedomResidualdevianee:1007.6on2135degreesoffreedomAlC:1725.7NumberofFisherScoringiterations:6得到的模型拟合结果为：In()23.1402.7700.640.400.710.370.21xi0.08x20.03x3在模型中，定性自变量的各个水平的单独效应是不可估计的，必须加上约束条件，这里的约束条件是每个定性变量第一个水平为00即效应1(hi)及1(factorl)按照R的默认约束条件都等于0结果分析：首先，各个变量都很显著，相比较而言factor3和factor5的显著性较差一些其次，当设定hiv1的效应为0时，hiv2对于死亡数的效应为正，且效应比hiv1的效应大；当设定factoid的效应为0时，factor的其余四个选项对于死亡数的效应均为负，且factor4的效应最大，factor5的效应最小；变量year对于死亡数的影响较大，其余两个变量对其影响较小。由模型拟合结果来分析实际情况，可知hiv为阳性时对血友病有较坏的影响，且影响较大；而在使用凝血因子制剂之后，对于病情均有改善，第二种和第四种制剂对于病情的改善效果较为明显；而变量year对于死亡数的影响明显比age和py的影响大，分析情况可能是因为医疗条件的进步，对于血友病的治疗有了明显的帮助。十折交叉验证对测试集的拟合结果:4、拟似然(quasi-likeIihood)模型对于所有的族，响应变量的方差依赖于均值并且拥有作为系数(multiplier)的尺度参数。方差对均值的依赖方式是响应分布的一个特性；例如对于poisson分布Var(y)=muo对于拟似然估计和推断，我们不是设定精确的响应分布而是设定关联函数和方差函数的形式。因为关联函数和方差函数都依赖于均#即拟似然模型为响应变量分布没有明确给定>ap=glm(deaths；famiIy=,quasi1,data=w)>summary(ap)Call:gIm(formula=deaths~famiIy=”quasi，”data=w)DevianceResiduaIs:Min1QMedian3QMax-0.8530-0.2895-0.08740.14475.2069Coefficients:0.1Signif.codes:0'***0.001'**'0.01l*0.05EstimateStd.ErrortvaIuePr(>|t|)(Intercep-3.37999530.2172690-15.557<2e~16***hiv0.37833770.023658215.992<2e-16factor-0.04356310.00863225047488e~--07year0.03406880.002465413.819<2eT6age0.01696870.00328605-1642-64e07py0.00667690.00067269.927<2e-16(Dispersionparameterforquasifamilytakentobe0.2946817)NuIIdevianee:795.31on2143degreesoffreedomResiduaIdevianee:630.03on2138degreesoffreedomAIC:NANumberofFisherScoringiterations:2由于没有明确的分布，这里并不区分分类变量的各个选项，只给出此变量的效应值，得到的模型拟合结果与Poisson对数线性模型基本一致，具体分析在这里不再赘述。十折交叉验证对测试集的拟合结果:5、多项logit模型模型:多项logit模型在类别上仍可归为广义线性模型，是二项分布的logistic回归向多项分布的推广，但是在R语言的glm()函数中只能进行二项分布的回归，而无法进行多项分布的回归。因此我们利用R语言mlogit包中的mlogit()函数来进行模型的拟合。函数mlogit的用法：mlogit(formula,data,subset,weights,na.action,start二NULL,alt.subset二NULL,refIeveI=NULL,nests=NULL,un.nest,eI=FALSE,unscaIed=FALSE,heterosc=FALSE,rpar=NULL,probit二FALSE,R=40,correIation=FALSE,haIton=NULL,random,nb=NULL,paneI=FALSE,estimate=TRUE,seed=10,...)mlogit.data(data,choice,shape=c("wide","Iong"),varying=NULL,sep=H.",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,如果每一行是表示一个选择，那么就应该选择longoalt.var:对于shape为long的数据，需要标明所有选择名称。a=mIogit(deaths~01hiv+factor+year+age+py,data=w1)summary(a)Call:mlogit(formula=deaths0|hiv+factor+year+age+pyTdata=w1,method="nr",print,level二0)Frequenciesofalternatives:01234560.854944030.098880600.028917910.013059700.002798510.000932840.00046642nrmethod21iterations,Oh:0m:9sg(-H)‘八Tg=7.61E-07gradientcIosetozeroCoefficients:EstimateStd.Errort一vaIue1:(intercept)-2.6616e+012.0816e+00-12.7863<2.2e-16***2:(intercept)-8.2781e+012.1635e+04-0.00380.99694713:(intercept)4:-8.9760e+011.8693e+04-0.00480.9961687(intercept)5:-1.1801e+021?9020e+04-0.00620.9950497(intercept)6:-1?3355e+021.2532e+04-0.01070.9914969(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.99825323hiv2.3319e+019.3463e+030.00250.99800934hiv2.1054e+019.5101e+030.00220.99823365hiv1.7966e+016.1673e+030.00290.99767576hiv2.2628e+017.3922e+030.00310.99755761facto-4.6267e〜5.8487e-02-0.07910.93694852facto-1.2710e-1.0385e-01-1.22390.22098193facto-4.5252e〜1.6078e-01-2.81450.00488594facto-7.3439e-3.7890e-01-1.93825facto2.2125e+03-0.00700.99441756facto6.6159e-011.8911e+000.34980.72645231year2.1558e-012.1647e-029.9590<2.2e-162year3.5500e-014.8365e-027.34002.136e-133year4.3400e-017.9469e-025.46124.728e-084year7.7129e-012.5454e-013.03020.00244425year1.1373e+007.1766e〜011.58480.11301946year8.7575e-016.3633e-011.37620.16874481age1.2110e-012.4454e-024.95207.344e-072age1.1098e-014.7636e-022.32990.0198138*3age1.3162e-017.4321e-021.77100.0765611.4age1.6933e-011.5747e-011.07530.28225025age2.5431e-013.5224e-010.72200.47030656age-4.8536e~9.4106e-01-0.51580.60602511py2.6013e-025.2107e-034.99235.966e-072py5.9033e-029.1636e-036.44211.179e-103py7.8478e-021.3123e-025.98012.230e-094py1.0475e-013.1720e-023.30240.00095885py1.8260e-011.0825e-011.68690.0916288'6py1.3230e-015.0361e-022.62710.00861270.010.05\10.1'11Signif.codes:0'***10.001Log-Likelihood:-841.11McFadden旷2:0.28468LikeIihoodratiotestchisq=669.49(p.vaIue=<2.22e〜16)从输出结果可以看出，对于不同的死亡数，自变量的系数不同。注：k多项Logit模型虽然好用，但从上面的叙述可以看出，多项Logit模型最大的限制在于各个类别必须是对等的，因此在可供选择的类别中，不可有主要类别和次要类别混杂在一起的情形。例如在研究旅游交通工具的选择时，可将交通工具的类别粗分为航空、火车、公用汽车、自用汽车四大类，但若将航空类别再依三家航空公司细分出三类而得到总共六个类别，则多项Logit模型就不适用，因为航空、火车、公用汽车'自用汽车均属同一等级的主要类别，而航空公司的区别则很明显的是较次要的类别，不应该混杂在一起。2、多项logit模型的因变量没有所列水平之外的可能，即当分类变量有两个以上的水平且这些水平为仅有的可能时，可以考虑多项logit模型。六、作为比较：用机器学习的算法模型拟合计数因变量数据仁随机森林拟合数据的十折交叉验证2、决策树拟合数据的十折交叉验证七、各种方法之间的比较各种方法关于测试集十折交叉验证的NMSE回归方法测试集NMSEPoisson对数线性模型0.7391105拟似然模型0.807315随机森林0.6519371决策树0.730649我们在这里所关注的是测试集的标准化均方误差(NM$E。显然，对于这个数据，经典的计数模型中Poisson对数线性模型较好；在算法模型中，随机森林较好。但总体来说，按照NMS从优到劣排序为：随机森林、决策树、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(iin1:Z){m=mm[[i]]a=glm(deaths:,family="poisson",data二w[-m,])y1=predict(a,w[m,],type二"response")NMSE[i]二mean((w[m,D]-y1厂2)/mean((w[m,D]-mean(w[m,D])厂2)}(MNMSE二mean(NMSE))#2.拟似然模型w二read,csv('hemophilia.csv')ap=glm(deaths~.,family二'quasi',date二w)summary(ap)#十折交叉验证NMSE=rep(0,Z)for(iin1:Z){m=mm[[i]]a=glm(deaths".,family="quasi",data—二w[-m,])y1=predict(a,w[m,],type^'response")NMSE[i]—二mean((w[m,D]-y1厂2)/mean((w[m,D]-mean(w[m,D])厂2)}(MNMSE二mean(NMSE))#3.多项logit模型instalI.packages("mlogit")Iibrary(mlogit)w二read,csv('hemophilia,csv')w1=mlogit.data(w,choice二"deaths",shape二"wide”)a二mlogit(deaths~0|hiv+factor+year+age+py,data=w1)summary(a)#十折交叉验证NMSE=rep(0,Z)for(iin1:Z){a二mlogit(deaths~0|hiv+factor+year+age+py,data=w1[-m,],na.action二F)y1二predict(a,w1[m,],type="response")NMSE[i]二mean((w1[m,D]-y1)"2)/mean((w1[m,D]-mean(w1[m,D])厂2)}(MNMSE=mean(NMSE))#4?随机森林拟合数据的十折交叉验证instalI.packages("randomForest")Iibrary(randomForest)FoIdefunction(Z=10,w,D,seed=7777){n=nrow(w)d=1:n;dd=list()e=levels(w[,D])T=length(e);set.seed(seed)for(iin1:T){dO=d[w[,D]=e[i]];j二Iength(dO)ZT=rep(1:Z,ceiling(j/Z))[1:j]id=cbind(sampIe(ZT,Iength(ZT)),dO)；dd[[i]]二id}mm二Iist();for(iin1:Z){u二NULL;for(jin1订)u=c(u,dd[[j]][dd[[j]][,1]==i,2])mm[[i]]=u}return(mm)}for(iin1:2)w[,i]=factor(w[,i])#将定性变量用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(iin1:Z){m=mm[[i]]A=randomForest(deaths~.,data二w[-m,],importance二TRUE,proximity二TRUE)y1=predict(A,w[m,])NMSE[i]二mean((w[m,D]-y1厂2)/mean((w[m,D]-mean(w[m,D])厂2)}(MNMSE=mean(NMSE))#5.决策树拟合数据的十折交叉验证Iibrary(rpart)NMSE二rep(0,Z)for(iin1:Z){m=mm[[i]]a=rpart(deathsA.,data=w[-m,])y1=predict(a,w[m,])NMSE[i]二mean((w[m,D]-y1厂2)/mean((w[m,D]-mean(w[m,D])厂2)}(MNMSE二mean(NMSE))