chapter 3 -回归试验

线性回归

library(MASS)#R中自带的有，未安装的时候可以有
library(ISLR)#加载ISLR数据包

#使用MASS库中的Boston房价数据集，查看变量名
names(Boston)
[1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"     "dis"     "rad"    
[10] "tax"     "ptratio" "black"   "lstat"   "medv"   

#简单线性回归

lm.fit <- lm(medv~lstat,data=Boston)
summary(lm.fit)
Call:
lm(formula = medv ~ lstat, data = Boston)

Residuals:
    Min      1Q  Median      3Q     Max 
-15.168  -3.990  -1.318   2.034  24.500 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 34.55384    0.56263   61.41   <2e-16 ***
lstat       -0.95005    0.03873  -24.53   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.216 on 504 degrees of freedom
Multiple R-squared:  0.5441,    Adjusted R-squared:  0.5432 
F-statistic: 601.6 on 1 and 504 DF,  p-value: < 2.2e-16

#提取系数估计
coef(lm.fit)
(Intercept)       lstat 
 34.5538409  -0.9500494 

#提取系数估计的置信区间
confint(lm.fit)
  2.5 %     97.5 %
(Intercept) 33.448457 35.6592247
lstat       -1.026148 -0.8739505


#计算置信区间和预测区间,newdata要写出数据框的形式，哪怕只有一个数值。
predict(lm.fit,newdata =data.frame(lstat=c(5,10,15)),interval="confidence") 
  fit      lwr      upr
1 29.80359 29.00741 30.59978
2 25.05335 24.47413 25.63256
3 20.30310 19.73159 20.87461

predict(lm.fit,newdata = data.frame(lstat=c(5,10,15)),interval="prediction")
     fit       lwr      upr
1 29.80359 17.565675 42.04151
2 25.05335 12.827626 37.27907
3 20.30310  8.077742 32.52846

正如预期的一样，置信区间和预测区间有相同的中心，但预测区间比置信区间要宽。

#绘制预测变量和响应变量的散点图和最小二乘线
attach(Boston)
plot(lstat,medv)
abline(lm.fit,col="red")

#根据散点图，似乎有二次关系。
#lm()的输出结果直接用plot()命令将自动生成四幅诊断图
#用par()函数可以直接查看
par(mfrow=c(2,2))
plot(lm.fit)

查看离群点

#绘制残差对拟合值的散点图
plot(predict(lm.fit),residuals(lm.fit))
#绘制学生化残差对拟合值的散点图
plot(predict(lm.fit),residuals(lm.fit))

#量化观测的杠杆作用
plot(hatvalues(lm.fit))

hatvalues()函数可以计算任意多个预测变量的杠杆统计量

#返回向量中最大杠杆的索引
which.max(hatvalues(lm.fit))

#方差膨胀因子，用car包中的vif函数
library(car)
#检验预测变量的多重共线性
vif(lm.fit)

交互项，语句lstat:black可以将lstat和black的交互项纳入模型，
语句lstat*age将lstat,age和交互项lstat:black同时作为预测变量，是一种简写形式。

summary(lm(medv~lstat*age,data=Boston))

#预测变量的非线性变换 ，加入平方项：I(X^2)
lm.fit2 <- lm(medv~lstat+I(lstat^2))
summary(lm.fit)

#模型二
lm.fit <- lm(medv~lstat)

#比较两个模型的优劣
anova(lm.fit,lm.fit2)

anova函数通过假设检验比较两个模型
Analysis of Variance Table

Model 1: medv ~ lstat
Model 2: medv ~ lstat + I(lstat^2)
Res.Df   RSS Df Sum of Sq     F    Pr(>F)    
1    504 19472                                 
2    503 15347  1    4125.1 135.2 < 2.2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
 模型二更显著

#查看诊断图
par(mfrow=c(2,2))
plot(lm.fit)


#向模型中加入高阶多项式,5次多项式
lm.fit5 <- lm(medv~poly(lstat,5))
summary(lm.fit5)
Call:
lm(formula = medv ~ poly(lstat, 5))

Residuals:
 Min       1Q   Median       3Q      Max 
-13.5433  -3.1039  -0.7052   2.0844  27.1153 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept)       22.5328     0.2318  97.197  < 2e-16 ***
poly(lstat, 5)1 -152.4595     5.2148 -29.236  < 2e-16 ***
poly(lstat, 5)2   64.2272     5.2148  12.316  < 2e-16 ***
poly(lstat, 5)3  -27.0511     5.2148  -5.187 3.10e-07 ***
poly(lstat, 5)4   25.4517     5.2148   4.881 1.42e-06 ***
poly(lstat, 5)5  -19.2524     5.2148  -3.692 0.000247 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 5.215 on 500 degrees of freedom
Multiple R-squared:  0.6817,    Adjusted R-squared:  0.6785 
F-statistic: 214.2 on 5 and 500 DF,  p-value: < 2.2e-16


#进行多项式变换
summary(lm(medv~log(rm),data=Boston))

#定性预测变量

names(Carseats)
#查看定性预测变量前几行
head(Carseats$ShelveLoc) 

#构建含有交互项的多元回归模型，给出定性变量，R将自动生成虚拟变量
lm.fit <- lm(Sales~.+Income:Advertising+Price:Age,data=Carseats)
summary(lm.fit)

#查看虚拟变量的编码
attach(Carseats)
contrasts(ShelveLoc)