统计学习导论(ISLR)

文章目录

  • 统计学习导论(ISLR)
  • 5 交叉验证和bootstrapR语言代码实战
  • 5.1 验证集方法
  • 5.2 留一法交叉验证（LOOCV）
  • 5.3 K折交叉验证
  • 5.4 Bootstrap
  • 5.5 bootstrap进行区间估计

5 交叉验证和bootstrapR语言代码实战

我们将探讨本章介绍的重采样技术的r语言代码实现

5.1 验证集方法

我们使用验证集方法对Auto数据集拟合各种线性模型来测试测试误差
首先，我们使用set.seed()函数设置随机种子来使得我们代码的复现，在包含随机性元素的交叉验证分析时，我们最好设置随机种子。使用sample()函数将样本随机分为两部分

# 导入相关库
library(ISLR2)
set.seed(1)#设置随机种子
# 将观测数据随机平均分成两部分，注意下述代码返回的是索引值
train <- sample(392,196)

我们通过lm()函数中的subset选择子集索引，在这里我们直接令subset=train

attach(Auto)
lm.fit <- lm(mpg~horsepower, data =Auto, subset = train)

我们可以通过[-train]来得到测试数据集,下面计算测试数据集上的MSE

mean((mpg-predict(lm.fit, Auto))[-train]^2)

23.2660086465003

可以看出线性回归的MSE为23.27，下面我们使用poly()函数来估计不同二次模型和三次模型的MSE

lm.fit2 <- lm(mpg ~ poly(horsepower, 2), data=Auto,
             subset = train)
mean((mpg-predict(lm.fit2,Auto))[-train]^2)

18.7164594933828

lm.fit3 <- lm(mpg ~ poly( horsepower, 3), data = Auto,
             subset = train)
mean((mpg-predict(lm.fit3,Auto))[-train]^2)

18.7940067973945

这些MSE分别为18.72和18.79。如果我们选择一个不同的训练集，那么我们将在验证集上获得一些不同的错误。

从测试集的MSE来看，二次拟合会比简单的线性拟合效果好，三次拟合并没有明显的提升拟合效果

5.2 留一法交叉验证（LOOCV）

loocv的结果可以直接使用glm()函数和cv.glm()函数得到。在第四章中，我们使用glm()函数拟合了logistic回归，还记得我们使用了family参数=binomal,如果我们不指定family值，则其默认是线性回归，和lm()函数一样。我们下面来验证这一结果

glm.fit <- glm(mpg~horsepower, data = Auto)
coef(glm.fit)


39.9358610211705 horsepower -0.157844733353654 

lm.fit <- lm(mpg~horsepower, data = Auto)
coef(lm.fit)


39.9358610211705 horsepower -0.157844733353654

从结果可以看出上述两种方法得到的结果一样。下面我们使用glm()和cv.glm()来实现LOOCV，首先导入boot包

library(boot)
glm.fit <- glm(mpg~horsepower, data = Auto)
cv.err <- cv.glm(Auto,glm.fit)
cv.err$delta
  1. 24.2315135179293
  2. 24.2311440937562

可以看出上述得到的两个MSE几乎是一致的，因为LOOCV几乎每一次训练的结果差不多。注意，在使用cv.glm()时，不需要再指定train，会自动实现LOOCV过程。下面我们来看一下不同次方的多项式结果

cv.error <- rep(0,10)
for (i in 1:10) {
    glm.fit <- glm(mpg~poly(horsepower,i),data= Auto)
    cv.error[i] <- cv.glm(Auto, glm.fit)$delta[1]
}
cv.error
  1. 24.2315135179292
  2. 19.2482131244897
  3. 19.334984064029
  4. 19.4244303104303
  5. 19.0332138547041
  6. 18.9786436582254
  7. 18.8330450653182
  8. 18.9611507120532
  9. 19.0686299814599
  10. 19.490932299334

正如正文中的图5.4所示，在二次拟合的时候错误率较低较大

loocv <- data.frame(x,cv=cv.error)
library(ggplot2)
ggplot(loocv, aes(x, cv.error)) +geom_point() + geom_line(lwd=1,col='blue')

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-agKiOEwt-1648890761147)(output_22_0.png)]

5.3 K折交叉验证

cv.glm()函数也可以实现k折交叉验证，下面我们使用k=10，也是我们常用的一个k值。我们重新设置随机种子数，在Auto数据集拟合模型。

# 7是我的幸运数字，我就这么设置了
set.seed(7)
# 初始化向量来保存10折交叉验证的结果
cv_error_10 <- rep(0,10)
for (i in 1:10){
    glm.fit <- glm(mpg~poly(horsepower,i), data = Auto)
    cv_error_10[i] <- cv.glm(Auto,glm.fit,K = 10)$delta[1]
}
cv_error_10
  1. 24.1463716629577
  2. 19.3130825829741
  3. 19.434897545051
  4. 19.5493689322887
  5. 19.0736379228708
  6. 18.7058531603005
  7. 19.2522869995751
  8. 18.8552270777634
  9. 18.9304332711781
  10. 20.4425474405408

可以得出和loocv一直的结果，在二次拟合正确率提升较高，再更高次的拟合没有明显的提升模型正确率，还记得在LOOCV中$delta的两个MSE结果基本一致，但是在k-折交叉验证中，这两个数的结果略有区别，第一个是标准的k-cv,第二个是修正了偏差后的k-cv。在这个数据集中，两个值依然十分接近

glm.fit <- glm(mpg~poly(horsepower,i), data = Auto)
cv.glm(Auto,glm.fit,K = 10)$delta
  1. 19.0008700163861
  2. 18.9443399638009

5.4 Bootstrap

bootstrap方法的一大优点是，它几乎可以应用于所有情况。不需要复杂的数学计算。在R中只需要两个步骤。

  • 首先，我们必须创建一个函数来计算我们要研究的统计量。
  • 其次，我们使用boot()函数，通过重复采样数据集中的观察值来进行bootstrap

ISLR2包中的Portfolio数据集是100对回报的模拟数据，按照正文5.2所述的方式生成。为了说明如何在这些数据上使用bootstrap，我们必须首先创建一个函数alpha()，需要两个参数，一个是数据集，另一个参数是索引，指定使用哪些数据来估计fromjson r语言 r语言smooth spline_交叉验证。然后，该函数根据选定的数据输出fromjson r语言 r语言smooth spline_交叉验证的估计值。根据我们在正文中给出的fromjson r语言 r语言smooth spline_交叉验证估计表达式，可以编写函数如下

alpha <- function(data,index){
    X <- data$X[index]
    Y <- data$Y[index]
    (var(Y)-cov(X,Y))/(var(Y)+var(X)-2*cov(X,Y))
}

该函数输出给定数量数据集的fromjson r语言 r语言smooth spline_统计学习_04,例如我们把元素数据集全部输入

alpha(Portfolio,1:100)

0.57583207459283

上述给出了fromjson r语言 r语言smooth spline_交叉验证的估计值，下面我们先使用sample()函数，从原始数据集中随机取100个观测值，replace设置为T，这样能选择重复的数据，默认值为F，不能选择重复的数据。下面相当于做了一个简单的bootstrap

set.seed(7)
alpha(Portfolio, sample(100, 100, replace = T))

0.538532591946793

我们需要做多次bootstrap，得到多个fromjson r语言 r语言smooth spline_交叉验证的估计值，最后取平均值得到fromjson r语言 r语言smooth spline_交叉验证，boot函数自动为我们计算了这些步骤，需要制定R参数表示次数，一般要大于1000次

alpha_hat <- boot(Portfolio,alpha,R=1000)
alpha_hat


ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Portfolio, statistic = alpha, R = 1000)


Bootstrap Statistics :
     original     bias    std. error
t1* 0.5758321 0.00105432  0.09045912

上述结果表明原始数据的fromjson r语言 r语言smooth spline_统计学习_08,并且bootstrap估计fromjson r语言 r语言smooth spline_fromjson r语言_09

bootstrap方法可用于评估统计学习方法中系数估计和预测的稳定性。我们使用bootstrap来衡量fromjson r语言 r语言smooth spline_bootstrap_10和fromjson r语言 r语言smooth spline_bootstrap_11的稳定性，
首先定义函数得到回归系数

boot.fn <- function(data, index)
  {coef(lm(mpg ~ horsepower, data = data, subset = index))}
boot.fn(Auto, 1:392)


39.9358610211705 horsepower -0.157844733353654 

set.seed(1)
boot.fn(Auto, sample(392, 392, replace = T))


40.3404516830189 horsepower -0.163486837689938 

# 使用boot()函数计算经过1000次bootstrap的标准差
boot(Auto, boot.fn, 1000)


ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
      original        bias    std. error
t1* 39.9358610  0.0549915227 0.841925746
t2* -0.1578447 -0.0006210818 0.007348956

通过bootstrap得到的参数的标准差估计值分别为0.84和0.0073，下面我们使用summary()函数得到系数的具体信息

summary(lm(mpg ~ horsepower, data = Auto))$coef

A matrix: 2 × 4 of type dbl

Estimate

Std. Error

t value

Pr(>|t|)

(Intercept)

39.9358610

0.717498656

55.65984

1.220362e-187

horsepower

-0.1578447

0.006445501

-24.48914

7.031989e-81

从上述结果来看，两个系数的标准差分别为0.72和0.0064，与bootstrap的结果略有不同，但这并不意味着bootstrap的结果有很大的偏差，相反，我们之前在线性回归中做了很多假设，但实际情况这些假设往往不满足，而bootstrap不需要任何假设，因此从这个角度来看，bootstrap得到的结果比summary()函数更准确

下面我们在使用bootstrap计算二次模型结果如下

boot.fn <- function(data, index)
  coef(
      lm(mpg ~ horsepower + I(horsepower^2), data = data, subset = index)
    )
set.seed(1)
boot(Auto, boot.fn, 1000)


ORDINARY NONPARAMETRIC BOOTSTRAP


Call:
boot(data = Auto, statistic = boot.fn, R = 1000)


Bootstrap Statistics :
        original        bias     std. error
t1* 56.900099702  3.511640e-02 2.0300222526
t2* -0.466189630 -7.080834e-04 0.0324241984
t3*  0.001230536  2.840324e-06 0.0001172164
summary(
    lm(mpg ~ horsepower + I(horsepower^2), data = Auto)
  )$coef

A matrix: 3 × 4 of type dbl

Estimate

Std. Error

t value

Pr(>|t|)

(Intercept)

56.900099702

1.8004268063

31.60367

1.740911e-109

horsepower

-0.466189630

0.0311246171

-14.97816

2.289429e-40

I(horsepower^2)

0.001230536

0.0001220759

10.08009

2.196340e-21

5.5 bootstrap进行区间估计

我们可以根据bootstrap得到的统计量估计值和标准差估计值进一步计算其置信区间，可以使用boot.ci()函数计算期间

result <- boot(Portfolio,alpha,1000)
boot.ci(result)
Warning message in boot.ci(result):
"bootstrap variances needed for studentized intervals"



BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 1000 bootstrap replicates

CALL : 
boot.ci(boot.out = result)

Intervals : 
Level      Normal              Basic         
95%   ( 0.3888,  0.7527 )   ( 0.3826,  0.7430 )  

Level     Percentile            BCa          
95%   ( 0.4086,  0.7690 )   ( 0.4006,  0.7565 )  
Calculations and Intervals on Original Scale