R之回歸分析廣義線性模型(Generalized Linear Model)glm

1. 介紹

廣義線性模型(Generalized  Linear Model)是一般線性模型的推廣，它使因變量的總體均值通過一個非線性連接函數而依賴于線性預測值，允許響應概率分布為指數分布族中的任何一員。許多廣泛應用的統計模型都屬于廣義線性模型，如常用于研究二元分類響應變量的Logistic回歸、Poisson回歸和負二項回歸模型等。一個廣義線性模型包含以下三個部分:

①隨機成分。

②線性成分。

③連接函數g。

各種常見的指數型分布及其主要參數

典型的連接函數及對應分布

廣義線性模型的參數估計一般不能用最小二乘估計，常用加權最小二乘法或最大似然法估計，各回歸系數β需用迭代方法求解。

2. 實現

R提供了擬合廣義線性模型的函數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, ...)

其中，

formula為擬合公式，與函數lm()中的參數formula用法相同;

family用于指定分布族，包括正態分布(gaussian)、二項分布(binomial)、泊松分布(poisson)和偽伽馬分布(Gamma)；

分布族還可以通過選項link來指定連接函數，默認值為family=gaussian (link=identity),二項分布默認值為family=binomial(link=logit);

data指定數據集;

offset指定線性函數的常數部分，通常反映已知信息;

control用于對待估參數的范圍進行設置。

例：

車險保單索賠次數分組數據

已知索賠次數服從泊松分布，相應的連接函數常用對數連接函數，模型可以寫為

下面用R實現，首先建立數據集，分類變量直接輸入定性的取值即可，glm()分析時會自動轉換成矩陣X,注意參數family的寫法。

> dat=data.frame(

y=c(42, 37, 10, 101, 73, 14),

n=c(500, 1200, 100, 400, 500, 300),

type=rep(c('小','中','大'),2),

gender=rep(c('男','女'),each=3)

)

> dat$logn=log(dat$n)  #風險暴露數取對數

#offset風險單位數事先已知

> dat.glm=glm(y~type+gender,offset=logn,data=dat,family=poisson(link=log))

> summary(dat.glm)  #glm的輸出結果

估計的回歸系數都是非常顯著的;Null deviance可以認為是模型的殘差，它的值越小說明模型擬合效果越好;模型的AIC統計量為61.68，它和deviance一起可以用來作為判斷標準，選取合適的分布族和鏈接函數。

下面通過作圖來觀察模型擬合的效果，首先提取模型的預測值，注意函數predict()提取的是線性部分的擬合值，在對數連接函數下，要得到Y的擬合值，應當再做一次指數變換。以實際觀測值為橫坐標，模型擬合值為縱坐標作圖，散點越接近直線y=x，說明模型的擬合效果越好。

> dat.pre=predict(dat.glm)

> layout(1)  #取消繪圖區域分割

> plot(y,exp(dat.pre),xlab='觀測值',ylab='擬合值',main="索賠次數的擬合效果",pch="*")

> abline(0,1)  #添加直線y=x，截距為0，斜率為1

若假設上例中的索賠次數服從負二項分布，在R中應輸入指令:

> library(MASS)

> attach(dat)

> dat.glmnb=glm.nb(y~type+gender+offset(logn))  #負二項回歸

> summary(dat.glmnb)  #輸出結果

負二項回歸擬合的模型AIC為60.45，殘差Null deviance為16.6831，小于泊松回歸擬合的殘差值，說明負二項分布的廣義線性模型更加穩定，但從回歸系數的顯著性上看，泊松回歸擬合的變量系數更加顯著。