Однофакторная и многофакторная линейные регрессии

Найти параметры \(a\) и \(b\) уравнения прямой \[y = a + b x,\]

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

Для определения коэффициентов a и b используется метод наименьших квадратов, при котором минимизируется сумма квадратов отклонений предсказанного значения \(\hat y_i = y(x_i)\) от соответствующего значения зависимой переменной \(y_i\) (Sum of Squared Errors - SSE): \[\text{SSE}(a,b)=\text{SS}_{res[iduals]}=\sum_{i=1}^N{\text{отклонение}_i}^2=\sum_{i=1}^N(y_i-f(x_i))^2=\sum_{i=1}^N(y_i-a-b\cdot x_i)^2\]

Решается задача \[\min_{a, b} \text{SSE}(a,b)\]

Пример. Рассмотрим зависимость объема продаж от затрат на рекламу по телевидению, не учитывая других факторов. Построим график статистической зависимости объема продаж в зависимости от затрат на ТВ рекламу.

Подключаем необходимые библиотеки:

import pandas as pd
import numpy as np
import scipy.stats 

import io
import requests

from matplotlib import pyplot as plt
import seaborn as sns
sns.set(font_scale=1.2, palette='Set2')

Загружаем данные:

url="https://classmech.ru/assets/files/advertising.csv"
s=requests.get(url).content
advert=pd.read_csv(io.StringIO(s.decode('utf-8')))
advert.head()

Вид таблицы с данными (первые пять строк):

	TV	Radio	Newspaper	Sales
230.1	37.8	69.2	22.1
44.5	39.3	45.1	10.4
17.2	45.9	69.3	12.0
151.5	41.3	58.5	16.5
180.8	10.8	58.4	17.9

Построим диаграмму рассеяния:

plt.scatter(advert['TV'], advert['Sales'])

Определим параметры линейной функции, аппроксимирующей статистические данные при помощи scipy.stats.linregress:

result = scipy.stats.linregress(advert['TV'], advert['Sales'])

Атрибуты результата (results):

result.slope = b - наклон прямой или коэффициент регрессии (Regression Coefficient)
result.intercept = a - начальное значение a = y(0)
result.rvalue = r коэффициент корреляции. \(R = r^2\) - коэффиицент детерминации, который показывает, в какой степени изменчивость одной переменной (y) обусловлена (детерминирована) влиянием другой переменной (x). Величина коэффициента детерминации равна 1, означает что функция идеально ложится на все точки — данные идеально скоррелированны.
result.pvalue - p-значение для проверки нулевой гипотезы о равенстве нулю наклона прямой \(H_0: b = 0\). По умолчанию альтернативной гипотезой является \(H_1: b \neq 0\).
result.stderr - стандартная ошибка наклона прямой b
result.intercept_stderr - стандартная ошибка начального значения a

plt.scatter(advert['TV'], advert['Sales'])

TV_x = np.linspace(min(advert['TV']),max(advert['TV']),100)
sales_predicted = result.intercept+result.slope*TV_x 

plt.plot(TV_x, sales_predicted, 'r--', lw = 3)
plt.show()

Найденное p-значение практически равно нулю

result.pvalue

>> 7.927911625322733e-74

поэтому нулевая гипотеза отвергается, следовательно между рекламой на ТВ и продажами есть статистическая связь.

Стандартные ошибки a и b позволяют оценить доверительный интервал для этих параметров \(a = \hat a \pm t_\alpha SE_a\) \[b = \hat b \pm t_\alpha SE_b\]

где \(\hat a\), \(\hat b\) – точечные оценки \(a\) и \(b\), \(SE_a\), \(SE_b\) - стандартные ошибки \(a\) и \(b\), \(t_{1-\alpha/2}\) - квантиль распределения Стьюдента уровня \(\alpha\) (двухсторонний) с n−2 степенями свободы.

# Используем распределение Стьюдента для оценки доверительных интервалов
# p - вероятность, df - число степеней свободы
tinv = lambda p, df: abs(scipy.stats.t.ppf(p/2, df))

ts = tinv(0.05, len(advert['Sales'])-2)
print('С уровнем доверия 95 %:')
print('a принадлежит интервалу от {:.3f} до {:.3f}'.format(result.intercept-ts*result.intercept_stderr, result.intercept+ts*result.intercept_stderr))
print('b принадлежит интервалу от {:.3f} до {:.3f}'.format(result.slope-ts*result.stderr, result.slope+ts*result.stderr))

Результат:

С уровнем доверия 95 %:
a принадлежит интервалу от 6.339 до 7.611
b принадлежит интервалу от 0.052 до 0.059

После построения доверительного интервала коэффициента корреляции, делается проверка на попадание нуля в этот интервал.

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

Если ноль не попал в доверительный интервал, то с высокой вероятностью в генеральной совокупности не может быть нулевого значения коэффициента корреляции, что означает: связь между переменными существует. В таком случае коэффициент корреляции является статистически значимым.

Мультилинейная регрессия

Множественной или мультилинейной называют линейную регрессию, в модели которой число независимых переменных две или более. Уравнение множественной линейной регрессии имеет вид: \[y = a_0 + b_1 x_1 + b_2 x_2 + \ldots + b_n x_n\]

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

import sklearn 
import sklearn.model_selection
import sklearn.linear_model 
import sklearn.metrics

Разобъем выборку на две части при помощи функции train_test_split:

выборка для определения коэффициентов модели (обучающая)
выборка для оценки работы модели (тестовая)

train, test = sklearn.model_selection.train_test_split(advert, test_size=0.2)
# Размер выборки для "обучения" модели 
# (определения коэффициентов мультилинейной регрессии)
print(train.shape)
# Размер выборки для последующего тестирования модели 
print(test.shape)

(160, 4)
(40, 4)

# Столбцы со свободными переменными - признаками
variables = ['TV', 'Radio', 'Newspaper']   
# Столбец со значениями функции - целевой признак
target    = 'Sales'  

# Значения независимых переменных обучающей выборки
X_train = train[variables]
# Значения независимых переменных тестовой выборки
X_test  = test[variables]

Для обзора данных построим отношения между парами переменных:

На диагоналях - плотности распределения
Выше диагонали - диаграммы рассеяния
Ниже диагонали - совместные плотности распределения

g = sns.PairGrid(advert[['TV',  'Radio', 'Newspaper', 'Sales']], diag_sharey=False, height=2)
g.map_lower(sns.kdeplot, alpha=0.5) 
g.map_upper(plt.scatter, alpha=0.5) 
g.map_diag(sns.kdeplot, lw=2, alpha=0.9, common_norm=True) 
g.add_legend()
plt.show()

Создаем модель

model = sklearn.linear_model.LinearRegression(fit_intercept=True)  

“Тренируем” ее на обучающей выборке (определяем коэффициенты)

model.fit(X_train, train[target])  

Свободный член мультилинейной функции

model.intercept_

>> 5.0293435246287075

Значения коэффициентов

model.coef_

>> array([ 0.05269666,  0.10909234, -0.00382442])

Оценим значения Y при помощи модели для значений независимых переменных тестового набора данных и сравним результаты с соответсвующими значениями тестовой набора

compareData  = pd.DataFrame({'Tested': test[target],'Predicted': model.predict(X_test)})
compareData = pd.concat([X_test, compareData],axis = 1)
# Добавим столбец с относительной погрешностью
compareData['Rel. error'] = np.round((compareData['Tested'] - compareData['Predicted'])*100/compareData['Tested'],1)
compareData.head()

TV	Radio	Newspaper	Tested	Predicted	Rel. error
4.1	    11.6	    5.7	    3.2	        6.489072	-102.8
75.1	35.0	    52.7	12.6	    12.603548	-0.0
168.4	7.1	        12.8	16.7	    14.629065	12.4
76.4	26.7	    22.3	11.8	    11.882850	-0.7
53.5	2.0	        21.4	8.1	        7.984957	1.4

Для оценки качества модели библиотека sklearn предлагает несколько метрик для оценки погрешностей построенной модели

Среднеквадратичная ошибка

\[MSE = \frac{1}{N} \sum_{i=1}^{N} (y_i - y(x_i))^2\]

sklearn.metrics.mean_squared_error(test[target], model.predict(X_test))

>> 3.0153747854144215

Коэффициент детерминации

sklearn.metrics.r2_score(test[target], model.predict(X_test))  

>> 0.9068127798062458

Коэффициент детерминации может быть определен, используя метод модели score

model.score(X_test,test[target])

>> 0.9068127798062458

Абсолютная средняя ошибка

\[\text{MAE}(y, \hat{y}) = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}}-1} \left| y_i - \hat{y}_i \right|.\]

sklearn.metrics.mean_absolute_error(test[target], model.predict(X_test))

>> 1.0890322497403002

Cредняя абсолютная процентная ошибка

\[\text{MAPE}(y, \hat{y}) = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}}-1} \frac{\left| y_i - \hat{y}_i \right|}{\max(\epsilon, \left| y_i \right|)}\]

sklearn.metrics.mean_absolute_percentage_error(test[target], model.predict(X_test))

>> 0.2162503623522544

Максимальная по модулю ошибка:

\[\text{Max Error}(y, y(x_i)) = \max(| y_i - y(x_i) |)\]

sklearn.metrics.max_error(test[target], model.predict(X_test))

>> 7.753015549113567

Средняя по модулю ошибка

\[\text{MAE}(y, \hat{y}) = \frac{1}{n_{\text{samples}}} \sum_{i=0}^{n_{\text{samples}}-1} \left| y_i - \hat{y}_i \right|.\]

Мультилинейная регрессия

Templates:

Error