Задание 1.4. Написать на R функцию stepCV, аналогичную stepAIC, но использующую кросс-валидацию для проверки значимости признака. Нужно следовать при этом принципу иерархии — нельзя выкидывать признаки более низкого порядка, если есть признаки более высокого. Функция должна работать со всеми методами, с которыми работает stepAIC.
Hint: Можно расковырять оригинальную функцию stepAIC, а кросс-валидацию взять из e1071.
library(MASS)
library(e1071)
source("stepCV.R")
Пример посложнее
Advertising <- read.csv("data/Advertising.csv")
l <- lm(y~. , epil)
stepCV(l, trace = TRUE)
## Start: CV.performance=61.23
## y ~ trt + base + age + V4 + subject + period + lbase + lage
##
## Df rank CV
## - base -42.25576 8 103.745
## - lbase -12.50889 8 73.998
## - trt -5.35505 8 66.845
## - period -3.42826 8 64.918
## - V4 -2.32123 8 63.811
## - subject -2.09985 8 63.589
## - age -0.85821 8 62.348
## - lage -0.79301 8 62.283
## <none> 9 61.490
##
## Step: CV.performance=100.78
## y ~ trt + age + V4 + subject + period + lbase + lage
##
## Df rank CV
## - lbase -56.66179 7 159.12
## - period -3.06349 7 105.53
## - subject -1.41415 7 103.88
## - lage -0.87568 7 103.34
## <none> 8 102.46
## - age 1.15964 7 101.30
## - V4 1.23309 7 101.23
## - trt 1.83984 7 100.62
##
## Step: CV.performance=156.77
## y ~ trt + age + V4 + subject + period + lage
##
## Df rank CV
## - age -5.68037 6 162.75
## - trt -3.31845 6 160.39
## - period -2.34772 6 159.42
## <none> 7 157.07
## - subject 0.17981 6 156.89
## - V4 0.56260 6 156.51
## - lage 3.13741 6 153.94
##
## Step: CV.performance=157.49
## y ~ trt + V4 + subject + period + lage
##
## Df rank CV
## - trt -3.569553 5 160.06
## - lage -2.899836 5 159.39
## - subject -0.202992 5 156.69
## - V4 -0.015523 5 156.50
## <none> 6 156.49
## - period 1.648713 5 154.84
##
## Step: CV.performance=158.62
## y ~ V4 + subject + period + lage
##
## Df rank CV
## <none> 5 160.13
## - subject 2.8453 4 157.29
## - period 3.6510 4 156.48
## - V4 4.6237 4 155.51
## - lage 5.0908 4 155.04
##
## Call:
## lm(formula = y ~ V4 + subject + period + lage, data = epil)
##
## Coefficients:
## (Intercept) V4 subject period lage
## 8.827070 -0.723164 0.009531 -0.271186 -2.825520
Можно сравнить результат с stepAIC:
stepAIC(l, trace = TRUE)
## Start: AIC=963.38
## y ~ trt + base + age + V4 + subject + period + lbase + lage
##
## Df Sum of Sq RSS AIC
## - period 1 8.7 12969 961.53
## - V4 1 9.3 12970 961.54
## - age 1 67.3 13028 962.60
## - trt 1 103.7 13064 963.26
## - lage 1 104.4 13065 963.27
## <none> 12960 963.38
## - subject 1 132.8 13093 963.78
## - lbase 1 1376.2 14337 985.19
## - base 1 9300.0 22260 1089.03
##
## Step: AIC=961.53
## y ~ trt + base + age + V4 + subject + lbase + lage
##
## Df Sum of Sq RSS AIC
## - age 1 67.3 13036 960.76
## - V4 1 70.9 13040 960.82
## - trt 1 103.7 13073 961.41
## - lage 1 104.4 13074 961.43
## <none> 12969 961.53
## - subject 1 132.8 13102 961.94
## - lbase 1 1376.2 14345 983.34
## - base 1 9300.0 22269 1087.12
##
## Step: AIC=960.76
## y ~ trt + base + V4 + subject + lbase + lage
##
## Df Sum of Sq RSS AIC
## - V4 1 70.9 13107 960.04
## <none> 13036 960.76
## - trt 1 111.8 13148 960.77
## - subject 1 130.0 13166 961.10
## - lage 1 329.4 13366 964.64
## - lbase 1 1431.0 14468 983.34
## - base 1 9368.6 22405 1086.56
##
## Step: AIC=960.04
## y ~ trt + base + subject + lbase + lage
##
## Df Sum of Sq RSS AIC
## <none> 13107 960.04
## - trt 1 111.8 13219 960.04
## - subject 1 130.0 13237 960.37
## - lage 1 329.4 13437 963.89
## - lbase 1 1431.0 14538 982.49
## - base 1 9368.6 22476 1085.30
##
## Call:
## lm(formula = y ~ trt + base + subject + lbase + lage, data = epil)
##
## Coefficients:
## (Intercept) trtprogabide base subject lbase
## -10.2977 -2.7702 0.5570 0.0873 -7.7688
## lage
## 5.4846
Результаты получаются разные (Sales~X + Newspaper при вызове stepCV, при вызове stepAIC - другой результат)