回归分析

#广义线性模型-logistic 与泊松回归 Y值类型为0/1或者一些计数型的

#logistics-AER包，婚外情数据affairs，各种因素是不是跟婚外情有关系

data(Affairs,package = "AER")

#变量x的统计分析

summary(Affairs)

#相应变量y的统计数值

table(Affairs$affairs)

#这时候我们需要把婚外情变化成二值型数据0/1  只关心有没有  ynaffair可以任意取名字，新增的一列而已

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)

##此时因为Y变为二值型变量0/1，可以用logistic回归了

fit<-glm(ynaffair~gender+age+yearsmarried+children+religiousness+education+occupation+rating,data=Affairs,family=binomial())

summary(fit)

#去除拟合不好的再拟合

fit_reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

summary(fit_reduced)

##因为fit_reduced拟合是fit拟合的子集，可以用卡方检验的anova进行比较  p>0.05，没有差别，因此可以删去无关的数据

anova(fit_reduced,fit,test = "Chisq")

##可以用coef单独看看系数

coef(fit_reduced)

##方法一：用优势比解释系数

##因为响应变量是Y=1的log值，不好解释，可以e-指数化后恢复原来的值；保持年龄什么的不变，婚龄增加一年，婚外情

    ##优势比乘于1.106.。。因此大于1，上升，小于1下降优势比  若分析几年，就乘于几年就好了*n --n年

exp(coef(fit_reduced))

##系数置信区间

confint(fit_reduced)

##指数化的结果--一般用这个

exp(confint(fit_reduced))

##方法二：用概率解释系数

##创建感兴趣的数据集

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_reduced,newdata = testdata,type = "response")

testdata

##过度离势的判断，避免不准确的显著性检验

fit_reduced<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family=binomial())

fit.od<-glm(ynaffair~age+yearsmarried+religiousness+rating,data=Affairs,family = quasibinomial())

#检验 >0.34,没有过度离势，之前的检验显著可信；

pchisq(summary(fit.od)$dispersion*fit_reduced$df.residual,fit_reduced$df.residual,lower=F)

##常见的回归诊断可以用线性回归诊断也可以用下面的

##一般方法：

op<-par(mfrow=c(2,2))

plot(fit_reduced)

par(op)

###本章的方法  预测值与残差值

plot(predict(fit_reduced,type = "response"),residuals(fit_reduced,type="deviance"))

###异常值--hat value  学生化残差，cook距离

plot(hatvalues(fit_reduced))

plot(rstudent(fit_reduced))

plot(cooks.distance(fit_reduced))

###泊松回归分析---robust包

library(robust)

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

names(breslow.dat)

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

##基本图形描述

op<-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(op)

###回归拟合：

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

summary(fit_BS)

###系数分析

coef(fit_BS)

##指数化

exp(coef(fit_BS))

##过度离势分析

##残差偏差与残差自由度比值  远远大于1  说明过度离势

deviance(fit_BS)/df.residual(fit_BS)

#或用qcc包中进行过度离势分析  p<0.05 进一步说明过度离势

library(qcc)

qcc.overdispersion.test(breslow.dat$sumY,type="poisson")

##将poisson换成quasipoisson 重新拟合

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

summary(fit)

##########################################################

##第五题  回归分析

data_5<-read.table("final5/diabetes.txt",header = T,sep = "\t")

data_5<-as.data.frame(data_5)

boxplot(data_5)

library(car)

scatterplotMatrix(data_5,spread=F,main="Scatter Plot Matrix")

##拟合

fit_5<-lm(Y~X1+X2+X3+X4,data = data_5)

summary(fit_5)

#回归诊断--标准方法

op<-par(mfrow=c(2,2))

plot(fit_5)

par(op)

par(par(mfrow=c(2,2)))

##离群点

outlierTest(fit_5)

##逐步回归分析 P大于0.05可以删去x1 x2

fit_5_new<-lm(Y~X3+X4,data = data_5)

anova(fit_5,fit_5_new)

##交互作用

fit_5_ab<-lm(Y~X3+X4+X3:X4,data = data_5)

fit_5_ab

##交互结果展示

library(effects)

plot(effect("X3:X4",fit_5_ab),multiline=T)