Введение

Подключим библиотеки, которые будем использовать.

library(knitr)
library(car)
library(ggplot2)
library(leaps)

Файл с данными содержит \(22\) зависимые переменные, \(10\) предикторов. Тренировочная выборка содержит \(100\) объектов, также как и тестовая. Требуется определить наилучшую модель связи между десятой зависимой переменной и набором предикторов, а также спрогнозировать значения зависимой переменной по заданным значениям предикторов. По условиям задачи связь может быть четырёх типов:

  1. Линейная зависимость от нескольких предикторов: \[ y \sim \beta_{i_{n}} + \beta_{i_{1}} x_{i_{1}} + \beta_{i_{2}} x_{i_{2}} + \ldots + \beta_{i_{n}} x_{i_{n}}; \]
  2. Квадратичная зависимость от одного из предикторов: \[ y \sim \beta_{0} + \beta_{1}x_{i} + \beta_{2}x_{i}^2; \]
  3. Линейная зависимость от пары предикторов и их произведения: \[ y \sim \beta_{0} \beta_{1}x_{i} + \beta_{2}x_{j} + \beta_{3}x_{i} * x_{j}; \]
  4. Линейная зависимость логарифмов зависимой и независимой переменной: \[ ln(y) \sim ln(x_{i} + \beta_{0}). \]

Работа с данными

Считаем данные в переменную data. Затем разделим ее на \(3\) части:

data <- read.csv("~/Projects/r_prak/tasks/task_1/data.csv", header = TRUE)
cols = names(data)
train <- data[1:100, grepl("X", names(data))]
train$Y <- data$Y21[1:100]
test <- data[101:200, grepl("X", names(data))]
kable(head(train,9))
X1 X2 X3 X4 X5 X6 X7 X8 X9 X10 Y
5.543269 10.315853 female 17.520172 2.1174891 5.2764380 red 9.572415 25.802848 14.6671107 120.10926
-2.802719 8.190485 female 6.568124 1.6576938 1.0044001 blue 4.276859 3.319735 0.8894010 48.51310
17.751634 4.491375 male 29.118749 6.3340438 8.0580712 green 1.616824 -17.464989 -19.9837341 127.19103
1.873201 11.059216 female 2.901179 -9.5299775 12.2856157 red -9.481376 3.294523 14.1826723 -33.98842
11.425261 12.593616 male 19.837411 -5.9173843 0.6525432 red 16.705655 25.821052 10.6493785 165.05452
4.155261 6.061951 male 9.705372 -0.0105406 4.2573849 blue 2.070889 4.192606 1.5223514 50.72742
12.295066 14.083780 female 24.789746 -0.4408137 3.4601678 red -7.341028 -8.137743 0.2848794 64.00656
2.366797 16.957718 male 21.579876 0.6623643 7.3963745 green -20.433548 -8.949928 11.4883807 -13.02837
-3.653828 13.182149 male 6.976611 -3.8977056 -14.3627654 red 13.153703 11.983920 -1.5238825 95.42675

Посмотрим возможные зависимости между предиктором и зависимой переменной:

for (col in names(train)[-length(names(train))]) {
  tmp = data.frame(x = train[col], Y = train["Y"])
  colnames(tmp) <- c("X","Y")
  if ((col == "X3") || (col == "X1")) {
    print( qplot ( data = tmp, x = X, y = Y, geom=c("boxplot") ) + xlab(col) )
  } else {
    print( qplot ( data = tmp, x = X, y = Y, geom=c("smooth", "point") ) + xlab(col) )
  }
}
## Warning: Continuous x aesthetic -- did you forget aes(group=...)?

Поиск наилучшей модели в каждой категории

1. Линейная зависимость от нескольких предикторов

Используем функцию vif, чтобы избавиться от коррелированных предикторов.

vif(lm(Y~.,data=train))
##           GVIF Df GVIF^(1/(2*Df))
## X1   89.443154  1        9.457439
## X2   47.963672  1        6.925581
## X3    1.374745  1        1.172495
## X4  217.149243  1       14.735985
## X5   95.547202  1        9.774825
## X6    1.214118  1        1.101870
## X7    1.226791  2        1.052429
## X8  107.158424  1       10.351735
## X9  187.277476  1       13.684936
## X10  95.218716  1        9.758008

Так как \(X4\) имеет наибольшее значение, исключим его.

vif(lm(Y~.-X4,data=train))
##           GVIF Df GVIF^(1/(2*Df))
## X1    1.051438  1        1.025396
## X2    1.161395  1        1.077680
## X3    1.373970  1        1.172165
## X5    1.184769  1        1.088471
## X6    1.213296  1        1.101497
## X7    1.199761  2        1.046583
## X8  101.972521  1       10.098144
## X9  177.981588  1       13.340974
## X10  90.247435  1        9.499865

Удалим по той же причине предиктор \(X9\).

vif(lm(Y~.-X4-X9,data=train))
##         GVIF Df GVIF^(1/(2*Df))
## X1  1.049701  1        1.024549
## X2  1.160814  1        1.077411
## X3  1.367172  1        1.169261
## X5  1.171612  1        1.082410
## X6  1.212083  1        1.100947
## X7  1.159672  2        1.037729
## X8  1.065060  1        1.032017
## X10 1.180842  1        1.086665

Теперь мы имеем линейную регрессию с исключенным фактором коллинеарности. Сохранним соотвествующюю матрицу признаков в uncorrelated:

uncorrelated <- subset(train, select = -c(X4, X9))
summary(uncorrelated)
##        X1                 X2                 X3           X5          
##  Min.   :-20.2704   Min.   : 0.008959   female:49   Min.   :-24.0382  
##  1st Qu.: -7.5623   1st Qu.: 3.742684   male  :51   1st Qu.: -5.7687  
##  Median : -0.7927   Median : 6.639532               Median :  1.1429  
##  Mean   : -0.6183   Mean   : 8.901576               Mean   :  0.3478  
##  3rd Qu.:  5.5737   3rd Qu.:12.932392               3rd Qu.:  6.4084  
##  Max.   : 18.3616   Max.   :29.056886               Max.   : 23.7103  
##        X6               X7           X8                X10          
##  Min.   :-18.1154   blue :21   Min.   :-24.9308   Min.   :-25.8075  
##  1st Qu.: -0.9153   green:54   1st Qu.: -6.8219   1st Qu.: -5.1914  
##  Median :  3.7625   red  :25   Median : -0.1508   Median :  0.3031  
##  Mean   :  2.7234              Mean   : -0.3471   Mean   :  1.0193  
##  3rd Qu.:  7.9727              3rd Qu.:  7.6276   3rd Qu.:  7.8190  
##  Max.   : 17.6291              Max.   : 18.2632   Max.   : 24.4446  
##        Y           
##  Min.   :-169.132  
##  1st Qu.:  -2.437  
##  Median :  34.360  
##  Mean   :  34.847  
##  3rd Qu.:  82.217  
##  Max.   : 220.758

Теперь попытаемся рассмотрим всевозможные наборы предикторов и отберем лучшие.

rsb <- regsubsets(Y~., data = uncorrelated)
summary(rsb)
## Subset selection object
## Call: regsubsets.formula(Y ~ ., data = uncorrelated)
## 9 Variables  (and intercept)
##         Forced in Forced out
## X1          FALSE      FALSE
## X2          FALSE      FALSE
## X3male      FALSE      FALSE
## X5          FALSE      FALSE
## X6          FALSE      FALSE
## X7green     FALSE      FALSE
## X7red       FALSE      FALSE
## X8          FALSE      FALSE
## X10         FALSE      FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
##          X1  X2  X3male X5  X6  X7green X7red X8  X10
## 1  ( 1 ) " " " " " "    " " " " " "     " "   "*" " "
## 2  ( 1 ) "*" " " " "    " " " " " "     " "   "*" " "
## 3  ( 1 ) "*" " " " "    "*" " " " "     " "   "*" " "
## 4  ( 1 ) "*" "*" " "    "*" " " " "     " "   "*" " "
## 5  ( 1 ) "*" "*" " "    "*" " " " "     "*"   "*" " "
## 6  ( 1 ) "*" "*" " "    "*" " " " "     "*"   "*" "*"
## 7  ( 1 ) "*" "*" "*"    "*" " " " "     "*"   "*" "*"
## 8  ( 1 ) "*" "*" "*"    "*" "*" " "     "*"   "*" "*"
plot(rsb)

Закрашенные клетки отвечают использованным в регрессии предикторам. Соответственно наилучший результат(BIC) имеет регрессия с использованием предикторов \(X3\) и \(X8\): \[ y \sim \beta_{0} + \beta_{8} x_{8} + \beta_{5} x_{5}. \] Пострим данную регресиию:

best_model_1 <- lm(Y ~ X8+X5,data=train)
best_bic_1 <- BIC(best_model_1)

2. Квадратичная зависимость от одного из предикторов

Из предыдущих исследований мы можем предположить, что три наиболее информативных предиктора: \(X1, X8, X5\). Соответственно, мы будем рассматривать три возможных варианта регрессии:

  1. \(y \sim \beta_{0} + \beta_{1} x_{1} + {\beta_{1} x_{1}}^2;\)
  2. \(y \sim \beta_{0} + \beta_{5} x_{5} + {\beta_{5} x_{5}}^2;\)
  3. \(y \sim \beta_{0} + \beta_{8} x_{8} + {\beta_{8} x_{8}}^2.\)

Посчитаем BIC для каждой модели и построим график:

bics <- BIC(lm(Y ~ X1 + X1 * X1, data = train),
             lm(Y ~ X5 + X5 * X5, data = train),
             lm(Y ~ X8 + X8 * X8, data = train)
            )
qplot(1:3, bics$BIC, geom = c("path","point")) + xlab("Kind of regression") + ylab("BIC")

Как видим, наилучшая модель - третья. Построим данную регрессию:

best_model_2 <- lm(Y ~ X8 + X8 * X8,data=train)
best_bic_2 <- BIC(best_model_2)

3. Линейная зависимость от пары предикторов и их произведения

По аналогии с предыдущим рассуждением, мы будем рассматривать пары из трех наиболее информативных предикторов \(X1, X8, X5\):

  1. \(y \sim \beta_{0} + \beta_{1} x_{1} + \beta_{5} x_{5} + \beta_{1} x_{1} \beta_{5} x_{5};\)
  2. \(y \sim \beta_{0} + \beta_{1} x_{1} + \beta_{8} x_{8} + \beta_{1} x_{1} \beta_{8} x_{8};\)
  3. \(y \sim \beta_{0} + \beta_{5} x_{5} + \beta_{8} x_{8} + \beta_{5} x_{5} \beta_{8} x_{8}.\)

Посчитаем BIC для каждой модели и построим график:

bics <- BIC(lm(Y ~ X1 + X5 + X1 * X5, data = train),
             lm(Y ~ X1 + X8 + X1 * X8, data = train),
             lm(Y ~ X10 + X8 + X5 * X8, data = train)
            )
qplot(1:3, bics$BIC, geom = c("point","path")) + xlab("Kind of regression") + ylab("BIC")

Как видим, наилучшая модель - вторая. Построим данную регрессию:

best_model_3 <- lm(Y ~ X1 + X8 + X1 * X8, data = train)
best_bic_3 <- BIC(best_model_3)

Итоговый выбор модели

Построим график BIC для всех моделей, лучших для своих типов, чтобы выбрать итоговую модель:

qplot(1:3, c(best_bic_1, best_bic_2, best_bic_3), geom = c("point","path")) + xlab("Kind of regression") + ylab("BIC")

print(c(best_bic_1, best_bic_2, best_bic_3))
## [1] 1078.837 1111.929 1076.404
best_model <- best_model_1
best_bic <- best_bic_1

Таким образом наилучшая модель это: \[ y \sim \beta_{0} + \beta_{8} x_{8} + \beta_{5} x_{5}. \] Рассмотрим различные характеристики нашей лучшей модели.

summary(best_model)
## 
## Call:
## lm(formula = Y ~ X8 + X5, data = train)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -100.463  -32.203   -1.962   38.999   97.881 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   35.239      4.938   7.137 1.74e-10 ***
## X8             4.566      0.522   8.747 6.77e-14 ***
## X5             3.432      0.515   6.664 1.62e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 49.32 on 97 degrees of freedom
## Multiple R-squared:  0.5176, Adjusted R-squared:  0.5076 
## F-statistic: 52.04 on 2 and 97 DF,  p-value: 4.423e-16
plot(best_model)

Таким образом мы получили весьма неплохую модель с отличными показателями.

Построение доверительного и предиктивного интервалов

Построим \(95\%\)-доверительный интервал для нашей модели и построим график для наглядности.

predicted <- as.data.frame(predict(best_model, test, level = 0.95, interval = "confidence"))
predicted$num <- c(1:100)
ggplot(predicted, aes(num)) + geom_line(aes(y=fit), colour="darkgreen") +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, fill = "red")

На самом деле интервал построен, однако его не видно, так как он весьма мал. Вот непостредственно предсказание нашей модели:

kable(predicted[,names(predicted)!="x"])
fit lwr upr num
101 73.2246094 59.5812228 86.8679960 1
102 -66.3767498 -89.9018777 -42.8516219 2
103 34.8662974 25.0117524 44.7208425 3
104 84.2769385 63.6444333 104.9094437 4
105 13.8495957 -2.5748386 30.2740300 5
106 103.9651052 85.9823090 121.9479013 6
107 36.7540458 22.1751629 51.3329287 7
108 52.4788915 36.2866966 68.6710863 8
109 21.8494572 3.8429747 39.8559397 9
110 -81.6611350 -106.4261183 -56.8961516 10
111 -0.7435148 -16.8494387 15.3624090 11
112 74.6633723 59.3482104 89.9785343 12
113 49.6406095 31.5390459 67.7421732 13
114 75.4550836 60.8038817 90.1062854 14
115 85.3753387 70.9136038 99.8370737 15
116 -7.8361368 -20.6802231 5.0079495 16
117 -44.8171478 -63.1623619 -26.4719337 17
118 77.7994164 64.9177442 90.6810885 18
119 76.8089669 59.9723163 93.6456175 19
120 43.6353941 33.5324723 53.7383159 20
121 -26.1100809 -42.5136469 -9.7065148 21
122 -15.3534170 -30.2089212 -0.4979128 22
123 146.2905839 121.8880253 170.6931426 23
124 21.4482734 10.8219152 32.0746316 24
125 40.0886270 14.8560332 65.3212209 25
126 43.4571381 29.6099595 57.3043167 26
127 -73.9812940 -97.5849134 -50.3776745 27
128 24.5812007 14.2813730 34.8810284 28
129 99.2372609 82.9742885 115.5002332 29
130 -80.5950684 -106.1432175 -55.0469193 30
131 78.2638789 60.9684411 95.5593167 31
132 73.9735485 54.1103110 93.8367860 32
133 6.1315285 -6.2136911 18.4767481 33
134 31.5735593 13.2226997 49.9244188 34
135 66.6091085 51.8991384 81.3190787 35
136 26.1563797 11.4949073 40.8178520 36
137 -31.1096615 -47.2548282 -14.9644948 37
138 65.7735413 49.6118592 81.9352235 38
139 23.2357127 12.7027113 33.7687142 39
140 115.7431373 94.9528279 136.5334467 40
141 -104.6038044 -133.5522082 -75.6554005 41
142 -33.4943947 -51.4719865 -15.5168029 42
143 45.1414126 32.7181076 57.5647176 43
144 110.9877164 91.3934868 130.5819460 44
145 78.6240027 64.2270006 93.0210049 45
146 89.7231120 73.9040030 105.5422209 46
147 16.6000070 -0.4431253 33.6431393 47
148 25.4527677 14.9505798 35.9549557 48
149 -27.0708690 -45.5537530 -8.5879849 49
150 137.5184750 113.2547360 161.7822140 50
151 -26.5351671 -44.6623337 -8.4080005 51
152 -10.1740339 -23.9147894 3.5667216 52
153 48.8838662 38.5263291 59.2414033 53
154 147.3624829 122.6826584 172.0423074 54
155 5.8011050 -5.9397134 17.5419234 55
156 56.8830184 44.1812330 69.5848039 56
157 49.4111410 35.1944093 63.6278728 57
158 140.0086470 111.1691957 168.8480983 58
159 112.8005660 93.3991489 132.2019831 59
160 -60.6188772 -82.5303222 -38.7074322 60
161 25.7129350 15.4619112 35.9639589 61
162 -25.7325480 -46.2636538 -5.2014422 62
163 57.7155980 43.2549058 72.1762901 63
164 42.2311327 31.7340057 52.7282597 64
165 136.6207840 114.2842645 158.9573034 65
166 224.3371932 184.9748660 263.6995205 66
167 148.3497439 124.0397237 172.6597642 67
168 11.0726940 -0.1674054 22.3127934 68
169 -1.0223761 -15.5199424 13.4751902 69
170 22.5685737 12.4695210 32.6676264 70
171 87.5388702 72.8600422 102.2176983 71
172 53.1549939 39.0853979 67.2245900 72
173 5.4159238 -7.1966085 18.0284561 73
174 8.7942407 -2.5713529 20.1598342 74
175 36.3051967 25.2099074 47.4004861 75
176 -41.1991429 -59.5897243 -22.8085616 76
177 62.6370516 51.3478776 73.9262255 77
178 10.0604560 -1.3100519 21.4309638 78
179 4.0936952 -8.3068566 16.4942470 79
180 2.3779832 -10.5870523 15.3430186 80
181 107.2018605 90.0517990 124.3519221 81
182 58.4140982 45.2604119 71.5677845 82
183 90.0976148 72.5479652 107.6472643 83
184 -28.9685405 -45.1443579 -12.7927231 84
185 -19.2032652 -34.2558894 -4.1506409 85
186 -35.5826323 -53.4238979 -17.7413668 86
187 73.9650150 59.5893702 88.3406598 87
188 33.7327732 23.9163639 43.5491825 88
189 7.9659654 -3.4447439 19.3766747 89
190 37.9734476 21.2687271 54.6781681 90
191 104.5139081 82.2752198 126.7525965 91
192 -39.0968648 -56.5124388 -21.6812909 92
193 70.0455403 57.9905145 82.1005662 93
194 49.1607498 38.7866221 59.5348774 94
195 -98.5036708 -127.4652819 -69.5420596 95
196 65.3207235 53.6860259 76.9554212 96
197 -34.0258283 -61.3781872 -6.6734693 97
198 -5.3053829 -29.2419879 18.6312221 98
199 1.0708028 -20.1268177 22.2684234 99
200 45.3761807 31.1880702 59.5642912 100

Проделаем аналогичное для \(95\%\)-прогностического интервала:

predicted <- as.data.frame(predict(best_model, test, level = 0.95, interval = "prediction"))
predicted$num <- c(1:100)
ggplot(predicted, aes(num)) + geom_line(aes(y=fit), colour="red") +
  geom_ribbon(aes(ymin = lwr, ymax = upr), alpha = 0.3, fill = "green")

kable(predicted[,names(predicted)!="x"])
fit lwr upr num
101 73.2246094 -25.6121107 172.061329 1
102 -66.3767498 -167.0543883 34.300889 2
103 34.8662974 -63.5190028 133.251598 3
104 84.2769385 -15.7643306 184.318208 4
105 13.8495957 -85.4092441 113.108435 5
106 103.9651052 4.4365354 203.493675 6
107 36.7540458 -62.2161470 135.724239 7
108 52.4788915 -46.7417839 151.699567 8
109 21.8494572 -77.6833949 121.382309 9
110 -81.6611350 -182.6356845 19.313415 10
111 -0.7435148 -99.9501478 98.463118 11
112 74.6633723 -24.4179550 173.744700 12
113 49.6406095 -49.9094876 149.190707 13
114 75.4550836 -23.5257881 174.435955 14
115 85.3753387 -13.5776654 184.328343 15
116 -7.8361368 -106.5656957 90.893422 16
117 -44.8171478 -144.4118371 54.777541 17
118 77.7994164 -20.9350392 176.533872 18
119 76.8089669 -22.5189146 176.136848 19
120 43.6353941 -54.7750945 142.045883 20
121 -26.1100809 -125.3654698 73.145308 21
122 -15.3534170 -114.3647358 83.657902 22
123 146.2905839 45.4043108 247.176857 23
124 21.4482734 -77.0173284 119.913875 24
125 40.0886270 -61.0016249 141.178879 25
126 43.4571381 -55.4079194 142.322196 26
127 -73.9812940 -174.6773023 26.714714 27
128 24.5812007 -73.8496974 123.012099 28
129 99.2372609 0.0050105 198.469511 29
130 -80.5950684 -181.7645458 20.574409 30
131 78.2638789 -21.1427980 177.670556 31
132 73.9735485 -25.9119031 173.859000 32
133 6.1315285 -92.5343709 104.797428 33
134 31.5735593 -68.0221700 131.169288 34
135 66.6091085 -32.3804791 165.598696 35
136 26.1563797 -72.8260128 125.138772 36
137 -31.1096615 -130.3226730 68.103350 37
138 65.7735413 -33.4421591 164.989242 38
139 23.2357127 -75.2198578 121.691283 39
140 115.7431373 15.6692036 215.817071 40
141 -104.6038044 -206.6849746 -2.522634 41
142 -33.4943947 -133.0220243 66.033235 42
143 45.1414126 -53.5342873 143.817113 43
144 110.9877164 11.1554117 210.820021 44
145 78.6240027 -20.3195616 177.567567 45
146 89.7231120 -9.4373613 188.883585 46
147 16.6000070 -82.7630827 115.963097 47
148 25.4527677 -72.9995110 123.905046 48
149 -27.0708690 -126.6910087 72.549271 49
150 137.5184750 36.6656898 238.371260 50
151 -26.5351671 -126.0899230 73.019589 51
152 -10.1740339 -109.0242418 88.676174 52
153 48.8838662 -49.5530873 147.320820 53
154 147.3624829 46.4087857 248.316180 54
155 5.8011050 -92.7909944 104.393204 55
156 56.8830184 -41.8281289 155.594166 56
157 49.4111410 -49.5063534 148.328635 57
158 140.0086470 37.9583203 242.058974 58
159 112.8005660 13.0059257 212.595206 59
160 -60.6188772 -160.9317211 39.693967 60
161 25.7129350 -72.7128682 124.138738 61
162 -25.7325480 -125.7529536 74.287858 62
163 57.7155980 -41.2372537 156.668450 63
164 42.2311327 -56.2206063 140.682872 64
165 136.6207840 36.2142337 237.027334 65
166 224.3371932 118.8291459 329.845240 66
167 148.3497439 47.4858142 249.213674 67
168 11.0726940 -87.4610316 109.606420 68
169 -1.0223761 -99.9806232 97.935871 69
170 22.5685737 -75.8415178 120.978665 70
171 87.5388702 -11.4460944 186.523835 71
172 53.1549939 -45.7414609 152.051449 72
173 5.4159238 -93.2837784 104.115626 73
174 8.7942407 -89.7538793 107.342361 74
175 36.3051967 -62.2121149 134.822508 75
176 -41.1991429 -140.8021987 58.403913 76
177 62.6370516 -35.9022842 161.176387 77
178 10.0604560 -88.4882309 108.609143 78
179 4.0936952 -94.5791426 102.766533 79
180 2.3779832 -96.3673832 101.123350 80
181 107.2018605 7.8203741 206.583347 81
182 58.4140982 -40.3562147 157.184411 82
183 90.0976148 -9.3536067 189.548836 83
184 -28.9685405 -128.1865445 70.249464 84
185 -19.2032652 -118.2443513 79.837821 85
186 -35.5826323 -135.0857278 63.920463 86
187 73.9650150 -24.9754439 172.905474 87
188 33.7327732 -64.6487145 132.114261 88
189 7.9659654 -90.5873680 106.519299 89
190 37.9734476 -61.3321561 137.279051 90
191 104.5139081 4.1290761 204.898740 91
192 -39.0968648 -138.5245142 60.330785 92
193 70.0455403 -28.5844698 168.675550 93
194 49.1607498 -49.2779508 147.599450 94
195 -98.5036708 -200.5885872 3.581246 95
196 65.3207235 -33.2587948 163.900242 96
197 -34.0258283 -135.6659128 67.614256 97
198 -5.3053829 -106.0799645 95.469199 98
199 1.0708028 -99.0885420 101.230148 99
200 45.3761807 -53.5372042 144.289566 100