《多元統(tǒng)計(jì)分析-基于R（第3版）》課件 第6章-廣義線性模型

2024/7/31主編：費(fèi)宇1主編：費(fèi)宇中國人民大學(xué)出版社2024/7/31主編：費(fèi)宇2第6章廣義線性模型6.1廣義線性模型的定義6.2Logistic模型6.3Probit模型6.4多項(xiàng)Logit模型6.5泊松對(duì)數(shù)線性模型6.6零膨脹計(jì)數(shù)模型6.7多項(xiàng)分布對(duì)數(shù)線性模型2024/7/31主編：費(fèi)宇36.1

廣義線性模型的定義第5章我們研究了多元線性模型,該模型的一個(gè)重要假定是因變量是連續(xù)型的變量(通常假定服從正態(tài)分布),但在許多情況下,這種假定并不合理,例如下面這兩種情況.(1)結(jié)果變量可能是類型變量.二值分類變量和多分類變量.(比如:是/否,差/一般/良好/優(yōu)秀等)顯然都不是連續(xù)型變量.2024/7/31主編：費(fèi)宇46.1廣義線性模型概述(2)結(jié)果變量可能是計(jì)數(shù)型變量(比如:一周交通事故的數(shù)目)這類變量都是非負(fù)的有限值,而且它們的均值和方差通常是相關(guān)的(一般線性模型假定因變量是正態(tài)變量,而且相互獨(dú)立).普通線性回歸模型(5.3)假定因變量y服從正態(tài)分布,其均值滿足關(guān)系式:μ=Xβ,這表明因變量的條件均值是自變量的線性組合.本章介紹六種常見的廣義線性模型：Logistic模型、Probit模型、多項(xiàng)Logit模型、泊松對(duì)數(shù)線性模型、零膨脹計(jì)數(shù)模型和多項(xiàng)分布對(duì)數(shù)線性模型.2024/7/31主編：費(fèi)宇56.1廣義線性模型的定義1.廣義線性模型的定義：(1)隨機(jī)成分:設(shè)y1,y2,…,yn是來自于指數(shù)分布族的隨機(jī)樣本,即yi的密度函數(shù)為其中ai(.),b(.),ci(.)是已知函數(shù),參數(shù)αi是典則參數(shù),?是散度參數(shù).2024/7/31主編：費(fèi)宇61.廣義線性模型的定義：(2)連接函數(shù):設(shè)yi的均值為μi而函數(shù)g(.)是單調(diào)可微的連接函數(shù),使得其中是協(xié)變量,是未知參數(shù)向量.指數(shù)分布族正態(tài)分布二項(xiàng)分布泊松分布2024/7/31主編：費(fèi)宇72024/7/31主編：費(fèi)宇82.正態(tài)線性回歸模型正態(tài)分布屬于指數(shù)分布族,其密度函數(shù)為與(6.1)對(duì)照可知2024/7/31主編：費(fèi)宇92.正態(tài)線性回歸模型只要取連接函數(shù)為,則正態(tài)線性回歸模型滿足廣義線性模型的定義.類似的,容易驗(yàn)證,二項(xiàng)分布和泊松分布都屬于指數(shù)分布族.下面介紹實(shí)際中應(yīng)用廣泛的Logistic模型、Probit模型、多項(xiàng)Logit模型、泊松對(duì)數(shù)線性模型、零膨脹計(jì)數(shù)模型和多項(xiàng)分布對(duì)數(shù)線性模型.2024/7/31主編：費(fèi)宇106.2Logistic模型1.模型定義 設(shè)yi服從參數(shù)為pi的二項(xiàng)分布,則μi=E(yi)=pi采用邏輯連接函數(shù),即這個(gè)廣義線性模型稱為Logistic模型.2024/7/31主編：費(fèi)宇11例6.1(數(shù)據(jù)文件為eg6.1)表6-1

某城市48個(gè)家庭的調(diào)查數(shù)據(jù)2024/7/31主編：費(fèi)宇122.模型的參數(shù)估計(jì)和檢驗(yàn)采用R軟件中的廣義線性模型過程glm()可以完成回歸系數(shù)的估計(jì),以及模型回歸系數(shù)的顯著性檢驗(yàn).程序如下：#例6.1廣義線性模型:Logistic模型setwd("C:/data")#設(shè)定工作路徑d6.1<-read.csv("exam6.1.csv",header=T)#將exam6.1.csv數(shù)據(jù)讀入到d6.1中g(shù)lm.logit<-glm(y~x1+x2,family=binomial(link=logit),data=d6.1)#建立y關(guān)于x1，x2的logistic回歸模型，數(shù)據(jù)為d6.1summary(glm.logit)#模型匯總注意邏輯連接函數(shù)是二項(xiàng)分布的典則連接函數(shù)，是默認(rèn)的連接函數(shù)，因此代碼中的(link=logit)可以省略.2024/7/31主編：費(fèi)宇13運(yùn)行以上程序可得如下結(jié)果：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主編：費(fèi)宇14注意到x2對(duì)應(yīng)的p值(0.948)比較大，即x2不顯著，所以考慮采用逐步回歸.glm.step<-step(glm.logit)#逐步回歸summary(glm.step)#給出模型回歸系數(shù)的估計(jì)和顯著性檢驗(yàn)等運(yùn)行以上程序可得如下結(jié)果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)#給出模型回歸系數(shù)的估計(jì)和顯著性檢驗(yàn)等2024/7/31主編：費(fèi)宇15Call: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模型預(yù)測如果要預(yù)測年收入為20萬元（x1=20）、家里有孩子（x2=1）的家庭有購買住房的可能性，可以采用以下命令：>yp<-predict(glm.step,data.frame(x1=20))>p.fit<-exp(yp)/(1+exp(yp));p.fit#估計(jì)x1=20時(shí)y=1的概率10.7728122容易看出，當(dāng)x1=20，x2=1時(shí)，估計(jì)y=1的概率約為0.77，即年收入為20萬元、家里有孩子的家庭有購買住房的可能性約為77%.2024/7/31主編：費(fèi)宇166.3Probit模型2024/7/31主編：費(fèi)宇176.3Probit模型:例6.1（續(xù)）2024/7/31主編：費(fèi)宇18運(yùn)行以上程序可得如下結(jié)果：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主編：費(fèi)宇192024/7/31主編：費(fèi)宇206.4多項(xiàng)Logit模型2024/7/31主編：費(fèi)宇21例6.2某城市48個(gè)家庭的調(diào)查數(shù)據(jù)2024/7/31主編：費(fèi)宇22問題：根據(jù)這個(gè)數(shù)據(jù)建立多項(xiàng)分布回歸模型并估計(jì)年收入為20萬元、家里有孩子的家庭有購買住房但還在還貸款的可能性.解:采用nnet程序包中的multinom()可以完成多項(xiàng)logit模型的擬合.#例6.2廣義線性模型:多項(xiàng)分布回歸模型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主編：費(fèi)宇23運(yùn)行以上程序可得如下結(jié)果：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對(duì)應(yīng)標(biāo)準(zhǔn)誤相對(duì)于x2的系數(shù)比較大，所以估計(jì)x2可能不顯著，采用step()函數(shù)對(duì)模型進(jìn)行逐步回歸.2024/7/31主編：費(fèi)宇24mlog.s<-step(mlog)#對(duì)mlog進(jìn)行逐步回歸summary(mlog.s)#查看所擬合的模型運(yùn)行以上程序可得如下結(jié)果：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主編：費(fèi)宇252024/7/31主編：費(fèi)宇262024/7/31主編：費(fèi)宇272024/7/31主編：費(fèi)宇28估計(jì)48個(gè)家庭最可能屬于3類家庭中的哪一類？2024/7/31主編：費(fèi)宇296.5泊松對(duì)數(shù)線性模型1.模型的定義設(shè)y服從參數(shù)λ為的泊松分布,則μ=E(y)=λ,采用對(duì)數(shù)連接函數(shù),即這個(gè)廣義線性模型稱為泊松對(duì)數(shù)線性模型.2024/7/31主編：費(fèi)宇30例6.3(數(shù)據(jù)文件為eg6.3)表6-3

Breslow癲癇數(shù)據(jù)2024/7/31主編：費(fèi)宇31例6.3(數(shù)據(jù)文件為eg6.3)這個(gè)數(shù)據(jù)是robust包中的Breslow癲癇數(shù)據(jù)(Breslow,1993).我們討論在治療初期的八周內(nèi),癲癇藥物對(duì)癲癇發(fā)病數(shù)的影響,響應(yīng)變量為八周內(nèi)癲癇發(fā)病數(shù)(y),預(yù)測變量為前八周內(nèi)的基礎(chǔ)發(fā)病次數(shù)(x1),年齡(x2)和治療條件(x3),其中治療條件是二值變量,x3=0表示服用安慰劑,x3=1表示服用藥物.根據(jù)這個(gè)數(shù)據(jù)建立泊松對(duì)數(shù)線性模型并對(duì)模型的系數(shù)進(jìn)行顯著性檢驗(yàn).2024/7/31主編：費(fèi)宇322.模型的參數(shù)估計(jì)和檢驗(yàn)采用R軟件中的廣義線性模型過程glm()來建立泊松對(duì)數(shù)線性模型并對(duì)模型的系數(shù)進(jìn)行顯著性檢驗(yàn).程序如下：#例6.3廣義線性模型:泊松對(duì)數(shù)線性模型setwd("C:/data")d6.3<-read.csv("exam6.3.csv",header=T)#將exam6.3.scv數(shù)據(jù)讀入到d6.3中g(shù)lm.ln<-glm(y~x1+x2+x3,family=poisson(link=log),data=d6.3)#建立y關(guān)于x1,x2,x3的泊松對(duì)數(shù)線性模型summary(glm.ln)#模型匯總,給出模型回歸系數(shù)的估計(jì)和顯著性檢驗(yàn)等泊松分布的默認(rèn)連接函數(shù)是對(duì)數(shù)連接函數(shù)，因此代碼中的(link=log)可以省略.2024/7/31主編：費(fèi)宇33運(yùn)行以上程序可得如下結(jié)果：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主編：費(fèi)宇342.模型的參數(shù)估計(jì)和檢驗(yàn)于是得回歸模型:從檢驗(yàn)結(jié)果可以看出:x1和x2的系數(shù)都顯著,說明基礎(chǔ)發(fā)病次數(shù)(x1),年齡(x2)和治療條件(x3)對(duì)八周內(nèi)癲癇發(fā)病數(shù)(y)重要影響.年齡(x2)的回歸系數(shù)為0.0227,表明保持其他預(yù)測變量不變,年齡增加1歲,癲癇發(fā)病數(shù)的對(duì)數(shù)均值將相應(yīng)的增加0.0227.2024/7/31主編：費(fèi)宇352.模型的參數(shù)估計(jì)和檢驗(yàn)在因變量的初始尺度(癲癇發(fā)病數(shù),而不是癲癇發(fā)病數(shù)的對(duì)數(shù))上解釋回歸系數(shù)比較容易,因此,指數(shù)化系數(shù):可以看出:保持其他預(yù)測變量不變,年齡增加1歲,癲癇發(fā)病數(shù)將乘以1.023;治療條件變化一個(gè)單位(即從安慰劑到藥物),癲癇發(fā)病數(shù)將乘以0.86,換言之,保持基礎(chǔ)癲癇發(fā)病數(shù)和年齡不變,服藥相對(duì)于安慰劑組癲癇發(fā)病數(shù)降低了14%.>exp(coef(glm.ln))(Intercept)x1x2x3

7.02044031.02291021.02300070.85838646.6零膨脹計(jì)數(shù)模型2024/7/31主編：費(fèi)宇366.6零膨脹計(jì)數(shù)模型2024/7/31主編：費(fèi)宇37例6.4(數(shù)據(jù)文件為eg6.4)2024/7/31主編：費(fèi)宇38例6.4：血友病數(shù)據(jù)解：先讀入數(shù)據(jù)并查看變量deaths死亡人數(shù)的分布：>#例6.4血友病數(shù)據(jù):先讀入數(shù)據(jù)并查看變量deaths>setwd("C:/data")>d6.4<-read.csv("exam6.4.csv",header=T)#將exam6.4.xls數(shù)據(jù)讀入到d6.4中>table(d6.4$deaths)#查看變量deaths

012345618332126228621>barplot(table(d6.4$deaths))#畫條形圖2024/7/31主編：費(fèi)宇392024/7/31主編：費(fèi)宇40如果不考慮零膨脹問題，直接采用泊松對(duì)數(shù)線性模型來擬合數(shù)據(jù)，程序如下：#例6.4血友病數(shù)據(jù):直接采用泊松對(duì)數(shù)模型setwd("C:/data")d6.4<-read.csv("exam6.4.csv",header=T)#將exam6.4.csv數(shù)據(jù)讀入到d6.4中hiv<-factor(d6.4$hiv)#將變量hiv因子化fac<-factor(d6.4$factor)#將變量factor因子化a1<-glm(deaths~hiv+fac+age+py,family=poisson(link=log),data=d6.4)#建立deaths關(guān)于hiv,fac,age,py的泊松對(duì)數(shù)線性模型summary(a1)#模型匯總,給出模型回歸系數(shù)的估計(jì)和顯著性檢驗(yàn)等2024/7/31主編：費(fèi)宇412024/7/31主編：費(fèi)宇422024/7/31