文章目录
1. 原始数据的分析2. 回归模型的拟合(参数估计和检验)3. 变量选择4. 回归诊断5. 改进措施6. 岭回归8. 其他基本假设：正态性、独立性、线性、同方差性
常用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)%*%Ye <- Y-X%*%Beta #<or> resid <- resid(lm) #残差的计算SSe <- t(e)%*%esigma2 <- 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)%*%YRSS_H <- n*mean(Y)^2SS_He <- t(Y)%*%Y-RSS_HRSS <- t(Beta)%*%t(X)%*%Y || RSS <- RSS_H + RSS_ISSe <- t(Y)%*%Y-RSSF <- (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 / TSSAdjusted_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*qBIC <- (-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是replacetest=(!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)]%*%coefival.errors[i]=mean((Hitters$Salary[test]-pred)^2)}val.errorswhich.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=10set.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)#这里可以简化为predictcv.errors[j,i]=mean( (dataset$y[folds==j]-pred)^2)}}mean.cv.errors=apply(cv.errors,2,mean)mean.cv.errorspar(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:]));Xkappa(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)] lambday_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)^2r2cv <- cor(y,results$cv.fit)^2cat("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$vectorsev <- svd$valuesdelta <- diag(sqrt(ev))lambda <- evec %*% delta %*% t(evec)lambdasq <- lambda^2beta <- solve(lambda) %*% rxyrsquare <- colSums(beta^2)rawwgt <- lambdasq %*% beta^2import <- (rawwgt / rsquare)*100import <- 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)
如果觉得《R语言与线性回归分析》对你有帮助，请点赞、收藏，并留下你的观点哦！
