Линейная регрессия — это контролируемый метод моделирования непрерывных данных. Модель соответствует линии, наиболее близкой ко всем наблюдениям в наборе данных. Основное предположение здесь состоит в том, что функциональная форма — это линия, и можно подобрать линию, которая будет наиболее близка ко всем наблюдениям в наборе данных. Обратите внимание, что если базовое предположение о линейности модели не соответствует действительности, то в модели обязательно будет ошибка (смещение в сторону линейности), однако лучше всего пытаться подогнать модель.
Цель линейной регрессии состоит в том, чтобы найти математическое уравнение для непрерывной переменной отклика Y как функции одной или нескольких переменных X. Так что вы можете использовать эту регрессионную модель для прогнозирования Y, когда известно только X.
Это математическое уравнение можно обобщить следующим образом:
Y=β1+β2 X+ϵ
где β1 — точка пересечения, а β2 — наклон.
В совокупности они называются коэффициентами регрессии, а ϵ — это член ошибки, часть Y, которую регрессионная модель не может объяснить.
Как применить линейную регрессию
Коэффициент линейной регрессии рассчитывается на основе выборочных данных. Основное предположение здесь состоит в том, что выборка не является предвзятой. Это предположение гарантирует, что выборка не обязательно всегда завышает или занижает коэффициенты. Идея состоит в том, что конкретная выборка может завышать или занижать оценку, но если взять несколько выборок и попытаться оценить коэффициент несколько раз, то среднее значение коэффициента из нескольких выборок будет точным.
Проблема:
Для этого поста давайте возьмем набор данных Boston, который является частью библиотеки MASS в R Studio. Ниже приведены функции, доступные в наборе данных Boston. Постановка задачи состоит в том, чтобы предсказать «медв» на основе набора входных признаков.
#Extract the data and create the training and testing sample library(MASS) library(ggplot2) attach(Boston) > names(Boston) [1] "crim" "zn" "indus" "chas" "nox" "rm" [7] "age" "dis" "rad" "tax" "ptratio" "black" [13] "lstat" "medv"
Графический анализ
Давайте попробуем понять эти переменные графически.
Как правило, для каждого из предикторов визуализировать закономерности помогают следующие графики:
- График рассеивания. Визуализируйте линейную зависимость между предиктором и ответом.
- График плотности: чтобы увидеть распределение переменной-предиктора. В идеале предпочтительнее распределение, близкое к нормальному (колоколообразная кривая), без перекоса влево или вправо.
Давайте посмотрим, как сделать каждый из них.
Использование точечной диаграммы для визуализации отношений
Диаграммы рассеивания могут помочь визуализировать линейные отношения между откликом и предикторными переменными.
В идеале, если у вас много переменных-предикторов, для каждой из них строится точечная диаграмма в зависимости от ответа, а также линия наилучшего соответствия, как показано ниже. (Ниже скриншот для переменной «crim» против «medv»)
Использование BoxPlot для проверки выбросов
Как правило, выбросом является любая точка данных, которая находится за пределами 1,5 * межквартильный диапазон (IQR).
IQR рассчитывается как расстояние между значениями 25-го процентиля и 75-го процентиля для этой переменной.
Как правило, выбросом является любая точка данных, которая находится за пределами 1,5 * межквартильный диапазон (IQR).
IQR рассчитывается как расстояние между значениями 25-го процентиля и 75-го процентиля для этой переменной.
Многомерный модельный подход
Объявление наблюдения выбросом, основанным только на одном (весьма неважном) признаке, может привести к нереалистичным выводам. Когда вам нужно решить, является ли отдельный объект (представленный строкой или наблюдением) экстремальным значением или нет, лучше коллективно рассмотреть важные функции (X). Это когда Дистанция Кука играет роль.
Расстояние повара
Расстояние Кука — это мера, вычисляемая по отношению к данной регрессионной модели, и поэтому на нее влияют только X-переменные, включенные в модель. Но что означает дистанция повара? Он вычисляет влияние каждой точки данных (строки) на прогнозируемый результат.
Расстояние повара для каждого наблюдения i измеряет изменение Y^Y^ (подобранное Y) для всех наблюдений с присутствием наблюдения i и без него, поэтому мы знаем, насколько наблюдение i повлияло на подобранные значения. Математически расстояние Кука Di для наблюдения i вычисляется как:
где,
- Y^jY^j — это значение j-го подходящего ответа, когда включены все наблюдения.
- Y^j(i)Y^j(i) — это значение j-го подходящего ответа, где подбор не включает наблюдение i.
- MSE – среднеквадратическая ошибка.
- p — количество коэффициентов в регрессионной модели.
mod <- lm(ozone_reading ~ ., data=ozone)
cooksd <- cooks.distance(model1) #model1 develop using the baseline approachМеры влияния
В общем, те наблюдения, которые имеют расстояние повара, превышающее среднее значение более чем в 4 раза, могут быть классифицированы как влиятельные. Это не жесткая граница.
plot(cooksd, pch="*", cex=2, main="Influential Obs by Cooks distance") # plot cook's distance
abline(h = 4*mean(cooksd, na.rm=T), col="red") # add cutoff line
text(x=1:length(cooksd)+1, y=cooksd, labels=ifelse(cooksd>4*mean(cooksd, na.rm=T),names(cooksd),""), col="red") # add labelsЕсли вы извлечете и изучите каждую влиятельную строку 1 на 1 (из вывода ниже), вы сможете понять, почему эта строка оказалась влиятельной. Вполне вероятно, что одна из переменных X, включенных в модель, имела экстремальные значения.
Использование графика плотности, чтобы проверить, близка ли переменная отклика к нормальной
Давайте проверим распределение переменной ответа «medv». На следующем рисунке показаны три распределения исходного «medv», логарифмического преобразования и преобразования квадратного корня. Мы видим, что и «log», и «sqrt» делают достойную работу по преобразованию распределения «medv» ближе к нормальному. В следующей модели я выбрал преобразование «log», но также можно попробовать преобразование «sqrt».
Корреляционный анализ
Корреляционный анализ изучает силу связи между двумя непрерывными переменными. Он включает в себя вычисление коэффициента корреляции между двумя переменными.
Так что же такое корреляция? И как это полезно в линейной регрессии?
Корреляция — это статистическая мера, показывающая степень линейной зависимости между двумя переменными.
Чтобы вычислить корреляцию, две переменные должны встречаться парами.
Корреляция может принимать значения от -1 до +1.
Если одна переменная постоянно увеличивается с увеличением значения другой, то они имеют сильную положительную корреляцию (значение близкое к +1).
Точно так же, если один последовательно уменьшается, когда другой увеличивается, они имеют сильную отрицательную корреляцию (значение, близкое к -1).
Значение ближе к 0 предполагает слабую связь между переменными.
Низкая корреляция (-0,2 ‹ x ‹ 0,2), вероятно, предполагает, что большая часть вариации переменной отклика (Y) не объясняется предиктором (X). В этом случае вам, вероятно, следует искать лучшие объясняющие переменные.
Если вы наблюдаете за набором данных Boston в консоли R, для каждого случая, когда rm увеличивается, medv также увеличивается вместе с ним.
Это означает, что между ними существует сильная положительная связь. Значит, корреляция между ними будет ближе к 1.
Однако корреляция не означает причинно-следственной связи.
Другими словами, если две переменные имеют высокую корреляцию, это не означает, что одна переменная «вызывает» увеличение значения другой переменной.
Корреляция - это только помощь для понимания отношений. Вы можете полагаться только на логику и деловые рассуждения, чтобы сделать это суждение.
Итак, как вычислить корреляцию в R?
Просто используйте функцию cor() с двумя числовыми переменными в качестве аргументов.
cor(Boston$rm, Boston$medv) # calculate correlation between rm and medv
#> [1] 0.6953599
Разделите образцы данных и создайте модель
Продолжим построение модели. Разделите входные данные на обучающий и оценочный наборы и создайте модель для обучающего набора данных. Можно видеть, что набор обучающих данных содержит 404 наблюдения, а набор тестовых данных содержит 102 наблюдения на основе разделения 80–20.
##Sample the dataset. set.seed(1) row.number <- sample(1:nrow(Boston), 0.8*nrow(Boston)) train = Boston[row.number,] test = Boston[-row.number,] > dim(train) [1] 404 14 > dim(test) [1] 102 14
Построение модели — Модель 1
Теперь в качестве первого шага мы подгоним модели множественной регрессии. Мы начнем с того, что возьмем все входные переменные в множественной регрессии.
#Let’s make default model.#Let’s make default model. model1 = lm(log(medv)~., data=train) summary(model1) par(mfrow=c(2,2)) for(i in 1:4) plot(model1, which = i) Call: lm(formula = log(medv) ~ ., data = train) Residuals: Min 1Q Median 3Q Max -0.71392 -0.10435 -0.00913 0.10259 0.83290 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 4.5387176 0.2331496 19.467 < 2e-16 *** crim -0.0110546 0.0015143 -7.300 1.63e-12 *** zn 0.0014176 0.0006281 2.257 0.024574 * indus 0.0020512 0.0028308 0.725 0.469120 chas 0.0853159 0.0402646 2.119 0.034732 * nox -0.9285807 0.1693534 -5.483 7.52e-08 *** rm 0.0589055 0.0189839 3.103 0.002056 ** age 0.0002373 0.0006075 0.391 0.696247 dis -0.0598220 0.0091621 -6.529 2.06e-10 *** rad 0.0152004 0.0030216 5.031 7.47e-07 *** tax -0.0005681 0.0001709 -3.325 0.000968 *** ptratio -0.0427382 0.0059698 -7.159 4.07e-12 *** black 0.0003423 0.0001209 2.831 0.004885 ** lstat -0.0319466 0.0022703 -14.071 < 2e-16 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 0.193 on 390 degrees of freedom Multiple R-squared: 0.7933, Adjusted R-squared: 0.7864 F-statistic: 115.1 on 13 and 390 DF, p-value: < 2.2e-16
Дает этот сюжет:
Наблюдение из сводки (модель 1)
Есть ли связь между предиктором и переменными ответа?
Мы можем ответить на этот вопрос, используя F-статистику. Это определяет коллективный эффект всех переменных-предикторов на переменную отклика. В этой модели F = 102,3 намного больше 1, поэтому можно сделать вывод, что существует связь между предиктором и переменной отклика.
Какие из переменных-предикторов являются значимыми?
Основываясь на «p-значении», мы можем сделать вывод об этом. Чем меньше значение «p», тем более значимой является переменная. Из «суммарного» дампа видно, что «zn», «возраст» и «инд» являются менее значимыми признаками, так как значение «p» для них велико. В следующей модели мы можем удалить эти переменные из модели.
Подходит ли эта модель?
Мы можем ответить на этот вопрос, основываясь на значении R2 (множественный R-квадрат), поскольку оно показывает, насколько вариации учитываются моделью. R2 ближе к 1 указывает на то, что модель объясняет большое значение дисперсии модели и, следовательно, хорошо подходит. В данном случае значение равно 0,7933 (ближе к 1), и, следовательно, модель подходит.
Наблюдение за сюжетом
График сопоставления и остатка
Графики остатка должны быть случайными по своей природе, и на графике не должно быть закономерностей. Среднее значение остаточного графика должно быть близко к нулю. На приведенном выше графике мы видим, что красная линия тренда находится почти на нуле, за исключением начального местоположения.
Нормальный график Q-Q
График Q-Q показывает, нормально ли распределены остатки. В идеале сюжет должен быть на пунктирной линии. Если график Q-Q не находится на прямой, то модели необходимо переработать, чтобы сделать невязку нормальной. На приведенном выше графике мы видим, что большинство графиков находятся на линии, кроме конца.
Масштаб-местоположение
Показывает, как распределяются остатки и имеют ли остатки одинаковую дисперсию или нет.
Остатки против кредитного плеча
График помогает найти важные наблюдения. Здесь нам нужно проверить точки, которые находятся за пределами пунктирной линии. Точка за пределами пунктирной линии будет влиятельной точкой, и ее удаление повлияет на коэффициенты регрессии.
Построение модели — Модель 2
В качестве следующего шага мы можем удалить четыре менее значимых признака («zn», «возраст» и «indus») и снова проверить модель.
# remove the less significant feature
model2 = update(model1, ~.-zn-indus-age)
summary(model2)
par(mfrow=c(2,2))
plot(model2)
Call:
lm(formula = log(medv) ~ crim + chas + nox + rm + dis + rad +
tax + ptratio + black + lstat, data = train)
Residuals:
Min 1Q Median 3Q Max
-0.71008 -0.10962 -0.01296 0.10330 0.84030
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 4.5489757 0.2327067 19.548 < 2e-16 ***
crim -0.0107965 0.0015120 -7.140 4.54e-12 ***
chas 0.0854786 0.0400433 2.135 0.033407 *
nox -0.9117698 0.1566753 -5.819 1.23e-08 ***
rm 0.0644686 0.0184026 3.503 0.000513 ***
dis -0.0523327 0.0072616 -7.207 2.95e-12 ***
rad 0.0143808 0.0028974 4.963 1.03e-06 ***
tax -0.0004624 0.0001505 -3.073 0.002263 **
ptratio -0.0464620 0.0055623 -8.353 1.16e-15 ***
black 0.0003416 0.0001211 2.821 0.005026 **
lstat -0.0314967 0.0021598 -14.583 < 2e-16 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1935 on 393 degrees of freedom
Multiple R-squared: 0.7905, Adjusted R-squared: 0.7851
F-statistic: 148.3 on 10 and 393 DF, p-value: < 2.2e-16
Дает этот сюжет:
Наблюдение из сводки (модель 1)
Есть ли связь между предиктором и переменной ответа?
F=131,2 намного больше 1, и это значение больше, чем значение F предыдущей модели. Можно сделать вывод, что существует связь между предиктором и переменной отклика.
Какая из переменных является значимой?
Теперь в этой модели все предикторы являются значимыми.
Подходит ли эта модель?
R2 = 0,7905 ближе к 1, поэтому эта модель подходит. Обратите внимание, что это значение немного уменьшилось по сравнению с первой моделью, но это должно быть нормально, так как удаление трех предикторов привело к падению с 0,7933 до 0,7905, и это небольшое падение. Другими словами, вклад трех предикторов в объяснение дисперсии составляет лишь небольшое значение (0,0028), и, следовательно, предиктор лучше исключить.
Наблюдение за графиком
Все четыре графика похожи на предыдущую модель, и мы не видим какого-либо существенного эффекта.
Построение модели — Модель 3 и Модель 4
Теперь мы можем улучшить модель, добавив квадратный член для проверки нелинейности. Сначала мы можем попробовать модель 3, введя квадратные члены для всех функций (из модели 2). А на следующей итерации мы можем удалить незначительную особенность из модели.
#Lets make default model and add square term in the model.
model3 = lm(log(medv)~crim+chas+nox+rm+dis+rad+tax+ptratio+
black+lstat+ I(crim^2)+ I(chas^2)+I(nox^2)+ I(rm^2)+ I(dis^2)+
I(rad^2)+ I(tax^2)+ I(ptratio^2)+ I(black^2)+ I(lstat^2), data=train)
summary(model3)
Call:
lm(formula = log(medv) ~ crim + chas + nox + rm + dis + rad +
tax + ptratio + black + lstat + I(crim^2) + I(chas^2) + I(nox^2) +
I(rm^2) + I(dis^2) + I(rad^2) + I(tax^2) + I(ptratio^2) +
I(black^2) + I(lstat^2), data = train)
Residuals:
Min 1Q Median 3Q Max
-0.74190 -0.08032 -0.00534 0.09105 0.75405
Coefficients: (1 not defined because of singularities)
Estimate Std. Error t value Pr(>|t|)
(Intercept) 8.266e+00 9.571e-01 8.636 < 2e-16 ***
crim -2.754e-02 4.094e-03 -6.726 6.33e-11 ***
chas 7.821e-02 3.670e-02 2.131 0.033686 *
nox -1.723e+00 1.099e+00 -1.567 0.117826
rm -6.355e-01 1.406e-01 -4.518 8.31e-06 ***
dis -1.267e-01 2.534e-02 -5.001 8.71e-07 ***
rad 2.150e-02 1.051e-02 2.046 0.041442 *
tax -1.140e-03 5.835e-04 -1.954 0.051428 .
ptratio -1.543e-01 7.778e-02 -1.984 0.047922 *
black 1.622e-03 5.486e-04 2.957 0.003299 **
lstat -4.756e-02 6.012e-03 -7.910 2.76e-14 ***
I(crim^2) 2.052e-04 5.299e-05 3.873 0.000126 ***
I(chas^2) NA NA NA NA
I(nox^2) 4.512e-01 8.146e-01 0.554 0.579939
I(rm^2) 5.402e-02 1.105e-02 4.888 1.50e-06 ***
I(dis^2) 6.531e-03 2.017e-03 3.238 0.001307 **
I(rad^2) -1.448e-04 4.148e-04 -0.349 0.727159
I(tax^2) 7.655e-07 7.322e-07 1.045 0.296450
I(ptratio^2) 3.347e-03 2.216e-03 1.510 0.131741
I(black^2) -3.056e-06 1.203e-06 -2.540 0.011487 *
I(lstat^2) 5.130e-04 1.684e-04 3.046 0.002477 **
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1741 on 384 degrees of freedom
Multiple R-squared: 0.8343, Adjusted R-squared: 0.8261
F-statistic: 101.7 on 19 and 384 DF, p-value: < 2.2e-16
##Removing the insignificant variables.
model4=update(model3, ~.-nox-rad-tax-I(crim^2)-I(chas^2)-I(rad^2)-
I(tax^2)-I(ptratio^2)-I(black^2))
summary(model4)
par(mfrow=c(2,2))
plot(model4)
Call:
lm(formula = log(medv) ~ crim + chas + rm + dis + ptratio + black +
lstat + I(nox^2) + I(rm^2) + I(dis^2) + I(lstat^2), data = train)
Residuals:
Min 1Q Median 3Q Max
-0.73918 -0.09787 -0.00723 0.08868 0.82585
Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) 6.6689322 0.4622249 14.428 < 2e-16 ***
crim -0.0104797 0.0013390 -7.827 4.72e-14 ***
chas 0.0909580 0.0379316 2.398 0.016954 *
rm -0.7616673 0.1445966 -5.268 2.29e-07 ***
dis -0.0918453 0.0239087 -3.842 0.000143 ***
ptratio -0.0308017 0.0049005 -6.285 8.72e-10 ***
black 0.0002632 0.0001134 2.321 0.020824 *
lstat -0.0440988 0.0060377 -7.304 1.57e-12 ***
I(nox^2) -0.6629324 0.1148901 -5.770 1.61e-08 ***
I(rm^2) 0.0654052 0.0113152 5.780 1.52e-08 ***
I(dis^2) 0.0045419 0.0020195 2.249 0.025063 *
I(lstat^2) 0.0003587 0.0001657 2.165 0.030997 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 0.1834 on 392 degrees of freedom
Multiple R-squared: 0.8123, Adjusted R-squared: 0.807
F-statistic: 154.2 on 11 and 392 DF, p-value: < 2.2e-16
Наблюдение из сводки (модель 4)
Существует ли связь между предиктором и переменными отклика?
F-Stat равен 154,2 и намного больше 1. Таким образом, между предиктором и откликом существует взаимосвязь.
Какие из переменных-предикторов являются значимыми?
Все переменные-предикторы являются значимыми.
Подходит ли эта модель?
R2 равно 0,8123, и это больше (и лучше), чем у первой и второй модели.
Прогноз
До сих пор мы проверяли ошибку обучения, но настоящая цель модели — уменьшить ошибку тестирования. Поскольку мы уже разделили набор данных образца на набор данных для обучения и тестирования, мы будем использовать набор тестовых данных для оценки модели, к которой мы пришли. Мы сделаем прогноз на основе «Модели 4» и оценим модель. В качестве последнего шага мы предскажем «тестовое» наблюдение и увидим сравнение между предсказанным откликом и фактическим значением отклика. RMSE объясняет в среднем, какая часть прогнозируемого значения будет отличаться от фактического значения. Основываясь на RMSE = 3,682, мы можем сделать вывод, что среднее прогнозируемое значение будет на 3,682 отличаться от фактического значения.
pred1 <- predict(model4, newdata = test) rmse <- sqrt(sum((exp(pred1) - test$medv)^2)/length(test$medv)) c(RMSE = rmse, R2=summary(model4)$r.squared) RMSE R2 3.6825564 0.8122965 par(mfrow=c(1,1)) plot(test$medv, exp(pred1))
Дает этот сюжет:
Как узнать, какая модель регрессии лучше всего подходит для данных?
Наиболее распространенными показателями, на которые следует обратить внимание при выборе модели, являются:
Вывод
Пример показывает, как подойти к моделированию линейной регрессии. Созданная модель все еще имеет возможности для улучшения, поскольку мы можем применять такие методы, как удаление выбросов с использованием методов вменения, обнаружение корреляции для дальнейшего повышения точности более точного прогноза. Существуют и другие методы проверки возможности дальнейшего повышения точности модели, такие как метод случайного леса и бустинга.
