[Rcode]线性回归

#简单线性回归： ##常用绘图： fit<-lm(weight~height,data=women) summary(fit) plot(women$height,women$weight,xlab="Height (in inches)",ylab="Weight (in pounds)") abline(fit) fit2<-lm(mpg~wt+I(wt^2),data=mtcars) summary(fit2)  plot(mtcars$wt,mtcars$mpg)  lines(mtcars$wt,fitted(fit2))   ##最好用car包中的scatterplot（）函数，显示的内容更全面 library(car) scatterplot(weight~height,data=women,spread=F,smoother.args=list(lty=2),pch=19,main="Women Age 30 -39",xlab="Height (in inches)",ylab="Weight (in pounds)")  ###该图提供了散点图、线性拟合图和平滑拟合（loess）曲线，还展示了两个变量的箱线图。spread=F删除了残差正负均方根在平滑曲线的展开和非对称信息。smoother.args设置loess拟合曲线为虚线。pch=19设置为实心圆。 #多元线性回归 ##lm函数需要数据框，这里把数据转化为数据框 ##多元回归分析第一步总是检查变量之间的相关性： states<-as.data.frame(state.x77[,c("Murder","Population","Illiteracy","Income","Frost")]) ##检查变量相关性,用car包的函数可视化比较强 cor(states) library(car) scatterplotMatrix(states,spread=F,smoother.args=list(lty=2),main="Scatter Plot Matrix")  ##线性回归 fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states) summary(fit)  ##考虑交互项： ###交互项显著说明：相应变量与其中一个预测变量的关系依赖于另一个变量的水平 fit<-lm(mpg~hp+wt+hp:wt,data=mtcars) ###展示交互项时可用effects包中的effect函数 install.packages("effects") library(effects)  plot(effect("hp:wt",fit,,list(wt=c(2.2,3.2,4.2))),multiline=T)  #回归诊断： ##可用于查看置信区间 fit<-lm(Murder~Population+Illiteracy+Income+Frost,data=states) confint(fit)   ##标准方法：可查看fit的四幅图形 par(mfrow=c(2,2)) plot(fit)  ###解决方法：二次拟合或者删除异常点    ##改进方法： ###（1）正态性： library(car) qqPlot(fit,labels=row.names(states),id.method="identify",simulate=T,main="Q-Q Plot")  ###id.method可以交互式绘图，simulat=T时，95%置信区间会用参数自助法生成。 ###（2）误差的独立性：car包中的函数可做Duibin-Watson检验 library(car) durbinWatsonTest(fit) ###（3）线性:通过成分残差图，查看自变量与因变量是否呈现非线性关系 library(car) crPlots(fit)  ###（4）同方差性： library(car) ncvTest(fit)  spreadLevelPlot(fit)  ###若代码建议幂次转换为0.5，则用sqrt（y）代替y，如果幂次转换建议为0，则使用对数变换，如果接近1，说明异方差性不是很明显，不需要做变换。  ###（5）线性模型假设的综合验证 ####多重共线性:会导致模型参数的置信区间过大，使单个系数解释起来很困难 library(car) vif(fit)  sqrt(vif(fit))>2  ####一般情况下，sqrt（vif）大于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 Value")   abline(h=c(2,3)*p/n,col="red",lty=2)   identify(1:n,hatvalues(fit),names(hatvalues(fit))) } hat.plot(fit)  ####强影响点：对模型参数估计值影响有些比例失衡的点。用Cook距离，即D统计量。 cutoff<-4/(nrow(states)-length(fit$coefficients)-2) plot(fit,which=4,cook.levels=cutoff)  abline(h=cutoff,lty=2,col="red")  #####这部分没怎么看懂 #改进措施： ##删除观测点 ##变量变换：  ###当模型违反正态假设时，可以对响应变量尝试某种变换。  library(car)  summary(powerTransform(states$Murder))  ###违背线性假设时，对预测变量进行变换常常比较有用。  library(car)  boxTidwell(Murder~Population+Illiteracy,data=states) ##增删变量  #选择最佳的回归模型： ##模型比较 fit1<-lm(Murder~Population+Illiteracy,data="states") fit2<-lm(Murder~Population+Illiteracy+Income+Frost,data="states") anova(fit1,fit2) ##也可以用AIC原则 AIC(fit1,fit2)  ##逐步回归： library(MASS) stepAIC(fit,direction = "backward") ##全子集回归可以用leaps包中的regsubsets（）函数实现 #深层次分析：交叉验证和相对重要性（以后再研究）