读书笔记_第十三章

#广义线性模型

#广义线性模型扩展了线性模型的框架，它包含了非正态因变量的分析

#Logistic回归
#当通过一系列连续型、类别性预测变量来预测二值型结果变量时，logistic回归是一个非常有用的工具

#范例一
install.packages("AER")
data(Affairs,package="AER")  #进行包中指定的数据加载,Affairs变成全局变量
summary(Affairs)
#对数据进行5数法概括

#数值列：min,25百分位数，50百分位数，平均数，75百分位数，最大值
#affairs       婚外私通的次数
#age           年龄
#yearsmarried  婚龄
#religiousness 宗教信仰程度
#education     学历
#occupation    职业
#rating        对婚姻的自我评分

#类别列：各个类别值下的频数统计
#gender        参与者的性别
#children      是否有小孩

#基于因变量affairs，分组统计频数
mytable <- table(Affairs$affairs)
prop.table(mytable)  #基于频数表，生成百分比表格

#将因变量affairs,基于是否有过婚外情，改造成2值类型,未有过的为0
Affairs$ynaffairs[Affairs$affairs >0 ] <- 1  #根据条件，创造新列及值,有过的为1
Affairs$ynaffairs[Affairs$affairs == 0] <- 0 #根据条件，创造新列及值,未有过的为0

class(Affairs$ynaffairs) #新列类型numeric
Affairs$ynaffairs <- factor(Affairs$ynaffairs,levels=c(0,1),labels=c("No","Yes"))
class(Affairs$ynaffairs) #新列类型factor

table(Affairs$ynaffairs) #基于新列，生成频数表

#glm(formula,family=family(link=function),data=)
#分布函数             默认的连接函数      
#binomial             link="logit"
fit.full <- glm(ynaffairs~gender+age+yearsmarried+children+
                    religiousness+education+occupation+rating,
                    data=Affairs,family=binomial())

summary(fit.full)
#从回归系数的p值可以看到，gendermale性别,
#childrenyes是否有孩子,education学历,occupation职业
#对方差的贡献都不显著，落入拒绝域，接受原假设：回归系数为0
#去除这些变量，重新拟合模型

fit.reduced <- glm(ynaffairs~age+yearsmarried+religiousness+
                       rating,data=Affairs,family=binomial())
summary(fit.reduced)
#新模型的每个回归系数都非常显著

#anova,比较广义线性回归
anova(fit.reduced,fit.full,test="Chisq")
#Pr(>Chi)=0.2108,未落入拒绝域，表明减少后的模型和之前的完整模型拟合程度一样好

#解释模型参数
#回归系数的含义是当其他预测变量不变时，一单位预测变量的变化可引起的响应变量对数优势比的变化
coef(fit.reduced)

#指数化
#保持年龄，宗教信仰，和婚姻评定不变，婚龄增加一年，婚外情的优势比将乘以1.106
#假如婚龄增加10年，优势比将乘以1.06^10，10次方，它反映的信息更为重要
#保持婚龄，宗教信仰，和婚姻评定不变，年龄增加一岁，婚外情的优势比将乘以0.965
#因为预测变量不能等于0，截距项在此处没有什么特定含义
exp(coef(fit.reduced))

#在优势比的尺度上得到系数95%的置信区间
exp(confint(fit.reduced))

#评估预测变量对结果概率的影响

#范例一

#创建虚拟数据集
#基于列向量，生成data.frame
#向量长度不一致时，参考最长的，进行重复值填充
testdata <- data.frame(rating=c(1,2,3,4,5),age=mean(Affairs$age),
                       yearsmarried=mean(Affairs$yearsmarried),
                       religiousness=mean(Affairs$religiousness)
                       )

#预测
#predict() 用拟合模型对新数据集进行预测
testdata$prob <- predict(fit.reduced,newdata=testdata,type="response")
testdata
#结论：
#假定年龄，婚龄，宗教信仰不变，对婚姻的评分从1~5，婚外情的概率从0.53降低为0.15

#范例二

#创建虚拟数据集
#基于列向量，生成data.frame
#向量长度不一致时，参考最长的，进行重复值填充
testdata <- data.frame(rating=mean(Affairs$rating),age=seq(17,57,10),
                       yearsmarried=mean(Affairs$yearsmarried),
                       religiousness=mean(Affairs$religiousness)
)

#预测
#predict() 用拟合模型对新数据集进行预测
testdata$prob <- predict(fit.reduced,newdata=testdata,type="response")
testdata
#结论：
#假定对婚姻的评分，婚龄，宗教信仰不变，对年龄从17~57，婚外情的概率从0.33降低为0.10

#过度离势
#所谓过度离势，即观测到响应变量的方差大于期望的二项分布的方差
#过度离势会导致奇异的标准误检验和不精确的显著性检验

#范例一

#计算二项分布模型的残差偏差与残差自由度的比值，然后与1比较，大1很多，则认为存在过度离差
deviance(fit.reduced)/df.residual(fit.reduced)
#deviance,返回拟合模型对象的偏差,等价于fit.reduced对象的Residual Deviance: 615.4
#df.residual，返回模型对象的残差自由度，等价于fit.reduced对象的Degrees of Freedom: 600 Total (i.e. Null);  596 Residual
#返回1.03248，不比1大很多，表明没有过度离差

fit <- glm(ynaffairs ~ age + yearsmarried + religiousness +rating,
           family=binomial(),data=Affairs) #第一次拟合，二项分布

fit.od <- glm(ynaffairs ~ age + yearsmarried + religiousness + rating,
              family=quasibinomial(),data=Affairs) #第二次拟合，类二项分布

#pchisq函数有很多不懂
#pchisq(q, df, ncp = 0, lower.tail = TRUE, log.p = FALSE)
pchisq(summary(fit.od)$dispersion * fit$df.residual,
       fit$df.residual,lower=F)
#结论：
#p=0.340122,未落入拒绝域，接受原假设：不存在过度离势，比值=1

#泊松回归
#当通过一系列连续型或类别性预测变量来预测计数型结果变量时，采用泊松回归

#范例一
install.packages("robust")
data(breslow.dat,package="robust") #加载数据到全局环境变量
names(breslow.dat) #返回列名，也对应ID序号

#对选定列，进行summary五数汇总
#响应变量sumY,随机化后八周内癫痫发病数
#预测变量：
#Trt,   治疗条件
#Age，  年龄
#Base,  前八周内基础癫痫发病数
summary(breslow.dat[c(6,7,8,10)]) 

opar <- par(no.readonly = TRUE)
par(mfrow=c(1,2)) #图片布局，一行2列
attach(breslow.dat)
hist(sumY,breaks=20,xlab="Seizure Count",
     main="Distribution of Seizures") #画sumY，连续型变量的直方图
boxplot(sumY~Trt,xlab="Treatment",main="Group Comparisons") #画箱线图
detach(breslow.dat)

fit <- glm(sumY~Base+Age+Trt,data=breslow.dat,family=poisson())
summary(fit)
#结果：
#Estimate      回归系数
#Std. Error    标准误
#z value       偏差
#Pr(>|z|)      参数为0的检验

#解释模型参数
#(Intercept) =1.94882593，当其他预测变量都为0时，癫痫发病数的对数均值，但不可能为0，所以截距项没有意义
#Base=0.02265174
#Age=0.02274013，表示保持其他预测变量不变，年龄增加一岁，癫痫发病数的对数均值则相应增加0.03
#Trtprogabide=-0.15270095，
coef(fit)

#通常在因变量的初始尺度（癫痫发病数，而非发病数的对数）上解释回归系数比较容易
#Age=1.0230007，保持其他变量不变，年龄增加一岁，期望的癫痫发病数将乘以1.023
#这意味着年龄的增加与较高的癫痫发病数相关联

#Trtprogabide=0.8583864，保持其他变量不变，1单位Trt的变化（即从安慰剂到治疗组），
#期望的癫痫发病数将乘以0.86，也就是说，保持基础癫痫发病数和年龄不变，
#服药组相对于安慰剂组癫痫发病数将降低了20%

#与logistic回归中的指数化参数相似，泊松模型中的指数化参数对响应变量的影响都是成倍增加的，
#而不是线性相加

#过度离势

#泊松分布的方差和均值相等
#当响应变量观测的方差比依据泊松分布预测的方差大时，泊松回归可能发生过度离势

#判断是否过度离势
#如果存在过度离势，在模型中你无法进行解释，那么可能会得到很小的标准误和置信区间
#也就是说，你将会发现并不真实存在的效应
#与logistic回归类似，此处如果残差偏差与残差自由度的比例远远大于1，那么表明存在过度离势

deviance(fit)/df.residual(fit) #返回 10.1717，远大于1，存在过度离势

#可能产生过度离势的原因
#遗漏了某个重要的预测变量
#可能因为事件相关
#在纵向数据分析中，重复测量的数据由于内在群聚特性可导致过度离势

#范例一
install.packages("qcc")
library("qcc")
qcc.overdispersion.test(breslow.dat$sumY,type="poisson")

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