R语言实战-第八章 R in action-chapter8

#第八章回归

#简单线性回归
#用到基础包中的women数据集，研究身高与体重的关系
head(women)
fit <- lm(weight~height,data=women)
summary(fit)
fitted(fit)#列出拟合模型的预测值
residuals(fit)#列出拟合模型的残差值
plot(women$height,women$weight)
abline(fit)
#或者lines(women$height,fitted(fit)) 两种方法的差异在右上角
#拟合后方差解释率为99.1%
print(fit)

#多项式回归观察到women散点图有曲线的趋势，尝试用多项式回归看回归结果怎样
head(women)
fit1 <- lm(weight~height+I(height^2),data=women)
summary(fit1)
fitted(fit1)#列出拟合模型的预测值
residuals(fit1)#列出拟合模型的残差值
plot(women$height,women$weight)
lines(women$height,fitted(fit1))#注意这里和简单线性回归绘制的方式不一样
#拟合后方差解释率为99.9%

#多元线性回归
#以state.x77数据为例研究人口文盲率等对谋杀率的贡献作用(未考虑预测变量的交互项)
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
#先检查变量之间的相关性
cor(states)
library(car)
car::scatterplotMatrix(states,spread= FALSE,smoother.args = list(lty=2),
main = "Scatterplot matrix")
#非对角线绘制变量间的散点图，并添加平滑和线性拟合曲线
#对角线区域绘制每个变量的密度图和轴须图

#使用lm()进行多元线性回归模型
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
summary(fit)
#forest/income不与murder呈线性关系，总体看来，
#所有预测变量解释了各州谋杀率57%的方差

#有交互项的多元线性回归
#以mtcars数据集为例，探讨马力hp和车重wt对每英里耗油量mpg的影响
fit <- lm(mpg~hp+wt+hp:wt,data = mtcars)
summary(fit)
#交互项也很显著，表明响应变量与其中一个预测变量的关系依赖于另一个预测变量的水平
#effects::effect()展示交互项的结果
#plot(effect(term,mod,,xlevels),multiline=TRUE)
library(effects)
plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=TRUE)

#回归诊断
#基本的方法 confint
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
summary(fit)
confint(fit)
#标准方法 plot函数检验正态性独立性线性和同方差性
fit <- lm(weight~height,data = women)
par(mfrow=c(2,2))
plot(fit)

fit2 <- lm(weight~height+I(height^2),data = women)
par(mfrow=c(2,2))
plot(fit2)

#删除点后(删除数据需谨慎，本应该是模型去匹配数据，而不是反过来)
newfit2 <- lm(weight~height+I(height^2),data = women[-c(13,15),])
par(mfrow=c(2,2))
plot(newfit2)

states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit2 <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
par(mfrow=c(2,2))
plot(fit2)

#改进的方法
#1 正态性：car包qqPlot函数
##car包已经更新原来的id.method="identify"已经不可用了，变成了id=list()
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
#用下面的新代码
car::qqPlot(fit,main="Q-Q Plot",id=list(method="identify",
labels=row.names(states)),simulate=TRUE,pch=16)

#这里出现的“警告: 没有点在0.25英尺内”是Rstudio的问题，
#放入R本来的环境运行不存在该问题。建议进入原环境运行

#car::qqPlot(fit,labels=row.names(states),id.method="identify",
# simulate=TRUE,main="Q-Q plot")无法交互作用

#Nevada在置信区间外，关注一下它：
#观察实际谋杀率和模拟的谋杀率差别
states["Nevada",]
fitted(fit)["Nevada"]
residuals(fit)["Nevada"]
rstudent(fit)["Nevada"]

#利用residplot可视化误差，生成学生化残差柱状图
residplot <- function(fit,nbreaks=10){
z <- rstudent(fit)
hist(z,breaks=nbreaks,freq=FALSE,
xlab="studentized residual",
main="distribution of errors")
rug(jitter(z),col = "brown")
curve(dnorm(x,mean = mean(z),sd=sd(z)),
add = TRUE,col="blue",lwd=2)
lines(density(z)$x,density(z)$y,
col="red",lwd=2,lty=2)
legend("topright",
legend = c("normal curve","kernel density curve"),
lty=1:2,col = c("blue","red"),cex=0.7)
}
residplot(fit)

#2 误差的独立性最好的方法是依据数据收集方式的先验知识
#car::durbinWatsonTest也可以检测相关性
library(car)
car::durbinWatsonTest(fit) #如果加上simulate=TRUE则运行结果与没有时有些许不同
#p值不显著说明无自相关性

#3 线性通过成分残差图观察
car::crPlots(fit)
#若图像存在非线性则说明预测变量的函数形式可能建模不够充分，可能需要添加一些曲线成分
#比如多项式或对一个或多个变量进行变换

#4 同方差性 car包中的ncvTest和spreadLevelPlot函数
car::ncvTest(fit)
#p值不显著说明满足同方差性
car::spreadLevelPlot(fit)
#会显示建议的变换，此处异方差性很不明显因此建议幂次接近1，不需要进行变换；
#注意建议幂变换为0则使用对数变换

#线性模型假设的综合验证
#gvlma::gvlma会给出综合验证，并进行偏斜度峰度异方差性的评价
library(gvlma)
gvmodel <- gvlma::gvlma(fit)
summary(gvmodel)

#多重共线性
#car::vif 一般原则下根号下vif>2表明存在多重共线性问题
car::vif(fit)
sqrt(car::vif(fit))>2
#均为FALSE 说明不存在多重共线性问题

#异常值观测
#1 离群点
#前面用Q-Qplot判断
#也可以用标准化残差值绝对值>2认为是离群点粗略判断
#更好的方法 car::outlierTest
car::outlierTest(fit)
#Nevada的Bonferroni p被认为是离群点删除后还需要在检验是否有其他离群点存在

#高杠杆值点 :由许多异常的预测变量值组合起来的与响应变量值没有多大关系
#使用帽子统计量判断帽子均值为p/n 一般认为观测点的帽子值大于帽子均值的2或3倍就可以认为是高杠杆值点
hat.plot <- function(fit){
p <- length(coefficients(fit)) #模型估计的参数数目(包括截距)
n <- length(fitted(fit)) #样本量
plot(hatvalues(fit),main = "Index plot of hat values")
abline(h=c(2,3)*p/n,col="red",lty=2)
identify(1:n,hatvalues(fit),names(hatvalues(fit)))
}
hat.plot(fit)
#高杠杆值点可能是强影响点也可能不是这要看他们是否是离群点
#警告: 没有点在0.25英尺内是Rstudio的问题，放入R本来的环境运行不存在该问题
#点击finish结束吧

#强影响点
#Cook's D值大于4/(n-k-1)表示有强影响点
#可以通过变量添加图判断

cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)#length(fit$coefficients)包括了截距所以这里-2
plot(fit,which = 4,cook,levels=cutoff)
abline(h=cutoff,lty=2,col="red")

car::avPlots(fit,ask=FALSE,id.method="identify")

#利用car::influencePlot将离群点杠杆值和强影响点信息整合到一幅图
car::influencePlot(fit,id.method="identify")
#纵坐标在±2之外可被认为是离群点；水平轴超过0.2或0.3具有高杠杆值；
#圆圈大小与影响成比例，圆圈很大的点可能是模型参数的估计造成的不成比例影响的强影响点

#违背回归假设的数据能改进的措施
#删除观测点变量变换添加或删除变量使用其他回归方法
#删除观测点持谨慎态度

#变量变换
library(car)
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
#用下面的新代码
car::qqPlot(fit,main="Q-Q Plot",id=list(method="identify",
labels=row.names(states)),simulate=TRUE,pch=16)
#1 当模型违反正态假设时通常可以对响应变量尝试某种变换 car::powerTransform
#2 当模型违反线性假设时通常可以对预测变量尝试某种变换 car::boxTidwell

#1 响应变量变换 Box-Cox正态变换 car::powerTransform
library(car)
summary(powerTransform(states$Murder))
#结果建议用y^0.6进行变换 λ=0.6接近0.5 可以尝试平方根转换
#但LR test, lambda = (1) p=0.14512 即λ=1时也无法拒绝证明本例可以不进行变换

#2 预测变量变换 car::boxTidwell
#用人口和文盲率预测谋杀率
car::boxTidwell(Murder~Population+Illiteracy,data=states)
#结果建议population^0.87和illiteracy^1.36可以改善线性关系
#但计分检验p值表明不需要进行变换
#变量变换应该谨慎因为变换有意义才好解释

#增删变量
#删除某个存在多重共线性的变量（某个变量平方根vif >2）
#或者使用岭回归——多元回归的变体专门用来处理多重共线性问题

#尝试其它回归方法
#存在离群点或强影响点:OLS回归→稳健回归模型
#违背正态性假设:使用非参数回归模型
#存在显著的非线性:尝试非线性回归模型
#违背误差独立性假设:利用专门研究误差结构的模型如时间序列模型或多层次回归模型
#转向广泛应用的广义线性模型它适用于许多OLS回归假设不成立的情况
#多重共线性问题使用岭回归——多元回归的变体

#选择"最佳"的回归模型
#模型比较基础包anova() AIC()函数
#anova()需要嵌套模型的模型才可以作比较
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit1 <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
fit2 <- lm(Murder~Population+Illiteracy,data = states)
anova(fit2,fit1)
#检验不显著表明可以删除Income+Frost

#AIC()函数不需要嵌套模型也可以比较
AIC(fit2,fit1)
#AIC值较小的模型要优先选择

#变量选择
#逐步回归模型(向前向后向前向后) 全自己回归
#逐步回归模型 MASS::stepAIC()
library(MASS)
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
MASS::stepAIC(fit,direction="backward")
#逐步回归存在争议因为不是每一个可能的模型都被评价了

#全子集回归
#leaps::regsubsets()展示结果
library(leaps)
leapss <- leaps::regsubsets(Murder~Population+Illiteracy+Income+Frost,
data = states,nbest=4)
plot(leapss,scale="adjr2")
#双预测变量模型(最上面的0.55)population和illiteracy是最佳模型

#car::subsets()展示结果
car::subsets(leapss,statistic="cp",
main="Cp plot for all subsets regression")
abline(1,1,tyl=2,col="red")
#越好的模型里截距和斜率均为1的直线越近

#注大部分情况中全子集回归都优于逐步回归，但是全子集回归较慢
#我们需要认识到拟合效果家而没有意义的模型对我们毫无帮助，我们应该借助自己对知识背景的理解才能获得最理想的模型

#深层次分析
#交叉验证评价回归方程的泛化能力 bootstrap::crossval
library(bootstrap)
#生成shrinkage函数用以评价泛化能力
shrinkage <- function(fit,k=10){
require(bootstrap)

theta.fit <- function(x,y){lsfit(x,y)}
theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef}

x <- fit$model[,2:ncol(fit$model)]
y <- fit$model[,1]

results <- crossval(x,y,theta.fit,theta.predict,ngroup=k)
r2 <- cor(y,fit$fitted.values)^2
r2cv <- cor(y,results$cv.fit)^2
cat("Original R-square =",r2,"\n")
cat(k,"Fold Cross-Validated R-square =",r2cv,"\n")
cat("Chang =",r2-r2cv,"\n")
}
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
fit <- lm(Murder~Population+Illiteracy+Income+Frost,data = states)
shrinkage(fit)
fit2 <- lm(Murder~Population+Illiteracy,data = states)
shrinkage(fit2)
#可见 fit2的模型具有更好的泛化能力，Chang值更小(因为观测值被随机分配到k个群组，
#所以每次运shrinkage函数结果有少许不同)

#相对重要性可以比较标准化回归系数或相对权重
#比较标准化回归系数 scale(注scale函数返回的是矩阵而lm函数要求数据框需要进行转换)
states <- as.data.frame(state.x77[,c("Murder","Population",
"Illiteracy","Income","Frost")])
zstates <- as.data.frame(scale(states))
zfit <- lm(Murder~Population+Illiteracy+Income+Frost,data = zstates)
coef(zfit)
# Illiteracy 标准化回归系数最大因此解释率最高

#相对权重
#生成relweights函数用以预测变量的相对权重
relweights <- function(fit,...){
R <- cor(fit$model)
nvar <- ncol(R)
rxx <- R[2:nvar,2:nvar]
rxy <- R[2:nvar,1]
svd <- eigen(rxx)
evec <- svd$vectors
ev <- svd$values
delta <- diag(sqrt(ev))
lambda <- evec %*% delta %*% t(evec)
lambdasq <- lambda^2
beta <- solve(lambda) %*%rxy
rsquare <- colSums(beta^2)
rawwgt <- lambdasq %*% beta^2
import <- (rawwgt/rsquare)*100
import <- as.data.frame(import)
row.names(import) <- names(fit$model[2:nvar])
names(import) <- "Weights"
import <- import[order(import),1,drop=FALSE]
dotchart(import$Weights,labels=row.names(import),
xlab="% of R-Square",pch=19,#决定系数
main="Relative importance of predictor variables",#预测变量的相对重要性
sub=paste("Total R-Square=",round(rsquare,digits = 3)),#总决定系数
...)
return(import)
}
relweights(fit,col="blue")

R语言实战-第八章 R in action-chapter8相关推荐

R语言实战应用-基于R语言的对应分析
一.基本概念和原理对应分析(Correspondence Analysis)是在因子分析的基础上发展起来的,对应分析是多变量统计分析中有用的分析方法.对应分析把R型和Q型因子统一起来,通过R型因子分 ...
R语言实战-第九章 R in action-chapter9
#第九章方差分析 #需要的packages:car gplots HH rrcov multicomp effects MASS mvotlier #单因素方差分析 #数据集来源multcomp包的c ...
R语言实战笔记--第八章 OLS回归分析
R语言实战笔记–第八章 OLS回归分析标签(空格分隔): R语言回归分析首先,是之前的文章,数理统计里面的简单回归分析,这里简单回顾一下: 简单回归分析的原理:最小二乘法,即使回归函数与实际值之 ...
备受推崇的《R语言实战》真的值得如此好评吗？
作者:刘洋溢 R语言中文社区专栏作者知乎ID:https://zhuanlan.zhihu.com/p/51396601 阅前思考: R语言入门必看的<R语言实战>真的是很好的入门书籍 ...
R语言实战笔记--第十五章处理缺失数据
R语言实战笔记–第十五章处理缺失数据标签(空格分隔): R语言处理缺失数据 VIM mice 缺失值(NA),是导致我们计算错误的一大来源,处理缺失数据在实际的应用中有着较为重要的作用. 基本方 ...
R 语言实战-Part 4 笔记
R 语言实战(第二版) ## part 4 高级方法 -------------第13章广义线性模型------------------ #前面分析了线性模型中的回归和方差分析,前提都是假设因变量服 ...
R语言pacman包管理R编程语言需要的包实战：使用p_load函数安装和加载多个R包、使用p_unload函数卸载多个R包、使用p_update函数更新过期的R包
R语言pacman包管理R编程语言需要的包实战:使用p_load函数安装和加载多个R包.使用p_unload函数卸载多个R包.使用p_update函数更新过期的R包目录
R语言实战（七）图形进阶
本文对应<R语言实战>第11章:中级绘图:第16章:高级图形进阶基础图形一章,侧重展示单类别型或连续型变量的分布情况:中级绘图一章,侧重展示双变量间关系(二元关系)和多变量间关系(多元关 ...
R语言必看推荐：R语言入门经典版（中文版）+R语言实战第二版（中文完整版）
R语言入门经典(中文版)R for beginners R语言经典教材第二版适合初学者作者:Emmanuel Paradis R 语言实战第二版(中文完整版) R语言实战(第2版)注重实用性,是 ...

R语言实战-第八章 R in action-chapter8

R语言实战-第八章 R in action-chapter8相关推荐

最新文章

热门文章