R软件学习笔记-8（回归分析）

CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)

pairs(~MPG+weight+displacement+horsepower,data=CarData)                    ##生成一个配对的散点图矩阵，矩阵由X中的每列的列变量对其他各列列变量的散点图组成，得到的矩阵中每个散点图行、列长度都是固定的。

##########建立 线性回归经验方程
Fit<-lm(MPG~weight+displacement+horsepower,data=CarData)
#############浏览回归分析结果
coefficients(Fit)
summary(Fit)
confint(Fit)


##########重新建立二元回归经验方程
Fit<-lm(MPG~weight+horsepower,data=CarData)
summary(Fit)
FitMPG<-predict(Fit,CarData,type="response")    ##预测
plot(CarData$weight,CarData$MPG,pch=1,xlab="自重",ylab="MPG")
points(CarData$weight,FitMPG,pch=10,col=2)
legend("topright",c("实际值","拟合值"),pch=c(1,10),col=c(1,2))

############拟合值与残差值
fitted(Fit)
residuals(Fit)
#######绘制残差图
par(mfrow=c(2,2))
plot(Fit)

##########重新建立对数-水平的二元回归经验方程
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
summary(Fit)
par(mfrow=c(2,2))
plot(Fit)
##########残差满足高斯-马尔科夫假定的其他诊断方法
library("car")
durbinWatsonTest(Fit)    ##独立性检验
summary(powerTransform(CarData$MPG))   #非正态性处理
spreadLevelPlot(Fit)   #等方差性假定图形
ncvTest(Fit)   #等方差性假定检验
crPlots(Fit)   ####被解释变量和解释变量的线性相关性判断

###########定义绘制杠杆值的用户自定义函数
LeveragePlot<-function(fit){
 Np<-length(coefficients(fit))-1
 N<-length(fitted(fit))
 plot(hatvalues(fit),main="观测点的杠杆值序列图",ylab="杠杆值",xlab="观测编号")
 abline(2*(Np+1)/N,0,col="red",lty=2)
 abline(3*(Np+1)/N,0,col="red",lty=2)
 identify(1:N,hatvalues(fit),names(hatvalues(fit)))
 }
LeveragePlot(Fit)

###############探测离群点
library("car")
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
rstudent(Fit)
outlierTest(Fit)
Fit<-lm(log(MPG)~weight+horsepower,data=CarData[-388,])
outlierTest(Fit)

#######################探测强影响点
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
par(mfrow=c(2,1))
plot(cooks.distance(Fit),main="Cook's distance",cex=0.5)      #获得Cook距离
Np<-length(coefficients(Fit))-1
N<-length(fitted(Fit))
CutLevel<-4/(N-Np-1)
plot(Fit,which=4)
abline(CutLevel,0,lty=2,col="red")

library("car")
avPlots(Fit,ask=FALSE,onepage=TRUE,id.method="identify")
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
Fit0<-lm(log(MPG)~weight+horsepower,data=CarData[-117,])

##############异常观测点的综合展示
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
influencePlot(Fit,id.method="identify",main="异常观测点的可视化")


##########多重共线性检验
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
library("car")
vif(Fit)

###########建模策略：拟合优度
Fit1<-lm(log(MPG)~weight+horsepower,data=CarData)
Fit2<-lm(log(MPG)~weight+horsepower+displacement,data=CarData)
summary(Fit1)
summary(Fit2)
anova(Fit1,Fit2)

###########建模策略：AIC
Fit1<-lm(log(MPG)~weight+horsepower,data=CarData)
Fit2<-lm(log(MPG)~weight+horsepower+displacement,data=CarData)
AIC(Fit1,Fit2)
BIC(Fit1,Fit2)

#######不同解释变量筛选策略下的AIC
Fit<-lm(log(MPG)~weight+horsepower+displacement+cylinders+acceleration,data=CarData)
library("MASS")
stepAIC(Fit,direction="backward")


###############全子集回归
library("leaps")
leapsFit<-regsubsets(log(MPG)~weight+horsepower+displacement+cylinders+acceleration,data=CarData,nbest=2)
summary(leapsFit)
coef(leapsFit, c(3,5))

####全子集回归的可视化评价
plot(leapsFit,scale="bic")
plot(leapsFit,scale="adjr2")

library("car")
subsets(leapsFit,statistic="cp",main="全子集回归模型评价图")
abline(1,1,lty=2,col="red")

############N折交叉验证
install.packages("bootstrap")
library("bootstrap")

MyNcross<-function(fit,k){
X<-as.matrix(fit$model[,2:ncol(fit$model)])
Y<-fit$model[,1]
NcrossR<-crossval(x=X,y=Y,theta.fit=lsfit(X,Y),theta.predict=cbind(1,X)%*%lsfit(X,Y)$coef,ngroup=k)
}
Fit<-lm(log(MPG)~weight+horsepower,data=CarData)
set.seed(12345)
Result<-MyNcross(Fit,10)
(R2<-cor(Fit$model[,1],Fit$fitted.values)^2)
(R2CV<-cor(Fit$model[,1],Result$cv.fit)^2)

###############自举法验证
BootLm<-function(DataSet,indices,formula){
 ReSample<-DataSet[indices,]
 fit<-lm(formula,data=ReSample)
 R2<-cor(fit$model[,1],fit$fitted.values)^2
 return(R2)
}
library("boot")
set.seed(12345)
CarData<-read.table(file="CarData.txt",header=TRUE)
BootObject<-boot(data=CarData,statistic=BootLm,R=100,formula=log(MPG)~weight+horsepower)
BootObject$t0
print(BootObject)
mean(BootObject$t)
plot(BootObject)
boot.ci(BootObject,conf=0.95,type=c("perc"),index=1)

###################回归方程的置换检验
library("lmPerm")
set.seed(12345)
Fit<-lmp(log(MPG)~weight+horsepower,data=CarData)
summary(Fit)

#####################带虚拟变量的线性回归
CarData<-read.table(file="CarData.txt",header=TRUE)
CarData$ModelYear<-as.factor(CarData$ModelYear)
Fit<-lm(log(MPG)~weight+horsepower+ModelYear,data=CarData)
summary(Fit)