[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（）函数实现
 
#深层次分析：交叉验证和相对重要性（以后再研究）