ISLR_LinearRegression

• Simple linear regression
  • Estimating the coefficients
  • Assessing the accuracy of the coeff
    • conference interval
    • hypothesis test
  • Assessing the accuracy of model
    • RSE
    • R2
• Multiple linear regression
  • add simple linear regressiom
  • Estimating the regression coefficients
  • Some import questions
• Other considerations in regression model
  • Qualitative predictors
  • Extensions of LR
    • additive assumption
    • non-linear relationship
  • Potential problems
    • non-linearity
    • correlation of error terms
    • non-constant variance of error terms
    • outliers
    • high-leverage points
    • collinearity
• Marketing plan
• Comparison of ls with knn
• Lab
  • libraries
  • simple linear regression
  • multiple linear regression
  • interaction terms
  • non-linear transformations of predictor
  • qualitative predictor
  • writing functions
• Exercise
• 参考资料

Simple linear regression

考虑Statistical Learning中的销售与三种媒介广告的关系。 提炼出以下几个问题：

1. 销售与媒介广告之间有关系么
2. 如果有关系的话，关系有多强，是不是只是比随机猜测好一点
3. 三种媒介广告，哪种对销售的促进作用更强
4. 有多大的精度，预测每个媒介广告对销售的促进作用
5. 有多大的精度，预测三种媒介广告对销售的促进作用
6. 销售与媒介广告之间的关系是线性的么
7. 在媒介广告之间是否存在synergy，或者说在在统计上是否存在interaction effect现象

Estimating the coefficients

线性回归首先应该估计线性回归的系数。

估计系数有两种方法：

1. normal equation
2. gradient descent

Assessing the accuracy of the coeff

conference interval

根据Y=f(x)+ϵY=f(x)+\epsilon，结合线性模型，得到：

Y=β0+β1X+ϵ

Y = \beta_0 + \beta_1 X + \epsilon

其中， ϵ\epsilon是mean-zero的误差项，代表非线性的因素与其他测量误差等，并且一般假设ϵ\epsilon和XX是独立的。

上图中，红色代表真实的Y=2+3XY=2+3X，点根据分布Y=2+3X+ϵY=2+3X+\epsilon产生，蓝色代表根据不同的数据点利用least squares拟合出的直线。可以看到，不同的数据点，拟合出的蓝线和红线有微小差异，但是总体很接近。 总结来说，就是先假设数据是符合线性关系的。然后利用样本的参数去估计群体的线性回归的参数。

具体来说： 先利用normal equation计算出参数的点估计值

然后各个参数的standard error为

其中，σ2\sigma^2代表Var(ϵ)Var(\epsilon)，这个值可以用RSS去估计

RSE=RSS/(n−2)‾‾‾‾‾‾‾‾‾‾‾√

RSE = \sqrt{RSS/(n-2)}

接下来，就可以计算置信区间了，95%的置信区间为

hypothesis test

同样地，根据standard error也可以进行hypothesis test

先计算t统计量

然后，根据t的大小决定是否推翻null hypothesis

Assessing the accuracy of model

一旦确定了线性关系的存在，那么很自然的一个问题是：模型拟合数据的程度是多少？ 评估的依据有两种统计量：

1. RSE(residual standard error)
2. R2R^2

RSE

模型的X和Y永远不可能完全拟合，即使知道了真实的ff。 原因是，有误差项ϵ\epsilon的存在。

RSE，反映的是预测值偏离真实值的平均程度。根据自由度做了保守估计。

例如，RSE=3.26，意味着每个市场的预测销售平均偏离真正的model大约3260单位。

RSE一般用来衡量模型的fit情况。如果RSE相对预测值很小，那么表明平均下来，预测值和真实值很接近，那么这个模型fit地就很好。

R2

RSE和Y的scale有很大关系，如果Y的scale很大，那么RSE也很大，解决办法一种是求RSE/mean(Y)，另一种是计算R2R^2统计量。

R2R^2的含义是：Y有多少的variance被X成功解释了？

R2R^2的公式是：

其中，TSS表示的是原有的YY的variance，RSS表示没有被解释的variance，TSS-RSS表示已经被解释的variance。R2越大，表示X解释Y的成分越多。 没有被解释的原因，可能是模型的错误，也可能是σ2=Var(ϵ)\sigma^2=Var(\epsilon)很大。

R2R^2好坏的标准是：根据实际应用定，精确的模型需要精确的R2；market等粗略的忽略了很多其他因素的模型需要很小的R2就足够。

R2R^2和相关系数的关系是： r2=R2r^2=R^2（只在简单线性回归合适，多元回归相关系数不可用），也就是说相关系数的平方代表了决定系数，表示一个变量能被另一个变量解释的比例。

Multiple linear regression

add simple linear regressiom

如果根据每个predictor都实行一个简单线性回归，会遇到以下问题：

1. 无法将多个predictor在一个模型中混合
2. 每个单独的模型中忽略了其他predictor的作用，这可能会造成很大的misunderstanding。

Estimating the regression coefficients

转化为最小值的优化问题

多元线性回归中，会遇到一个问题：Xi的系数和简单线性回归差别很大，甚至原来统计显著的参变量会变得统计不显著。比如newspaper。 原因主要在： 简单线性回归中，忽视了其他的predictor。多元线性回归中，假设其他predictor不变。

那么，newspaper到底和sale有关系么，关系多强呢？ 答案是：即使简单线性回归说两者有关系，但是要根据多元线性回归的结果，两者并没有关系。 解释是：sale和newspaper并没有关系。但是newspaper和ratio的相关系数很高，因此ratio提高了，newspaper提高了，sale提高了。 就像shark attacks和ice cream的关系，简单线性回归两者是有关系的，但是加入temperature后两者之间的关系就不显著了。

Some import questions

Q1: 是否至少有一个predictor是有用的，线性关系的存在性？

Q2: 所有的X都有用还是只有subset部分有用，X的效用性？

Q3: 模型fit的程度如何，模型的评价？

Q4: 具体的模型，以及模型的精度？

Q1: 是否至少有一个predictor是有用的，线性关系的存在性？ A1: 使用假设检验的F统计量。

如果线性关系存在的话，

如果线性关系不存在的话，

也就是说，如果线性关系不存在，那么F约等于1，线性关系存在F大于1。

当H0为真，ϵ\epsilon呈现正态分布的情况下，F统计量遵循F分布（即使ϵ\epsilon不呈现正态分布，如果sample size n足够大的话依然满足F分布）

上面的H0中，是所有X的系数都为0，有的时候只想要检验一部分subset大小为q的X的系数为0。

In this case we fit a second model that uses all the variables except those last q. Suppose that the residual sum of squares for that model is RSS0.

下图中newspaper的系数t就表示：t的平方就是缺少了nespaper对应的F统计量。p值很大表示不足以推翻只有newpaper系数为0的H0假设。

有一个问题，按照上图的逻辑：每一个predictor都有对应的t和p因此我们知道了该predictor对模型的重要程度，那么为什么还需要F？或者如果有一个predictor的p值显著的话，是不是就能够断定线性关系存在呢？ 答案是否定的。假定共有100个X，每一个X和Y都没有任何关系，但是由于随机误差都有5%的可能性得到p-value<0.05。这时候，如果只是根据每一个X的p来判断线性模型是否成立，那么成立的可能性高达1-0.95^100 = 0.994。所以不能用这种方法，尤其p很大的时候。

F检定一般适用于p < n的情况。当p>n的时候，不能使用least squares去fit线性模型。可以考虑forward selection。

Q2: 所有的X都有用还是只有subset部分有用，X的效用性？ A2: 首先用F统计量判断线性模型的有效性，然后开始寻找重要的参数。 遍历的方法一共有2n2^n种情况，一般来说采用启发式的方法：

1. forward selection: 先从空模型开始，然后逐次往模型中加入predictor，加入的准则是RSS最小，加入的终止条件是RSS变化小于一定范围。
2. backward selection: 从包含所有predictor的模型开始，移去p值最高的predictor，直到所有的p-value < threshold。
3. mixed selection: 将上面两种方法混合，从空模型开始，不断加入predictor，但是加入的同时注意每个predictor的p值，保留小于threshold的p值，去除高于threshold的p值，直到遍历了每一个predictor。

Q3: 模型fit的程度如何，模型的评价？

R2 简单线性回归中，R2=r2R^2=r^2，多元线性回归中，R2=Cor(Y,Ŷ )2R^2 = Cor(Y, \hat{Y})^2。多元回归一个重要的特征就是最大化了Cor(Y,Ŷ ) Cor(Y, \hat{Y})。 R2如果变化大的话那么可以考虑加入predictor，反之不需要加。因为对应的predictor如果只是获得了很小的R2提升，那么很有可能是对数据过度拟合造成的。这时候可以查看p-value以及RSE做决定。

RSE 随着predictor个数的增加，可能会出现RSS变小但是RSE变大的情况。

PIC 有的时候，画图是可以发现潜在的不易被数学描述的协同关系的。

Q4: 具体的模型，以及模型的精度？

这里包括三种不确定性：

1. 系数β\beta的不确定性
2. model bias，假设线性模型
3. 随机误差ϵ\epsilon。注意predict interval往往比confidence interval更大，原因是既包括了对f(x)估计的误差，也包括了ϵ\epsilon的误差。同时，predict interval可以应用于每一个对象，confidence interval往往用于多个对象的平均值。

Other considerations in regression model

Qualitative predictors

如果只有两个level的话 可以一个设置为基础0，一个设置为偏置。 根据TABLE3.7，男性的均值是509，女性的均值是509+19。同时，注意到gender的p值是不显著的，表明均值上并无男女的性别差异。

同时，也可以改变编码值，对结果无影响，但是对于模型的解释有影响。

如果有多个level的话 选择上设置多个dummy variable，并且其个数永远比level少一个，也就是需要一个baseline。baseline的设置不影响最后的每个level的值但是影响p值。 因此，不能只是根据单个level的p值判断线性关系是否存在，应该使用F检定。

Extensions of LR

线性回归有两个很重要的假设（现实中经常不符合）：

1. additive(means that the effect of changes in Xi on the Y is independent of values of other predictors)
2. linear(means that the change in the response Y due to one-unit change in X is constant regardless of the value of X)

additive assumption

Removing the additive assumption

上图中，X1每增加一个单位对Y的影响，和X2的值是有关系的，也就是两者之间有商业上的synergy，统计上的interaction现象。

有的时候，会出现interaction很显著但是main effect不显著的情况，这时候仍然需要包含main effect。

如果是factor数据的intersection

non-linear relationship

Potential problems

non-linearity

如果数据本身是非线性的，那么讨论的结果的意义就值得怀疑。 这时候一般使用residual−xresidual-x或者residual−ŷ residual-\hat{y}图来看模型是否符合线性的假设。 如果不符合的话，一个简单的方法是进行变换，比如log(X),X‾‾√,X2log(X), \sqrt{X},X^2等。

correlation of error terms

error term相关的坏处 低估了standard error，进而对CI,HT造成影响。

errors不相关的含义是，the fact that εi is positive provides little or no information about the sign of εi+1，即一个error term的情况并不会对其他的造成影响。 如果相关，会造成standard error变大的结果，对confidence interval和hypothesis test造成影响。p值会比实际的偏低，这样会造成误判统计显著。模型的置信得不到保障。

The standard errors that are computed for the estimated regression coefficients or the fitted values are based on the assumption of uncorrelated error terms. If in fact there is correlation among the error terms, then the estimated standard errors will tend to underestimate the true standard errors.

error term的相关一般在time series中出现的较多。 同时，在比较身高与体重的关系中，如果调查的对象是一家人或者同一个环境中的人的话，也会出现error的相关。

non-constant variance of error terms

that the error terms have a constant variance, Var(εi) = σ2. The standard errors, confidence intervals, and hypothesis tests associated with the linear model rely upon this assumption.

outliers

outliers可能对model fit没有影响，但是会影响CI，HT。 比如包含与不包含outlier的RSE分别是1.09和0.77（RSE用来计算CI等），R2分别是0.805和0.892。 一般使用studentized residuals图去侦测outliers。

high-leverage points

outliers指的是有一个不寻常的Y，high-leverage是有一个不寻常的X。 41是一个很危险的点，因为既是outlier又是high-leverage。

high-leverage的计算如下

collinearity

含义:

Collinearity refers to the situation in which two or more predictor variables are closely related to one another.

问题：

The presence of collinearity can pose problems in the regression context, since it can be difficult to separate out the individual effects of collinear variables on the response. In other words, since limit and rating tend to increase or decrease together, it can be difficult to determine how each one separately is associated with the response, balance. 无法估计每个X单独对Y的作用 减少了模型系数估计的精度，造成了standard error的上升，t统计量的下降，因此拒绝H0的能力下降，侦测出非0系数的能力下降。

下图中，M1无共线性，M2有共线性。limit实际上很重要，但是因为共线性造成它的p值很低，误判不重要。

侦测共线性，一方面可以做相关系数矩阵。 但是也可能出现共线性存在但是相关系数矩阵都很小的情况，这叫做multi-collinearity。 解决的办法是计算VIF(variance inflation factor)，这个值是1意味着没有一点共线性，大于5到10说明共线性问题很显著。

解决途径

1. 扔掉对应的variable，因为共线性很高，这个variable解释的东西也可以被其他的解释掉，所以不用担心。
2. 将共线性的变量结合成一个新的变量。

Marketing plan

回答前面提出的问题

1. 销售与媒介广告之间有关系么 先拟合多元线性回归，然后用F统计量去坚定H0
2. 如果有关系的话，关系有多强，是不是只是比随机猜测好一点

可以用RSE和R2去判断 RSE估计了Y偏离拟合直线的标准距离，本题是1681，大约占了meanY-14022的12%。 R2为90，表明了解释了很高的Y variation。

3 三种媒介广告，哪种对销售的促进作用更强

可以查看多元拟合的p值，newspaper的p值很低，表明在有TV和ratio的情况下，可以将其去除。 同时，也可以采用forward selection和bacword selection等进一步去判断重要的变量。

4 有多大的精度，预测每个媒介广告对销售的促进作用

根据多元回归的每个系数的点估计以及对应的se算出置信区间CI，CI包括0说明这个系数不是统计显著的。 同时，共线性会对se有很大影响，造成se偏大。因此，对tv, ratio, newspaper三者的VIF做检定，分别是1.005, 1.145, 1.145，都很小接近1，远小于10的阈值，因此没有共线性。

5 有多大的精度，预测未来的销售

预测的点估计可以直接采用线性回归，预测的区间考虑每个回归系数的CI。 同时需要注意，对总体平均的预测f(X)f(X)采用confidence interval，对单个群体的预测采用Y=f(x)+ϵY=f(x)+\epsilon，即prediction interval，其值考虑了单个obes的误差，较大。

6 销售与媒介广告之间的关系是线性的么

可以使用residual plot，如果没有明显模式存在的话就是线性的。

7 在媒介广告之间是否存在synergy，或者说在在统计上是否存在interaction effect现象

直观的方法是画图 或者可以做出intersection的模型，然后看对应的p值是不是统计显著，RSE有没有下降，R2有没有显著的上升。

Comparison of ls with knn

ls是参数方法，knn是非参数方法。

Q：什么时候ls比knn好么 A：当关系真实是线性或者很接近线性的时候，但也只是好一点而已

Q：什么时候knn比ls好么 A：当关系远非线性，这时候会好很多

Q：按照上面的讨论，线性的话knn稍差，非线性的话knn超好，那么是否现实生活中（大部分问题都是非线性的）直接用knn就可以了 A：现实生活中，ls一般比knn效果好，尤其是高维情况下。原因是数据量太少，比如100个数据，在二维空间中可能很密集，但是在100维空间中，即使距离最近的点可能也很远，这叫做curse of dimensionality，因此一般线性的方法会更好一些。

Lab

libraries

library(MASS)
library(ISLR)

simple linear regression

# Simple Linear Regression

lm.fit=lm(medv~lstat,data=Boston)
attach(Boston)
lm.fit=lm(medv~lstat)
lm.fit
summary(lm.fit)
names(lm.fit)
coef(lm.fit)
confint(lm.fit)
predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval="confidence")
predict(lm.fit,data.frame(lstat=(c(5,10,15))), interval="prediction")

plot(lstat,medv)
abline(lm.fit)
abline(lm.fit,lwd=3)
abline(lm.fit,lwd=3,col="red")
plot(lstat,medv,col="red")
plot(lstat,medv,pch=20)
plot(lstat,medv,pch="+")
plot(1:20,1:20,pch=1:20)
par(mfrow=c(2,2))
plot(lm.fit)
plot(predict(lm.fit), residuals(lm.fit))
plot(predict(lm.fit), rstudent(lm.fit))
#find leverage points location
plot(hatvalues(lm.fit))
which.max(hatvalues(lm.fit))

multiple linear regression

lm.fit=lm(medv~lstat+age,data=Boston)
summary(lm.fit)
lm.fit=lm(medv~.,data=Boston)
summary(lm.fit)
?summary.lm
summary(lm.fit)$r.sq  #r2
summary(lm.fit)$sigma  #rse


library(car)
vif(lm.fit)
lm.fit1=lm(medv~.-age,data=Boston)
summary(lm.fit1)
lm.fit1=update(lm.fit, ~.-age)

interaction terms

summary(lm(medv~lstat*age,data=Boston)) #包含main effect
summary(lm(medv~lstat:age,data=Boston)) #不包含main effect

non-linear transformations of predictor

lm.fit2=lm(medv~lstat+I(lstat^2))
summary(lm.fit2)
lm.fit=lm(medv~lstat)
anova(lm.fit,lm.fit2)
par(mfrow=c(2,2))
plot(lm.fit2)
plot(lm.fit)
#ANOVA检定的F远大于1，说明了lm.fit2更优，同时比较residual plot，可以看出更优的原因是更符合线性假设。
lm.fit5=lm(medv~poly(lstat,5))
summary(lm.fit5)
summary(lm(medv~log(rm),data=Boston))

qualitative predictor

fix(Carseats)
names(Carseats)
lm.fit=lm(Sales~.+Income:Advertising+Price:Age,data=Carseats)
summary(lm.fit)
attach(Carseats)
contrasts(ShelveLoc)
?contrasts

writing functions

LoadLibraries
LoadLibraries()
LoadLibraries=function(){
 library(ISLR)
 library(MASS)
 print("The libraries have been loaded.")
 }
LoadLibraries
LoadLibraries()

Exercise

Q8:

Auto = read.csv("../data/Auto.csv", header=T, na.strings="?")
Auto = na.omit(Auto)
summary(Auto)
attach(Auto)
lm.fit = lm(mpg ~ horsepower)
summary(lm.fit)

#RSE and R2
mean(mpg, na.rm=T)
summary(lm.fit)$sigma
summary(lm.fit)$sigma/mean(mpg, na.rm=T)
summary(lm.fit)$r.sq
?summary.lm

#confidence interval and prediction interval
predict(lm.fit, data.frame(horsepower=c(98)), interval="confidence")
predict(lm.fit, data.frame(horsepower=c(98)), interval="prediction")

# linear fit 
plot(horsepower, mpg)
abline(lm.fit)

# diagnose plot, evidence for non-linearity.
par(mfrow=c(2,2))
plot(lm.fit)

Q9:

pairs(Auto)
cor(subset(Auto, select=-name))
lm.fit1 = lm(mpg~.-name, data=Auto)
summary(lm.fit1)

par(mfrow=c(2,2))
plot(lm.fit1)
plot(predict(lm.fit1), rstudent(lm.fit1))

lm.fit2 = lm(mpg~cylinders*displacement+displacement*weight)
summary(lm.fit2)

# 诊断模型
# 线性假设，误差方差齐性假设，误差正态分布假设
# 1
lm.fit3 = lm(mpg~log(weight)+sqrt(horsepower)+acceleration+I(acceleration^2))
summary(lm.fit3)
par(mfrow=c(2,2))
plot(lm.fit3)
plot(predict(lm.fit3), rstudent(lm.fit3))
# 2
lm.fit2<-lm(log(mpg)~cylinders+displacement+horsepower+weight+acceleration+year+origin,data=Auto)
summary(lm.fit2)
par(mfrow=c(2,2)) 
plot(lm.fit2)
plot(predict(lm.fit2),rstudent(lm.fit2))

Q10：

library(ISLR)
summary(Carseats)
attach(Carseats)
lm.fit = lm(Sales~Price+Urban+US)
summary(lm.fit)

lm.fit2 = lm(Sales ~ Price + US)
summary(lm.fit2)

confint(lm.fit2)

plot(predict(lm.fit2), rstudent(lm.fit2))
# 没有outlier
par(mfrow=c(2,2))
plot(lm.fit2)
# leverage点的标准是 (p+1)/n

Q11:

set.seed(1)
x = rnorm(100)
y = 2*x + rnorm(100)

lm.fit = lm(y~x+0)
summary(lm.fit)

lm.fit = lm(x~y+0)
summary(lm.fit)

lm.fit = lm(y~x)
lm.fit2 = lm(x~y)
summary(lm.fit)
summary(lm.fit2)

Q12：

set.seed(1)
x = rnorm(100)
y = 2*x
lm.fit = lm(y~x+0)
lm.fit2 = lm(x~y+0)
summary(lm.fit)
summary(lm.fit2)

set.seed(1)
x <- rnorm(100)
y <- -sample(x, 100)
sum(x^2)
sum(y^2)
lm.fit <- lm(y~x+0)
lm.fit2 <- lm(x~y+0)
summary(lm.fit)
summary(lm.fit2)

Q13:

set.seed(1)
x = rnorm(100)
eps = rnorm(100, 0, sqrt(0.25))
y = -1 + 0.5*x + eps

plot(x, y)

lm.fit = lm(y~x)
summary(lm.fit)

plot(x, y)
abline(lm.fit, lwd=3, col=2)
abline(-1, 0.5, lwd=3, col=3)
legend(-1, legend = c("model fit", "pop. regression"), col=2:3, lwd=3)

lm.fit_sq = lm(y~x+I(x^2))
summary(lm.fit_sq)

set.seed(1)
eps1 = rnorm(100, 0, 0.125)
x1 = rnorm(100)
y1 = -1 + 0.5*x1 + eps1
plot(x1, y1)
lm.fit1 = lm(y1~x1)
summary(lm.fit1)
abline(lm.fit1, lwd=3, col=2)
abline(-1, 0.5, lwd=3, col=3)
legend(-1, legend = c("model fit", "pop. regression"), col=2:3, lwd=3)

set.seed(1)
eps2 = rnorm(100, 0, 0.5)
x2 = rnorm(100)
y2 = -1 + 0.5*x2 + eps2
plot(x2, y2)
lm.fit2 = lm(y2~x2)
summary(lm.fit2)
abline(lm.fit2, lwd=3, col=2)
abline(-1, 0.5, lwd=3, col=3)
legend(-1, legend = c("model fit", "pop. regression"), col=2:3, lwd=3)

confint(lm.fit)
confint(lm.fit1)
confint(lm.fit2)

Q14:

set.seed(1)
x1 = runif(100)
x2 = 0.5 * x1 + rnorm(100)/10
y = 2 + 2*x1 + 0.3*x2 + rnorm(100)

cor(x1, x2)
plot(x1, x2)

lm.fit = lm(y~x1+x2)
summary(lm.fit)

lm.fit = lm(y~x1)
summary(lm.fit)
lm.fit = lm(y~x2)
summary(lm.fit)

x1 = c(x1, 0.1)
x2 = c(x2, 0.8)
y = c(y, 6)
lm.fit1 = lm(y~x1+x2)
summary(lm.fit1)
lm.fit2 = lm(y~x1)
summary(lm.fit2)
lm.fit3 = lm(y~x2)
summary(lm.fit3)

Q15:

library(MASS)
summary(Boston)
Boston$chas <- factor(Boston$chas, labels = c("N","Y"))
summary(Boston)
attach(Boston)
lm.zn = lm(crim~zn)
summary(lm.zn) # yes
lm.indus = lm(crim~indus)
summary(lm.indus) # yes
lm.chas = lm(crim~chas) 
summary(lm.chas) # no
lm.nox = lm(crim~nox)
summary(lm.nox) # yes
lm.rm = lm(crim~rm)
summary(lm.rm) # yes
lm.age = lm(crim~age)
summary(lm.age) # yes
lm.dis = lm(crim~dis)
summary(lm.dis) # yes
lm.rad = lm(crim~rad)
summary(lm.rad) # yes
lm.tax = lm(crim~tax)
summary(lm.tax) # yes
lm.ptratio = lm(crim~ptratio)
summary(lm.ptratio) # yes
lm.black = lm(crim~black)
summary(lm.black) # yes
lm.lstat = lm(crim~lstat)
summary(lm.lstat) # yes
lm.medv = lm(crim~medv)
summary(lm.medv) # yes

lm.all = lm(crim~., data=Boston)
summary(lm.all)

x = c(coefficients(lm.zn)[2],
      coefficients(lm.indus)[2],
      coefficients(lm.chas)[2],
      coefficients(lm.nox)[2],
      coefficients(lm.rm)[2],
      coefficients(lm.age)[2],
      coefficients(lm.dis)[2],
      coefficients(lm.rad)[2],
      coefficients(lm.tax)[2],
      coefficients(lm.ptratio)[2],
      coefficients(lm.black)[2],
      coefficients(lm.lstat)[2],
      coefficients(lm.medv)[2])
y = coefficients(lm.all)[2:14]
plot(x, y)

lm.zn = lm(crim~poly(zn,3))
summary(lm.zn) # 1, 2
lm.indus = lm(crim~poly(indus,3))
summary(lm.indus) # 1, 2, 3
# lm.chas = lm(crim~poly(chas,3)) : qualitative predictor
lm.nox = lm(crim~poly(nox,3))
summary(lm.nox) # 1, 2, 3
lm.rm = lm(crim~poly(rm,3))
summary(lm.rm) # 1, 2
lm.age = lm(crim~poly(age,3))
summary(lm.age) # 1, 2, 3
lm.dis = lm(crim~poly(dis,3))
summary(lm.dis) # 1, 2, 3
lm.rad = lm(crim~poly(rad,3))
summary(lm.rad) # 1, 2
lm.tax = lm(crim~poly(tax,3))
summary(lm.tax) # 1, 2
lm.ptratio = lm(crim~poly(ptratio,3))
summary(lm.ptratio) # 1, 2, 3
lm.black = lm(crim~poly(black,3))
summary(lm.black) # 1
lm.lstat = lm(crim~poly(lstat,3))
summary(lm.lstat) # 1, 2
lm.medv = lm(crim~poly(medv,3))
summary(lm.medv) # 1, 2, 3

参考资料

1. solutions