Классический линейный регрессионный анализ

2. КЛАССИЧЕСКИЙ ЛИНЕЙНЫЙ РЕГРЕССИОННЫЙ АНАЛИЗ 2.1. Оценки модели линейной регрессии 2.1.1. Построение модели Слайд 2 Пусть n-вектор Y, связан с q-мерной неслучайной векторной переменной x = [x1, …, xq]′. Значения Yi, i = 1, …, n, полученные в эксперименте при заданных xi = [xi1, …, xiq]′, случайным образом изменяются около некоторого неизвестного истинного уровня (xi). Тогда можем записать Yi = (xi) + i, i = 1, …, n, (2.1) где i – случайная ошибка, которая объясняет отклонение Yi от величины (xi). При этом  может быть случайной компонентой, присущей величине (x), и представлять случайную ошибку измерения значений Y или влияние различных неучтенных факторов. Предположим, что (x) можно описать линейной моделью первого порядка по xj с q переменными η(x, β)  β 0  β1x1    β q x q . , где  = [0, 1, …, q]′  вектор неизвестных параметров (коэффициентов), подлежащий оцениванию. Тогда получим Y  β 0  β1x1    β q x q  ε . (2.2) Если в формуле (2.2) M = 0, то условное математическое ожидание случайного вектора Y при заданных переменных xj, j = 1, …, q равно M[Y / x]  β 0  β1x1    β q x q . (2.3) Уравнение (2.3), в котором x играет роль «независимой» переменной, называется уравнением регрессии или просто регрессией. Термин «регрессия» впервые был введен Ф. Гальтоном1 (1886) в теории наследственности для обозначения явления «возврата к среднему состоянию» (regression to mediocrity), состоящего в том, что дети тех родителей, рост которых превышает среднее значение на a единиц, имеют в среднем рост, превышающий среднее значение меньше чем на a единиц. В дальнейшем переменные1 Y, xj, j = 1, …, q и  будем называть откликом, регрессорами и остатком (используются другие названия этих переменных: выход, зависимая или эндогенная переменная; факторы, предикторы, входные, экзогенные или независимые переменные; ошибка, помеха, невязка). Гальтон (Galton) Фрэнсис (16.02.1822 – 17.01.1911) – английский психолог и антрополог. Следуя большинству книг по регрессионному анализу, случайный вектор-отклик и матрицы будем обозначать полужирными прописными буквами; векторы-регрессоры и вектор-остаток строчными полужирными буквами. 1 1 Когда используется уравнение (2.2) при анализе совокупности данных и оценивается вектор параметров , то предполагается, что элементы этой совокупности однородны в смысле подчинения одному и тому же причинному закону. Это означает, что параметры j приемлемы для каждого отдельно взятого наблюдения. Слайд 3 Пример 2.1. Идентификация статических характеристик сложного объекта, выходы которого, измеряемые со случайными ошибками, является функциями многих входных переменных. ε1  ε q Необходимо по наблюдениям входов и выходов z1 Y1 определить эти функции. В   F(z, u, ) общем случае совокупzp Ym ность переменных, опредеu1  u r ляющих текущее состояние сложного объекта, можно описать следующими груп- Рис.2.1. Модель сложного пами входных и выходных объекта переменных (рис. 2.1). 1. Контролируемые неуправляемые переменные z = (z1, …, zp), значения этих переменных можно измерить, но нельзя произвольно изменить. 2. Контролируемые управляющие переменные u = (u1, …, ur), значения, которых в любой момент времени можно изменить в пределах допустимого диапазона. 3. Неконтролируемые неуправляемые переменные  = (1, …, q), которые характеризуют множество реально существующих факторов, влияющих на текущее состояние объекта, но недоступных контролю и управлению. 4. Контролируемые управляемые переменные Y = (Y1, …, Ym), которые характеризуют результат функционирования объекта. Входные переменные z, u,  могут рассматриваться как причины, оказывающие влияние на каждую из выходных переменных Yi. При общем рассмотрении нет необходимости разделять контролируемые переменные (z,u,) поэтому объединим их в одну группу и обозначим X. Далее будем полагать, что xj при j = 1, …, q – неслучайные контролируемые независимые переменные;  – случайная неконтролируемая переменная (остаток, помеха, ошибка). Так как каждая из выходных переменных Yi полностью определяется в вероятностном смысле группой входных переменных X и остатком , то достаточно рассмотреть схему с одной выходной переменной (откликом). Будем полагать, что случайная остаток  аддитивно приложен к выходной переменной Y, т. е. Y=  + . Тогда физическую модель, характеризующую зависимость Y от Х можно выразить уравнением (2.1). Слайд 4 Структурная схема объекта, соответствующая этой модели, приведена на рис. 2.2. В общем случае функε ция (x,) нелинейна отноx1 сительно вектора Y  ( x, β ) параметров . Простейшим  и важнейшим для практики xq является случай линейной зависимости (x,β) от . Рис. 2.2. Структурная схема Линейную регрессионную объекта модель можно получить, разложив (x,) в ряд Тейлора в точке x0 = 0. Тогда q  ∂η Yi  η(0, , 0)     j 1 ∂x j   2  xij  1   ∂ η 0 2 j k  ∂x j ∂xk    xij xik    ε . i 0  Обозначив постоянные β 0  (0, , 0) , β j  (∂η / ∂x j ) 0 , β jk  (1/ 2)(∂ 2 η / ∂x j ∂xk ) 0 , получим Yi  β 0   qj 1β j xij   j  k β jk xij xik    εi . Ограничимся рассмотрением в этом уравнении только первых двух членов, случайные ошибки и ошибки за счет неучтенных членов ряда отнесем к остатку . При этом будем полагать, что неучтенные члены не коррелированы с учтенными. Тогда уравнение можно переписать в виде модели линейной регрессии (2.2).  Модель вида (2.2) является весьма общей и очень широко используется. Частными случаями ее являются, например, полиномиальная модель q-го порядка одной переменной Yi  β 0  β1 xi  β 2 xi2    β q xiq  ε i , Основное свойство модели вида (2.2) заключается в ее линейности по отношению к неизвестному вектору коэффициентов . По сравнению с ней, например, модель Yi  β 0  β1e β 2 xi  ε i нелинейная по параметру 2. Слайд 5 Рассмотрим оценки β̂ вектора коэффициентов  регрессионной модели (2.3). При этом будем различать два типа оценок. Первый – точечные оценки, получаемые на основании наблюдаемых данных регрессоров и отклика. Второй тип оценок связан с построением доверительных областей (интервалов) в пространстве оценок, которые с заданной вероятностью «накрывают» неизвестное истинное значение. Анализ уравнения (2.3) и оценку его коэффициентов будем проводить с использованием матричной алгебры. Применение матриц упрощает расчеты и придает им наглядность. 2.1.2. Оценивание параметров. Свойства оценок Рассмотрим схему, изображенную на рис. 2.2, где пунктиром выделена ненаблюдаемая часть. Пусть отклик Y связан с входами полиномом вида (2.2). Записывая эти n уравнений в матричной форме, получаем  Y1   x10 x11 Y   x x  2    20 21        Yn   xn 0 xn1     x1q  β 0   ε1  x2 q   β1  ε 2               xnq  β q  ε n  или Y  Xβ  ε, (2.4) где x10 = x20 = …= xn0 = 1. Матрица X типа «объект-признак» (см. п. 1.2.1) размера n×(q + 1) называется регрессионной матрицей, столбцами которой являются регрессоры xj, j = 1, …, q, а строками – n объектов или опытов; Y и  – n-векторы отклика и остатков,   подлежащий оцениванию (q + 1)-вектор неизвестных коэффициентов. В активных экспериментах элементы матрицы X выбираются равными только нулю и единице и в этом случае X называется матрицей плана. Необходимо по наблюдениям (xi1, …, xiq, Yi), i = 1, …, n найти наилучшую оценку ˆβ  (βˆ , βˆ , , βˆ ) вектора коэффициентов  = 0, 1, …, q)＇уравнения регрессии 1 q (x,). Одним из самых распространенных методов оценки вектора коэффициентов регрессии  является метод наименьших квадратов (МНК) (least squares method). Для обеспечения эффективности МНК-оценок должны соблюдаться следующие постулаты1: 1 В реальной жизни все эти постулаты редко соблюдаются. Влияние нарушений постулатов рассматривается в главе 7. Слайд 6 1. Число регрессоров q полинома (2.3) априори известно точно. 2. Все регрессоры измеряются без ошибок, а вычисления проводятся абсолютно точно. 3. Остаток  является независимой нормально распределенной случайной величиной с нулевым средним M[] = 0 и неизвестной постоянной дисперсией 2 при всех i = 1, …, n. 4. Дисперсия отклика Yi постоянна, или является известной функцией номера наблюдения i = 1, …, n. 5. Распределение Yi одинаково при всех i = 1, …, n. 6. Число опытов n существенно больше числа регрессоров q. Слайд 7 Если постулаты (1 – 6) соблюдаются и X′X обратимая матрица, то согласно фундаментальной теореме Гаусса-Маркова наилучшей оценкой вектора коэффициентов  является оценка β̂ , доставляющая минимум суммы квадратов остатков (невязок, ошибок, помех): Q(β̂)   in1(Yi   qk 0 xik β̂ k )2  min . Заметим, что этот остаток не может равняться нулю, так как число наблюдений n превосходит число неизвестных параметров q. Если Q(βˆ ) имеет производные по β̂ m , то необходимым условием минимума являются уравнения ∂Q / ∂β̂ m  2 in1(Yi   qk 0 xik β̂ k ) xim  0 , или n q n  i 1 xim  k 0 xik β̂ k   i 1 ximYi , m  0, 1, , q . (2.5) Система уравнений (2.5) называется системой нормальных уравнений (СНУ) МНК. Слово «нормальных» не связано с нормальным распределением вероятностей, а только подчеркивает, что уравнения, как правило, имеют такой «нормальный» вид. Слайд 8 Обозначим  = X, θˆ  Xβˆ . Будем минимизировать величину ′ = ||Y  ||2 по отношению к , где  – подпространство оценок θ̂ . Если изменять значения вектора  в пределах , то квадрат длины вектора ||Y  ||2 достигнет минимума при значении θ  θˆ , которое является проекцией вектора Y на подпространство . Тогда справедливо (Y  θˆ )  θˆ и, следовательно, (Y  θˆ )  X (рис. 2.3). Отсюда для скалярного произведения (Y  θˆ ) и X получаем (Y  θˆ , X)  0 или Xθˆ  XY . (2.6) Если столбцы матрицы X линейно независимы, Y  θˆ то существует единА ственный вектор параθ̂  метров , для которого Рис. 2.3. Нахождение точки θˆ  Xβˆ . Тогда система А, (2.6) выразится в виде для которой норма Y  θˆ миСНУ XX βˆ  XY . нимальна Y Слайд 9 Обозначим C = X′X, ,  = X′Y, , тогда СНУ запишется в виде (2.7) Cβ  Ψ . Свойства матрицы С: a) так как регрессоры x1, …, xq линейно независимы, то матрица С невырождена; b) С положительно определена и ранг ее в точности равен q; c) С–симметричная матрица, С = С′, т. е. является эрмитовой1. Отсюда следует, что СНУ (2.6) имеет единственное решение (2.8) βˆ  C1 Ψ  ( XX)1 XY . Оценка β̂ называется оценкой метода наименьших квадратов (МНК-оценкой). Так как в решении (2.8) матрица A = (X′X)1X′ неслучайная, то βˆ  AY является линейной комбинацией наблюдений Y. В соответствии с теоремой Гаусса-Маркова МНКоценка β̂ имеет наименьшую дисперсию среди всех возможных несмещенных линейных оценок. Слайд 10 Свойства оценок МНК-оценка β̂ является случайной величиной. Найдем математическое ожидание оценки β̂ . Используя решение (2.8) и, учитывая, что матрица X является детерминированной, получаем M[βˆ ]  M[(XX) 1 XY]  M[(XX) 1 X( Xβ  ε)]   ( XX) 1 XXβ  ( XX) 1 M[ε]  β. 1 Эрмит (Hermite) Шарль (24.12.1822 – 14.01.1901) – французский математик. Таким образом, математическое ожидание оценки вектора β̂ равно истинному значению β , т. е. β̂ является несмещенной оценкой. Другими словами, если эксперимент снова и снова повторяется при неизменной матрице Х, среднее значение β̂ будет равно β. Распределения До сих пор единственное предположение относительно  состояло в том, что M[] = 0, D[] = 2In. Если дополнительно предположить, что остаток  нормально распределен с параметрами 0, 2In, или, при кратком обозначении,   Nn(0, 2In), то Y  Nn (X, 2In). Отсюда получается целый ряд результатов, связанных с распределениями. Если Y  Nn (X, 2In), то: 1) βˆ ~ N n (β, σ ε2 ( XX) 1 ) , 2) (βˆ  β)XX(βˆ  β) / σ ε2 ~ χ (2q 1) , 3) β̂ не зависит от S e2 , 4) SSE / σε2  (n  q 1)Se2 / σε2 ~ χ 2nq1 . Предположение нормальности распределения остатков позволяет создать целостную систему статистической обработки, которая включает точечные, интервальные оценки и проверки статистических гипотез. Однако на практике распространенный миф нормальности распределения не всегда выполняется, а в случае малых выборок гипотезу нормальности распределения ошибок трудно проверить. Отклонение от нормальности может быть вызвано и засорением наблюдений чужеродными элементами. В этом случае для обнаружения и удаления этих элементов нужно применить методы, изложенные в предыдущей лекции. Другой подход связан с применением вместо МНК метода наименьших модулей (МНМ) [13]. Близким к МНМ является непараметрический регрессионный анализ, например, знаковый регрессионный анализ, который позволяет получать хорошие оценки и при сильно засоренных выборках [2, 8]. И, наконец, для таких данных можно использовать робастную регрессию [15, 21] или решать задачу регрессии с помощью нейронных сетей [10]. Слайд 10 2.1.3. Проверка статистических гипотез Поверка гипотез относительно коэффициентов регрессии В регрессионном анализе проверяются две нулевые гипотезы относительно коэффициентов уравнения. 1. Гипотеза равенства истинного коэффициента k значению k 0 H0: k = k 0, против альтернативной гипотезы H1: k  k 0. 2. Гипотеза одновременного равенства нулю всех, кроме 0, регрессионных коэф- фициентов H 0 : β1    β q  0 . Проверяем конкуренцию двух моделей Y   k βk xk  ε и Y  β 0  ε или, что то же самое, проверяем эффект от введения переменных x1, …, xq в модель регрессии. Первая гипотеза позволяет проверить значимость отдельно взятого коэффициента. Если β̂ k  Nn (k, 2 C2), то статистика β̂ k  β k 0 β̂ k  β k 0  t 1 Sβ̂ S C e kk k распределена по Стьюденту с числом с.с. v = n – q – 1, а квадрат ее имеет Fраспределение (распределение Снедокора-Фишера) с числом с.с. v1 = 1, v2 = n – q –1: Слайд 11 F (β̂ k  β 0k ) 2 . S e2 (C 1 ) kk (2.14) Если окажется, что вычисленное значение F связано с табличным для заданного уровня значимости  неравенством F  F (α, ν1 , ν 2 ) , ν1  1 , ν 2  n  q  1 , то гипотеза k = k 0 отвергается. Обычно значения k 0 неизвестны, поэтому проверяют гипотезу H0: k = 0. В этом случае получаем F β̂ k 2 (2.15) S e 2 (C 1 ) kk с числом с.с. v1 = 1, v2 = n – q – 1. Проверка этой гипотезы имеет важное значение, так как позволяет ответить на вопрос: можно ли считать k = 0? Если ответ положительный, то k-й регрессор можно удалить из рассматриваемой модели. Слайд 12 Для проверки второй гипотезы сначала выражение (2.14) перепишем в виде взвешенного расстояния между β̂ k  0 и k 0: 1 d (β̂ k ,β k 0 )  (β̂ k  β k 0 )(Sβ2 ) (β̂ k  β k 0 )  F ( ν1 , ν 2 ) . K Рассмотрим случай с центрированными данными ~ ~x  x  x , Y  Y Y . j j j Для центрированных данных оценка вектора коэффициентов βˆ   (β̂1, β̂2 , , β̂q ) не содержит свободного члена β̂ 0 , Оценку β̂ 0 находим по формуле (см. п. 2.4.4) β̂ 0  Y   qj1β̂ j x j . Тогда d (βˆ  ,0)  (βˆ  )(côv(βˆ  ))1 βˆ   Fур  q , ~~ где côv(βˆ  )  C1Se2  (XX)1 Se2 , а d (βˆ ,0) – расстояние Махаланобиса. Из этого уравнения находим статистику Fур ~ ~  ~ˆ ~ˆ (βˆ  ) X Xβˆ Y Y SSR / q , (2.16) Fyp    2 2 Se q S e q SSE /(n  q  1) где SSR   in1 (Yˆi  Yi ) 2 – сумма квадратов, объясненная уравнением регрессии (Sum of Squares about Regression), SSE   in1 ei2 – сумма квадратов остатков. Слайд 13 Формула (2.16) определяет отношение дисперсии оценки модели к дисперсии остатка Fур = S2/Se2. Статистика Fур имеет F-распределение с числом с.с. v1 = q, v2 = n – q – 1. Если Fур ≥ Fкр(, v1, v2), то нулевая гипотеза отвергается. Всякая сумма квадратов связана с числом с.с. Например, для SSE число с.с. равно числу опытов n минус (q + 1) коэффициентов регрессии. Используя формулы (2.12) и (2.16) мы можем построить таблицу дисперсионного анализа (ANOVA) (табл. 2.1). «Средний квадрат» получается при делении каждой суммы квадратов на соответствующее ей число с.с. Таблица дисперсионного анализа Т а б л и ц а 2 . 1. Источник Сумма диспер- квадратов сии Модели SSR Остатки SSE Полная SST ч.с.с. Средний квадрат νq MSE  ν  n  q  1 MSE  S e2  ν  n 1  SSR q SSE n  q 1 F-отношение F SSR / q SSE /( n  q  1)   Слайд 14 2.1.4. Адекватность модели Оценка постулируемой модели M[Y / x]  β 0  β1x1    β k x q , найденная по экспериментальным данным, равна ˆ  β̂0  β̂1x1    β̂ k x q . Y Адекватность модели означает, что постулируемая модель не противоречит наблюдениям. Она оценивается с помощью множественного коэффициента детерминации (coefficient of multiple determination), который равен квадрату коэффициента множественной корреляции R между Y и Ŷ R ˆ  Y) (Y  Y)(Y . ˆ  Y)(Y ˆ  Y)]1 / 2 [(Y  Y)(Y  Y)(Y ~ ~ Введем обозначения для центрированных данных: Y  Y  Y и Yˆ  Yˆ  Y . Тогда коэффициент детерминации запишется в виде R2 ~ ~ˆ 2 (YY )  ~~ ~~ . ˆ Y ˆ) (YY)(Y Сделаем следующие преобразования: ~ ~ˆ ~ ~ˆ ~ˆ ~ˆ ~ ~ˆ ~ˆ ~ˆ ~ˆ ~ˆ ~ˆ YY  (Y  Y  Y)Y  (Y  Y )Y  YY  YY  SSR , ~ ~ ~ˆ ~ˆ ~ ~ где (Y  Y)Y  0 в силу ортогональности (Y  Yˆ ) и Yˆ ; ~ ~ ~ ~ˆ ~ˆ ~ ~ˆ ~ˆ ~ ~ˆ ~ ~ˆ ~ˆ ~ˆ YY  (Y  Y  Y)(Y  Y  Y)  (Y  Y )(Y  Y)  YY  ~~ ˆ Y  SSE  SSR  SST .  ee  Y Здесь SST  in1 (Yi  Y )2 – полная сумма квадратов (Total Sum of Squares). В результате получаем R2 ~~ ˆ Y ˆ Y SSR SSR . ~~  YY SSR  SSE SST (2.18) Отсюда видно, что адекватность модели определяется отношением доли дисперсии, объясненной уравнением регрессии вариации откликов SSR к общей вариации SST. Если окажется, что величина R2 не меньше табличного R2(, v1, v2), где v1 = q, v2 = n – q – 1, то постулируемая модель адекватна. Так как R есть мера взаимосвязи, то значение коэффициента детерминации всегда находится между нулем и единицей 0 ≤ R2 ≤ 1. Равенство R2 единице свидетельствует о том, что модель полная и полнее не может быть (e′e = 0), тогда, как R2 = 0 означает, что модель абсолютно не способна объяснить вариацию наблюдаемых данных. Коэффициент детерминации, определяемый выражением (2.30), обладает одним существенным недостатком. При равенстве числа регрессоров q числу наблюдений n величина R2 равна 1. Кроме того, по мере добавления регрессоров в уравнение, значение R2 неизбежно возрастает. Это ведет к неоправданному предпочтению моделей с большим числом регрессоров. Отсюда следует, что необходима поправка к R2, которая бы учитывала число регрессоров и число наблюдений. В результате получаем скорректированный коэффициент детерминации (adjusted coefficient of multiple determination) R 2 (см. [6], кн.1, с. 128) R 2 1 n 1 (1  R 2 ) . n  q 1 Если в уравнение регрессии добавить некоторую переменную, то R 2 увеличится тогда и только тогда, когда F-статистика для соответствующего коэффициента переменной будет больше единицы. Коэффициент детерминации R2связан с Fур соотношением Fyp  R2 q SSR q ( SSR / SST ) / q .   SSE (n  q 1) (1  SSR / SST ) /(n  q 1) (1  R 2 ) /(n  q 1)