2.2. Множественная регрессия

Обобщением линейной регрессионной модели с двумя переменными является многомерная регрессионная модель (или модель множественной регрессии). Пусть n раз измерены значения факторов x₁ , x₂ , ..., x_k и соответствующие значения переменной y; предполагается, что

y_i = _o + ₁x_i1 + ... + _k x_ik+ _i , i = 1, ..., n, (12)

(второй индекс у х относится к номеру фактора, а первый - к номеру наблюдения); предполагается также, что

M_i = 0, M =  ²,

M(_i _j) = 0, i  j, (12a)

т.е. _i - некоррелированные случайные величины . Соотношения (12) удобно записывать в матричной форме:

Y = X +  , (13)

где Y = (y₁, ..., y_k)^T - вектор-столбец значений зависимой переменной, Т - символ транспонирования,  = (₀, ₁, ..., _k)^T - вектор-столбец (размерности k) неизвестных коэффициентов регрессии,  = (₁ , ..., _n)^T - вектор случайных отклонений,

-матрица n (k + 1); в i - й строке (1, x_i1, ...,x_ik) находятся значения независимых переменных в i-м наблюдении первая переменная - константа, равная 1.

Оценка коэффициентов регрессии. Построим оценку для вектора  так, чтобы вектор оценок = Х зависимой переменной минимально (в смысле квадрата нормы разности) отличался от вектора Y заданных значений:

по .

Решением является (если ранг матрицы Х равен k +1) оценка

= (X^TX)^-1 X^TY (14)

Нетрудно проверить, что она несмещенная. Ковариационная (дисперсионная) матрица равна

D = (  ) (  )^T =  ² (X^TX)^1 =  ² Z , (15)

где обозначено Z = (X^TX)^1.

Справедлива

теорема Гаусса - Маркова. В условиях (12а) оценка (14) является наилучшей (в смысле минимума дисперсии) оценкой в классе линейных несмещенных оценок.

Оценка дисперсии  ² ошибок. Обозначим

e = Y  = Y  Х = [I  X (X^TX)^1 X^T] Y = BY (16)

вектор остатков (или невязок); B = I  X (X^TX)^1 X^T - матрица; можно проверить, что B² = B. Для остаточной суммы квадратов справедливо соотношение

M = M (n - k -1)  ² ,

откуда следует, что несмещенной оценкой для  ² является

s² = . (17)

Если предположить, что _i в (12) нормально распределены, то справедливы следующие свойства оценок:

1) (n - k - 1) имеет распределение хи квадрат с n-k-1 степенями свободы;

2. оценки и s² независимы.

Как и в случае простой регрессии, справедливо соотношение:

или

T_ss = E_ss + R_ss , (18)

в векторном виде:

,

где = ( . Поделив обе части на полную вариацию игреков

T_ss = , получим коэффициент детерминации

R² = (19)

Коэффициент R² показывает качество подгонки регрессионной модели к наблюдённым значениям y_i. Если R² = 0, то регрессия Y на x₁ , ..., x_k не улучшает качество предсказания y_i по сравнению с тривиальным предсказанием . Другой крайний случай R² = 1 означает точную подгонку: все e_i = 0, т.е. все точки наблюдений лежат на регрессионной плоскости. Однако, значение R² возрастает с ростом числа переменных (регрессоров) в регрессии, что не означает улучшения качества предсказания, и потому вводится скорректированный (adjusted) коэффициент детерминации

(20)

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

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

s_j = , j = 0, 1, ..., k, (21)

где z_jj- диагональный элемент матрицы Z. Если ошибки _i распределены нормально, то, в силу свойств 1) и 2), приведенных выше, статистика

(22)

распределена по закону Стьюдента с (n - k - 1) степенями свободы, и потому неравенство

 t_p s_j , (23)

где t_p - квантиль уровня (1 + P_Д) / 2 этого распределения, задает доверительный интервал для _j с уровнем доверия Р_Д.

Проверка гипотезы о нулевых значениях коэффициентов регрессии. Для проверки гипотезы Н₀ об отсутствии какой бы то ни было линейной связи между y и совокупностью факторов, Н₀: ₁ = ₂ = ... = _k = 0, т.е. об одновременном равенстве нулю всех коэффициентов, кроме коэффициента ₀ при константе, используется статистика

F = = = , (24)

распределенная, если Н₀ верна, по закону Фишера с k и n - k - 1 степенями свободы. Н₀ отклоняется, если

F > F_ (k, n - k - 1), (25)

где F_ - квантиль уровня 1 - .

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

Пример [5]. Исследуется зависимость урожайности y зерновых культур ( ц/га ) от ряда факторов (переменных) сельскохозяйственного производства, а именно,

х₁ - число тракторов на 100 га;

х₂ - число зерноуборочных комбайнов на 100 га;

х₃ - число орудий поверхностной обработки почвы на 100 га;

х₄ - количество удобрений, расходуемых на гектар (т/га);

х₅ - количество химических средств защиты растений, расходуемых на гектар (ц/га).

Исходные данные для 20 районов области приведены в табл. 2.

Таблица 2

	y	x1	x 2	x 3	x 4	x 5
1	9.7	1.59	.26	2.05	.32	.14
2	8.4	.34	.28	.46	.59	.66
3	9.0	2.53	.31	2.46	.30	.31
4	9.9	4.63	.40	6.44	.43	.59
5	9.6	2.16	.26	2.16	.39	.16
6	8.6	2.16	.30	2.69	.32	.17
7	12.5	.68	.29	.73	.42	.23
8	7.6	.35	.26	.42	.21	.08
9	6.9	.52	.24	.49	.20	.08
10	13.5	3.42	.31	3.02	1.37	.73
11	9.7	1.78	.30	3.19	.73	.17
12	10.7	2.40	.32	3.30	.25	.14
13	12.1	9.36	.40	11.51	.39	.38
14	9.7	1.72	.28	2.26	.82	.17
15	7.0	.59	.29	.60	.13	.35
16	7.2	.28	.26	.30	.09	.15
17	8.2	1.64	.29	1.44	.20	.08
18	8.4	.09	.22	.05	.43	.20
19	13.1	.08	.25	.03	.73	.20
20	8.7	1.36	.26	.17	.99	.42

Здесь мы располагаем выборкой объема n = 20; число независимых переменных (факторов) k = 5. Матрица Х должна содержать 6 столбцов размерности 20; первый столбец состоит из единиц, а столбцы со 2-го по 6-й представлены соответственно столбцами 37 таблицы (файл Harvest 2. sta.). Специальный анализ (здесь не приводимый) технологии сбора исходных данных показал, что допущения (12а) могут быть приняты в качестве рабочей гипотезы, поэтому можем записать уравнения статистической связи между y_i и X_i = (x_i1, x_i2, ..., x_i5), i = 1, ..., n в виде (13).

Выполнение в пакете STATISTICA

Работаем в модуле Multiple Regression (множественная регрессия).

Ввод данных. Образуем таблицу 6v  20c с 6 столбцами (variables - переменными) и 20 строками (cases). Столбцы назовем y, x₁, x₂ , ..., x₅ . Введем в таблицу исходные данные.

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

Graphs - Stats 2D Graphs - Scatterplots - Variables - X: x1, Y: y, Graph Type: Regular, Fit (подбор): Linear - OK.

Наблюдаем диаграмму рассеяния с подобранной прямой парной регрессии, параметры которой отражены в заголовке. Повторим это еще 4 раза, заменяя х₁ на другие факторы: х₂ , ..., х₅ . Иногда такой просмотр позволяет увидеть основную зависимость. В нашем примере этого нет.

Выполнение регрессионного анализа:

Analysis - Startup Panel - кнопка Variables: - отбираем зависимую переменную Dependent var: y и независимые переменные Independent var: x₁  x₅ (при нажатой клавише Ctrl) - OK - Input file (входной файл): Raw Data (необработанные файлы) - OK - в окне Model Definition (уточнения) Metod: Standart, Intercept: Include in model (постоянную составляющую включить в модель) - ОК..

В окне Mult. Regr. Results имеем основные результаты: коэффициент детерминации (19) R² = 0.517; для проверки гипотезы Н₀ об отсутствии какой бы то ни было линейной связи между переменной y и совокупностью факторов определена статистика (24) F = 3.00; это значение соответствует уровню значимости р = 0.048 (эквивалент (25) согласно распределению F (5,14) Фишера с df = 5 и 14 степенями свободы. поскольку значение р весьма мало, гипотеза Н₀ отклоняется.

Кнопка Regression summary - имеем таблицу результатов:

Regression Summary for Dependent Variable: Y
R = .71923865 RІ = .51730424 Adjusted RІ = .34491290 F(5,14) = 3.0008 p<.04787 Std. Error of estimate: 1.5990
	B	St. Err of B	t(14)	p-level
Intercpt	3.51460	5.41853	.648625	.527078
X1	-.00613	.93167	-.006580	.994843
X2	15.54246	21.50311	.722800	.481704
X3	.10990	.83254	.132004	.896859
X4	4.47458	1.54345	2.899065	.011664
X5	-2.93251	3.08833	-.949546	.358448

В ее заголовке повторены результаты предыдущего окна; в столбце В указаны оценки неизвестных коэффициентов по (14). Таким образом, оценка (x) неизвестной функции регрессии f (x) в данном случае:

(x) = 3.51  0.06 x₁ + 15.5 x₂ + 0.11 x₃ + 4.47 x₄  2.93 x₅ (26)

В столбце St. Err. of B указаны стандартные ошибки s_j оценок коэффициентов (по (21)); видно, что стандартные ошибки в оценке всех коэффициентов, кроме ₄ , превышают значения самих коэффициентов, что говорит о статистической ненадежности последних. В столбце t(14) -значение статистики Стьюдента (22) для проверки гипотезы о нулевом значении соответствующих коэффициентов; в столбце p-level -уровень значимости отклонения этой гипотезы; достаточно малым (0.01) этот уровень является только для коэффициента при x₄ . Только переменная x₄ - количество удобрений, подтвердила свое право на включение в модель. В то же время проверка гипотезы об отсутствии какой бы то ни было линейной связи между y и (х₁ , ..., х₅) с помощью статистики (24) (об этом сказано выше)

F = 3.00 , p = 0.048 ,

говорит о том, что следует продолжить изучение линейной связи между y и (х₁ , ..., х₅), анализируя как их содержательный смысл, так и матрицу парных корреляций, которая определяется так:

возврат в окно Multi. Regr. Results - кнопка Correlations and desc. Stats - Correlations. Из матрицы видно, что х₁ , х₂ и х₃ (оснащенность техникой)

Correlations (harvest2.sta)
	X1	X2	X3	X4	X5	Y
X1	1.000	.854	.978	.110	.341	.430
X2	.854	1.000	.882	.027	.460	.374
X3	.978	.882	1.000	.030	.278	.403
X4	.110	.027	.030	1.000	.571	.577
X5	.341	.460	.278	.571	1.000	.332
Y	.430	.374	.403	.577	.332	1.000

сильно коррелированы (парные коэффициенты корреляции 0.854, 0.882 и 0.978), т.е. имеет место дублирование информации, и потому, по-видимому, есть возможность перехода от исходного числа признаков (переменных) к меньшему.

Сравнение различных регрессий. Пошаговый отбор переменных.

На 1-м шаге (k = 1) найдем один наиболее информативную переменную. При k = 1 величина R² совпадает с квадратом обычного (парного) коэффициента корреляции

R² = r² (y, x) ,

из матрицы корреляций находим:

r² (y, x_j) = r² (y, x₄) = (0.577)² = 0.333

Так что в классе однофакторных регрессионных моделей наиболее информативным предиктором (предсказателем) является x₄ - количество удобрений. Вычисление скорректированного (adjusted) коэффициента детерминации по (20) дает

R²_adj (1) = 0.296.

Это значение получаем возвратом в окно Select dep. And indep. Var. Lists: Dep. Var: y, Indep. Var.: x4 -OK - OK.

2-й шаг (k = 2). Среди всевозможных пар (х₄ , х_j ), j = 1, 2, 3, 5, выбирается наиболее информативная (в смысле R² или, что то же самое, в смысле R²_adj ) пара:

возврат в окно Select dep. And indep. Var. и перебор различных пар; результат:

(х₄ , х₁) = 0.406, (х₄ , х₂) = 0.399,

(х₄ , х₃ ) = 0.421, (х₄ , х₅) = 0.255,

откуда видно, что наиболее информативной парой является (х₄ , х₃ ), которая дает

(2) = (х₄ , х_j) = 0.421

Оценка уравнения регрессии урожайности по факторам х₃ и х₄ имеет вид (х₃ , х₄) = 7.29 + 0.28 х₃ + 3.47 х₄ (27)

(0.66) (0.13) (1.07)

Внизу в скобках указаны стандартные ошибки, взятые из столбца Std. Err. of B таблицы Regression Results для варианта независимых переменных (х₃ , х₄) Все три коэффициента статистически значимо отличаются от нуля при уровне значимости  = 0.05, что видно из столбца p-level той же таблицы.

3-й шаг (k = 3). Среди всевозможных троек (х₄ , х₃ ,х_j), j = 1, 2, 5, выбираем аналогично наиболее информативную: (х₄ , х₃ ,х₅), которая дает (3) = 0.404,

что меньше, чем (2) = 0.421; это означает, что третью переменную в модель включать нецелесообразно, т.к. она не повышает значение (более того, уменьшает). Итак, результатом анализа является (28).