廣義線性模型擴展了線性模型的框架，它包含了非正態的因變量分析
廣義線性模型擬合形式：
$$g(\mu_\lambda) = \beta_0 + \sum_{j=1}^m\beta_jX_j$$
$g(\mu_\lambda)為連接函數$. 假設響應變量服從指數分布族中某個分布(不僅僅是正態分布)，極大擴展了標准線性模型，模型參數估計的推導依據是極大似然估計，而非最小二乘法.
可以放松Y為正態分布的假設，改為Y服從指數分布族中的一種分布即可
glm()函數：glm(formula,family=family(link=function), data = )
Logistic regression：響應變量為二值(0,1), 模型假設Y服從二項分布，線性模型擬合形式：
可用以下代碼擬合Logistic回歸：glm(Y~X1+X2+X3, family = binomial(link='logit',data=mydata)
#通過婚外情數據來預測婚外情情況
#每個參與者身上有9個變量：性別，年齡，婚齡，是否有小孩，宗教信仰程度，
#學歷，職業，婚姻自我評分
library(AER)
data(Affairs,package = 'AER')
#查看描述性統計信息
summary(Affairs)
#將affairs轉化Wie二值因子ynaaffair
Affairs$ynaffair[Affairs$affairs > 0] <-1
Affairs$ynaffair[Affairs$affairs == 0] <-0
Affairs$ynaaffair <- factor(Affairs$ynaffair,levels=c(0,1),labels=c("NO","Yes"))
table(Affairs$ynaffair)
#因子化之后的值可以作為Logistic回歸的結果變量
fit.full <- glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family=binomial())
#描述模型
summary(fit.full)
#由P值得出性別，孩子，學歷，職業對方程的貢獻不顯著，去除這些變量重新擬合
fit.reduced <- glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())
summary(fit.reduced)
#由結果可以看出這次個每個回歸系數都很顯著
#由於兩個模型嵌套，可以使用anova()對他們進行比較
#卡方值p=0.21，表明四個預測變量的新模型與九個預測變量的模擬擬合程度一樣好
anova(fit.reduced,fit.full,test='Chisq')
#解釋模型系數：對數優勢
coef(fit.reduced)
#指數優勢
exp(coef(fit.reduced))
#評價預測變量對結果概率的影響
#婚姻評分對婚外情概率的影響
#創建虛擬數據集，年齡，婚齡，宗教信仰均為均值，婚姻評分為1-5
testdata <- data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),yearsmarried=mean(Affairs$yearsmarried),religiousness=mean(Affairs$religiousness))
testdata$prob = predict(fit.reduced,newdata = testdata,type="response")
testdata
Logistic 回歸變種：
- 穩健Logistic regression: robust包中glmRob()函數可以擬合文件的廣義線性模型，當擬合回歸模型出現離群點的強影響點時，穩健logistic regression便可以派生用場.
- 多項布回歸，若響應變量包含兩個以上無序類比(已婚，寡居，離婚)，便可以使用mlogit包中的mlogit()函數擬合多項Logistic回歸
- 序數Logistic回歸，若響應變量是一組有序類別(好，中，差)，便可以使用rms()包中的mlogit()函數擬合多項logistic回歸
Poisson regression：響應變量為計數型, 模型假設Y服從泊松分布, 線性模型擬合形式：
許多分析標准線性模型lm()連用的函數在glm()中都有對應的形式：
#使用robust包中的癲癇數據Breslow，來討論癲癇數據對癲癇發病率的影響 library(robust) data(breslow.dat,package = "robust") names(breslow.dat) #我們只關注，Trt治療條件，Age:年齡，基礎癲癇發病數Base, 響應變量為sumY隨機化后八周內癲癇發病數 summary(breslow.dat[c(6,7,8,10)]) opar <- par(no.readonly = TRUE) par(mfrow=c(1,2)) attach(breslow.dat) #圖中可以看到因變量的偏倚特性和可能的離群點 hist(sumY,breaks = 20,xlab = "Seizure count",main="Distribution of Seizures") boxplot(sumY~Trt,xlab="Treatment",main="Group Comarisons") par(opar) #擬合泊松回歸 fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,family=poisson()) summary(fit) #獲取模型系數 coef(fit) exp(coef(fit))
泊松回歸的的變種：
- 時間段變化泊松回歸
- 零膨脹泊松回歸：logisitic regression + poisson regression
模型擬合和回歸診斷：
1. 初始響應變量的預測值和殘差的圖形
2. 檢驗模型是否過度離勢