Предыдущий пост см. здесь.

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

Матрицы

Матрица, — это двумерный массив чисел. Размерность матрицы выражается числом строк и столбцов.

Например, A — это матрица с четырьмя строками и двумя столбцами:

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

Массив numpy можно сконструировать из набора данных при помощи функции pandas df.values:

def ex_3_16():
    '''Конвертирование в массив (матрицу) numpy 
       таблицы данных роста и веса'''
    df = swimmer_data()[['Рост, см', 'Вес']]
    return df.values

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

array([[166.,  68.],
       [192.,  89.],
       [173.,  65.],
       ...,
       [188.,  79.],
       [187.,  78.],
       [183.,  70.]])

Можно также воспользоваться функцией из библиотеки numpy np.array, которая принимает последовательность скалярных величин либо последовательность последовательностей и, в случае если это возможно, конвертирует их в одномерный массив (в формате numpy.ndarray):

def ex_3_17():
    '''Конвертирование в массив (матрицу) numpy 
       данных числового ряда с данными о росте'''
    return swimmer_data()['Рост, см'].head(20).values # первые 20

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

array([166., 192., 173., 179., 201., 190., 175., 160., 202., 173., 175.,
       205., 185., 175., 185., 170., 165., 179., 165., 170.])

Матрицы часто могут вырастать до больших размеров, поэтому, чтобы не переполнять информацией окно интерпретатора, можно ограничить число выводимых элементов, воспользовавшись функциями pandas head и tail либо функционалом индексации библиотеки numpy (result_array[:5]); в обоих случаях будет выведено заданное число элементов.

^{Библиотека pandas взаимодействует с функциями библиотеки numpy напрямую, получая выгоду от векторизованных операций, таких как log, exp, sqrt и др., на массивах/матрицах. С другой стороны, различные функции numpy могут без проблем использоваться с кадрами данных DataFrame (и рядами Series) библиотеки pandas при условии, что содержащиеся внутри данные являются числовыми. Например, np.exp(df), np.asarray(df), df.T.dot(df)).}

Размерность

Элемент в i-ой строке и j-ом столбце обозначается как A_ij. И поэтому в приведенном выше примере индексация будет такой:

$A_{31}=2$

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

Векторы

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

Здесь y — это 4-мерный вектор; его i-й элемент обозначается как y_i. Векторы в математической литературе индексируются, начиная с единицы, если не указано иное.

Так, обозначает первый элемент, не второй. В уравнениях векторы в основном закрепляются за переменными, обозначаемыми строчными буквами.

^{Программный интерфейс библиотеки numpy в отличие от библиотеки pandas (где Series и DataFrame могут служить соответственно для представления векторов и матриц) не делает различие между векторами и одностолбцовыми матрицами, и мы можем создать вектор, передав в функцию np.array единственную последовательность.}

Сборка

Как мы уже убедились, матрицы можно собирать из последовательностей Python и наборов данных pandas. Кроме того, матрицы можно собирать из более мелких конструктивных составляющих, при условии совпадения размерностей, надстраивая столбцы бок о бок и добавляя строки. В простейшем случае мы можем добавить столбец единиц в начало или конец таблицы и затем преобразовать в матрицу следующим образом:

'''Добавление столбца в таблицу данных (массив)'''
df = pd.DataFrame({'x':[2, 3, 6, 7],'y':[8, 7, 4, 3]})
df['константа'] = 1
df

	x	y	константа
0	2	8	1
1	3	7	1
2	6	4	1
3	7	3	1

На самом деле, нам потребуется это делать для члена смещения. Напомним, что β₁ выражает константное значение, поэтому мы должны обеспечить, чтобы соответствующий x₁ тоже был константой. Без смещения переменная y равнялось бы нулю, в случае когда значения x равны нулю.

Сложение и скалярное произведение

Скаляр — это название для обыкновенного числа. Когда мы прибавляем скаляр в матрицу, мы на самом деле прибавляем это число отдельно в каждый элемент матрицы.

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

Имплементация на Python с использованием pandas:

df1 = pd.DataFrame([[1,0],[2,5],[3,1]])
df2 = pd.DataFrame([[4,0.5],[2,5],[0,1]])
df1 + df2

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

Имплементация на Python с использованием pandas:

df1 * 3

Матрично-векторное умножение

В функции dot с применением сложного алгоритма матричного умножения имплементирован стандартный способ умножения матриц. Например, результатом умножения матрицы размера 3 × 2 на матрицу размера 2 × 1 является матрица размера 3 × 1, при этом число столбцов слева должно совпадать с числом строк справа:

Для получения Ax надо помножить каждую строку A поэлементно с соответствующим элементом матрицы x и сложить результаты. Например, первая строка матрицы A содержит элементы 1 и 3. Они попарно умножаются на элементы в векторе x: 1 и 5. Затем, произведения чисел складываются и в результате получаем 16. Эта операция называется точечным произведением, или скалярным произведением, для чего, собственно, и предназначено матричное умножение.

Имплементация на Python с использованием pandas:

df3 = pd.DataFrame([[1,3],[0,4],[2,1]])
vec = [1,5]
df3.dot(vec)

0    16
1    20
2     7
dtype: int64

Матрично-матричное умножение

Выполнение матрично-матричного умножения матриц очень походит на матрично-векторное умножение. Из соответствующих элементов матриц A и B попарно, строка за строкой и столбец за столбцом, берется сумма произведений.

Как и ранее, мы можем умножать матрицы между собой, только когда число столбцов в первой матрице равно числу строк во второй. Если первая матрица A имеет размерность m_A × n_A, а вторая матрица B — размерность m_B × n_B, то для того чтобы их перемножить, n_A и m_B должны быть эквивалентными.

В приведенном выше наглядном примере:

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

df3 = pd.DataFrame([[1,3],[0,4],[2,1]])
df4 = pd.DataFrame([[1,0],[5,6]])     
df3.dot(df4)

Или то же самое в библиотеке numpy:

np.matmul(df3,np.asarray(df4))

Транспонирование

Транспонировать матрицу означает перевернуть матрицу по главной диагонали, проходящей из верхнего левого угла в нижний правый угол. Транспонирование матрицы A обозначается как A^T:

Столбцы и строки изменились таким образом, что:

Следовательно, если:

то:

Имплементация на Python с использованием pandas:

df3.T

	0	1	2
0	1	0	2
1	3	4	1

Подробнее об операциях с матрицами и векторами см. в моем репо на Gitgub для этой серии постов.

Нейтральная матрица

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

Единичная матрица — это нейтральная по умножению матрица (или нейтральный элемент для умножения). Как и при скалярном умножении на число 1, матричное умножение на нейтральную матрицу не имеет никакого эффекта.

Для получения нейтральной матрицы нам придется обратиться к функции библиотеке numpy np.identity, которая конструирует нейтральные матрицы. Учитывая, что они всегда квадратные, мы передаем всего один аргумент, который обозначает ширину и высоту одновременно.

Имплементация на Python с использованием pandas:

df = pd.DataFrame(np.identity(5))
df

	  0	  1	  2	  3	  4
0	1.0	0.0	0.0	0.0	0.0
1	0.0	1.0	0.0	0.0	0.0
2	0.0	0.0	1.0	0.0	0.0
3	0.0	0.0	0.0	1.0	0.0
4	0.0	0.0	0.0	0.0	1.0

Обратная матрица

Если мы имеем квадратную матрицу A, то обратная для A матрица обозначается как A^-¹, и она будет иметь следующие свойства, где I — это нейтральная матрица:

Нейтральная матрица является обратной самой себе. Не все матрицы обратимы. Необратимые матрицы также называются сингулярными или вырожденными. Обратная матрица вычисляется посредством функции np.linalg.pinv из модуля линейной алгебры библиотеки numpy.

Имплементация на Python с использованием pandas:

df5 = pd.DataFrame(np.random.rand(3, 3), list('abc'), list('xyz'))
print(df5)
df_inv = pd.DataFrame(np.linalg.pinv(df5.values), df5.columns, df5.index)
print(df_inv)

          x         y         z
a  0.625754  0.385261  0.462726
b  0.615084  0.111360  0.255420
c  0.723909  0.270869  0.221620
          a         b         c
x -1.451613  1.303231  1.528861
y  1.584699 -6.402303  4.070011
z  2.804750  3.568103 -5.456161

Нормальное уравнение

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

Мы читаем «для отыскания β надо умножить инверсию произведения транспонированной X и X, на произведение транспонированной X и y», где X — это матрица независимых переменных (включая член пересечения) и y — вектор, содержащий зависимые переменные нашей выборки. Результат β содержит вычисленные коэффициенты. Нормальное уравнение относительно легко выводится из уравнения множественной регрессии, применяя правила матричного умножения, однако соответствующие математические выкладки лежат за пределами объема данного поста.

Указанное нормальное уравнение можно имплементировать на Python, используя для этого только те функции, с которыми мы только что познакомились:

def normal_equation(x, y):
    '''Имплементация нормального уравнения'''
    # numpy.linalg.inv(A) фактически вызывает numpy.linalg.solve(A,I), 
    # где I - это нейтральная матрица, и находит решение разложением 
    # LU матрицы средствами динамической библиотеки lapack
    xtx  = np.matmul(x.T.values, x.values) 
    # вычислить мультипликативную инверсию матрицы
    xtxi = np.matmul(np.linalg.inv(np.matmul(xtx.T,xtx)),xtx.T)
    xty  = np.matmul(x.T.values, y.values) 
    return np.matmul(xtxi, xty)

Нормальное уравнение выражает математику регрессии методом обычных наименьших квадратов в очень сжатой форме. Его можно использовать следующим образом (помня о добавлении смещения):

def ex_3_18():
    '''Решение нормального уравнения 
       на примере данных роста и веса'''
    df = swimmer_data()
    X = df[['Рост, см']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log)
    return normal_equation(X, y)

В результате выполнения этого примера получим следующую матрицу:

array([ 1.69103131,  0.01429648])

Это значения представляют коэффициенты β₁ и β₂, которые соответствуют параметрам пересечения и наклона. К счастью, они согласуются со значениями, которые мы вычислили ранее.

Дополнительные признаки

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

^{В машинном усвоении закономерностей (Да-да, именно так. См. мои посты «}^{Никто никого не обучает}^{» и «}^{Что такое machine learning?}^{») в качестве синонима для независимой переменной широко используется понятие «признак», англ. feature. Другими синонимами являются «предсказатель», «предиктор», «регрессор», «объяснительная переменная», либо просто «входная переменная».}

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

def ex_3_19():
    '''Пример создания матрицы признаков NumPy
       на примере данных роста и возраста'''
    X = swimmer_data()[['Рост, см', 'Возраст']]
    X.insert(0, 'константа', 1)
    return X.values

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

array([[  1., 166.,  23.],
       [  1., 192.,  22.],
       [  1., 173.,  20.],
       ...,
       [  1., 188.,  24.],
       [  1., 187.,  19.],
       [  1., 183.,  22.]])

Наша функция нормального уравнения примет эту новую матрицу без какого-либо дальнейшего изменения:

def ex_3_20():
    '''Решение нормального уравнения 
       для данных роста и возраста в качестве независимых и 
       веса в качестве зависимой переменной'''
    df = swimmer_data()
    X = df[['Рост, см', 'Возраст']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log)
    return normal_equation(X, y)

В результате получим нижеследующие коэффициенты:

array([1.69002036, 0.01395437, 0.00279859])

Эти три числа соответствуют соответственно пересечению, наклону (угловому коэффициенту) для роста (0.013954) и наклону для возраста (0.002799). В целях установления факта улучшения нашей модели за счет этих новых данных можно рассчитать значение R² нашей новой модели и сравнить его с представленным ранее.

Множественный R-квадрат

При расчете R² в ранее рассмотренном случае мы увидели, каким образом он выражает объем дисперсии, объясненной моделью:

Учитывая, что дисперсия — это средневзвешенная квадратичная ошибка, мы можем умножить оба члена var(ε) и var(y) на размер выборки и прийти к приведенному ниже альтернативному уравнению для R²:

Это попросту сумма квадратичных остатков на сумме квадратичных отклонений от среднего. При помощи функций библиотеки pandas dot функция суммы квадратов имплементируется элементарно, что во многом упрощает имплементацию матричного R-квадрата в исходном коде:

def matrix_r_squared(coefs, x, y):
    '''Вычислить матричный R-квадрат'''
    fitted      = x.dot(coefs) 
    residuals   = y - fitted 
    difference  = y - y.mean()  
    rss         = residuals.dot(residuals)  # сумма квадратов
    ess         = difference.dot(difference)
    return 1 - (rss / ess)

Переменная rss обозначает остаточную сумму квадратов, от англ. residual sum of squares (RSS), переменная ess — объясненную сумму квадратов, от англ. explained sum of squares (ESS). Мы можем вычислить матричный R² для нашей новой модели следующим образом:

def ex_3_21():
    '''Вычислить матричный R-квадрат 
       на данных роста и возраста в качестве независимых и 
       веса в качестве зависимой переменной'''
    df = swimmer_data()
    X = df[['Рост, см', 'Возраст']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log)
    beta = normal_equation(X, y) 
    return matrix_r_squared(beta, X, y)

0.7568466547183842

В результате выполнения этого примера получим значение 0.757. Значение R² увеличилось на небольшую величину за счет включения возраста. Учитывая, что мы использовали несколько независимых переменных, R² теперь называется коэффициентом множественной детерминации.

Скорректированный матричный R-квадрат

Мотивом для добавления в регрессию больше независимых переменных может быть то, что наше значение R² всегда растет. Добавление новых независимых переменных не сделает предсказание зависимой переменной сложнее — если новая объяснительная переменная не имеет объяснительной силы, то ее коэффициент просто будет равен 0, и R² останется таким же, каким он был без этой независимой переменной.

Вместе с тем, добавление новой переменной не говорит нам об улучшении или ухудшении модели. Если мы хотим узнать о том, помогает ли новая переменная на самом деле генерировать более качественную подгонку к данным, мы можем воспользоваться скорректированным , нередко записываемым как R̅². В отличие от R², скорректированный R̅² будет расти, только если новая независимая переменная увеличивает R² больше, чем ожидается вследствие случайности:

def matrix_adj_r_squared(coefs, x, y):
    '''Вычислить скорректированный матричный R-квадрат'''
    r_squared = matrix_r_squared(coefs, x, y) 
    n = y.shape[0]  # строки
    p = coefs.shape[0]
    dec = lambda x: x-1
    return 1 - (1 - r_squared) * (dec(n) / dec(n-p))

Скорректированный R̅² зависит от двух дополнительных параметров, n и p, которые относятся соответственно к размеру выборки и числу модельных параметров:

def ex_3_22():
    '''Вычислить скорректированный матричный R-квадрат 
       на данных роста и возраста в качестве независимых и 
       веса в качестве зависимой переменной'''
    df = swimmer_data()
    X = df[['Рост, см', 'Возраст']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log)
    beta = normal_equation(X, y) 
    return matrix_adj_r_squared(beta, X, y)

0.7559934850858171

Этот пример возвращает значение 0.756. Оно по-прежнему крупнее изначальной модели, поэтому возраст определенно несет некую объяснительную силу.

Линейная модель в numpy и scipy

Хотя имплементация нашей собственной версии нормального уравнения и R̅² предоставляет ценную возможность познакомиться с матричной алгеброй, важно отметить, что библиотеки numpy и scipy предлагают соответственно функции np.linalg.lstsq и stats.linregress, которые делают все то, что мы рассмотрели и даже больше. В прилагаемых к данной серии постов примерах исходного кода имеется пример, в котором продемонстрирована работа этих функций.

К примеру, функция numpy np.linalg.lstsq ожидает вызова с аргументами x и y (в виде последовательностей либо матриц, в случае множественной регрессии). Указанная функция вернет коллекцию x, содержащую решение методом обычных наименьших квадратов, остатки residuals, эффективный ранг матрицы rank и сингулярные значения s. Мы воспользуемся этой функцией для написания простой обертки, которую будем использовать вместо собственной имплементации нормального уравнения. Наша имплементация функции линейной модели будет возвращать только коэффициенты модели и будет использоваться для расчета беты, в частности, в приведенном ниже F-тесте:

def linear_model(x, y):
    '''Обертка вокруг библиотечной функции 
       линейной регрессии наименьшими квадратами, 
       вместо собственной имплементации нормального 
       уравнения normal_equation'''
    return np.linalg.lstsq(x,y,rcond=-1)[0]

F-тест значимости модели

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

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

Здесь j — это некий индекс в векторе параметров за исключением пересечения, т.е. свободного члена. Вычисляемая нами F-статистика представляет собой отношение объясненной дисперсии на необъясненной (остаточной) дисперсии. Она может быть выражена через отношение средневзвешенного квадрата регрессионной модели, от англ. mean squared model (MSM) на средневзвешенной квадратичной ошибке, от англ. mean square error (MSE):

Средневзвешенный квадрат регрессионной модели (MSM) равен объясненной сумме квадратов (ESS) деленной на модельную степень свободы, где модельная степень свободы — это число параметров в модели за исключением свободного члена. Средневзвешенная квадратичная ошибка (MSE) равна остаточной сумме квадратов (RSS) деленной на остаточную степень свободы, где остаточная степень свободы — это размер выборки минус число модельных параметров.

После расчета F-статистики мы отыскиваем ее в F-распределении, параметризованном теми же двумя степенями свободы:

def f_test(fitted, x, y):
    '''F-тест коэффициентов регрессии'''
    difference = fitted - y.mean() 
    residuals  = y - fitted
    ess        = difference.dot(difference) # сумма квадратов
    rss        = residuals.dot(residuals)
    p          = x.shape[1]    # столбцы
    n          = y.shape[0]    # строки
    df1        = p - 1
    df2        = n - p
    msm        = ess / df1
    mse        = rss / df2
    f_stat     = msm / mse     # mse модели / mse остатков
    f_test     = 1-stats.f.cdf(f_stat, df1, df2) 
    return f_test

def ex_3_23():
    '''Проверка значимости модели на основе F-теста
       на примере данных роста, возраста и веса'''
    df = swimmer_data()
    X = df[['Рост, см', 'Возраст']]
    X.insert(0, 'константа', 1.0)
    y = df['Вес'].apply(np.log)
    beta = linear_model(X, y)    
    fittedvalues = np.dot(X,beta) 

    # проверка коэффициентов модели
    return ('F-тест', f_test(fittedvalues, X, y))

('F-тест', 1.1102230246251565e-16)

В результате проверки будет получено число 1.11x10e-16. Это ничтожно малое число, и, как следствие, можно быть уверенными в том, что модель значима.

Отметим, что при малых выборках F-тест количественно измеряет увеличивающуюся неопределенность, что линейная модель допустима. В условиях случайной выборки из пяти элементов, например, данные иногда едва показывают какую-либо линейную связь вообще, и F-тест трактует данные, как незначимые даже при 50%-ом интервале уверенности.

Категориальные и фиктивные переменные

Теперь мы могли бы попытаться включить в регрессионный анализ «Пол» в качестве признака, однако мы столкнемся с проблемой. Входные данные выражены не числом, а как «М» или «Ж». Это пример категориальной переменной: переменной, которая может принимать одно из конечного множества неупорядоченных и (обычно) нечисловых значений. Другими примерами категориальных переменных является вид спорта, в котором спортсмен специализируется, или отдельно взятое состязание, в котором он наиболее квалифицирован.

Обычные наименьшие квадраты опираются на числовое значение минимизируемого остаточного расстояния. Каким может быть расстояние между плаванием и легкой атлетикой в числовом выражении? Из этого может следовать, что включить категориальные переменные в наше уравнение регрессии невозможно.

^{Категориальные или номинальные переменные стоят в отдельном ряду по сравнению с непрерывными переменными, т.к. они не лежат на числовой оси. Иногда категории представлены цифрами, как например, почтовые индексы, при этом мы не должны исходить из того, что числовые категории с неизбежностью упорядочены либо что интервал между категориями равный.}

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

Когда такая категория, как вид спорта, принимает более двух значений, для каждого вида спорта можно внести независимую переменную. При этом создают переменную для плавания и еще одну для тяжелой атлетики и т.д. Значение для плавания будет равно 1 для пловцов и 0 — для остальных.

Поскольку половая принадлежность может оказаться для нашей регрессионной модели полезной объяснительной переменной, давайте преобразуем женский пол в 0 и мужской — в 1 и добавим производный столбец, который будет содержать фиктивную переменную.

Давайте рассчитаем значение R̅², чтобы засвидетельствовать наличие или отсутствие улучшения в этом скорректированном показателе:

def ex_3_25():
    '''Обработка категориальных признаков 
       (создание двоичной переменной)'''
    df = swimmer_data()
    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int) # строковое -> числовое

    X = df[['Рост, см', 'Возраст', 'бин_Пол']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log)  
    
    beta = linear_model(X, y) 
    return matrix_adj_r_squared(beta, X, y)

0.8082954905432824

В результате выполнения этого примера получим 0.809. С участием таких признаков, как рост, возраст и пол, мы успешно объяснили более 80% дисперсии в весе наших олимпийских пловцов.

Относительная мощность

На этом этапе, возможно, было бы целесообразным поинтересоваться, какой признак играет самую важную роль в объяснении наблюдавшегося веса: возраст, пол или рост? Мы могли бы воспользоваться нашим скорректированным R² и взглянуть, насколько его значение изменяется, но это потребовало бы от нас повторно рассчитывать регрессию для каждой переменной, которую мы хотим проверить.

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

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

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

def beta_weight(coefs, x, y):
    '''Вычисление относительного вклада каждого признака'''
    sdx = x.std()
    sdy = y.std()
    return [x / sdy * c for x,c in zip(sdx,coefs)]

def ex_3_26():
    '''Относительный вклад каждого признака в предсказании веса
       на примере данных роста, возраста и пола'''
    df = swimmer_data()
    # получить двоичное поле
    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int)
    X = df[['Рост, см', 'Возраст', 'бин_Пол']] 
    X.insert(0, 'константа', 1)
    y = df['Вес'].apply(np.log) 
    beta = linear_model(X, y) 
    res = beta_weight(beta, X, y)
    return res

В результате мы получим (округлено до трех десятичных знаков):

[0.0, 0.6501469135033348, 0.05842998157513067, 0.30387262631851747]

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

Коллинеарность

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

Например, в нашем распоряжении также имеется столбец «Дата рождения», и может возникнуть соблазн попытаться внести и его. Это дата, но мы легко могли бы конвертировать ее в число, подходящее для использования в регрессии. Это можно сделать, попросту взяв год из даты рождения, воспользовавшись для этого библиотечной функцией pandas pd.to_datetime:

'''Служебная функция приведения строкового 
   представления даты к типу DateTime и извлечение года'''
str_to_year = lambda x: pd.to_datetime(x).year

def ex_3_27():
    '''Относительный вклад признаков в предсказании веса
       с участием признака с датой (год)'''
    df = swimmer_data()
    df['бин_Пол'] = df['Пол'].map({'М': 1, 'Ж': 0}).astype(int) 
    df['Год рождения'] = df['Дата рождения'].map(str_to_year)
    X = df[['Рост, см', 'Возраст', 'бин_Пол', 'Год рождения']] 
    X.insert(0, 'константа', 1.0)
    y = df['Вес'].apply(np.log) 
    
    beta = linear_model(X, y) 
    return beta_weight(beta, X, y)

[-0.0,
 0.650070475196164,
 0.09580282723307212,
 0.3041431115029873,
 0.03769748899125406]

Новый признак «Год рождения» имеет бета-коэффициент всего 0.038, меньше веса признака «Возраст», который мы вычислили ранее. Однако, вес признака «Возраст» теперь показывает значение 0.096. Его относительная важность увеличилась более чем на 65%, поскольку мы добавили признак «Год рождения». Тот факт, что добавление нового признака изменило важность существующего признака, указывает на то, что имеется проблема.

Включив дополнительный параметр «Год рождения», мы непреднамеренно нарушили правило регрессионного оценивания. Посмотрим почему:

def ex_3_28():
    '''График коллинеарности возраста спортсменов и даты их рождения'''
    df = swimmer_data()
    df['Год рождения'] = df['Дата рождения'].map(str_to_year)
    xs = df['Возраст'].apply(jitter(0.5))
    ys = df['Год рождения']
    pd.DataFrame(np.array([xs,ys]).T).plot.scatter(0, 1, s=3, grid=True)
    plt.xlabel('Возраст')
    plt.ylabel('Год рождения')
    #saveplot('ex_3_28.png')
    plt.show()

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

Мультиколлинеарность

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

^{На практике независимые переменные часто в определенной мере коллинеарны. Обратите внимание, например, что возраст и рост или пол и рост сами коррелируют друг с другом. Серьезные ошибки в коэффициентах могут возникать, только когда это условие принимает крайние значения.}

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

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

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

Самый верный метод оценить степень мультиколлинеарности состоит в оценке регрессии каждой независимой переменной на всех других независимых переменных. Если любой R² из этих уравнений близок к 1.0, то имеется высокая мультиколлинеарность. На деле самый крупный из этих R² служит в качестве индикатора степени существующей мультиколлинеарности.

После выявления мультиколлинеарности имеется несколько способов решить эту проблему:

Увеличить размер выборки. Больше данных могут дать более точные оценки параметров с меньшими стандартными ошибками.
Совместить признаки в один признак. Если у вас несколько признаков, которые в сущности измеряют тот же атрибут, то следует найти способ их унификации в один признак.
Отбросить проблемную переменную (или переменные).
Ограничить уравнение предсказания. Коллинеарность влияет на коэффициенты модели, но результат может по-прежнему хорошо вписываться в данные.

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

«Возраст» R²= 0.1049, тогда как «Год рождения» R² = 0.1050.

Как и ожидалось, между двумя переменными практически нет разницы, и каждая из них объясняет порядка 10% дисперсии в весе. Поскольку год рождения объясняет чуть больше дисперсии, мы его оставим и отбросим признак «Возраст».

Примеры исходного кода для этого поста находятся в моем репо на Github. Все исходные данные взяты в репозитории автора книги.

В следующем коротком посте, посте №4, будет рассмотрен процесс предсказания.

Python, корреляция и регрессия: часть 3