R语言广义线性模型

Posted 2021-04-18 一只数据狗的成长之路

#广义线性模型和glm（）函数

#Logistic回归

library(AER)

data(Affairs,package = "AER")

summary(Affairs)

table(Affairs$affairs)#偷腥次数和频率

#展示二值型变量

Affairs$ynaffair[Affairs$affairs > 0] <- 1

Affairs$ynaffair[Affairs$affairs == 0] <- 0

Affairs$ynaffair <- 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)

#剔除性别，是否有孩子，学历和职业

fit.reduce <- glm(ynaffair ~ age + yearsmarried + religiousness

                + rating, data = Affairs,family=binomial())

summary(fit.reduce)

#新模型的每个回归系数都非常显著，用anova函数对两模型嵌套比较，可用卡方检验

anova(fit.reduce,fit.full,test = "Chisq")

#卡方p=0.2108表示不显著，表明四个预测变量的拟合程度和九个变量的拟合程度一样好，可以用更简单的

#解释模型参数

coef(fit.reduce)#对数模型解释

exp(coef(fit.reduce))#指数模型解释

exp(confint(fit.reduce))#获得系数的置信区间

testdata <- data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),

                       yearsmarried=mean(Affairs$yearsmarried),

                       religiousness=mean(Affairs$religiousness))#评价预测变量对结果概率的影响

testdata

testdata$prob <- predict(fit.reduce,newdata = testdata,type = "response")

testdata#使用测试数据集测试预测相应的概率

#假定年龄婚龄和宗教信仰不变，再看看年龄的影响

testdata <- data.frame(rating=mean(Affairs$rating),

                       age=seq(17,57,10),

                       yearsmarried=mean(Affairs$yearsmarried),

                       religiousness=mean(Affairs$religiousness))

testdata#年龄的影响

testdata$prob <- predict(fit.reduce,newdata = testdata,type = "response")

testdata#当年龄从17增加到57时婚外情概率从0.34降到了0.109

#泊松回归

library(robust)

data(breslow.dat,package = "robust")

names(breslow.dat)#数据集统计汇总信息

summary(breslow.dat[c(6,7,8,10)])

#画图

opar <- par(no.readonly = T)

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 comparisons")

par(opar)

#拟合泊松回归

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

summary(fit)

#解释模型参数

coef(fit)

#指数化系数

exp(coef(fit))