R语言-广义线性模型

Posted 2020-10-25 月上贺兰

使用场景:结果变量是类别型,二值变量和多分类变量,不满足正态分布

　　　　  结果变量是计数型,并且他们的均值和方差都是相关的

解决方法:使用广义线性模型,它包含费正太因变量的分析

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,所以并没有充分的证据表明药物治疗相对于使用安慰剂能够降低癫痫的发病次数