《统计学习导论》R语言代码整理

这学期的统计课学习了《统计学习导论：基于R应用》，觉得这本书非常好，是机器学习也是R语言的入门教材。和纯统计理论的课比起来，学这本书更符合当下“学习”的热潮。

英文原版教材：
An Introduction to Statistical Learning: with Applications in R

电子版PDF下载、原始数据集和原书代码请看：
http://faculty.marshall.usc.edu/gareth-james/ISL/index.html

中文版：
《统计学习导论：基于R应用》机械工业出版社

在复习期末考期间，整理了一下这本书中关键的R代码。以后如果还要用R来做实证分析或者画图，也方便查找。

一、特殊函数

lm(y~x, data) #简单线性回归
  lm(y~x+0) #不含截距的简单线性回归
  lm(y~x1+x2+x3, data) #多元线性回归
  lm(wage ~ poly(age, 3), data) #多项式回归
  lm(wage ~ cut(age, 8), data) #阶梯函数拟合
  lm(wage ~ bs(age, knots/df), data) #（三次）回归样条，bs()-splines包
  lm(wage ~ ns(age, knots/df), data) #自然样条，ns()-splines包
  lm(wage ~ ns(age, 5) + education, data) #自然样条来拟合广义可加模型
glm(y~x, data, family=binomial, subset) #逻辑斯蒂回归
  #如果glm()没有设定family参数，和lm()一样执行线性回归
lda(y~x, data, subset) #线性判别分析-MASS库
qda(y~x, data, subset) #二次判别分析-MASS库
knn(train.X, test.X, train.Direction, k) #K最近邻法-class库

regsubset(y~., data, nvmax, method) #最优预测变量子集的选择-leaps库
glmnet(x, y, alpha = 0, lambda) #岭回归-glmnet库
glmnet(x, y, alpha = 1, lambda) #LASSO-glmnet库
pcr(y~., data, scale=T, validation="CV") #主成分回归-pls库
plsr(y~., data, scale=T, validation="CV") #偏最小二乘-pls库

bs() #默认生成三次样条，产生针对给定结点的所有样条基函数的矩阵– splines包
ns() #拟合自然样条 – splines包
    bs(age, knots=c(25,40,60)), data)  #给定结点
    bs(age, df=6), data)   #给定自由度
    #ns()同
smooth.spline(y, data, df) #光滑样条拟合(,cv=TRUE) -splines包
loess(y ~ x, span, data) #局部回归-splines包
s() #拟合光滑样条来拟合GAM - gam库
gam(y ~ x1 + x2 + s(age, 4), data) #光滑样条来拟合广义可加模型-gam包

svm(y~., data, kernel=”linear”,cost, scale=FALSE) #支持向量分类器-e1071库
tune(svm, y~., data, kernel=”linear”, ranges=list(cost=0.1,1,10)) #实现交叉验证（默认十折），比较不同cost参数下支持向量分类器的表现（交叉验证错误率）

tree(y~., data, subset) #建立分类树和回归树-tree包
randomForest(y~., data, subset, mtry, ntree) #装袋法、随机森林-randomForest包
gbm(y~., data, distribution = "gaussian"/ "bernoulli", n.trees, interaction.depth, shrinkage) #提升法-gbm包

kmeans(x, k, nstart) #K均值聚类
hclust(dist(x), method = "complete") #系统聚类
cutree(hc.complete, 3) #根据谱系图的切割获得各个观测的类标签

cv.glm(data, glm.fit) #留一/k折交叉验证法-boot库
cv.glmnet(X, Y, alpha) #用交叉验证选择岭回归/Lasso调节参数λ
cv.tree(oj.tree, FUN = prune.tree) #用交叉验证确定最优的树规模
boot(data, boot.fn, 1000) #自助法估计标准误差-boot库

二、基本函数

c() #“连接” 建立一个数值向量
length() #得到向量的长度
ls() #查看所有的对象列表
rm() #去除不想要的对象  
    rm(list=ls())#消除所有的对象
is.na() #识别有缺失值的观测，返回一个与输入向量等长的向量，元素TRUE表示该位置的向量缺失
na.omit() #剔除缺失数据行
list() #创建一个列表
data.frame()#创建数据框,由行和列组成,与matrix不同的是,每个列可以是不同的数据类型
matrix() #建立一个数值矩阵 
    matrix(data,nrow,ncol) #数据、行数、列数
dim() #输出一个矩阵的行数和列数 
    dim()[1] #输出行数
cor() #计算相关系数 cor(x,y) cor(subset())
dist() #计算欧氏距离矩阵
median() #求中位数
subset() #筛选子集 
    subset(data,select=-x) #排除定性变量x
rownames() #批量为矩阵的行命名
colnames() #批量为矩阵的列命名
    colnames(x, do.NULL = FALSE, prefix = "x.") #把列命名为x.1 x.2 x.3
scale() #对变量进行标准化处理
rnorm() #生成正态分布随机数 
    rnorm(n,mean,sd.) #样本容量、均值、标准差
runif() #生成均匀分布随机数
    runif(n,min,max) #样本容量、最小值、最大值
seq() #生成一个规则序列 
    seq(1,500) #生成1-500的规则序列
poly(x,3，raw=T) #返回一个矩阵，其每一列都是正交多项式的基，即x,x^2,x^3
    #raw=TRUE,用poly函数直接估计x,x^2,x^3
hatvalues() #计算杠杆统计量 hatvalues(lm.fit)
which.max() #识别出向量中最大元素的索引
which.min() #识别出向量中最大元素的索引
range() #输出两个值:最小值和最大值
confint() #得到系数估计值的95%置信区间
sample() #随机选取样本组成子集 
    sample(200,100) #从200个观测中随机选取100个
coefficients()/coef() #获取拟合模型的系数  coef(glm.fit) coef(mod.fwd, id = 3)
predict() #预测概率 
    predict(glm.fit, type = "response") #输出概率P(Y=1|X)
    predict(oj.tree, OJ.test, type = "class")
    predict(ridge.mod, s = 50, type = "coefficients")
rep("Down", length) #创建一个向量，由length个"Down"组成
table() #产生混淆矩阵来判断有多少观测被正确/错误分类了 
    table(glm.fit, Direction) #拟合的预测值，真实值
mean(glm.fit==Direction) #计算正确预测的比例
contrasts() #创造一个哑变量，只有0-1
cbind() #列合并 
    cbind(x,y,z)[train,] #将xyz三个变量train部分合并成一个矩阵
as.matrix() #把一个vector或者data frame变成矩阵
    #as.matrix一般是将整个数据框(各列数据类型必须相同，否则进行强制转换)转换为维数相同的矩阵，目的是利用那些矩阵运算的函数。cbind可以按列将若干部分合并成矩阵或数据框，目的一般是为了整合数据。
model.matrix() #生成一个与n个预测变量相对应的矩阵，并自动将定性变量转化为哑变量
as.dist() #将任意一个对称方阵转换成距离矩阵的形式
apply(data, 1/2, mean/var) #对数据集的每一行或列使用同一函数，1表示行，2表示列
anova() #实现方差分析,前一个是后一个的子集 anova(fit.1, fit.2)
cut() #自动对变量进行分割点的选择 
    cut(age,4) #将age变量分割成四段，产生三个分割点

三、画图

一些函数

par(mfrow=c(2,2)) #将绘图区域分成2*2网格面板
par(mar,oma) #调整图的边界
plot(lm.fit) #生成诊断图（4个）
plot(x,y) #绘制散点图，如果x是定性变量，自动绘制箱线图
boxplot(x,y) #绘制箱线图
barplot() #绘制柱状图
abline(lm.fit) #绘制最小二乘回归直线
abline(a,b) #画一条截距为a，斜率为b的直线
pairs()  #作出数据集中所有变量的散点图矩阵
points() #在一张图上添加点，与plot()配合用
lines() #在一张图上添加线，与plot()配合用
hist() #绘制直方图
title() #产生整个图的标题（非子图）

一些参数

ylim = c(0, 100)  #y轴显示的范围
cex = 0.8    #绘图符号相对于默认大小的缩放倍数
lwd = 7      #线条宽度相对于默认宽度的倍数
col = "red"  #绘图符号颜色
type= "p"       #在图形中数据显示为点（详见下方）
pch = 4         #数据点类型为×（详见下方）
lty = "dashed"  #线条类型为虚线
lty = 2         #线条类型为虚线（详见下方）

type

type= "l" #在图形中数据显示为线
type= "b" #在图形中数据显示为点和连接线
type= "o" #在图形中数据点覆盖在线上
type= "h" #在图形中数据显示为从点到x轴的垂直线
type= "s" #在图形中数据显示为阶梯图

pch （plotting character）

lty（line types）

特定问题里的画图

回归的残差分析

plot(predict(lm.fit), residuals(lm.fit))  #生成残差图
plot(predict(lm.fit), rstudent(lm.fit))  #生成学生化残差图

岭回归/lasso交叉验证

mod.lasso = cv.glmnet(Xmatrix, Y, alpha = 1) 
plot(mod.lasso)  #lasso交叉验证选择参数λ，作出交叉验证误差MSE~Log(λ)的图像

主成分分析交叉验证

validationplot(pcr.fit, val.type="MSEP")#作每个可能的主成分个数M所对应的的交叉验证得分图像

广义可加模型拟合图

gam1 = lm(wage ~ ns(age, 5) + education, data) #自然样条来拟合广义可加模型
plot.gam(gam1, se=TRUE, col="red") #画拟合结果 
gam3 = gam(wage ~ s(age, 5)+ education, data) #光滑样条来拟合广义可加模型-gam包 
plot(gam3, se=TRUE, col="red") #画拟合结果

plot(gam.fit, se = T, col = "blue") #画GAM拟合图（有几项就有几个图，p197图7-11）

决策树结构图

plot(tree.carseats) #画决策树的结构
text(tree.carseats, pretty = 0, cex = 0.7) #显示决策树上结点标记

聚类谱系图

hc.complete = hclust(dist(USArrests), method="complete")  #用最长距离法进行系统聚类
plot(hc.complete) #绘出聚类谱系图

四、几类问题的关键代码

第3章是回归问题，非常简单，应用前面列举的特殊函数即可完成。
下面从第4章开始贴一下几个重要方法的关键代码（来自习题）

第4章 分类问题

Logistis Regression - 逻辑斯蒂回归- p112

glm.fit = glm(Direction ~ Lag2, data = Weekly, family = binomial, subset = train) 
glm.probs = predict(glm.fit, Weekly.0910, type = "response") 
glm.pred = rep("Down", length(glm.probs))
glm.pred[glm.probs > 0.5] = "Up"
table(glm.pred, Direction.0910)
mean(glm.pred == Direction.0910)

LDA/QDA -线性判别分析/二次判别分析- p113

library(MASS)
lda.fit = lda(Direction ~ Lag2, data = Weekly, subset = train)
lda.pred = predict(lda.fit, Weekly.0910)
lda.class = lda.pred$class        #取预测的列表中第一个元素class
table(lda.class, Direction.0910)
mean(lda.class == Direction.0910)

KNN - K邻近法- p114

library(class)
train.X = as.matrix(Lag2[train])#与训练数据相关的预测变量矩阵,TRUE部分
test.X = as.matrix(Lag2[!train])#与测试数据相关的预测变量矩阵,FALSE部分
train.Direction = Direction[train]#训练观测类标签的向量
set.seed(1)  #设置一个随机种子
knn.pred = knn(train.X, test.X, train.Direction, k = 1)  
table(knn.pred, Direction.0910)
mean(knn.pred == Direction.0910)

第5章 重抽样方法

LOOCV - 留一交叉验证法– p133

library(boot)
Data = data.frame(x, y)
set.seed(1)
glm.fit = glm(y ~ x)
cv.err=cv.glm(Data, glm.fit)
cv.err$delta   
#cv.glm()会生成一个有几个组成部分的列表，其中一个组成部分delta向量中两个数字为交叉验证的结果

k-fold CV（k折交叉验证法）与LOOCV相似，区别仅在于

cv.err=cv.glm(Data, glm.fit, K=10)  #十折交叉验证

bootstrap - 自助法 – p135

boot.fn = function(data, index) 
  + return(coef(glm(default ~ income + balance, data = data, family = binomial, subset = index)))  #创建boot.fn函数，返回截距和斜率的估计
library(boot)
boot(Default, boot.fn, 50)   #产生50个系数的自助法估计的标准误差

第6章 线性模型选择与正则化

最优子集选择 – p168

library(leaps)
regfit.full = regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "forward")
regfit.full = regsubsets(Salary ~ ., data = Hitters, nvmax = 19, method = "backward") 
#nvmax设置用户所需的预测变量个数，method设置向前/向后选择

Ridge regression岭回归 / LASSO – p173-176

library(glmnet)
set.seed(1)
x = model.matrix(Salary ~ ., data = Hitters.train) #构造x矩阵
y = Hitters.train$Salary #构造y向量
x.test = model.matrix(Salary ~ ., data = Hitters.test)
#best.lambda = mod.lasso$lambda.min #选择使得交叉验证误差最小的λ（lambda.min）
#ridge.fit = glmnet(x, y, alpha = 0)  #岭回归
lasso.fit = glmnet(x, y, alpha = 1)  #Lasso
#glmnet()函数在使用时，必须要输入一个x矩阵和一个y向量，但不会用到y~x这个句法
lasso.pred = predict(lasso.fit, s = 0.01, newx = x.test)
mean((Hitters.test$Salary - lasso.pred)^2)

PCR - 主成分回归 – p177

library(pls)
pcr.fit = pcr(Apps~., data = College.train, scale = TRUE, validation = "CV") 
summary(pls.fit)
  # scale=TRUE：在生成主成分之前标准化每个预测变量
  # validation="CV" 使pcr()函数使用十折交叉验证计算每个可能的主成分个数M所对应的的交叉验证误差
validationplot(pcr.fit, val.type="MSEP")  #作出每个可能的主成分个数M所对应的的交叉验证得分图像

pcr.pred = predict(pcr.fit, College.test, ncomp = 10)   #主成分个数M=10，看adjCV极小值
mean((College.test[, "Apps"] - data.frame(pcr.pred))^2)

PLS - 偏最小二乘 – p179

library(pls)
pls.fit = plsr(Apps~., data=College.train, scale=T, validation="CV")
summary(pls.fit)
validationplot(pls.fit, val.type="MSEP")
pls.pred = predict(pls.fit, College.test, ncomp = 11)
mean((College.test[, "Apps"] - pls.pred)^2)

第7章 非线性模型

Polynomial regression - 多项式回归 – p200

lm.fit = lm(wage ~ poly(age, 3), data = Wage)   #3次多项式回归
agelims = range(Wage$age)  #范围
age.grid = seq(from = agelims[1], to = agelims[2])
lm.pred = predict(lm.fit, data.frame(age = age.grid))  #预测
plot(wage ~ age, data = Wage, col="darkgrey") #画点
lines(age.grid, lm.pred, col="blue", lwd=2) #画线

Step function - 阶梯函数 – p203

lm.fit = lm(wage ~ cut(age, 8), data = Wage)
agelims = range(Wage$age)
age.grid = seq(from = agelims[1], to = agelims[2])
lm.pred = predict(lm.fit, data.frame(age = age.grid))
plot(wage~age, data=Wage, col="darkgrey")
lines(age.grid, lm.pred, col="red", lwd=2)

Spline - 回归样条 – p203

library(splines)
attr(bs(dis, df = 4), "knots") #选择结点 50%-3.20745 
sp.fit = lm(nox ~ bs(dis, df = 4, knots = 3.20745), data = Boston)
summary(sp.fit)
#dislim = range(dis)
#dis.grid = seq(from = dislim[1], to = dislim[2], by = 0.1)  #x轴
sp.pred = predict(sp.fit, list(dis = dis.grid))  #y轴
plot(nox ~ dis, data = Boston, col = "darkgrey")
lines(dis.grid, sp.pred, col = "red", lwd = 2)

GAM - 广义可加模型 – p205

library(gam)
fit = gam(wage ~ maritl + jobclass + s(age, 4), data = Wage)
deviance(fit)

第9章 支持向量机

SVM - 支持向量机 - p248

library(e1071)

支持向量分类器：线性核函数

set.seed(1)
tune.out = tune(svm, mpglevel ~ ., data = Auto, kernel = "linear", ranges = list(cost = c(0.01,0.1, 1, 5, 10, 100)))
summary(tune.out) #求交叉验证错误率

支持向量机：非线性核函数-多项式核函数

set.seed(2)
tune.out = tune(svm, mpglevel ~ .,data = Auto, kernel = "polynomial", ranges = list(cost = c(0.1, 1, 5, 10), degree = c(2, 3, 4)))
summary(tune.out)

支持向量机：非线性核函数-径向核函数

set.seed(3)
tune.out = tune(svm, mpglevel ~ ., data = Auto, kernel = "radial", ranges = list(cost = c(0.1,1, 5, 10), gamma = c(0.01, 0.1, 1, 5, 10, 100)))
summary(tune.out)

第8章 基于树的方法

Tree - 决策树– p225

library(tree)
tree.carseats = tree(Sales ~ ., data = Carseats.train)

#回归树 p227

pred.carseats = predict(tree.carseats, Carseats.test)
mean((Carseats.test$Sales - pred.carseats)^2) #测试均方误差
#用交叉验证确定最优的树复杂性
cv.carseats = cv.tree(tree.carseats, FUN = prune.tree)
#最优的树复杂性时，终端节点数size = 6
pruned.carseats = prune.tree(tree.carseats, best = 6) #剪枝

#分类树 p225

oj.pred = predict(oj.tree, OJ.test, type = "class")
table(OJ.test$Purchase, oj.pred)
mean(OJ.test$Purchase!=oj.pred)

RandomForest / Bagging - 随机森林/装袋法 – p228

#回归树

library(randomForest)
rf.carseats = randomForest(Sales ~ ., data = Carseats.train, mtry = 5, ntree = 500, importance = T)
rf.pred = predict(rf.carseats, Carseats.test)
mean((Carseats.test$Sales - rf.pred)^2) #测试MSE
importance(rf.carseats) #看重要指标

Boosting - 提升法 – p230

library(gbm)

#回归树

boost.hitters = gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas)
summary(boost.hitters)
train.pred = predict(boost.hitters, Hitters.train, n.trees = 1000)
test.pred = predict(boost.hitters, Hitters.test, n.trees = 1000)
train.errors  = mean((Hitters.train$Salary - train.pred)^2) #训练错误率
test.errors = mean((Hitters.test$Salary - test.pred)^2) #测试错误率

#分类树

boost.caravan = gbm(Purchase ~ ., data = Caravan.train, distribution = "bernoulli", n.trees = 1000, shrinkage = 0.01)
boost.prob = predict(boost.caravan, Caravan.test, n.trees = 1000, type = "response")
boost.pred = ifelse(boost.prob > 0.2, 1, 0)
table(Caravan.test$Purchase, boost.pred)

第10章 无指导学习

K-means clustering - K均值聚类 – p280

km.out = kmeans(x, 3, nstart=20)
km.out$cluster  #输出分类情况
table(km.out$cluster, c(rep(1,20), rep(2,20), rep(3,20)))

Hierarchical clustering - 系统聚类 – p281

hc.complete = hclust(dist(USArrests), method="complete")  #用最长距离法进行系统聚类
plot(hc.complete) #画谱系图
cutree(hc.complete, 3) #根据谱系图的切割，将数据分成3个不同的类
table(cutree(hc.complete, 3)) #看每一类有几个“州”

附：全书思维导图（来自网络）