Линейную регрессию можно рассматривать как один из основных алгоритмов, который долгое время использовался на бумаге, прежде чем его начали кодировать на компьютерах. Многие студенты, изучающие эконометрику, сталкиваются с этим алгоритмом в своих исследованиях. Это пример алгоритма Контролируемое машинное обучение.

АЛГОРИТМ ОБУЧЕНИЯ С УПРАВЛЕНИЕМ
Алгоритмы обучения с учителем — это тип алгоритмов, в которых у нас есть функции, а также результат набора данных поезда, и мы используем функции набора тестовых данных для создания переменных результата. для тестовых наборов данных.

Вас могут заинтересовать несколько задач, которые можно решить с помощью линейной регрессии:

  • Как далеко поднимется мяч, если вы бросите мяч с начальной скоростью v м/с (может быть, это требовалось до Ньютона)
  • Какая страховая премия может потребоваться в зависимости от возраста клиента
  • Какая прибыль будет получена, если вы добавите функцию № 1, 2, 3 … n для конкретного приложения

МАТЕМАТИКА РЕГРЕССИИ
Мы перейдем к шагам, необходимым для режима регрессии, но перед этим давайте разберемся с некоторыми математическими принципами регрессионных моделей.

Обычно у нас есть n связанных точек данных, для которых мы хотим предсказать. (Помните, что мы пытаемся оценить значения этих параметров β0 и β1 статистически), а не пытаемся решить эти уравнения для точных значений (помните, что есть 2 переменных и n уравнений, и очень мала вероятность того, что эти уравнения согласуются с одним решение).

Поэтому у нас есть:

Если мы перепишем это уравнение, мы придем к уравнению:

Теперь попытаемся найти значения β1 и β0, минимизирующие значение параметра $\epsilon$ вообще (статистически).
Для того же мы вычисляем значение Squared Sum Errors (или суммы квадратов ошибок), также известного как SSE.

Одно и то же уравнение можно экстраполировать на несколько переменных. Как

СБОР ДАННЫХ
Первым шагом модели линейной регрессии и любой модели в этом отношении является сбор данных. Здесь мы будем использовать набор данных о стоимости страхования, который доступен на Kaggle. Скачать можно здесь.



import pandas as pd
insurance_cost = pd.read_csv("Insurance_Cost.csv")
insurance_cost.info()
<class 'pandas.core.frame.DataFrame'>
RangeIndex: 1338 entries, 0 to 1337
Data columns (total 7 columns):
 #   Column    Non-Null Count  Dtype  
---  ------    --------------  -----  
 0   age       1338 non-null   int64  
 1   sex       1338 non-null   object 
 2   bmi       1338 non-null   float64
 3   children  1338 non-null   int64  
 4   smoker    1338 non-null   object 
 5   region    1338 non-null   object 
 6   charges   1338 non-null   float64
dtypes: float64(2), int64(2), object(3)
memory usage: 73.3+ KB

Из вышеприведенной таблицы можно сделать следующие выводы:

  • Все точки данных заполнены (иначе хотя бы одна из переменных будет иметь другой ненулевой счетчик)
  • Есть 4 числовые переменные: возраст, ИМТ, дети, расходы.
  • Есть 3 категориальные переменные: Пол, Курильщик, Регион.
    Мы можем выполнить некоторую предварительную обработку данных, к чему мы сейчас и подойдем.
    Кроме того, для этого набора данных требуется подробный EDA, который мы рассмотрим отдельно.

РАБОТА С КАТЕГОРИЧЕСКИМИ ПЕРЕМЕННЫМИ
Как вы видели выше, уравнение линейной регрессии не может принимать категориальные переменные и ожидает только числовую переменную. Однако мы можем преобразовать категориальные переменные в числовые, закодировав их.
Кодирование — это не что иное, как определение 0 или 1 для каждого уровня категориальных переменных.
Возьмем примеры из наших наборов данных:

  • Пол будет иметь два уровня — M и F, мы можем закодировать это как двоичную переменную 0 для мужчины и 1 для женщины.
  • Smoker снова будет иметь два уровня — Y или N, которые можно преобразовать в двоичную переменную.
  • Регион немного сложнее. Для региона мы развернем 0/1 для каждого региона
    Давайте разберемся с этим немного подробнее.
insurance_cost['smoker'].unique() ##Giving two levels of unique points
array(['yes', 'no'], dtype=object)
insurance_cost['sex'].unique() ##Giving two levels of unique points
array(['female', 'male'], dtype=object)
insurance_cost['region'].unique() ##Giving 4 levels of unique points
array(['southwest', 'southeast', 'northwest', 'northeast'], dtype=object)
pd.get_dummies(insurance_cost['region'])

Для других переменных есть намного более простой (и более быстрый метод). Мы можем использовать функцию np.where.

import numpy as np
insurance_cost['smoker']=np.where(insurance_cost['smoker']== 'yes',1,0)
insurance_cost['sex'] = np.where(insurance_cost['sex']=='female',1,0)
#Works very similar to the if-else condition in excel

Давайте перечислим переменные-предикторы в наших данных (X-переменные), чтобы понять список переменных, которые будут влиять на наш результат.

insurance_cost_revised= pd.get_dummies(insurance_cost,columns = ['region'],drop_first = True)

Здесь мы используем drop_first = True, потому что для K переменных мы хотим создать фиктивные переменные K-1.
Таким образом, по сути, мы уменьшаем количество категориальных переменных от K = K-1
Потому что для модели понятно, что если ни одна из этих переменных не истинна, то другое значение будет истинным. Одно из фундаментальных требований модели линейной регрессии заключается в том, что переменные должны быть независимы друг от друга.

Мы также можем покрыть доказательство моделей линейной регрессии, если это необходимо. Перед этим давайте сначала извлечем функции (переменные-предикторы или переменные x).

Charges = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])

Мы создали набор данных переменной Y, теперь нам нужно выбрать набор данных переменной X.

insurance_cost_revised.columns
X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
##You might have to install statsmodels package using !pip install statsmodels
import statsmodels.api as sm
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
type(Y_var)
pandas.core.frame.DataFrame
X_variables.head()

Y_var

Теперь давайте разделим наши данные на тестовый и обучающий наборы, для которых мы импортируем пакет train_test_split из библиотеки sklearn.model_selection.
sklearn, пожалуй, самый известный пакет для машинного обучения.
Мы также добавляем константу для .... CHECK WHY DO YOU ADD CONSTANT

from sklearn.model_selection import train_test_split
import statsmodels.api as sm
X_variables = sm.add_constant(X_variables)
train_X, test_X, train_y, test_y = train_test_split(X_variables,Y_var,
                                                   train_size = 0.80,
                                                   random_state = 121)
model_1 = sm.OLS(train_y,train_X).fit()
model_1.summary2()

ПЕРВЫЕ ШАГИ В ПРОВЕРКЕ МОДЕЛИ
Мы создали нашу первую модель. Теперь нам нужно понять результаты, давайте посмотрим, что означает каждый из этих результатов и каковы будут наши первые шаги.

Наше основное внимание здесь будет сосредоточено на средней таблице (элементы фиолетового цвета). Остальные элементы следующие:

  • Желтые и синие элементы: это описание модели.
  • Красные элементы: наиболее важные и используются для улучшения качества модели. Насколько хорошо модель предсказывает
  • Черные элементы: используются для проверки того, насколько хорошо модель соответствует статистическим предположениям OLS (Обычный метод наименьших квадратов).
  • Фиолетовые элементы: используются для определения независимых переменных, используемых в модели.

Выбор наиболее важных функций: при выборе функций моделей мы выполняем простой тест гипотезы:

Если вы помните предпосылку проверки гипотезы, когда p-значение ‹ 0,05, то мы отвергаем нулевую гипотезу на уровне значимости 5%. Теперь это означает, что если какая-либо переменная имеет P-Value ‹ 0,05, это означает, что в соответствии с приведенным выше тестом ее коэффициент НЕ РАВЕН нулю и, следовательно, должен присутствовать в нашей системе тестирования. Точно так же, если какая-либо переменная имеет P-значение > 0,05, мы можем сделать вывод, что коэффициент для этой переменной равен нулю и, следовательно, может быть исключен.

Мы выбираем переменную с наивысшим P-значением, удаляем ее, а затем видим влияние на модель, а затем выбираем следующую переменную с наивысшим P-значением (в новой модели) и продолжаем делать это до тех пор, пока все переменные не будут иметь P-значение ‹ 0,05.

Это называется ОБРАТНОЕ ИСКЛЮЧЕНИЕ.

Здесь мы видим, что P-значение, связанное с region North-west, является самым высоким, мы удаляем эту переменную и снова запускаем модель.

insurance_cost_revised = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'region_northwest']
insurance_cost_revised.head()

Пройдем испытание:

  • Разделение X_Variables и Y_variables (Вы можете избежать этого шага, выполнив следующие шаги для переменных X и Y, но я намеренно не делаю этого, скоро объясню то же самое)
  • Добавление постоянного члена
  • Разделение данных в Test и Train
  • Запуск модели OLS

СНОВА !!!

X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
X_variables = sm.add_constant(X_variables)
print(Y_var.head())
print(X_variables.head())
charges
0  16884.92400
1   1725.55230
2   4449.46200
3  21984.47061
4   3866.85520
   const  age  sex     bmi  children  smoker  region_southeast  \
0    1.0   19    1  27.900         0       1                 0   
1    1.0   18    0  33.770         1       0                 1   
2    1.0   28    0  33.000         3       0                 1   
3    1.0   33    0  22.705         0       0                 0   
4    1.0   32    0  28.880         0       0                 0   

   region_southwest  
0                 1  
1                 0  
2                 0  
3                 0  
4                 0
train_X, test_X, train_y, test_y = train_test_split(X_variables,Y_var,
                                                   train_size = 0.80,
                                                   random_state = 121)
model_2 = sm.OLS(train_y,train_X).fit()
model_2.summary2()

В модели очень мало изменений, теперь давайте запустим модель с удалением переменной sex, так как она имеет самое высокое значение P. Мы сделаем все за один раз, а не отдельные сегменты кодов. Помните, что у нашего insurance_cost_revised больше нет region_northwest, так как мы удалили его на первом шаге.

#Removing the sex variable
insurance_cost_revised = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'sex']
print(insurance_cost_revised.head())

#Splitting of X and Y variables
X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
X_variables = sm.add_constant(X_variables)
print(Y_var.head())
print(X_variables.head())

#Splitting of test and train cases
train_X, test_X, train_y, test_y = train_test_split(X_variables,Y_var,
                                                   train_size = 0.80,
                                                   random_state = 121)
#Running the model
model_3 = sm.OLS(train_y,train_X).fit()
model_3.summary2()

OUTPUT
age     bmi  children  smoker      charges  region_southeast  \
0   19  27.900         0       1  16884.92400                 0   
1   18  33.770         1       0   1725.55230                 1   
2   28  33.000         3       0   4449.46200                 1   
3   33  22.705         0       0  21984.47061                 0   
4   32  28.880         0       0   3866.85520                 0   

   region_southwest  
0                 1  
1                 0  
2                 0  
3                 0  
4                 0  
       charges
0  16884.92400
1   1725.55230
2   4449.46200
3  21984.47061
4   3866.85520
   const  age     bmi  children  smoker  region_southeast  region_southwest
0    1.0   19  27.900         0       1                 0                 1
1    1.0   18  33.770         1       0                 1                 0
2    1.0   28  33.000         3       0                 1                 0
3    1.0   33  22.705         0       0                 0                 0
4    1.0   32  28.880         0       0                 0                 0

Мы видим небольшое улучшение adjusted R-squared с 0,745 до 0,746. Давайте повторим описанные выше шаги для последней переменной с P > 0,05, которая равна region_southwest, и запустим модель.

#Removing the region_southwest variable
insurance_cost_revised = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'region_southwest']
print(insurance_cost_revised.head())

#Splitting of X and Y variables
X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
from sklearn.model_selection import train_test_split
import statsmodels.api as sm
X_variables = sm.add_constant(X_variables)
print(Y_var.head())
print(X_variables.head())

#Splitting of test and train cases
train_X, test_X, train_y, test_y = train_test_split(X_variables,Y_var,
                                                   train_size = 0.80,
                                                   random_state = 121)
#Running the model
model_4 = sm.OLS(train_y,train_X).fit()
model_4.summary2()
OUTPUT
age     bmi  children  smoker      charges  region_southeast
0   19  27.900         0       1  16884.92400                 0
1   18  33.770         1       0   1725.55230                 1
2   28  33.000         3       0   4449.46200                 1
3   33  22.705         0       0  21984.47061                 0
4   32  28.880         0       0   3866.85520                 0
       charges
0  16884.92400
1   1725.55230
2   4449.46200
3  21984.47061
4   3866.85520
   const  age     bmi  children  smoker  region_southeast
0    1.0   19  27.900         0       1                 0
1    1.0   18  33.770         1       0                 1
2    1.0   28  33.000         3       0                 1
3    1.0   33  22.705         0       0                 0
4    1.0   32  28.880         0       0                 0

Теперь мы видим интересную проблему: наша adjusted R-squared отпала, и мы получаем еще одну переменную region_southeast как переменную с P > 0,05. Давайте проделаем все шаги еще раз. После удаления переменной region_southeast

insurance_cost_revised = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'region_southeast']
print(insurance_cost_revised.head())


X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
X_variables = sm.add_constant(X_variables)
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
print(Y_var.head())
print(X_variables.head())

train_X,test_X,train_y,test_y = train_test_split(X_variables,
                                                Y_var,
                                                train_size = 0.80,
                                                random_state = 121)

model_5 = sm.OLS(train_y,train_X).fit()
model_5.summary2()
OUTPUT
age     bmi  children  smoker      charges
0   19  27.900         0       1  16884.92400
1   18  33.770         1       0   1725.55230
2   28  33.000         3       0   4449.46200
3   33  22.705         0       0  21984.47061
4   32  28.880         0       0   3866.85520
       charges
0  16884.92400
1   1725.55230
2   4449.46200
3  21984.47061
4   3866.85520
   const  age     bmi  children  smoker
0    1.0   19  27.900         0       1
1    1.0   18  33.770         1       0
2    1.0   28  33.000         3       0
3    1.0   33  22.705         0       0
4    1.0   32  28.880         0       0

Теперь, когда мы решили все проблемы, связанные с P-значениями модели, и выбрали наиболее важные функции:

  • Возраст
  • ИМТ
  • Дети
  • Курильщик

Оценка и улучшение модели с точки зрения статистических предположений

Нам необходимо оценить нашу модель с точки зрения точности математических предположений метода OLS (Обычный метод наименьших квадратов). Вот некоторые из проверок, которые мы проведем:

  • Нормальное распределение остатков (остатки, т. е. прогнозируемое значение — фактические значения должны быть распределены нормально)
  • Анализ выбросов (насколько изменяется зависимая переменная при удалении определенного наблюдения)
  • Автокорреляция (корреляция между независимыми переменными отсутствует)
  • Гомоскедастичность (означает, что дисперсия членов ошибок одинакова)

Мы проведем некоторые тесты и изменим нашу модель в соответствии с результатами теста. Это будет последний шаг в тонкой настройке нашей модели.

Нормальное распределение остатков

В структуре OLS предполагается, что остатки нормально распределены. Это вытекает из теории центральной предельной теоремы и проверки гипотез. Если остатки не распределены нормально, это означает, что взятая выборка является НЕ случайной, и поэтому мы не можем использовать структуру проверки гипотез, чтобы установить, какие переменные важны, а какие нет.
Однако это допущение о ненормальности не обязательно должно быть очень сильным, приближенное нормальное распределение остаточных переменных работает в реальном мире.

Подробнее о проверке гипотез читайте здесь:



Чтобы проверить нормальность, мы делаем график PP наших остатков по сравнению с идеальным нормальным распределением, если остатки лежат вдоль линии идеального нормального графика, мы можем сказать, что наши остатки являются нормальными (или почти нормальными).

import matplotlib.pyplot as plt
import seaborn as sns
insurance_data_resid = model_5.resid
probplot = sm.ProbPlot(insurance_data_resid)
plt.figure(figsize = (8,6))
probplot.ppplot(line = '45')
plt.title('Normal PP Plot of Regression Standardized Residuals')
plt.show()

insurance_cost_revised['charges'].hist()
plt.show()

Мы видим, что переменная charge сильно искажена, поэтому мы можем преобразовать переменную в log, а затем попробовать еще раз.

Y_var = np.log(Y_var)
Y_var.hist()
plt.show()

train_X,test_X,train_y,test_y = train_test_split(X_variables,
                                                Y_var,
                                                train_size = 0.80,
                                                random_state = 121)

model_6 = sm.OLS(train_y,train_X).fit()
model_6.summary2()

Этот прогон не только устранил проблему ненормального остатка, но и дал следующие преимущества:

  • Уменьшены значения AIC и BIC, что означает улучшение прогностической способности модели.
  • Увеличен скорректированный R-Square с 0,745 в первой модели до 0,758.
  • Имеет уменьшенное значение числа условия, что снижает коллинеарность (хотя есть и другие лучшие тесты для того же)
import matplotlib.pyplot as plt
import seaborn as sns
insurance_data_resid = model_6.resid
probplot = sm.ProbPlot(insurance_data_resid)
plt.figure(figsize = (8,6))
probplot.ppplot(line = '45')
plt.title('Normal PP Plot of Regression Standardized Residuals')
plt.show()

У нас это намного лучше, чем в исходной модели, остатки распределены почти нормально, а скорректированный R-квадрат улучшился до 0,758 с 0,745.

plt.figure(figsize = (10,7))
sns.heatmap(insurance_cost_revised.corr(),annot = True)

На приведенной выше диаграмме показана корреляция со сборами и другими числовыми переменными. Как и предсказывает наша модель (самое высокое значение $\beta$), курение имеет наибольшую корреляцию со страховыми взносами.

pred = model_6.predict(test_X)
from sklearn.metrics import r2_score
print("R2 Score",r2_score(test_y,pred))
OUTPUT
R2 Score 0.7735003467233176

Тест гомоскедастичности:

Гомоскедастичность измеряет распределение остаточных переменных, что означает, что дисперсия остаточных переменных не сильно увеличивается с увеличением значения зависимой переменной. Это можно сделать, проверив график остатков с подобранной переменной, если форма представляет собой перевернутую воронку, это означает, что дисперсия увеличивается (гетроскедастичность).
Другим тестом для того же является Durbin Watson, который должен быть между 1 и 2, мы можем видеть то же самое в результате сводки нашей модели, его значение равно 1,945.

Мы можем удалить гетероскедастичность (в отличие от гомоскедастичности), изменив независимые переменные (взяв логарифм или стандартизировав переменные). Мы бы хотели, чтобы Durbin Watson был близок к нулю и чтобы не было перевернутой воронки.

Еще одним тестом для выполнения теста на гомоскедастичность является тест white's.

def standardized(vals):
    return((vals-vals.mean())/vals.std())

plt.scatter(standardized(model_6.fittedvalues),standardized(model_6.resid))
plt.xlabel('Fitted Values')
plt.ylabel('Residual')
plt.show()

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

  • Преобразование переменной X
  • Преобразовать переменную Y (это мы уже сделали, когда брали бревно)

Подтвердим наше предположение о гетероскедастичности с помощью теста Уайта.

from statsmodels.stats.diagnostic import het_white
white_test = het_white(model_6.resid,model_6.model.exog)
labels = ['Test Statistics','Test Statistics p-Value','F-Statistics','F-Test p-value']
#print(dict(zip(labels,white_test)))
print(pd.DataFrame(dict(zip(labels,white_test)).items(),columns = ['Test Label','Test Value']))

OUTPUT
Test Label    Test Value
0          Test Statistics  1.149814e+02
1  Test Statistics p-Value  1.965648e-18
2             F-Statistics  9.779939e+00
3           F-Test p-value  1.485640e-19

В приведенном выше тесте оба значения P равны ‹ 0,05, что означает, что мы должны отклонить нулевую гипотезу о гомоскедастичности и, следовательно, наличие гетероскедастичности. Мы можем попытаться передать входные переменные и посмотреть, принесет ли это пользу.

Тем временем мы также построим наш прогноз против фактических значений Y, чтобы увидеть распределение. Кроме того, мы можем видеть, что значение R-Square из прогноза составляет 0,773, что близко к фактической модели, что означает, что наша модель также хорошо справляется с фактическим прогнозом.

pred = model_6.predict(test_X)
from sklearn.metrics import r2_score
print("R2 Score",r2_score(test_y,pred))
plt.scatter(test_y,pred)
plt.xlabel('Y_Test')
plt.ylabel('Y_Predicted')
plt.show()
OUTPUT
R2 Score 0.7735003467233176

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

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

insurance_cost_revised[['bmi','age']].hist()

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

X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
X_variables = sm.add_constant(X_variables)
X_variables['bmi'] = np.log(X_variables['bmi'])
X_variables['age'] = np.log(X_variables['age'])
X_variables['children'] = standardized(X_variables['children'])
X_variables['smoker'] = standardized(X_variables['smoker'])

Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
Y_var = np.log(Y_var)
print(Y_var.head())
print(X_variables.head())

train_X,test_X,train_y,test_y = train_test_split(X_variables,
                                                Y_var,
                                                train_size = 0.80,
                                                random_state = 121)

model_7 = sm.OLS(train_y,train_X).fit()
model_7.summary2()
OUTPUT
charges
0  9.734176
1  7.453302
2  8.400538
3  9.998092
4  8.260197
   const       age       bmi  children    smoker
0    1.0  2.944439  3.328627 -0.908274  1.969850
1    1.0  2.890372  3.519573 -0.078738 -0.507273
2    1.0  3.332205  3.496508  1.580335 -0.507273
3    1.0  3.496508  3.122585 -0.908274 -0.507273
4    1.0  3.465736  3.363149 -0.908274 -0.507273

Мы достигли наивысшего значения R-квадрата, но давайте еще раз проведем тест на гомоскадистичность и нормальность остатков.
Хотя гомоскедастичность еще предстоит оценить, мы достигли более высокого значения R-квадрата, равного 0,761, что является хорошим знаком.

import matplotlib.pyplot as plt
import seaborn as sns
insurance_data_resid = model_7.resid
probplot = sm.ProbPlot(insurance_data_resid)
plt.figure(figsize = (8,6))
probplot.ppplot(line = '45')
plt.title('Normal PP Plot of Regression Standardized Residuals')
plt.show()

plt.scatter(standardized(model_7.fittedvalues),standardized(model_7.resid))
plt.xlabel('Fitted Values')
plt.ylabel('Residual')
plt.show()

from statsmodels.stats.diagnostic import het_white
white_test = het_white(model_7.resid,model_7.model.exog)
labels = ['Test Statistics','Test Statistics p-Value','F-Statistics','F-Test p-value']
#print(dict(zip(labels,white_test)))
print(pd.DataFrame(dict(zip(labels,white_test)).items(),columns = ['Test Label','Test Value']))
OUTPUT
Test Label    Test Value
0          Test Statistics  1.048361e+02
1  Test Statistics p-Value  1.906163e-16
2             F-Statistics  8.823286e+00
3           F-Test p-value  2.391576e-17

Как мы видим в приведенной выше таблице, значения p меньше 0,05, что означает, что мы должны отклонить нулевую гипотезу о наличии гомоскедастичности. Что означает наличие гетероскедастичности. Таким образом, мы не можем удалить гетероскедастичность из данных. Поэтому мы будем работать с Model_7 с улучшенным значением R-Squared.

Давайте проверим значения предсказанных значений Y, чтобы понять точность модели.

pred = model_7.predict(test_X)
from sklearn.metrics import r2_score
print("R2 Score",r2_score(test_y,pred))
plt.scatter(test_y,pred)
plt.xlabel('Y_Test')
plt.ylabel('Y_Predicted')
plt.show()
OUTPUT
R2 Score 0.7682519851617251

Это уменьшило значение R-Square для значений прогноза, поэтому мы можем использовать Model_6 в качестве окончательной модели. Мы запустим исправленную модель, чтобы переменные Train и Test приняли правильную форму.

insurance_cost_revised = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'region_southeast']
print(insurance_cost_revised.head())


X_variables = insurance_cost_revised.loc[:,insurance_cost_revised.columns != 'charges']
X_variables = sm.add_constant(X_variables)
Y_var = pd.DataFrame(insurance_cost_revised['charges'],columns = ['charges'])
print(Y_var.head())
print(X_variables.head())
train_X,test_X,train_y,test_y = train_test_split(X_variables,
                                                Y_var,
                                                train_size = 0.80,
                                                random_state = 121)

model_6 = sm.OLS(train_y,train_X).fit()
model_6.summary2()
age     bmi  children  smoker      charges
0   19  27.900         0       1  16884.92400
1   18  33.770         1       0   1725.55230
2   28  33.000         3       0   4449.46200
3   33  22.705         0       0  21984.47061
4   32  28.880         0       0   3866.85520
       charges
0  16884.92400
1   1725.55230
2   4449.46200
3  21984.47061
4   3866.85520
   const  age     bmi  children  smoker
0    1.0   19  27.900         0       1
1    1.0   18  33.770         1       0
2    1.0   28  33.000         3       0
3    1.0   33  22.705         0       0
4    1.0   32  28.880         0       0

Тестирование мультиколинеарности

Чтобы проверить мультиколиарность, мы можем проверить VIF, чтобы проверить мультиколинеарность. Значение VIF › 4 означает мультиколлинеарность.
Кроме того, номер условия показывает в модели 6 выше, что мультиколинеарность уменьшилась по сравнению с номером условия = 316.

from statsmodels.stats.outliers_influence import variance_inflation_factor
def get_vif_factors(X):
    X_matrix = np.matrix(X_variables)
    vif = [variance_inflation_factor(X_matrix,i) for i in range(X_matrix.shape[1])]
    vif_factors = pd.DataFrame()
    vif_factors['column'] = X_variables.columns
    vif_factors['VIF'] = vif
    return vif_factors
vif_factors  = get_vif_factors(X_variables)
vif_factors

Идентификация любой критической позиции

Если есть какое-либо очень влиятельное наблюдение, наша модель может быть искажена этой точкой данных. Проверяем это по:

  • Кредитное плечо
  • Расстояние повара

Очень большие круги на графике кредитного плеча означают очень важную точку данных (и нам, возможно, придется их удалить). Точно так же расстояние повара › 1 означает наличие важных точек данных. Мы делаем оба этих теста в следующих сегментах.

from statsmodels.graphics.regressionplots import influence_plot
fig,ax = plt.subplots(figsize = (10,10))
influence_plot(model_7,ax = ax)
plt.title("Figure - Leverage Values vs Residual")
plt.show()

Вычисление расстояния повара

insurance_influence = model_7.get_influence()
(c,p) = insurance_influence.cooks_distance
plt.stem(np.arange(len(train_X)),
        np.round(c,3),
        markerfmt = ',')
plt.title('Cooks Distance for Insurance Data')
plt.xlabel('Row Index')
plt.ylabel('Cooks Distance')
plt.show()

Подведение итогов

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

  • Скачать данные
  • Выполните EDA для данных (не рассматривается в этой записной книжке)
  • Подготовьте данные
  • Преобразование категориальных переменных в фиктивные переменные (K-1 для K уровней категорий в каждой категориальной переменной)
  • Запустите модель и попробуйте удалить ненужные переменные с помощью Backward Elimination.
  • Протестируйте модели на основе R-квадрата модели на тестовых переменных и переменных обучения.
  • Проверка статистической модели
  • Распределение остатков
  • Гомоскедастичность
  • Мультиколинеарность
  • Выявление и устранение критических замечаний

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

Это мой второй пост на этой платформе, ищу ваши честные отзывы о том же.

Надеюсь, это поможет, если в приведенных выше шагах/наблюдениях есть какие-либо ошибки, сообщите мне об этом в комментариях. Изображения взяты из данных Python, а остальные изображения созданы вручную.