Практические руководства
Модель линейной регрессии с Python
Как вы можете построить и проверить качество вашей регрессионной модели с помощью графических и числовых выходных данных
Регрессионные модели - это широко используемые инструменты машинного обучения, позволяющие делать прогнозы на основе данных, изучая взаимосвязь между функциями и непрерывными результатами. Проверка допущений модели и понимание их соответствия так же важна, как и проверка точности и достоверности модели.
Из этой статьи вы узнаете:
- Как построить модель линейной регрессии
- Как оценить модель по точности прогноза и R-квадрат
- Как проверить диагностику модели путем тестирования линейности, нормальности, гомоскедастичности, автокорреляции и мультиколлинеарности
- Как обнаружить выбросы с помощью расстояния Кука
- Как найти влияние каждого наблюдения с помощью стьюдентизированных остатков и H-плеча.
Следя за статьей, я рекомендую вам проверить Jupyter Notebook на моем GitHub для полного анализа и кода.
Вступление
Регрессионный анализ используется для моделирования взаимосвязи между одной зависимой переменной Y (также известной как ответ, цель или результат) и одной или несколькими независимыми переменными X (также известными как предиктор или характеристика). .
Когда у нас есть один предиктор, это «простая» линейная регрессия, а когда у нас более одного предиктора, это «множественная» линейная регрессия.
Как правило, данные не попадают точно на линейную линию, поэтому уравнение регрессии должно включать член ошибки (ϵ).
Подходящие (прогнозируемые) значения, обычно обозначаемые Y-шляпой и коэффициентами β -hat, выбираются моделью, поскольку они минимизируют квадрат расстояния (также называемого остатками) между фактическими Значение Y и наилучшее соответствие (Y-образная шляпа), минимизирующее следующее выражение:
Метод минимизации суммы квадратов остатков называется регрессией обыкновенных наименьших квадратов (OLS).
Модель линейной регрессии
Мы будем строить модель множественной линейной регрессии на наборе данных о жилищном строительстве в Бостоне конца 1970-х годов. Данные состоят из 506 случаев с 14 атрибутами. Давай посмотрим!
1. Чтение данных
names = ['CRIM', 'ZN', 'INDUS', 'CHAS', 'NOX', 'RM', 'AGE', 'DIS', 'RAD', 'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']
data = pd.read_csv("./housing.csv", delim_whitespace=True, names=names)
data.head()
2. Определите предикторные и конечные переменные
Мы будем прогнозировать цены на жилье («MEDV»), используя все остальные переменные.
predictors = ["CRIM", "ZN", "INDUS", "CHAS", "NOX", "RM", "AGE", "DIS", "RAD", "TAX", "PTRATIO", "B", "LSTAT"] outcome = "MEDV" -> Create X and y datasets X = data[predictors] y = data["MEDV"]
3. Создайте модель линейной регрессии.
Для линейной регрессии воспользуемся библиотекой Statsmodels. (Scikit-learn также можно использовать в качестве альтернативы, но здесь я предпочел статистические модели для более подробного анализа регрессионной модели). При этом мы по-прежнему будем использовать Scikit-learn для разделения поездов и тестов. Подробности ниже!
-> Because we are using statsmodels, we need to add the constant term to the X value. import statsmodels.api as sm X = sm.add_constant(X)
Мы разделим наш набор данных на обучающие и тестовые наборы (80% для обучения и 20% для тестирования). Модель регрессии будет учиться на данных обучения, где выходные данные известны, а позже мы обобщим модель на тестовом наборе. Мы будем прогнозировать выходные значения y тестового набора и сравнивать эти прогнозы с фактическими значениями. Это делается для выявления проблем переоснащения или недообучения.
-> Train-test split from sklearn.model_selection import train_test_split train_X, test_X, train_y, test_y = train_test_split(X, y, train_size = 0.8, random_state = 42) -> Linear regression model model = sm.OLS(train_y, train_X) model = model.fit() print(model.summary2())
- Нет. Наблюдений - 404 (80% от общего набора данных).
- R-квадрат (степень точности модели от 0 до 1) составляет 0,75
- p-значение ›0,05 для некоторых предикторов:« ZN »,« INDUS »,« AGE »,« TAX »означает, что мы рассмотрим возможность удаления этих переменных из списка предикторов! Это потому, что p-значения определяют, существуют ли отношения, которые мы наблюдаем в этой выборке, в большей совокупности.
5. Прогноз на тестовой выборке.
Мы будем использовать model.predict в test_X и сравним прогнозы с фактическими значениями test_y.
predictions = model.predict(test_X)
df_results = pd.DataFrame({‘Actual’: test_y, ‘Predicted’: predictions})
6. Оцените модель
Один из важнейших показателей - среднеквадратичная ошибка; квадратный корень из средней квадратичной ошибки в предсказанных значениях y (y-hat). Он измеряет общую точность модели.
Еще одна важная метрика - R-квадрат в диапазоне от 0 до 1; измеряет долю вариации данных, которая учитывается в регрессионной модели. Он показывает, насколько хорошо модель соответствует данным, а R-квадрат = 1 указывает, что прогнозы регрессии идеально соответствуют данным.
Давайте найдем RMSE и R-квадрат для прогнозов.
from sklearn.metrics import r2_score, mean_squared_error RMSE = np.sqrt(mean_squared_error(test_y, predictions)) r2 = r2_score(test_y, predictions) print(RMSE, r2)
RMSE = 4,92
R-квадрат = 0,66
Как мы видим, производительность нашей модели упала с 0,75 (на обучающих данных) до 0,66 (на тестовых данных), и мы ожидаем, что при следующих прогнозах с использованием этой модели мы будем отставать от 4,92.
7. Диагностика модели.
Прежде чем строить модель линейной регрессии, мы делаем следующие предположения:
- Линейность. Связь между X и средним значением Y линейна.
- Нормальность: остатки подчиняются нормальному распределению, а ожидаемое среднее значение остатков равно нулю.
- Гомоскедастичность: дисперсия остатка постоянна для всех значений X.
- Без автокорреляции: остатки не зависят друг от друга.
- Отсутствие мультиколлинеарности. Между переменными-предикторами нет высокой корреляции.
Прежде чем запутаться в производительности модели, важно проверить, полностью ли выполняются предположения. Может быть трудно полностью выполнить все предположения с реальными данными, однако важно понимать, насколько правильно интерпретируется наша модель.
Линейность
Мы предположили, что существует линейная связь между предикторами и результатом, но это может быть неверно. Когда мы делаем моделирование, мы пытаемся подогнать наши данные к функции, которая объясняет наши данные, но линейная регрессия может не быть подходящей оценкой в этом случае.
Мы проверим линейность с помощью диаграммы рассеяния, чтобы увидеть прогнозируемые значения в сравнении с фактическими значениями. В идеале точки данных должны лежать вокруг диагональной линии на графике.
-> Plot the actual vs predicted results sns.lmplot(x='Actual', y='Predicted', data=df_results, fit_reg=False) -> Plot the diagonal line d_line= np.arange(df_results.min().min(), df_results.max().max()) plt.plot(d_line, d_line, color='red', linestyle='--') plt.show()
В наборе данных Бостона мы видим, что идеальной взаимосвязи нет. Что касается нижнего и верхнего пределов значений, наши прогнозы необъективны.
Нормальность
Мы предположили, что члены ошибки (невязки) модели распределены нормально. Это предположение может быть нарушено из-за несоблюдения предположения о линейности или из-за выбросов в наборе данных.
Мы проверим нормальность, проверив график квантиль-квантиль, график Q-Q.
from statsmodels.graphics.gofplots import qqplot fig=qqplot(model.resid_pearson,line=’45',fit=’True’) plt.xlabel(“Theoretical quantiles”) plt.ylabel(“Sample quantiles”) plt.show()
Идеальное совпадение слева показывает, что распределение данных следует красной линии от нижнего левого угла до верхнего правого. В наборе данных Бостона мы видим, что наша модель склонна к заниженной оценке.
Вы можете пойти дальше и построить гистограмму остатков или выполнить тесты нормальности, такие как тест Андерсона-Дарлинга или Шапиро-Уилка.
Гомоскедастичность
Мы предположили, что дисперсия членов ошибки (остатков) постоянна для всех значений X. Однако гетероскедастичность, нарушение гомоскедастичности происходит, когда у нас нет одинаковой дисперсии между членами ошибки.
Это нарушение может произойти, когда модель присваивает разные веса подмножествам данных.
Мы проверим гомоскедастичность, построив график остатков и проверив, однородна ли дисперсия.
fig, ax = plt.subplots(figsize=(5, 5))
sns.regplot(model.fittedvalues,model.resid, scatter_kws={'alpha': 0.25}, line_kws={'color': 'C2', 'lw': 2}, ax=ax)
ax.set_xlabel('predicted')
ax.set_ylabel('residuals')
plt.tight_layout()
plt.show()
Идеальный график слева показывает равномерное распределение остатков, однако мы не можем сказать то же самое для Boston Dataset. Мы должны быть осторожны с гетероскедастичностью, которая указывает на то, что ошибки предсказания различаются для разных диапазонов предсказанных значений, а это означает, что наша модель может быть неполной.
Нет автокорреляции
Автокорреляция возникает, когда остатки не независимы друг от друга, другими словами, когда значение y_t не зависит от y_t-1. Обычно это касается данных временных рядов, таких как цены на акции. Чтобы узнать больше об этом, вы можете прочитать мою статью о прогнозировании временных рядов.
Мы выполним тест Дарбина-Ватсона, чтобы определить, выполняется ли предположение об отсутствии автокорреляции.
from statsmodels.stats.stattools import durbin_watson durbinWatson = durbin_watson(model.resid) print(durbinWatson)
››› Дурбин-Уотсон = 2.11
Статистика Дарбина-Ватсона находится в диапазоне от 0 до 4. Значение 2,0 означает, что автокорреляция отсутствует. Значения от 0 до 2 указывают на положительную автокорреляцию, а значения от 2 до 4 указывают на отрицательную автокорреляцию. В нашем случае статистика Дарбина-Ватсона очень близка к 2,0, поэтому мы можем сказать, что никакое предположение автокорреляции не нарушено.
Нет мультиколлинеарности
Мы предположили, что переменные-предикторы в регрессии не коррелируют друг с другом. Нарушение этого предположения не делает нашу модель непригодной для использования, но все же важно понимать и избавляться от некоторых ненужных предикторов, если они сильно коррелированы с другими.
Мы построим тепловую карту, чтобы визуально увидеть корреляцию между предикторами, а также рассчитаем коэффициент увеличения дисперсии (VIF), который измеряет мультиколлинеарность между независимыми переменными.
mask = np.zeros_like(train_X.drop("const", axis=1).corr())
mask[np.triu_indices_from(mask)] = True
sns.heatmap(train_X.drop("const", axis=1).corr().round(2), annot=True, mask=mask, cmap="cividis")
Как мы видим, между некоторыми переменными у нас есть высокая (отрицательная / положительная) корреляция в Boston Dataset.
И рассчитываем VIF:
from statsmodels.stats.outliers_influence import variance_inflation_factor as vif
for i in range(len(X.columns)):
v=vif(np.matrix(X),i)
print("Variance inflation factor for {}: {}".format(X.columns[i],round(v,2)))
Коэффициент инфляции дисперсии для CRIM: 1,79
Коэффициент инфляции дисперсии для ZN: 2,3
Коэффициент инфляции дисперсии для INDUS: 3,99
Коэффициент инфляции дисперсии для CHAS: 1,07
Коэффициент инфляции дисперсии для NOX: 4,39
Коэффициент увеличения дисперсии для RM: 1,93
Коэффициент увеличения дисперсии для AGE: 3,1
Коэффициент увеличения дисперсии для DIS: 3,96
Коэффициент увеличения дисперсии для RAD: 7,48
Дисперсия коэффициент инфляции для НАЛОГА: 9,01
Коэффициент инфляции дисперсии для PTRATIO: 1,8
Коэффициент инфляции дисперсии для B: 1,35
Коэффициент инфляции дисперсии для LSTAT: 2,94
VIF, равный 1, указывает, что две переменные не коррелированы, VIF между 1 и 5 указывает на умеренную корреляцию, а VIF выше 5 указывает на высокую корреляцию.
Как показывает результат, у нас есть 2 сильно коррелированных и 5 умеренно коррелированных переменных. Лучше, если мы избавимся от некоторых из них.
8. БОНУС: выбросы
Обнаружение выбросов с помощью графика расстояний Кука
Расстояние Кука определяет эффект удаления данного наблюдения из набора данных. Мы можем использовать дистанцию Кука при исследовании того, является ли наблюдение потенциальным выбросом или влиятельной переменной. Вот как построить расстояние Кука.
from statsmodels.stats.outliers_influence import OLSInfluence as influence
inf=influence(model)
(i, d) = inf.cooks_distance
plt.title("Cook's distance plot")
plt.stem(np.arange(len(i)), i, markerfmt=",")
plt.show()
Из сюжета Кука мы можем понять, на какие наблюдения нам нужно обратить больше внимания, и решить, отбросить их или нет. (Как правило, наблюдение имеет большое влияние, если расстояние Кука больше 4 / N-k-1 (N = количество наблюдений, k = количество предикторов, желтая горизонтальная линия на графике)
9. БОНУС: графики влияния.
Графики влияния показывают стьюдентизированные остатки в сравнении с левереджем каждого наблюдения, измеренными с помощью матрицы шляп.
from statsmodels.stats.outliers_influence import OLSInfluence influence = OLSInfluence(model) fig, ax = plt.subplots() ax.axhline(-2.5, linestyle='-', color='C1') ax.axhline(2.5, linestyle='-', color='C1') ax.scatter(influence.hat_matrix_diag, influence.resid_studentized_internal, s=1000 * np.sqrt(influence.cooks_distance[0]), alpha=0.5) ax.set_xlabel(‘H Leverage’) ax.set_ylabel(‘Studentized Residuals’) ax.set_title(“BOSTON DATASET \nInfluence Plot”) plt.tight_layout() plt.show()
В отличие от идеального графика, в Boston Dataset мы можем видеть, что некоторые точки данных имеют низкое кредитное плечо, но большие остатки. Эти очень влиятельные данные нуждаются в деликатной оценке перед включением в модель.
Заключение
В ходе этого анализа мы обнаружили, что некоторые из предположений, которые мы сделали перед построением модели линейной регрессии, были нарушены и вызвали значительные проблемы в производительности модели и прогнозировании результатов. Несмотря на то, что встреча с предположениями с реальными данными звучит немного поверхностно, нам все же необходимо проверить и узнать недостатки нашей модели и, если возможно, исправить или улучшить их, прежде чем принимать важные решения на основе результатов.
Если у вас есть такие результаты, как здесь, с набором данных Boston, я бы посоветовал вам попробовать моделировать с помощью Random Forest Regressor или XGBoost, чтобы понять, можете ли вы улучшить производительность модели. Вы можете обогатиться с помощью анализа главных компонентов - PCA для извлечения признаков, чтобы отбросить наименее важные предикторы, сохранив при этом наиболее ценные части всех переменных. Кроме того, вы также можете выполнить перекрестную проверку (CV) в K-кратном порядке, разделив весь набор данных на k складок и убедившись, что каждая кратность используется в качестве набора для тестирования.
Надеюсь, вам понравилось читать эту статью и вы сочтете ее полезной для анализа!
Если вам понравилась эта статья, вы можете прочитать другие мои статьи здесь и подписаться на меня на Medium. Дайте мне знать, если у вас есть вопросы или предложения.
Понравилась статья? Станьте участником и получите больше!
