使用場景:結果變量是類別型,二值變量和多分類變量,不滿足正態分布

　　　　 結果變量是計數型,並且他們的均值和方差都是相關的

解決方法:使用廣義線性模型,它包含費正太因變量的分析

1.Logistics回歸(因變量為類別型)

　　案例:匹配出發生婚外情的模型

　　1.查看數據集的統計信息

2 library(AER)
3 data(Affairs,package = ‘AER‘)
4 summary(Affairs)
5 table(Affairs$affairs)

　　結果:該數據從601位參與者收集了,婚外情次數,性別,年齡,結婚年限,是否有孩子,宗教信仰,教育背景,職業,婚姻的自我評價這9個變量

　　　　 結果變量是婚外情發生的次數72%的夫妻沒有婚外情,最多的是一年中每月都有婚外情占6%

　　2.將結果值轉換為二值型因子

1 Affairs$ynaffair[Affairs$affairs > 0] <- 1
2 Affairs$ynaffair[Affairs$affairs == 0] <- 0
3 Affairs$ynaffair <- factor(Affairs$ynaffair, 
4                            levels=c(0,1),
5                            labels=c("No"
 
,"Yes"))
6 table(Affairs$ynaffair)

　　3.將該因子作為二值型變量的結果變量

1 fit.full <- glm(ynaffair ~ gender + age + yearsmarried + children + 
2                   religiousness + education + occupation +rating,
3                 data=Affairs,family=binomial())
4 summary(fit.full)

　　結果:性別,是否有孩子,學歷和職業對模型不顯著,去除後進行分析

1 fit.reduced <- glm(ynaffair ~ age + yearsmarried + religiousness + 
2                      rating, data=Affairs, family=binomial())
3 summary(fit.reduced)

　　3.使用卡方檢驗來判斷比較

1 anova(fit.reduced,fit.full,test = ‘Chisq‘)

　　 結果:p=0.21,表示新模型的擬合更好

　　4.解釋模型參數

1 coef(fit.reduced)
2 exp(coef(fit.reduced))

　　結果:婚齡每增加1歲,婚外情發生的可能性將乘以1.106,相反年齡增加1歲,婚外情發生的可能性乘以0.9652

　　5.評價婚姻評分對婚外情的影響

1 # 1.手動生成數據集
2 # 2.使用predict函數來進行預測
3 testdata <- data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),
4                        yearsmarried=mean(Affairs$yearsmarried),
5                        religiousness=mean(Affairs$religiousness))
6 testdata
7 testdata$prob <- predict(fit.reduced,newdata = testdata,type=‘response‘)
8 testdata

　　結果:當婚姻評分從1(很不幸)變成5(很幸福)的時候,婚外情發生的概率從0.53降低到0.15

　　6.評價年齡對婚外情的影響

1 testdata <- data.frame(rating=mean(Affairs$rating),
2                        age=seq(17,57,10),
3                        yearsmarried=mean(Affairs$yearsmarried),
4                        religiousness=mean(Affairs$religiousness))
5 testdata$prob <- predict(fit.reduced,newdata = testdata,type=‘response‘)
6 testdata

　　結果:當其他變量不變時,年齡從17到57歲,婚外情的概率從0.34降低到0.11

　　7.判斷是否過度離勢

　　　　過度離勢會導致標準誤檢驗和不精確的顯著性檢驗,此時任然可以使用gml()擬合擬合Logistics回歸,但是把二項分布改為類二項分布

1 # 如果結果接近1,表示沒有過度離勢
2 deviance(fit.reduced)/df.residual(fit.reduced)

　　結果:沒有過度離勢

2.泊松回歸(因變量為計數型)

　　使用場景:通過一系列連續型或類別型預測變量來預測計數型結果變量時采用泊松分布

　　案例:藥物治療是否能減小癲癇的發病數

　　1.查看數據集

1 data(breslow.dat,package = ‘robust‘)
2 names(breslow.dat)
3 summary(breslow.dat[c(6,7,8,10)])

　　結果:我們分析年齡,治療條件,前八周的發病次數和隨機化後八周內的發病次數的關系,所以只采用4個變量

　　2.圖形

1 opar <- par(no.readonly = T)
2 par(mfrow=c(1,2))
3 attach(breslow.dat)
4 hist(sumY,breaks = 20,xlab = ‘Seizure Count‘,main = ‘Distribution of Sizeture‘)
5 boxplot(sumY~Trt,xlab=‘Treatment‘,main=‘Group Comparisons‘)
6 par(opar)

　　結果:可以看出使用藥物的組,癲癇的發病率有所減少

　　3.擬合泊松回歸

1 fit <- glm(sumY~Base+Age+Trt,data = breslow.dat,family = poisson())
2 summary(fit)

　　結果:偏差,回歸參數,標準誤差和參數為0的檢驗

　　4.解釋模型參數

1 coef(fit)
2 exp(coef(fit))

　　結果:年齡每增加1歲,癲癇的發病數將乘以1.023,如果從安慰劑組調到藥物組,則發病率會減少14%

　　5.判斷是否過度離勢

1 deviance(fit)/df.residual(fit)

　　結果:大於1,存在過度離勢

　　6.調整模型

1 fit.new <- glm(sumY~Base+Age+Trt,data = breslow.dat,family = quasipoisson())
2 summary(fit.new)

　　結果:標準誤差和第一次模型相比,大了許多,同時標準誤差越大會導致Trt的p值大於0.05,所以並沒有充分的證據表明藥物治療相對於使用安慰劑能夠降低癲癇的發病次數

R語言-廣義線性模型