R语言与线性回归分析

基本假设：正态性、独立性、线性、同方差性
常用R包：car(carData)、MASS、leap…

1. 原始数据的分析

plot(x,y); abline(lm(y~x))
library(carData);library(car)
scatterplot(x,y)  #生成拟合的二元关系图（散点+箱线）
scatterplotMatrix()  #生成散点图矩阵

cor()  #计算相关系数

2. 回归模型的拟合(参数估计和检验)

# 数据准备
n <- ; p <- 
Y <- as.matrix(dataset[,1])
X1 <- as.matrix(dataset[,2:p])
X <- cbind(1,X1)
XX <- t(X)%*%X

# 参数估计
Beta <- solve(XX)%*%t(X)%*%Y
e <- Y-X%*%Beta #<or> resid <- resid(lm) #残差的计算
SSe <- t(e)%*%e
sigma2 <- SSe/(n-p) || sigma2 <- sd(e)^2*(n-1)/(n-p)  #σ标准差的估计
covBeta <- sigma2*solve(XX)  #Beta的协方差阵
r <- matrix(nrow=p,ncol=p)
for (i in 1:p)  {
	for (j in 1:p)
		r[i,j] <- covBeta[i,j]/(sqrt(covBeta[i,i])*sqrt(covBeta[j,j]))
		}   #Beta的相关系数阵

# 模型检验
Beta1 <- beta[-1]
Xc <- (X-matrix(1,nrow=18)%*%colMeans(X))[,-1]
RSS_I <- t(Beta1)%*%t(Xc)%*%Y
RSS_H <- n*mean(Y)^2
SS_He <- t(Y)%*%Y-RSS_H
RSS <- t(Beta)%*%t(X)%*%Y  ||  RSS <- RSS_H + RSS_I
SSe <- t(Y)%*%Y-RSS
F <- (RSS-RSS_H)/(p-1)/(SSe/(n-p))
pf(F,p-1,n-p,lower.tail = T)

# 系数检验
t <- matrix(1,ncol=p)
for(j in 1:4) t[j]<-beta[j]/sqrt(sigma2*solve(t(X)%*%X)[j,j])
p-value <- 2*pt(abs(t),n-p,lower.tail = FALSE)

# 复相关系数
TSS<-t(Y-mean(Y))%*%(Y-mean(Y))
R2 <- RSS_I  / TSS
Adjusted_R2 <- 1-(1-R2)*(n-1)/(n-p)

myfit <- lm(y~x1+x2+x3,dataframe)

summary() 展示拟合模型的详细结果
coefficients()  列出拟合模型的模型参数（截距项和斜率）
confint()  列出模型参数的置信区间（95%）
fitted()  列出拟合模型的预测值
anova()  生成拟合模型的方差分析表
vcov()  列出模型参数的协方差矩阵
AIC()  输出赤池信息统计量
plot()  生成评价拟合模型的诊断图
predict()  用拟合模型对新的数据集预测响应变量值

3. 变量选择

变量选择准则

## 平均残差平方和准则 (RMSq)
RMSq <- SSeq/(n-q)
## Cp准则
Cp <- SSeq/sigma2-(n-2*q)
## 赤池信息准则 (AIC/BIC)
AIC <- (-2)*lnL+2*q
BIC <- (-2)*lnL+q*log(n)

# （下面代码假设有3个自变量）
e <- list()
SSeq<-matrix(1,nrow=7); RMS<-matrix(1,nrow=7); Cp<-matrix(1,nrow=7); 
lnL<-matrix(1,nrow=7); AIC<-matrix(1,nrow=7); BIC<-matrix(1,nrow=7);
q<-matrix(c(2,2,2,3,3,3,4),nrow=1)
e[[7]] <- resid(lmlist[[7]])
for(i in 1:7)
{
  e[[i]] <- resid(lmlist[[i]])
  SSeq[i]<- t(e[[i]])%*%e[[i]]
  RMS[i]<-SSeq[i]/(n-q[i])
  Cp[i] <-SSeq[i]/(sd(e[[7]])^2*(n-1)/(n-p))-n+2*q[i]
  lnL[i]<-(-n/2)*(log(SSeq[i])+log(2*pi/n)+1)
  AIC[i]<-(-2)*lnL[i]+2*q[i]   ||  AIC(lmlist[[i]])   # 两者的结果有一点差异
  BIC[i]<-(-2)*lnL[i]+q[i]*log(n)  ||  BIC(lmlist[[i]])   # 两者的结果有一点差异
} 
COM<-data.frame(row.names=c('x1','x2','x3','x1 x2','x1 x3','x2 x3','x1 x2 x3'),RMS,Cp,AIC)
COM

逐步回归

library(MASS)
exps <- stepAIC(lm,scope= ,direction=<"both"/"forward"/"backward">...)
summary(exps)

全子集回归

library(leaps)
leaps <- regsubsets(y~.,data= ,nbest=<n>)   # nbest表示展示n个最佳的子集回归模型
plot(leaps,scale="adjr2") #绘图展示子模型和调整R平方的结果
coef(leaps,6)

sleaps <- summary(leaps)
names(sleaps)  # "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj" 
plot(reg.summary$adjr2,xlab="Number of Variables",ylab="Adjusted RSq",type="l")
i <- which.max(reg.summary$adjr2)
points(i,reg.summary$adjr2[i], col="red",cex=2,pch=20)

res <- data.frame(sleaps$outmat,调整R平方=sleaps$adjr2)
res   #表格展示子模型和调整R平方的结果

library(car)
subsets(leaps,statistic="cp",main="Cp Plot for All subsets Regression")   #基于Cp统计量的不同子模型的比较
abline(1,1,lty=2,col="red")

向前/向后回归

regfit.fwd=regsubsets(y~.,data=,nvmax=,method="forward")#向前逐步选择
summary(regfit.fwd)
regfit.bwd=regsubsets(y~.,data=,nvmax=,method="backward") #向后逐步选择
summary(regfit.bwd)
coef(regfit.fwd,7)
coef(regfit.bwd,7)

模型选择

set.seed(1)
train=sample(c(TRUE,FALSE), nrow(dataset),rep=TRUE) #rep是replace
test=(!train)
regfit.best=regsubsets(y~.,data=dataset[train,],nvmax=19)
test.mat=model.matrix(y~.,data=dataset[test,])#此函数用于生成回归设计矩阵
val.errors=rep(NA,19)
for(i in 1:19){
  coefi=coef(regfit.best,id=i)
  pred=test.mat[,names(coefi)]%*%coefi
  val.errors[i]=mean((Hitters$Salary[test]-pred)^2)
}
val.errors
which.min(val.errors)
coef(regfit.best,10)

#以上程序有点冗余，可以编写预测函数
predict.regsubsets=function(object,newdata,id,...){
  form=as.formula(object$call[[2]])
  mat=model.matrix(form,newdata)
  coefi=coef(object,id=id)
  xvars=names(coefi)
  mat[,xvars]%*%coefi
}
regfit.best=regsubsets(y~.,data=dataset,nvmax=19)
coef(regfit.best,10)

# k折交叉验证
k=10
set.seed(1)
folds=sample(1:k,nrow(dataset),replace=TRUE)
cv.errors=matrix(NA,k,19, dimnames=list(NULL, paste(1:19)))
for(j in 1:k){
 best.fit=regsubsets(y~.,data=dataset[folds!=j,],nvmax=19)
  for(i in 1:19){
    pred=predict.regsubsets(best.fit,dataset[folds==j,],id=i)#这里可以简化为predict
    cv.errors[j,i]=mean( (dataset$y[folds==j]-pred)^2)
  }
}
mean.cv.errors=apply(cv.errors,2,mean)
mean.cv.errors
par(mfrow=c(1,1))
plot(mean.cv.errors,type='b')
reg.best=regsubsets(y~.,data=dataset, nvmax=19)
coef(reg.best,11)

4. 回归诊断

综合诊断

# 残差图、QQ图、位置尺度图、cook统计量图
plot(fit)
par(mfrow=c(2,2)); plot(fit)  # 生成在一张图中

# 模型假设通过与否的综合检验
library(gvlma)
gvmodel <- gvlma(fit)
summary(gvmodel)

对误差项假设的诊断

## 残差的计算
resid() or residuals()   # 拟合模型的残差值(SSe)
resid()/sd(resid())  #标准化残差
rstudent()  # 学生化残差(ri)

## 正态性诊断
# ①正态QQ图
library(car)
qqPlot()
# ②学生化残差柱状图
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(lm)

## 误差独立性诊断
# Durbin-Watson检验 (p-value>0.05 → 无自相关性，误差项之间独立)
durbinWatsonTest(lm)

## 线性诊断
# 成分残差图（偏残差图）
library(car)
crPlots(lm)

## 同方差性诊断
ncvTest(lm) #生成计分检验，零假设为误差方差不变
spreadLevelPlot(lm) #生成最佳拟合曲线的散点图

异常观测值的影响分析

## 异常值点(离群点——模型预测效果不佳的观测点，有很大的残差值)
# ①QQ图——落在置信区间以外的点
# ②标准化残差值大于2或小于-2的点
# ③根据单个最大残差值的显著性判断模型是否存在离群点
outlierTest()  # p-value<0.05, 判定为离群点

## 强影响点（对模型参数估计值影响比例失衡的点）
# ①Cook距离（D统计量 > 4/(n-p-1)）
cutoff <- 4/(nrow(data)-length(fit$coefficients)-2)
plot(fit,which=4,cook.levels=cutoff)
abline(h=cutoff,lty=2,col="red")
# ②变量添加图（提供关于强影响点如何影响模型的信息）
library(car)
avPlots(fit,ask=FALSE)

## 高杠杆值点（与其他预测变量有关的离群点，由许多异常的预测变量值组合起来，与响应变量没有关系）
# 帽子统计量（hatvalue > 2or3*p/n）
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)
	idenfify(1:n,hatvalues(fit),names(hatvalues(fit)))
	}  #identify选项可以在图中交互式标注离群点
hat.plot(fit)

## 综合
library(car)
influencePlot(fit,main="Influence Plot",sub="Circle size is proportional to Cook's distance")
# 纵坐标超过2或-2的可能是离群点，横坐标超过0.2或0.3的可能是高杠杆值点，圆圈较大的可能是强影响点

多重共线性诊断

## 方差膨胀因子 vif=1/(1-R^2), vif越大，变量之间的线性相关程度越大
#（一般vif>10表明存在严重多重共线性，剔除该变量;或计vif的均值>1）
#当样本量不算小，而R2接近1时，可以认为存在严重的多重共线性
library(car)
vif(lm) 

## 条件数 k（lambda(1)/lambda(p-1)） 
# (<=100:共线性的程度较小 100~1000：存在较强的多重共线性   >1000：存在严重的多重共线性)
X <- scale(as.matrix(data[2:]));X
kappa(t(X)%*%X,exact=TRUE)

## 另：计算矩阵的特征值
lambda <- eigen(t(X)%*%X)
lambda$values[1]/lambda$values[p-1]

5. 改进措施

1.删除观测点
2.增删变量（删除有多重共线性的变量）
3.Box-Cox变换（变换Y）

bc<-boxcox(y~ ,lambda=seq(-2,2,0.01))

lambda<-bc$x[which.max(bc$y)] 
lambda

y_bc<-(y^lambda-1)/lambda 
lm_bc<-lm(y_bc~I(x1^2)+x2)
summary(lm_bc)

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

4.Box-Cox变换、Box-Tidwell变换（变换X）

# Box-Cox变换（适用于模型违反正态假设的情况）
library(car)
summary(powerTranform(x))
# Box-Tidwell变换（适用于模型违反线性假设的情况）
library(car)
boxTidwell(y~x1+x2)

5.尝试其他方法（岭回归、稳健回归、非参数回归、非线性回归、时间序列、多层次回归、广义线性回归）

6. 岭回归

8. 其他

有交互项的线性回归

lm<-lm(y~x1+x2+x1:x2)
library(effects)
plot(effect("x1:x2",lm,,list(x1=c()),multiline=TRUE)

深层次分析——交叉验证（刀切法）

shrinkage <- function(fit,k=10){
	  require(bootstrap)
	  x <- fit$model[,2:ncol(fit$model)]
	  y <- fit$model[,1]
	  theta.fit <- function(x,y){lsfit(x,y)}
	  theta.predict <- function(fit,x){as.matrix(cbind(1,x),nrow=18)%*%(as.matrix(fit$coef))}
	  results <- crossval(t(as.matrix(x,nrow=18)),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")
}
shrinkage(lm)
# 代码有点问题

深层次分析——变量的相对重要性

# 标准化（先将数据标准化后再做回归分析）
zdata <- scale(data) 

# 相对权重（对所有可能子模型添加一个预测变量引起的R平方平均增加量的近似值，结果显示各变量对R平方的解释比例）
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)