При заданном наборе переменных-предикторов p и переменной отклика множественная линейная регрессия использует метод, известный как метод наименьших квадратов, для минимизации суммы квадратов остатков (RSS):
RSS = Σ(y i – ŷ i ) 2
куда:
- Σ : греческий символ, означающий сумму
- y i : Фактическое значение отклика для i -го наблюдения
- ŷ i : прогнозируемое значение отклика на основе модели множественной линейной регрессии.
Однако, когда переменные-предикторы сильно коррелированы, мультиколлинеарность может стать проблемой. Это может привести к тому, что оценки коэффициентов модели будут ненадежными и будут иметь высокую дисперсию.
Один из способов избежать этой проблемы — вместо этого использовать регрессию основных компонентов , которая находит M линейных комбинаций (известных как «главные компоненты») исходных p предикторов, а затем использует метод наименьших квадратов для подбора модели линейной регрессии, используя главные компоненты в качестве предикторов.
В этом руководстве представлен пошаговый пример того, как выполнить регрессию основных компонентов в Python.
Шаг 1: Импортируйте необходимые пакеты
Во-первых, мы импортируем необходимые пакеты для выполнения регрессии основных компонентов (PCR) в Python:
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn. preprocessing import scale
from sklearn import model_selection
from sklearn. model_selection import RepeatedKFold
from sklearn.model_selection import train_test_split
from sklearn. decomposition import PCA
from sklearn. linear_model import LinearRegression
from sklearn. metrics import mean_squared_error
Шаг 2: Загрузите данные
В этом примере мы будем использовать набор данных под названием mtcars , который содержит информацию о 33 различных автомобилях. Мы будем использовать hp в качестве переменной ответа и следующие переменные в качестве предикторов:
- миль на галлон
- дисп
- дрянь
- вес
- qсек
Следующий код показывает, как загрузить и просмотреть этот набор данных:
#define URL where data is located
url = "https://raw.githubusercontent.com/Statology/Python-Guides/main/mtcars.csv"
#read in data
data_full = pd.read_csv (url)
#select subset of data
data = data_full[["mpg", "disp", "drat", "wt", "qsec", "hp"]]
#view first six rows of data
data[0:6]
mpg disp drat wt qsec hp
0 21.0 160.0 3.90 2.620 16.46 110
1 21.0 160.0 3.90 2.875 17.02 110
2 22.8 108.0 3.85 2.320 18.61 93
3 21.4 258.0 3.08 3.215 19.44 110
4 18.7 360.0 3.15 3.440 17.02 175
5 18.1 225.0 2.76 3.460 20.22 105
Шаг 3: Подгонка модели ПЦР
Следующий код показывает, как подогнать модель PCR к этим данным. Обратите внимание на следующее:
- pca.fit_transform(scale(X)) : это сообщает Python, что каждая из переменных-предикторов должна быть масштабирована, чтобы иметь среднее значение 0 и стандартное отклонение 1. Это гарантирует, что никакая переменная-предиктор не будет чрезмерно влиять на модель, если это произойдет. измерять в разных единицах.
- cv = RepeatedKFold() : это говорит Python использовать перекрестную проверку k-fold для оценки производительности модели. Для этого примера мы выбираем k = 10 кратностей, повторенных 3 раза.
#define predictor and response variables
X = data[["mpg", "disp", "drat", "wt", "qsec"]]
y = data[["hp"]]
#scale predictor variables
pca = PCA()
X_reduced = pca. fit_transform ( scale (X))
#define cross validation method
cv = RepeatedKFold(n_splits= 10 , n_repeats= 3 , random_state= 1 )
regr = LinearRegression()
mse = []
# Calculate MSE with only the intercept
score = -1\*model_selection. cross_val_score (regr,
np.ones((len(X_reduced),1)), y, cv=cv,
scoring='neg_mean_squared_error').mean()
mse.append(score)
# Calculate MSE using cross-validation, adding one component at a time
for i in np.arange (1, 6):
score = -1\*model_selection. cross_val_score (regr,
X_reduced[:,:i], y, cv=cv, scoring='neg_mean_squared_error').mean()
mse.append(score)
# Plot cross-validation results
plt.plot (mse)
plt.xlabel('Number of Principal Components')
plt.ylabel('MSE')
plt.title('hp')
Регрессия основных компонентов в Python
На графике отображается количество основных компонентов по оси x и тестовая MSE (среднеквадратичная ошибка) по оси y.
Из графика видно, что тестовая MSE уменьшается при добавлении двух основных компонентов, но начинает увеличиваться при добавлении более двух основных компонентов.
Таким образом, оптимальная модель включает только первые две главные компоненты.
Мы также можем использовать следующий код для расчета процента дисперсии в переменной ответа, объясненной добавлением каждого основного компонента в модель:
np.cumsum (np.round (pca. explained_variance_ratio_ , decimals= 4 )\* 100 )
array([69.83, 89.35, 95.88, 98.95, 99.99])
Мы можем видеть следующее:
- Используя только первый главный компонент, мы можем объяснить 69,83% вариации переменной отклика.
- Добавляя второй главный компонент, мы можем объяснить 89,35% вариации переменной отклика.
Обратите внимание, что мы всегда сможем объяснить большую дисперсию, используя больше основных компонентов, но мы можем видеть, что добавление более двух основных компонентов фактически не увеличивает процент объясненной дисперсии намного.
Шаг 4: Используйте окончательную модель для прогнозирования
Мы можем использовать окончательную модель PCR с двумя основными компонентами, чтобы делать прогнозы по новым наблюдениям.
В следующем коде показано, как разделить исходный набор данных на набор для обучения и тестирования и использовать модель PCR с двумя основными компонентами для прогнозирования набора для тестирования.
#split the dataset into training (70%) and testing (30%) sets
X_train,X_test,y_train,y_test = train_test_split (X,y,test_size= 0.3 ,random_state= 0 )
#scale the training and testing data
X_reduced_train = pca. fit_transform ( scale (X_train))
X_reduced_test = pca. transform ( scale (X_test))[:,:1]
#train PCR model on training data
regr = LinearRegression()
regr. fit (X_reduced_train[:,:1], y_train)
#calculate RMSE
pred = regr. predict (X_reduced_test)
np.sqrt ( mean_squared_error (y_test, pred))
40.2096
Мы видим, что тестовый RMSE оказался равным 40,2096.Это среднее отклонение между прогнозируемым значением hp и наблюдаемым значением hp для наблюдений в тестовом наборе.
Полный код Python, используемый в этом примере, можно найти здесь .