Понимание множественной линейной регрессии.

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

Множественная линейная регрессия необходима, когда одной переменной недостаточно для создания хорошей модели и точных прогнозов.

Давайте начнем понимать это с набора данных о жилье….

Постановка задачи.

Рассмотрим компанию по недвижимости, у которой есть набор данных, содержащий цены на недвижимость в регионе Дели. Он хочет использовать данные для оптимизации продажных цен на недвижимость на основе таких важных факторов, как площадь, спальни, парковка и т. Д.

По сути, компания хочет -

• Чтобы определить переменные, влияющие на цены на жилье, например площадь, количество комнат, санузлов и т. д.
• Чтобы создать линейную модель, которая количественно связывает цены на жилье с такими переменными, как количество комнат, площадь, количество ванных комнат и т. Д.
• Чтобы знать точность модели, то есть насколько хорошо эти переменные могут предсказать цены на жилье.

Приступим к кодированию ...

Шаг 1. Чтение и понимание данных

Давайте сначала импортируем NumPy и Pandas и прочитаем набор данных жилья.

Импорт необходимых библиотек

# importing required libraries
import numpy as np
import pandas as pd

Чтение набора данных

housing = pd.read_csv("Housing.csv")

Отображение набора данных

# Check the head of the dataset
housing.head()

Шаг 2: Визуализация данных

А теперь давайте займемся самым важным шагом - пониманием данных.

• Если есть какая-то очевидная мультиколлинеарность, это первое место, где ее можно уловить.
• Здесь вы также сможете определить, имеют ли некоторые предикторы прямую связь с переменной результата.

Мы будем визуализировать наши данные с помощью matplotlib и seaborn.

import matplotlib.pyplot as plt
import seaborn as sns

Визуализация числовых переменных

Давайте сделаем парный график всех числовых переменных.

sns.pairplot(housing)
plt.show()

Визуализация категориальных переменных

Как вы могли заметить, есть также несколько категориальных переменных. Давайте построим диаграмму для некоторых из этих переменных.

plt.figure(figsize=(20, 12))
plt.subplot(2,3,1)
sns.boxplot(x = 'mainroad', y = 'price', data = housing)
plt.subplot(2,3,2)
sns.boxplot(x = 'guestroom', y = 'price', data = housing)
plt.subplot(2,3,3)
sns.boxplot(x = 'basement', y = 'price', data = housing)
plt.subplot(2,3,4)
sns.boxplot(x = 'hotwaterheating', y = 'price', data = housing)
plt.subplot(2,3,5)
sns.boxplot(x = 'airconditioning', y = 'price', data = housing)
plt.subplot(2,3,6)
sns.boxplot(x = 'furnishingstatus', y = 'price', data = housing)
plt.show()

Мы также можем визуализировать некоторые из этих категориальных характеристик параллельно, используя аргумент hue. Ниже приведен график для furnishingstatus с airconditioning в качестве оттенка.

plt.figure(figsize = (10, 5))
sns.boxplot(x = 'furnishingstatus', y = 'price', hue=airconditioning', data = housing)
plt.show()

Шаг 3: Подготовка данных

• Как видите, в нашем наборе данных много столбцов со значениями «Да» или «Нет».
• Но для того, чтобы соответствовать линии регрессии, нам потребуются числовые значения, а не строка. Следовательно, нам нужно преобразовать их в единицы и нули, где 1 - это «Да», а 0 - «Нет».

# List of variables to map
varlist =  ['mainroad', 'guestroom', 'basement', 'hotwaterheating', 'airconditioning', 'prefarea']
# Defining the map function
def binary_map(x):
    return x.map({'yes': 1, "no": 0})
# Applying the function to the housing list
housing[varlist] = housing[varlist].apply(binary_map)
# Check the housing dataframe now
housing.head()

Фиктивные переменные

Переменная furnishingstatus имеет три уровня. Нам также нужно преобразовать эти уровни в целые числа.

Для этого мы будем использовать что-то под названием dummy variables.

# Get the dummy variables for the feature 'furnishingstatus' and store it in a new variable - 'status'
status = pd.get_dummies(housing['furnishingstatus'])

Здесь мы создаем фиктивные переменные для переменной «Furnishingstatus».

# Check what the dataset 'status' looks like
status.head()

Теперь вам не нужны три столбца. Вы можете опустить столбец furnished, так как тип меблировки можно определить только по последним двум столбцам, где -

• 00 будет соответствовать furnished
• 01 будет соответствовать unfurnished
• 10 будет соответствовать semi-furnished

Это сделано для того, чтобы избежать эффектов избыточности и мультиколлинеарности.

# Let's drop the first column from status df using 'drop_first = True'
status = pd.get_dummies(housing['furnishingstatus'], drop_first = True)

здесь мы отбрасываем переменную снабженная

# Add the results to the original housing dataframe
housing = pd.concat([housing, status], axis = 1)

объединение новых фиктивных переменных с набором данных

# Drop 'furnishingstatus' as we have created the dummies for it
housing.drop(['furnishingstatus'], axis = 1, inplace = True)

Отбрасываем Furnishingstatus, поскольку мы уже создали для него фиктивные переменные.

housing.head()

Шаг 4: Разделение данных на наборы для обучения и тестирования

Как вы знаете, первый базовый шаг к регрессии - это разделение на поезд и тест.

from sklearn.model_selection import train_test_split
# We specify this so that the train and test data set always have the same rows, respectively
np.random.seed(0)
df_train, df_test = train_test_split(housing, train_size = 0.7, test_size = 0.3, random_state = 100)

Изменение масштаба функций

Здесь мы видим, что, за исключением area, все столбцы имеют небольшие целочисленные значения. Поэтому чрезвычайно важно масштабировать переменные так, чтобы они имели сопоставимый масштаб.

Если у нас нет сопоставимых шкал, то некоторые коэффициенты, полученные путем подбора регрессионной модели, могут быть очень большими или очень маленькими по сравнению с другими коэффициентами.

Это может сильно раздражать во время оценки модели. Поэтому рекомендуется использовать стандартизацию или нормализацию, чтобы единицы полученных коэффициентов находились в одной шкале.

Как мы знаем, есть два распространенных способа масштабирования:

1. Мин-макс масштабирование
2. Стандартизация (среднее-0, сигма-1)

На этот раз мы будем использовать масштабирование MinMax.

from sklearn.preprocessing import MinMaxScaler

Примените scaler () ко всем столбцам, кроме переменных «да-нет» и «фиктивные».

scaler = MinMaxScaler()
# Apply scaler() to all the columns except the 'yes-no' and 'dummy' variables
num_vars = ['area', 'bedrooms', 'bathrooms', 'stories', 'parking','price']
df_train[num_vars] = scaler.fit_transform(df_train[num_vars])

Давайте проверим коэффициенты корреляции, чтобы увидеть, какие переменные сильно коррелированы.

# Let's check the correlation coefficients to see which variables are highly correlated
plt.figure(figsize = (16, 10))
sns.heatmap(df_train.corr(), annot = True, cmap="YlGnBu")
plt.show()

Как вы могли заметить, area кажется наиболее коррелированным с price. Давайте посмотрим на график пары area против price.

plt.figure(figsize=[6,6])
plt.scatter(df_train.area, df_train.price)
plt.show()

Итак, мы выбираем area в качестве первой переменной и попытаемся подогнать под нее линию регрессии.

Разделение на наборы X и Y для построения модели

y_train = df_train.pop('price')
X_train = df_train

Шаг 5: Построение линейной модели

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

Помните, что в statsmodels вам нужно явно подогнать константу с помощью sm.add_constant(X), потому что, если мы не выполним этот шаг, statsmodels по умолчанию соответствует линии регрессии, проходящей через начало координат.

import statsmodels.api as sm
# Add a constant
X_train_lm = sm.add_constant(X_train[['area']])
# Create a first fitted model
lr = sm.OLS(y_train, X_train_lm).fit()

lr.params возвращает коэффициенты переменной.

# Check the parameters obtained
lr.params

# Let's visualise the data with a scatter plot and the fitted regression line
plt.scatter(X_train_lm.iloc[:, 1], y_train)
plt.plot(X_train_lm.iloc[:, 1], 0.127 + 0.462*X_train_lm.iloc[:, 1], 'r')
plt.show()

Распечатаем сводку модели линейной регрессии.

# Print a summary of the linear regression model obtained
print(lr.summary())

Полученное значение R-квадрата составляет 0.283.

Добавление еще одной переменной

Поскольку у нас так много переменных, мы явно можем добиться большего. Итак, давайте продолжим и добавим вторую наиболее коррелированную переменную, то есть bathrooms.

# Assign all the feature variables to X
X_train_lm = X_train[['area', 'bathrooms']]

импорт библиотеки статистики и установка OLS

# Build a linear model
import statsmodels.api as sm
X_train_lm = sm.add_constant(X_train_lm)
lr = sm.OLS(y_train, X_train_lm).fit()
lr.params

# Check the summary
print(lr.summary())

Мы явно улучшили модель, поскольку значение скорректированного R-квадрата, поскольку его значение выросло до 0.477 с 0.281. Давайте продолжим и добавим еще одну переменную bedrooms.

# Assign all the feature variables to X
X_train_lm = X_train[['area', 'bathrooms','bedrooms']]
# Build a linear model
import statsmodels.api as sm
X_train_lm = sm.add_constant(X_train_lm)
lr = sm.OLS(y_train, X_train_lm).fit()
lr.params

# Print the summary of the model
print(lr.summary())

Мы снова улучшили скорректированный R-квадрат. Теперь давайте продолжим и добавим все переменные функции.

Добавление всех переменных в модель

# Check all the columns of the dataframe
housing.columns

позволяет построить линейную модель

#Build a linear model
import statsmodels.api as sm
X_train_lm = sm.add_constant(X_train)
lr_1 = sm.OLS(y_train, X_train_lm).fit()
lr_1.params

print(lr_1.summary())

Глядя на p-значения, кажется, что некоторые из переменных не являются действительно значимыми (в присутствии других переменных).

Может, мы могли бы уронить немного?

Мы могли бы просто отбросить переменную с наивысшим несущественным значением p. Лучше было бы дополнить это информацией VIF.

Проверка VIF

Коэффициент инфляции дисперсии, или VIF, дает базовое количественное представление о том, насколько переменные характеристик коррелированы друг с другом. Это чрезвычайно важный параметр для тестирования нашей линейной модели. Формула для расчета VIF:

VIF = 1/1-R².

# Check for the VIF values of the feature variables. 
from statsmodels.stats.outliers_influence import variance_inflation_factor
# Create a dataframe that will contain the names of all the feature variables and their respective VIFs
vif = pd.DataFrame()
vif['Features'] = X_train.columns
vif['VIF'] = [variance_inflation_factor(X_train.values, i) for i in range(X_train.shape[1])]
vif['VIF'] = round(vif['VIF'], 2)
vif = vif.sort_values(by = "VIF", ascending = False)
vif

Обычно мы хотим, чтобы VIF был меньше 5. Поэтому явно есть некоторые переменные, которые нам нужно отбросить.

Почему мы вообще рассматриваем цифру 5 ??? Я объясню…

Предположим, VIF = 5. т.е.

1/1-R² = 5

1-R² = 1/5

1-R² = 0.2

R² = 0.8

Это означает, что любая переменная с показателем VIF выше 5 объясняет более 80% дисперсии данных.

В этом случае мы можем столкнуться с проблемами мультиколлинеарности.

Отсюда и цифра 5.

Удаление переменной и обновление модели

Как видно из сводки и фрейма данных VIF, некоторые переменные по-прежнему не имеют значения. Одна из этих переменных - semi-furnished, поскольку она имеет очень высокое значение p, равное 0.938. Давайте продолжим и отбросим эти переменные.

# Dropping highly correlated variables and insignificant variables
X = X_train.drop('semi-furnished', 1,)
# Build a third fitted model
X_train_lm = sm.add_constant(X)
lr_2 = sm.OLS(y_train, X_train_lm).fit()
# Print the summary of the model
print(lr_2.summary())

Давайте снова посчитаем оценку vif.

# Calculate the VIFs again for the new model
vif = pd.DataFrame()
vif['Features'] = X.columns
vif['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif['VIF'] = round(vif['VIF'], 2)
vif = vif.sort_values(by = "VIF", ascending = False)
vif

Удаление переменной и обновление модели

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

Как вы могли заметить, переменная bedroom имеет значительно высокий VIF (6.6) и высокое значение p (0.206). Следовательно, эта переменная не очень полезна, и ее следует отбросить.

# Dropping highly correlated variables and insignificant variables
X = X.drop('bedrooms', 1)
# Build a second fitted model
X_train_lm = sm.add_constant(X)
lr_3 = sm.OLS(y_train, X_train_lm).fit()
# Print the summary of the model
print(lr_3.summary())

# Calculate the VIFs again for the new model
vif = pd.DataFrame()
vif['Features'] = X.columns
vif['VIF'] = [variance_inflation_factor(X.values, i) for i in range(X.shape[1])]
vif['VIF'] = round(vif['VIF'], 2)
vif = vif.sort_values(by = "VIF", ascending = False)
vif

Удаление переменной и обновление модели

Как вы могли заметить, при падении semi-furnised VIF также уменьшился до mainroad, так что теперь он меньше 5.

Но из резюме мы все еще можем видеть, что некоторые из них имеют высокое значение p. basement, например, имеет p-значение 0,03. Мы также должны отбросить эту переменную.

X = X.drop('basement', 1)
# Build a fourth fitted model
X_train_lm = sm.add_constant(X)
lr_4 = sm.OLS(y_train, X_train_lm).fit()
lr_4.params

Как видите, значения VIF и p находятся в допустимом диапазоне. Поэтому мы делаем наши прогнозы, используя только эту модель.

Окончательная оценка R² для нашего набора данных о поездах составляет 0,676.

Шаг 7: Остаточный анализ данных поезда

Итак, теперь, чтобы проверить, нормально ли распределены члены ошибок (что на самом деле является одним из основных предположений линейной регрессии),

давайте построим гистограмму ошибок и посмотрим, как она выглядит.

y_train_price = lr_4.predict(X_train_lm)
# Plot the histogram of the error terms
fig = plt.figure()
sns.distplot((y_train - y_train_price), bins = 20)
fig.suptitle('Error Terms', fontsize = 20)   
               
plt.xlabel('Errors', fontsize = 18)                         

Шаг 8: Прогнозы с использованием окончательной модели

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

Применение масштабирования к тестовым наборам

num_vars = ['area', 'bedrooms', 'bathrooms', 'stories', 'parking','price']
df_test[num_vars] = scaler.transform(df_test[num_vars])

Разделение на X_test и y_test

y_test = df_test.pop('price')
X_test = df_test
# Adding constant variable to test dataframe
X_test_m4 = sm.add_constant(X_test)
# Creating X_test_m4 dataframe by dropping variables from X_test_m4
X_test_m4 = X_test_m4.drop(["bedrooms", "semi-furnished", "basement"], axis = 1)
# Making predictions using the fourth model
y_pred_m4 = lr_4.predict(X_test_m4)

Вычисление оценки R² для тестового набора данных

from sklearn.metrics import r2_score
r2_score(y_true=y_test,y_pred=y_pred_m4)

Окончательная оценка R² для нашего набора данных о поездах составляет 0,676.

Для нашего тестового набора данных мы получили 0,66 балла R².

Это означает, что наша модель также хорошо работает на тестовом наборе данных.

Шаг 9: Оценка модели

Давайте теперь построим график зависимости фактических значений от прогнозируемых.

# Plotting y_test and y_pred to understand the spread
fig = plt.figure()
plt.scatter(y_test, y_pred_m4)
fig.suptitle('y_test vs y_pred', fontsize = 20)  
           
plt.xlabel('y_test', fontsize = 18)  
                        
plt.ylabel('y_pred', fontsize = 16)

Мы видим, что уравнение нашей наиболее подходящей линии выглядит следующим образом:

𝑝𝑟𝑖𝑐𝑒=0.236×𝑎𝑟𝑒𝑎+0.202×𝑏𝑎𝑡ℎ𝑟𝑜𝑜𝑚𝑠+0.11×𝑠𝑡𝑜𝑟𝑖𝑒𝑠+0.05×𝑚𝑎𝑖𝑛𝑟𝑜𝑎𝑑+0.04×𝑔𝑢𝑒𝑠𝑡𝑟𝑜𝑜𝑚+0.0876×ℎ𝑜𝑡𝑤𝑎𝑡𝑒𝑟ℎ𝑒𝑎𝑡𝑖𝑛𝑔+0.0682×𝑎𝑖𝑟𝑐𝑜𝑛𝑑𝑖𝑡𝑖𝑜𝑛𝑖𝑛𝑔+0.0629×𝑝𝑎𝑟𝑘𝑖𝑛𝑔+0.0637×𝑝𝑟𝑒𝑓𝑎𝑟𝑒𝑎−0.0337×𝑢𝑛𝑓𝑢𝑟𝑛𝑖𝑠ℎ𝑒𝑑.

Примечание. Чтобы отбросить переменные, которые имеют низкую корреляцию с целевой переменной, мы можем рекурсивное исключение функций, предоставляемое sklearn.

Вывод:

Надеюсь, читатели получат интуитивное представление о множественной линейной регрессии.