Pavel Polishchuk, 2014
Примеры построения моделей в R.
Сгенерируем набор данных для моделирования в соответствии с формулой
\[ Y = x_1 + x_2^2 + x_1 \times x_2 + noise \]
Для этого вначале инициализируем генератор случайных чисел, чтобы получать всегда воспроизводимый результат при генерации случайных чисел.
set.seed(50)
Сгенерируем случайные переменные \( x_1 \) и \( x_2 \) в соответствии с законом нормального распределения
x1 <- rnorm(n=200, mean=5, sd=2)
x2 <- rnorm(200, 2, 3)
Сгенерируем зависимую переменную в соответствии с вышеприведенной формулой
y <- x1 + x2^2 + x1*x2 + rnorm(200, 0, 0.5)
Объединим эти параметры в одном data.frame
df <- data.frame(y, x1, x2)
Рассмотрим созданный набор данных. Построим гистограмму распределению значений Y и посмотрим как рапределены значения переменных между собой попарно.
hist(df$y)
pairs(df)
Посчитаем матрицу взаимных корреляций всех переменных между собой и округлим результат до двух цифр после запятой
round(cor(df), 2)
y x1 x2
y 1.00 0.18 0.86
x1 0.18 1.00 -0.05
x2 0.86 -0.05 1.00
Помимо традиционного коэффицента корреляции Пирсона функцией cor можно посчитать ранговый коэффицент корреляции Спирмена.
round(cor(df, method="spearman"), 2)
y x1 x2
y 1.00 0.12 0.94
x1 0.12 1.00 -0.02
x2 0.94 -0.02 1.00
И из рисунка и из таблиц видно, что связь между \( y \) и \( x_2 \) довольно тесная.
Альтернативный способ представления данных - использование функции pairs.panels пакета psych, которая возвращает одновременно и диаграммы распределения данных и значения коэффициентов корреляции.
psych::pairs.panels(df)
Зададим индексы для объектов обучающей и тестовой выборки
train.ids <- 1:100
test.ids <- 101:200
Построим несколько моделей в соответствии с различными формулами, которые определяют вид искомой зависимости.
Попробуем сперва обычную линейную комбинацию исходных переменных.
m1 <- lm(y ~ x1 + x2, df[train.ids,])
summary(m1)
Call:
lm(formula = y ~ x1 + x2, data = df[train.ids, ])
Residuals:
Min 1Q Median 3Q Max
-15.62 -7.56 -3.43 1.95 54.96
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.379 3.493 -1.83 0.071 .
x1 3.377 0.625 5.40 4.8e-07 ***
x2 8.567 0.415 20.65 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.3 on 97 degrees of freedom
Multiple R-squared: 0.818, Adjusted R-squared: 0.814
F-statistic: 218 on 2 and 97 DF, p-value: <2e-16
В результате мы получили уравнение вида \[ Y = 3.377\pm0.625*x_1 + 8,567\pm0,415*x_2 - 6,379\pm3,493 \]
Все коэффициенты в уравнении статистически значимы с высоким уронеем достоверности, кроме свободного члена, который имеет низкий уровень статистической достоверности. Величина \( R^2 = 0.818 \) является статистически достоверной.
\[ R^2 = 1 - \frac{\sigma^2}{\sigma_y^2} = 1 - \frac{^{ESS}/_{n}}{^{TSS}/_{n}} = 1 - \frac{ESS}{TSS} = 1 - \frac{\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}{\sum_{i=1}^{n}(y_i - \bar{y})^2} \]
\[ R_{adj}^2 = 1 - \frac{s^2}{s_y^2} = 1 - \frac{^{ESS}/_{n-k}}{^{TSS}/_{n-1}} =1 - (1- R^2)\frac{n-1}{n-k} \leq R^2 \]
\( y_i \) - наблюдаемое значение \( y \)
\( \hat{y}_i \) - предсказанное по модели значение \( y \)
\( \bar{y} \) - среднее значение наблюдаемых значений \( y \)
\( n \) - количество наблюдений
\( k \) - количество параметров
Применим одну из простейших процедур отбора переменных - исключим незначащие члены из уравнения
m2 <- lm(y ~ x1 + x2 + 0, df[train.ids,])
summary(m2)
Call:
lm(formula = y ~ x1 + x2 + 0, data = df[train.ids, ])
Residuals:
Min 1Q Median 3Q Max
-17.73 -8.88 -4.61 2.30 49.03
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x1 2.361 0.290 8.15 1.2e-12 ***
x2 8.267 0.386 21.44 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 12.4 on 98 degrees of freedom
Multiple R-squared: 0.914, Adjusted R-squared: 0.913
F-statistic: 524 on 2 and 98 DF, p-value: <2e-16
Коэффициент детерминации существенно увеличился \( R^2 = 0.914 \)
Определим с помощью дисперсионного анализа является ли отличие моделей значимым или нет.
anova(m1, m2)
Analysis of Variance Table
Model 1: y ~ x1 + x2
Model 2: y ~ x1 + x2 + 0
Res.Df RSS Df Sum of Sq F Pr(>F)
1 97 14603
2 98 15105 -1 -502 3.33 0.071 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Величина 0.071 > 0.05, поэтому можно утверждать, что с вероятностью 95% отличие моделей не значимо и мы в праве выбрать любую модель.
Попытаемся построить более сложные модели.
Модель учитывающее произведение независимых переменных
m3 <- lm(y ~ x1*x2, df[train.ids,])
summary(m3)
Call:
lm(formula = y ~ x1 * x2, data = df[train.ids, ])
Residuals:
Min 1Q Median 3Q Max
-10.35 -7.51 -3.06 3.77 51.30
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 5.132 3.623 1.42 0.160
x1 1.160 0.664 1.75 0.084 .
x2 3.866 0.890 4.34 3.5e-05 ***
x1:x2 0.962 0.167 5.77 9.5e-08 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 10.6 on 96 degrees of freedom
Multiple R-squared: 0.865, Adjusted R-squared: 0.861
F-statistic: 205 on 3 and 96 DF, p-value: <2e-16
Модель второго порядка, учитывающая квадраты независимых переменных
m4 <- lm(y ~ poly(x1, 2) + poly(x2, 2), df[train.ids,])
summary(m4)
Call:
lm(formula = y ~ poly(x1, 2) + poly(x2, 2), data = df[train.ids,
])
Residuals:
Min 1Q Median 3Q Max
-15.569 -2.232 -0.475 1.683 22.654
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 31.056 0.594 52.27 < 2e-16 ***
poly(x1, 2)1 65.297 5.983 10.91 < 2e-16 ***
poly(x1, 2)2 -24.414 5.979 -4.08 9.3e-05 ***
poly(x2, 2)1 253.662 5.991 42.34 < 2e-16 ***
poly(x2, 2)2 105.653 5.971 17.69 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 5.94 on 95 degrees of freedom
Multiple R-squared: 0.958, Adjusted R-squared: 0.956
F-statistic: 544 on 4 and 95 DF, p-value: <2e-16
В последнем случае видно, что все коэффициенты в уравнении имеют высокую значимость и коэффициент детерминации также высокий \( R^2 = 0.958 \)
Сравнение моделей указывает на значимое отличие последней модели от предыдущей.
anova(m2, m4)
Analysis of Variance Table
Model 1: y ~ x1 + x2 + 0
Model 2: y ~ poly(x1, 2) + poly(x2, 2)
Res.Df RSS Df Sum of Sq F Pr(>F)
1 98 15105
2 95 3354 3 11751 111 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Для сравнения построим модель с использованием правильной функциональной зависимости
m5 <- lm(y ~ x1 + x1:x2 + I(x2^2), df[train.ids,])
summary(m5)
Call:
lm(formula = y ~ x1 + x1:x2 + I(x2^2), data = df[train.ids, ])
Residuals:
Min 1Q Median 3Q Max
-1.0893 -0.3216 0.0275 0.3209 1.0808
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.15588 0.15289 -1.02 0.31
x1 1.03061 0.02750 37.48 <2e-16 ***
I(x2^2) 1.00983 0.00459 219.78 <2e-16 ***
x1:x2 0.99062 0.00436 226.95 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.518 on 96 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 9.97e+04 on 3 and 96 DF, p-value: <2e-16
Свободный членя является не значимым, следовательно зависимость проходит через начало отсчета. Построим модель без свободного члена.
m6 <- lm(y ~ x1 + x1:x2 + I(x2^2) + 0, df[train.ids,])
summary(m6)
Call:
lm(formula = y ~ x1 + x1:x2 + I(x2^2) + 0, data = df[train.ids,
])
Residuals:
Min 1Q Median 3Q Max
-1.0777 -0.3465 0.0152 0.3057 1.1711
Coefficients:
Estimate Std. Error t value Pr(>|t|)
x1 1.00601 0.01319 76.2 <2e-16 ***
I(x2^2) 1.00769 0.00409 246.3 <2e-16 ***
x1:x2 0.99165 0.00425 233.5 <2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Residual standard error: 0.518 on 97 degrees of freedom
Multiple R-squared: 1, Adjusted R-squared: 1
F-statistic: 2.2e+05 on 3 and 97 DF, p-value: <2e-16
pred.m4 <- predict(m4, df[test.ids,])
Сопоставим качество прогноза с истинным значением \( y \)
cor(pred.m4, df$y[test.ids])
[1] 0.9673
cor(pred.m4, df$y[test.ids]) ^ 2
[1] 0.9357
Зададим вид вывода результатов и применим функцию plot к моделям m4 и m6.
# 2 x 2 pictures on one plot and square plotting region, independent of device size
par(mfrow = c(2, 2), pty = "s")
plot(m4)
plot(m6)
Вернем вид вывода результатов к виду по умолчанию
par(mfrow = c(1, 1), pty = "m")
require(pls)
pls.m1 <- plsr(y ~ ., data=df[train.ids, ], ncomp=2, validation="CV", segments=5, segment.type="random")
Воспользуемся функцией summary для вывода статистики модели
summary(pls.m1)
Data: X dimension: 100 2
Y dimension: 100 1
Fit method: kernelpls
Number of components considered: 2
VALIDATION: RMSEP
Cross-validated using 5 random segments.
(Intercept) 1 comps 2 comps
CV 28.6 13.87 14.00
adjCV 28.6 13.76 13.79
TRAINING: % variance explained
1 comps 2 comps
X 69.23 100.00
y 78.68 81.78
Чтобы посчитать коэффициент детерминации для обучающей выборки и для случая кросс-валидации воспользуемся следующими функциями. Следует отметить, что модель pls в качестве прогнозированных значений возвращает массив!
# train
cor(pls.m1$fitted.values[,,2], df$y[train.ids]) ^ 2
[1] 0.8178
# cross-validation
cor(pls.m1$validation$pred[,,2], df$y[train.ids]) ^ 2
[1] 0.7604
Для pls модель также возможно использовать функцию plot, которая возвращает диаграмму зависимости предсказанных значений от наблюдаемых
plot(pls.m1)
abline(0, 1) # добавляет линию выходящую из начала координат и идущей под углом 45 градусов
Спрогнозируем значения тестовой выборки и оценим качество прогноза
pred.pls.m1 <- predict(pls.m1, df[test.ids,])
cor(pred.pls.m1[,,2], df$y[test.ids]) ^ 2
[1] 0.775
Используем другой вид зависимости, который был успешно применен ранее в случае линейной регрессии
pls.m2 <- plsr(y ~ poly(x1, 2) + poly(x2, 2), data=df[train.ids, ], ncomp=2, validation="CV", segments=5, segment.type="random")
summary(pls.m2)
Data: X dimension: 100 4
Y dimension: 100 1
Fit method: kernelpls
Number of components considered: 2
VALIDATION: RMSEP
Cross-validated using 5 random segments.
(Intercept) 1 comps 2 comps
CV 28.6 8.976 6.848
adjCV 28.6 8.246 6.711
TRAINING: % variance explained
1 comps 2 comps
X 24.62 50.23
y 94.86 95.81
cor(pls.m2$validation$pred[,,2], df$y[train.ids]) ^ 2
[1] 0.9415
Как видно качество модели значительно уеличилось. Подтвердим это оценив прогнозирующую способность модели на тестовом наборе данных.
pred.pls.m2 <- predict(pls.m2, df[test.ids,])
cor(pred.pls.m2[,,2], df$y[test.ids]) ^ 2
[1] 0.9344
plot(pls.m2)
abline(0, 1)
require(randomForest)
rf.m <- randomForest(y ~ ., data=df[train.ids, ], do.trace=50, mtry=2)
| Out-of-bag |
Tree | MSE %Var(y) |
50 | 37.52 4.68 |
100 | 36.98 4.61 |
150 | 36.74 4.58 |
200 | 35.9 4.48 |
250 | 36.21 4.52 |
300 | 36.34 4.53 |
350 | 36.35 4.54 |
400 | 36.82 4.59 |
450 | 35.87 4.47 |
500 | 35.72 4.46 |
Следует отметить, что большинство функций имеют несколько способов вызова через определение функции как в вышеприведенном случае, либо через указание матрицы X и вектора Y. Так построение модели можно осуществить через
rf.m <- randomForest(x=df[train.ids, -1], y=df[train.ids, "y"], do.trace=50, mtry=2)
Вызовем статистику модели
rf.m
Call:
randomForest(formula = y ~ ., data = df[train.ids, ], do.trace = 50, mtry = 2)
Type of random forest: regression
Number of trees: 500
No. of variables tried at each split: 2
Mean of squared residuals: 35.72
% Var explained: 95.54
Рассчитаем значения коэффициента детерминации для out-of-bag выборки (внутренняя валидация моделей)
cor(rf.m$predicted, df$y[train.ids]) ^ 2
[1] 0.9556
Функция plot в случае моделей случайного леса возвращает зависимость величины среднеквадратичной ошибки прогноза out-of-bag выборки от количества деревьев в лесе
plot(rf.m)
Чтобы оценить качество модели предскажем тестовый набор данных
pred.rf <- predict(rf.m, df[test.ids,])
cor(pred.rf, df$y[test.ids]) ^ 2
[1] 0.9644
Построим график зависимости между предсказанными и наблюдаемыми значениями
plot(pred.rf, df$y[test.ids])
abline(0, 1)
Конвертируем значения y из df в два класса по заданному граничному значению
df$yclass <- factor(ifelse(df$y > 30, 1, 0))
sapply(df, class)
y x1 x2 yclass
"numeric" "numeric" "numeric" "factor"
Построение классификационной модели ничем не отличается от регресссионных моделей. Алгоритм сам выбирает тип модели (скласссификационная или регрессионная) в зависимости от типа данных передаваемого Y. Если Y представлен числовыми значениями, то по умолчанию будет построена регрессионная модель, если Y представлен как factor, то будет построена классификационная модель.
rf.m2 <- randomForest(yclass ~ . - y, data=df[train.ids, ], do.trace=50, mtry=2, importance=TRUE, nPerm=3)
ntree OOB 1 2
50: 4.00% 3.57% 4.55%
100: 4.00% 3.57% 4.55%
150: 4.00% 3.57% 4.55%
200: 4.00% 3.57% 4.55%
250: 4.00% 3.57% 4.55%
300: 4.00% 3.57% 4.55%
350: 4.00% 3.57% 4.55%
400: 4.00% 3.57% 4.55%
450: 4.00% 3.57% 4.55%
500: 4.00% 3.57% 4.55%
Обратите внимание на способ представления формулы. Знак минус указывает, что эти переменные должны быть исключены из расммотрения при построении модели.
Выведем статистику модели
rf.m2
Call:
randomForest(formula = yclass ~ . - y, data = df[train.ids, ], do.trace = 50, mtry = 2, importance = TRUE, nPerm = 3)
Type of random forest: classification
Number of trees: 500
No. of variables tried at each split: 2
OOB estimate of error rate: 4%
Confusion matrix:
0 1 class.error
0 54 2 0.03571
1 2 42 0.04545
Проверим прогнозирующую способность на тестовом наборе данных
pred.rf.m2 <- predict(rf.m2, df[test.ids,])
Следует помнить, что результатом прогноза будет factor
class(pred.rf.m2)
[1] "factor"
Выведем матрицу ошибок классификации
table(pred.rf.m2, df$yclass[test.ids])
pred.rf.m2 0 1
0 60 6
1 1 33
Выведем значения важности переменных в последней модели
rf.m2$importance
0 1 MeanDecreaseAccuracy MeanDecreaseGini
x1 0.02934 0.002637 0.01767 2.659
x2 0.40098 0.476636 0.42748 46.118
rf.m2$importance[,"MeanDecreaseGini"]
x1 x2
2.659 46.118
Следует помнить, что в R реализовано множество различных алгоритмов построения регрессионных и классификационных моделей. Выбор метода построения модели лежит полностью на исследователе и требует некоторого опыта и навыков. Познакомится с некоторыми другими методами машинного обучения, которые присутствуют в R можно по ссылке http://cran.r-project.org/web/views/MachineLearning.html
xlsx). Построить модели зависимости heating.load и Cooling.load от остальных переменных методами, описанными в уроке (можно также поискать другие методы по приведенной в конце занятия ссылке и использовать их).