《多元统计分析-基于R（第3版）》课件 第5、6章-多元线性模型、广义线性模型

5.1多元正态分布5.2多元线性模型5.3变量选择5.4回归诊断5.5回归预测第5章多元线性模型2主编：费宇2024/7/315.1.1多元正态分布的定义5.1多元正态分布2024/7/313主编：费宇5.1.1多元正态分布的定义5.1多元正态分布2024/7/314主编：费宇5.1.2多元正态分布的性质5.1多元正态分布2024/7/315主编：费宇5.1.2多元正态分布的性质5.1多元正态分布2024/7/316主编：费宇5.2多元线性模型5.2.1模型定义其中x1,…,xp是非随机的自变量,y是随机的因变量,β0是常数项,β1,…,βp是回归系数,ε是随机误差项.7主编：费宇2024/7/315.2.1.模型定义模型(5.1)的样本形式其中i=1,…,n,表示有n组观测值.8主编：费宇2024/7/315.2.1.模型定义模型(5.2)的矩阵形式其中9主编：费宇2024/7/31例5.1(数据文件为eg5.1)10主编：费宇表5-1抽样调查得到的36个学生的相关成绩2024/7/31yx1x2x3x4x5yx1x2x3x4x5858386909076456065608678909288879280768081758075787076738590888582868580807281829088828081868790………………………………928385908580627865608588788482739083878083857883模型（5.3）的参数的最小二乘估计的最小二乘估计

5.2.2模型的参数估计和检验2024/7/3111主编：费宇1.回归方程的显著性检验其中是回归平方和，而

是残差平方和，拒绝域为

5.2.2模型的参数估计和检验2024/7/3112主编：费宇5.2.2模型的参数估计和检验2024/7/3113主编：费宇例5.1续1(数据文件为eg5.1)14主编：费宇建立y关于x1、x2、x3、x4和x5的线性回归方程,并对方程和回归系数进行显著性检验.2024/7/31yx1x2x3x4x5yx1x2x3x4x5858386909076456065608678909288879280768081758075787076738590888582868580807281829088828081868790………………………………928385908580627865608588788482739083878083857883#例5.1回归分析：全变量回归setwd("C:/data")#设定工作路径d5.1<-read.csv("exam5.1.csv",header=T)#将exam5.1.csv数据读入到d5.1中lm.exam<-lm(y~x1+x2+x3+x4+x5,data=d5.1)#建立y关于x1,x2,x3,x4和x5的线性回归方程,数据为d5.1summary(lm.exam)#给出回归系数的估计和显著性检验等例5.1续1(数据文件为eg5.1)15主编：费宇R程序为:Call:lm(formula=y~x1+x2+x3+x4+x5,data=d5.1)Residuals:Min1QMedian3QMax-10.0696-1.7983-0.15352.93616.8726

Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)-32.7353415.35701-2.1320.0413*x10.162710.150311.0820.2877x20.227840.138351.6470.1100x30.881160.111087.9337.46e-09***x4-0.051360.15476-0.3320.7423x50.168870.143761.1750.2494回归分析结果为2024/7/3116主编：费宇---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandarderror:4.021on30degreesoffreedomMultipleR-squared:0.8945,AdjustedR-squared:0.877F-statistic:50.89on5and30DF,p-value:9.359e-14例5.1续(数据文件为eg5.1)2024/7/3117主编：费宇回归方程的F值为50.89,相应的p值为9.359

10-14,说明回归方程是显著的;但t检验对应的p值则显示:常数项和x3是显著的,而x1、x2、x4和x5不显著.最优模型一般满足2个条件（1）模型反映了变量间的真实关系

（2）模型包含的变量尽量少例5.1（续2）建立y关于x1、x2、x3、x4和x5的线性回归方程,并对方程和回归系数进行显著性检验（逐步回归建立“最优方程”）.

5.3变量选择2024/7/3118主编：费宇>#例5.1回归分析：逐步回归>#假设exam5.1.xls中的数据已经读入到d5.1中>lm.exam<-lm(y~x1+x2+x3+x4+x5,data=d5.1)#建立全变量回归方程>lm.step<-step(lm.exam,direction="both")#进行逐步回归5.3变量选择2024/7/3119主编：费宇例5.1(续2)逐步回归程序Start:AIC=105.63y~x1+x2+x3+x4+x5

DfSumofSqRSSAIC-x411.78486.83103.76-x1118.95503.99105.01-x5122.31507.36105.25<none>485.05105.63-x2143.85528.90106.74-x311017.441502.49144.335.3变量选择2024/7/3120主编：费宇例5.1(续)回归结果Step:AIC=103.76y~x1+x2+x3+x5

DfSumofSqRSSAIC-x1117.91504.73103.06-x5120.57507.40103.25<none>486.83103.76-x2142.99529.81104.80+x411.78485.05105.63-x311112.961599.79144.595.3变量选择2024/7/3121主编：费宇Step:AIC=103.06y~x2+x3+x5

DfSumofSqRSSAIC-x5117.40522.14102.28<none>504.73103.06+x1117.91486.83103.76+x410.74503.99105.01-x2170.76575.50105.78-x311848.492353.23156.485.3变量选择2024/7/3122主编：费宇Step:AIC=102.28y~x2+x3

DfSumofSqRSSAIC<none>522.14102.28+x5117.40504.73103.06+x1114.74507.40103.25+x410.25521.89104.26-x2166.64588.78104.60-x311953.302475.43156.305.3变量选择2024/7/3123主编：费宇Call:lm(formula=y~x2+x3,data=d5.1)Residuals:Min1QMedian3QMax-10.4395-2.5508-0.44592.73677.2345

Coefficients:EstimateStd.ErrortvaluePr(>|t|)(Intercept)-18.842907.58902-2.4830.0183*x20.249230.121442.0520.0481*x30.968040.0871311.1111.09e-12***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

Residualstandarderror:3.978on33degreesoffreedomMultipleR-squared:0.8865,AdjustedR-squared:0.8796F-statistic:128.8on2and33DF,p-value:2.566e-16回归模型汇总信息:summary(lm.step)2024/7/3124主编：费宇5.4.1残差分析和异常点探测

残差向量e是模型中随机误差项

的估计，残差分析可以诊断模型的基本假定是否成立。

5.4

回归诊断2024/7/3125主编：费宇例5.2

计算例5.1得到的逐步回归模型lm.step的普通残差和标准化残差,判断可能存在的异常点,画出相应的残差散点图,并直观判断模型的基本假定是否成立.5.4

回归诊断2024/7/3126主编：费宇#例5.2#假设由例5.1已经得到逐步回归模型lm.stepy.res<-residuals(lm.exam)#计算模型lm.exam的普通残差y.rst<-rstandard(lm.step)#计算回归模型lm.step的标准化残差print(y.rst)#输出回归模型lm.step的标准化残差y.rsty.fit<-predict(lm.step)#计算回归模型lm.step的预测值plot(y.res~y.fit)#绘制以普通残差为纵坐标，预测值为横坐标的散点图plot(y.rst~y.fit)#绘制以标准化残差为纵坐标，预测值为横坐标的散点图5.4

回归诊断2024/7/3127主编：费宇分别采用residuals(),rstandard()和rstudent()来计算普通残差,标准化残差和学生化残差.123456-1.226479490.701233481.85465439-0.18487397-0.731575470.14591132789101112-0.651653781.37662024-0.28171298-0.96473838-0.798622470.81284419131415161718-0.48393343-1.176685880.913377160.564389020.658766891.49006874192021222324-2.871217390.527102680.81076269-0.668013511.20184149-1.040201892526272829300.32282704-0.04616114-0.159120010.21602487-0.21306706-0.23026109313233343536-0.24302334-2.03567204-0.33183300-0.073548931.804380090.737029325.4

回归诊断2024/7/3128主编：费宇回归模型lm.step的标准化残差y.rst如下第19号点是异常点2024/7/3129主编：费宇图5-1例5.2中的普通残差图(左)和标准化残差图(右)例5.3

通过方差稳定变换来更新例5.1得到的逐步回归模型lm.step,并计算更新后模型的标准化残差,画出相应的残差散点图,并直观判断模型的基本假定是否成立.5.4

回归诊断2024/7/3130主编：费宇#例5.3#假设由例5.1已经得到逐步回归模型lm.steplm.step_new<-update(lm.step,log(.)~.)#对模型进行对数变换y.rst<-rstandard(lm.step_new)#计算lm.step_new的标准化残差y.fit<-predict(lm.step_new)#计算lm.step_new的预测值plot(y.rst~y.fit)#绘制以标准化残差为纵坐标,预测值为横坐标的散点图（见图5-2）采用对数变换来解决方差非齐问题2024/7/3131主编：费宇对数变换后：第19号点是异常点2024/7/3132主编：费宇图5-2例5.3中的标准化残差图lm.exam<-lm(log(y)~x1+x2+x3+x4,data=d5.1[-c(19),])#去掉第19号观测值再建立全变量回归方程lm.step<-step(lm.exam,direction="both")#用一切子集回归法来进行逐步回归y.rst<-rstandard(lm.step)#计算回归模型lm.step的标准化残差y.fit<-predict(lm.step)#计算回归模型lm.step的预测值plot(y.rst~y.fit)#绘制以标准化残差为纵坐标，预测值为横坐标的散点图去掉19号观测值再回归2024/7/3133主编：费宇残差几乎全部落在[-2，2]区域内2024/7/3134主编：费宇图5-3例5.3中的标准化残差图:去掉19号观测值5.4.2回归诊断：一般的方法残差分析无法分析模型的影响点,即探测哪些点对模型的推断有重要影响,本节给出的回归诊断方法,可以诊断模型的基本假定是否成立,哪些值是异常点,哪些点是强影响点.在R中,函数plot()和influence.measures()可以用来绘制诊断图和计算诊断统计量.5.4

回归诊断2024/7/3135主编：费宇例5.4对例5.3得到的逐步回归模型lm.step_new进行回归诊断分析.5.4回归诊断2024/7/3136主编：费宇#例5.4#假定由例5.3已经获得模型lm.step_newpar(mfrow=c(2,2))#在一个2×2网格中创建4个绘图区plot(lm.step_new)#绘制模型诊断图influence.measures(lm.step_new)#计算各个观测值的诊断统计量运行上述程进行序可得回归诊断图(图5-4)和如下36个观测值对应的诊断统计量的值.第11，19和33号观测值被诊断为强影响点.2024/7/3137主编：费宇Influencemeasuresoflm(formula=log(y)~x2+x3，data=d5.1)：dfb.1_dfb.x2dfb.x3dffitcov.rcook.dhatinf10.172353-0.052013-1.36e-01-0.291711.0522.82e-020.06622-0.0789410.0622346.53e-030.110451.1604.17e-030.069130.196029-0.049262-1.17e-010.373080.8364.31e-020.0383…110.014319-0.1934602.46e-01-0.270291.3462.49e-020.2065*12-0.0378780.0103993.07e-020.060361.1881.25e-030.0806…19-2.0235661.0016117.40e-01-2.218990.2329.50e-010.1645*20-0.0066640.071366-8.14e-020.147121.0937.33e-030.0419…33-0.0808070.140284-1.13e-01-0.152341.4647.96e-030.2556*第3,19和35号观测值可能是异常点和强影响点.2024/7/3138主编：费宇图5-4例5.4的回归诊断图回归预测分为点预测和区间预测两种,可以采用函数predict()来实现.例5.5

给定解释变量x2=80，x3=90，利用例5.1得到的回归模型对y进行点预测和区间预测(置信度为95%).5.5回归预测2024/7/3139主编：费宇#例5.5#假定由例5.1已经获得模型lm.steppreds<-data.frame(x2=80,x3=90)#给定解释变量x2和x3的值predict(lm.step,newdata=preds,interval="prediction",level=0.95)#进行点预测和区间预测5.5回归预测2024/7/3140主编：费宇点预测和区间预测的程序如下:>predict(lm.step,newdata=preds,interval="prediction",level=0.95)#区间预测fitlwrupr188.2190779.7974896.640675.5回归预测2024/7/3141主编：费宇运行上述程序可得y的点预测和区间预测的结果如下:程序中选项interval="prediction"表示要给出预测区间,选项level=0.95表示置信水平是95%.计算结果y的点预测为88.22,预测区间为[79.80,96.64].2024/7/31主编：费宇2024/7/31主编：费宇43

2024/7/31主编：费宇44第6章广义线性模型6.1广义线性模型的定义6.2Logistic模型6.3Probit模型6.4多项Logit模型6.5泊松对数线性模型6.6零膨胀计数模型6.7多项分布对数线性模型2024/7/31主编：费宇456.1

广义线性模型的定义第5章我们研究了多元线性模型,该模型的一个重要假定是因变量是连续型的变量(通常假定服从正态分布),但在许多情况下,这种假定并不合理,例如下面这两种情况.(1)结果变量可能是类型变量.二值分类变量和多分类变量.(比如:是/否,差/一般/良好/优秀等)显然都不是连续型变量.2024/7/31主编：费宇466.1广义线性模型概述(2)结果变量可能是计数型变量(比如:一周交通事故的数目)这类变量都是非负的有限值,而且它们的均值和方差通常是相关的(一般线性模型假定因变量是正态变量,而且相互独立).普通线性回归模型(5.3)假定因变量y服从正态分布,其均值满足关系式:μ=Xβ,这表明因变量的条件均值是自变量的线性组合.本章介绍六种常见的广义线性模型：Logistic模型、Probit模型、多项Logit模型、泊松对数线性模型、零膨胀计数模型和多项分布对数线性模型.2024/7/31主编：费宇476.1广义线性模型的定义1.广义线性模型的定义：(1)随机成分:设y1,y2,…,yn是来自于指数分布族的随机样本,即yi的密度函数为其中ai(.),b(.),ci(.)是已知函数,参数αi是典则参数,ϕ是散度参数.2024/7/31主编：费宇481.广义线性模型的定义：(2)连接函数:设yi的均值为μi而函数g(.)是单调可微的连接函数,使得其中是协变量,是未知参数向量.指数分布族正态分布二项分布泊松分布2024/7/31主编：费宇492024/7/31主编：费宇502.正态线性回归模型正态分布属于指数分布族,其密度函数为与(6.1)对照可知2024/7/31主编：费宇512.正态线性回归模型只要取连接函数为,则正态线性回归模型满足广义线性模型的定义.类似的,容易验证,二项分布和泊松分布都属于指数分布族.下面介绍实际中应用广泛的Logistic模型、Probit模型、多项Logit模型、泊松对数线性模型、零膨胀计数模型和多项分布对数线性模型.2024/7/31主编：费宇526.2Logistic模型1.模型定义 设yi服从参数为pi的二项分布,则μi=E(yi)=pi采用逻辑连接函数,即这个广义线性模型称为Logistic模型.2024/7/31主编：费宇53例6.1(数据文件为eg6.1)表6-1

某城市48个家庭的调查数据2024/7/31主编：费宇542.模型的参数估计和检验采用R软件中的广义线性模型过程glm()可以完成回归系数的估计,以及模型回归系数的显著性检验.程序如下：#例6.1广义线性模型:Logistic模型setwd("C:/data")#设定工作路径d6.1<-read.csv("exam6.1.csv",header=T)#将exam6.1.csv数据读入到d6.1中glm.logit<-glm(y~x1+x2,family=binomial(link=logit),data=d6.1)#建立y关于x1，x2的logistic回归模型，数据为d6.1summary(glm.logit)#模型汇总注意逻辑连接函数是二项分布的典则连接函数，是默认的连接函数，因此代码中的(link=logit)可以省略.2024/7/31主编：费宇55运行以上程序可得如下结果：Call:glm(formula=y~x1+x2，family=binomial(link=logit)，data=d6.1)DevianceResiduals:Min1QMedian3QMax-2.30297-0.198320.022830.202511.59258Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)-7.531152.56352-2.9380.00331**x10.439560.138643.1700.00152**x2-0.081031.24747-0.0650.94821---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1(Dispersionparameterforbinomialfamilytakentobe1)Nulldeviance:61.105on47degreesoffreedomResidualdeviance:17.643on45degreesoffreedomAIC:23.643NumberofFisherScoringiterations:82024/7/31主编：费宇56注意到x2对应的p值(0.948)比较大，即x2不显著，所以考虑采用逐步回归.glm.step<-step(glm.logit)#逐步回归summary(glm.step)#给出模型回归系数的估计和显著性检验等运行以上程序可得如下结果Start:AIC=23.64y~x1+x2

DfDevianceAIC-x2117.64721.647<none>17.64323.643-x1159.00863.008

Step:AIC=21.65y~x1

DfDevianceAIC<none>17.64721.647-x1161.10563.105>summary(glm.step)#给出模型回归系数的估计和显著性检验等2024/7/31主编：费宇57Call:glm(formula=y~x1,family=binomial(link=logit),data=d6.1)DevianceResiduals:Min1QMedian3QMax-2.28859-0.197030.022760.204001.60887

Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)-7.56822.5101-3.0150.00257**x10.43960.13873.1690.00153**---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

(Dispersionparameterforbinomialfamilytakentobe1)

Nulldeviance:61.105on47degreesoffreedomResidualdeviance:17.647on46degreesoffreedomAIC:21.647

NumberofFisherScoringiterations:8模型预测如果要预测年收入为20万元（x1=20）、家里有孩子（x2=1）的家庭有购买住房的可能性，可以采用以下命令：>yp<-predict(glm.step,data.frame(x1=20))>p.fit<-exp(yp)/(1+exp(yp));p.fit#估计x1=20时y=1的概率10.7728122容易看出，当x1=20，x2=1时，估计y=1的概率约为0.77，即年收入为20万元、家里有孩子的家庭有购买住房的可能性约为77%.2024/7/31主编：费宇586.3Probit模型2024/7/31主编：费宇596.3Probit模型:例6.1（续）2024/7/31主编：费宇60运行以上程序可得如下结果：Call:glm(formula=y~x1,family=binomial(link=probit),data=d6.1)

DevianceResiduals:Min1QMedian3QMax-2.2493-0.15220.00180.17681.6024

Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)-4.340281.27539-3.4030.000666***x10.249890.069443.5990.000320***---Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1

(Dispersionparameterforbinomialfamilytakentobe1)

Nulldeviance:61.105on47degreesoffreedomResidualdeviance:17.349on46degreesoffreedomAIC:21.349

NumberofFisherScoringiterations:92024/7/31主编：费宇612024/7/31主编：费宇626.4多项Logit模型2024/7/31主编：费宇63例6.2某城市48个家庭的调查数据2024/7/31主编：费宇64问题：根据这个数据建立多项分布回归模型并估计年收入为20万元、家里有孩子的家庭有购买住房但还在还贷款的可能性.解:采用nnet程序包中的multinom()可以完成多项logit模型的拟合.#例6.2广义线性模型:多项分布回归模型library(nnet)setwd("C:/data")d6.2<-read.csv("exam6.2.csv",header=T)d6.2$x2<-as.factor(d6.2$x2)#将x2这一列因子化mlog<-multinom(y~x1+x2,data=d6.2)#建立模型summary(mlog)#查看所拟合的模型2024/7/31主编：费宇65运行以上程序可得如下结果：Call:multinom(formula=y~x1+x2，data=d6.2)Coefficients:(Intercept)x1x22-7.4438920.4329375-0.067896533-17.3785220.7438569-0.57429520Std.Errors:(Intercept)x1x222.5703380.13962821.24601334.4477300.18612381.704516ResidualDeviance:37.79579AIC:49.79579注意到x2对应标准误相对于x2的系数比较大，所以估计x2可能不显著，采用step()函数对模型进行逐步回归.2024/7/31主编：费宇66mlog.s<-step(mlog)#对mlog进行逐步回归summary(mlog.s)#查看所拟合的模型运行以上程序可得如下结果：Call:multinom(formula=y~x1,data=d6.2)

Coefficients:(Intercept)x12-7.4794080.43324433-17.2933710.7313709

Std.Errors:(Intercept)x122.5180900.139753034.4241140.1834096

ResidualDeviance:37.98674AIC:45.986742024/7/31主编：费宇672024/7/31主编：费宇682024/7/31主编：费宇692024/7/31主编：费宇70估计48个家庭最可能属于3类家庭中的哪一类？2024/7/31主编：费宇716.5泊松对数线性模型1.模型的定义设y服从参数λ为的泊松分布,则μ=E(y)=λ,采用对数连接函数,即这个广义线性模型称为泊松对数线性模型.2024/7/31主编：费宇72例6.3(数据文件为eg6.3)表6-3

Breslow癫痫数据2024/7/31主编：费宇73例6.3(数据文件为eg6.3)这个数据是robust包中的Breslow癫痫数据(Breslow,1993).我们讨论在治疗初期的八周内,癫痫药物对癫痫发病数的影响,响应变量为八周内癫痫发病数(y),预测变量为前八周内的基础发病次数(x1),年龄(x2)和治疗条件(x3),其中治疗条件是二值变量,x3=0表示服用安慰剂,x3=1表示服用药物.根据这个数据建立泊松对数线性模型并对模型的系数进行显著性检验.2024/7/31主编：费宇742.模型的参数估计和检验采用R软件中的广义线性模型过程glm()来建立泊松对数线性模型并对模型的系数进行显著性检验.程序如下：#例6.3广义线性模型:泊松对数线性模型setwd("C:/data")d6.3<-read.csv("exam6.3.csv",header=T)#将exam6.3.scv数据读入到d6.3中glm.ln<-glm(y~x1+x2+x3,family=poisson(link=log),data=d6.3)#建立y关于x1,x2,x3的泊松对数线性模型summary(glm.ln)#模型汇总,给出模型回归系数的估计和显著性检验等泊松分布的默认连接函数是对数连接函数，因此代码中的(link=log)可以省略.2024/7/31主编：费宇75运行以上程序可得如下结果：Call:glm(formula=y~x1+x2+x3,family=poisson(link=log),data=data6.3)DevianceResiduals:Min1QMedian3QMax-6.0569-2.0433-0.93970.792911.0061Coefficients:EstimateStd.ErrorzvaluePr(>|z|)(Intercept)1.94882590.135619114.370<2e-16***x10.02265170.000509344.476<2e-16***x20.02274010.00402405.6511.59e-08***x3-0.15270090.0478051-3.1940.0014**Signif.codes:0‘***’0.001‘**’0.01‘*’0.05‘.’0.1‘’1(Dispersionparameterforpoissonfamilytakentobe1)Nulldeviance:2122.73on58degreesoffreedomResidualdeviance:559.44on55degreesoffreedomAIC:850.71NumberofFisherScoringiterations:52024/7/31主编：费宇762.模型的参数估计和检验于是得回归模型:从检验结果可以看出:x1和x2的系数都显著,说明基础发病次数(x1),年龄(x2)和治疗条件(x3)对八周内癫痫发病数(y)重要影响.年龄(x2)的回归系数为0.0227,表明保持其他预测变量不变,年龄增加1岁,癫痫发病数的对数均值将相应的增加0.0227.2024/7/31主编：费宇772.模型的参数估计和检验在因变量的初始尺度(癫痫发病数,而不是癫痫发病数的对数)上解释回归系数比较容易,因此,指数化系数:可以看出:保持其他预测变量不变,年龄增加1岁,癫痫发病数将乘以1.023;治疗条件变化一个单位(即从安慰剂到药物),癫痫发病数将乘以0.86,换言之,保持基础癫痫发病数和年龄不变,服药相