R语言实战-读书笔记 （第8章 回归）

***********************************

与导图结合的脚本文件：
创建脚本：文件——新建脚本程序，将以下代码复制粘贴至脚本内，选中右键运行当前或所选代码。

##回归##

#简单线性回归
head(women)
fit <- lm(weight ~ height, data=women)
summary(fit)
fitted(fit)
residuals(fit)
plot(women$height,women$weight,
xlab="Height (in inches)",
ylab="Weight (in pounds)")
abline(fit)
#多项式回归
fit2 <- lm(weight ~ height + I(height^2), data=women)
summary(fit2)
plot(women$height,women$weight,
xlab="Height (in inches)",
ylab="Weight (in lbs)")
lines(women$height,fitted(fit2))
#加强功能图形
library(car)
scatterplot(weight ~ height, data=women,
spread=TRUE, smoother.args=list(lty=2), pch=19,
main="Women Age 30-39",
xlab="Height (inches)",
ylab="Weight (lbs.)")
#spread=FALSE选项删除残差正负均方根在平滑曲线上的展开和非对称信息#
#多元线性回归
#探究一个州的犯罪率和其他因素的关系，包括人口、文盲率、平均收入和结霜天数（温度在冰点以下的平均天数）。
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
head(states)
cor(states)#相关#
library(car)
scatterplotMatrix(states, spread=FALSE, smoother.args=list(lty=2),main="Scatter Plot Matrix")
#在非对角线区域绘制变量间的散点图，并添加平滑和线性拟合曲线。对角线区域绘制每个变量的密度图和轴须图。
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,data=states)
summary(fit)
#当预测变量不止一个时，回归系数含义为：一个预测变量增加一个单位，其他预测变量保持不变时，因变量将要增加的数量。其他自变量不变，文盲率上升1%，谋杀率上升4.14%
#交互作用
head(mtcars)
fit <- lm(mpg ~ hp + wt + hp:wt, data=mtcars)
summary(fit)
#交互作用极显著，说明响应变量与其中一个预测变量的关系依赖于另一个预测变量的水平，
#即每加仑汽油行驶英里数与汽车马力的关系依车重不同而不同。
#图形展示交互项结果
library(effects)
plot(effect("hp:wt", fit,, list(wt=c(2.2,3.2,4.2))), multiline=TRUE)

#回归诊断
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)
newfit <- lm(weight~ height + I(height^2), data=women[-c(13,15),])#剔除点13和15

states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
par(mfrow=c(2,2))
plot(fit)
#正态性
library(car)
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
qqPlot(fit, labels=row.names(states), id.method="identify",
simulate=TRUE, main="Q-Q Plot")
states["Nevada",]#查看置信区间外的点
fitted(fit)["Nevada"]
residuals(fit)["Nevada"]
rstudent(fit)["Nevada"]
#绘制学生化残差图
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=.7)
}
residplot(fit)
#独立性
durbinWatsonTest(fit)
#线性
#library(car)
crPlots(fit)
#同方差性
#library(car)
ncvTest(fit)
spreadLevelPlot(fit)#异方差性很不明显，因此建议幂次接近1（不需要进行变换）。
#综合验证
install.packages("gvlma")
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)
#多重共线性
#library(car)
vif(fit)
sqrt(vif(fit)) > 2
#离群点
#library(car)
outlierTest(fit)
#高杠杆值点
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)
#强影响点
cutoff <- 4/(nrow(states)-length(fit$coefficients)-2)
plot(fit, which=4, cook.levels=cutoff)
abline(h=cutoff, lty=2, col="red")#Cook距离
avPlots(fit, ask=FALSE, id.method="identify")#变量添加图，图中的直线表示相应预测变量的实际回归系数
influencePlot(fit, id.method="identify", main="Influence Plot",
sub="Circle size is proportional to Cook's distance")
#变换
summary(powerTransform(states$Murder))#正态
#结果表明，你可以用Murder0.6来正态化变量Murder。由于0.6很接近0.5，你可以尝试用平方根变换来提高模型正态性的符合程度。
#但在本例中，λ= 1的假设也无法拒绝（p=0.145），因此没有强有力的证据表明本例需要变量变换，这与Q-Q图结果一致。
boxTidwell(Murder~Population+Illiteracy,data=states)#线性
ncvTest(fit)#同方差性
spreadLevelPlot(fit)#异方差性很不明显，因此建议幂次接近1（不需要进行变换）。

#模型比较
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)
#由于检验不显著（p=0.994），我们可以得出结论：不需要将这两个变量添加到线性模型中，可以将它们从模型中删除。
AIC(fit1,fit2)
#此处AIC值表明没有Income和Frost的模型更佳。注意，ANOVA需要嵌套模型，而AIC方法不需要。
#逐步回归
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost,
data=states)
library(MASS)
stepAIC(fit, direction="backward")
#变量选择
library(leaps)
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
leaps <-regsubsets(Murder ~ Population + Illiteracy + Income +
Frost, data=states, nbest=4)
plot(leaps, scale="adjr2")
library(car)
subsets(leaps, statistic="cp",
main="Cp Plot for All Subsets Regression")
abline(1,1,lty=2,col="red")
#k重交叉验证
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("Change =", r2-r2cv, "\n")
}
states <- as.data.frame(state.x77[, c("Murder", "Population", "Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data = states)
fit2<-lm(Murder~Population+Illiteracy,data=states)
shrinkage(fit)
shrinkage(fit2)#R平方减少得越少，预测则越精确。
#相对重要性
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
zstates <- as.data.frame(scale(states))
#scale()函数将数据标准化为均值为0、标准差为1的数据集，这样用R回归即可获得标准化的回归系数。
#（注意，scale()函数返回的是一个矩阵，而lm()函数要求一个数据框，你需要用一个中间步骤来转换一下。
zfit <- lm(Murder~Population + Income + Illiteracy + Frost, data=zstates)
coef(zfit)#Illiteracy是最重要的预测变量，而Frost是最不重要的
#相对重要性~相对权重法
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)
}
states <- as.data.frame(state.x77[,c("Murder", "Population",
"Illiteracy", "Income", "Frost")])
fit <- lm(Murder ~ Population + Illiteracy + Income + Frost, data=states)
relweights(fit, col="blue")#在这个模型中Illiteracy比Income相对更重要