R语言_回归分析

#lm拟合回归模型
#简单线性回归
fit = lm(weight~height,data=women)
summary(fit)
women$weight
fitted(fit)
residuals(fit)
plot(density(residuals(fit)))
plot(women$height,women$weight,
     xlab="hight",
     ylab="weight")
abline(fit)
#多项式回归
fit2 = lm(weight~height+I(height^2),data=women)
summary(fit2)
plot(women$height,women$weight,
     xlab="hight",
     ylab="weight")
lines(women$height,fitted(fit2))
#scatterplot函数
library(car)
scatterplot(weight~height,
            data=women,
            spread=FALSE,
            lty.smooth=2,
            pcj=19)
#多元线性回归
states = as.data.frame(state.x77[,c("Murder","Population",
                                    "Illiteracy","Income","Frost")])
cor(states)
scatterplotMatrix(states,
                  spread=FALSE,
                  lty.smooth=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)
fit = lm(mpg~hp*wt,data=mtcars)
summary(fit)
install.packages("effects")
library(effects)
plot(effect("hp:wt",fit,xlevels=list(wt=c(2.2,3.2,4.2))),multiline=TRUE)

#回归推断
fit = lm(Murder~Population+Illiteracy+Income+Frost,
         data=states)
confint(fit)
#标准方法
#检验正态性，独立性，线性，同方差性
fit = lm(weight~height,data=women)
par(mfrow=c(2,2))
plot(fit)
fit2 = lm(weight~height+I(height^2),data=women)
plot(fit2)
newfit = lm(weight~height+I(height^2),data=women[-c(13,15),])
plot(newfit)
#又一个例子
fit = lm(Murder~Population+Illiteracy+Income+Frost,
         data=states)
plot(fit)
#改进的方法，更加推荐
#正态性
library(car)
fit = lm(Murder~Population+Illiteracy+Income+Frost,
         data=states)
par(mfrow=c(1,1))
qqPlot(fit,labels=row.names(states),
       id.method="identity",
       simulate=TRUE,main="Q-Q Plot")
states["Nevada",]
fitted(fit)["Nevada"]
residuals(fit)["Nevada"]
rstudent(fit)["Nevada"]
#误差的独立性
#针对误差的类时间序列，滞后项是1
durbinWatsonTest(fit)
#线性
library(car)
crPlots(fit)
#同方差性
library(car)
ncvTest(fit)
spreadLevelPlot(fit)
#线性模型假设的综合验证
install.packages("gvlma")
library(gvlma)
gvmodel = gvlma(fit)
summary(gvmodel)
#多重共线性
library(car)
vif(fit)
sqrt(vif(fit))>2  #>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)
  identity(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")
library(car)
influencePlot(fit)