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

Обычно в сельском хозяйстве применяются комбинированные удобрения (NPK) с внесением азотных (N), фосфорных (P) и калийных (K) удобрений в определенной пропорции. Пусть мы располагаем данными об урожайности некоторой культуры от возрастающих доз комбинированного удобрения NPK. Возможно ли по таким данным выделить отдельно эффекты каждой составляющей N, P, K? Ну, конечно же, нельзя. Компоненты N, P, K в этой задаче связаны точной линейной зависимостью и при подключении в модель одной компоненты (любой) процесс подключения автоматически завершится. При этом весь совместный вклад трех видов удобрений будет целиком приписан только одной компоненте, которую метод подключения случайным образом изберет первой. Относительно оставшихся (не подключенных в модель показателей делается стандартное (и совершенно неверное) заключение типа: “... в данных почвенно-климатических условиях (почему-то) не проявилось действие некоторых видов удобрений...”.

где A = X^TX .

Если матрица A - невырожденная, то для нее существует обратная матрица C=A^-1 и искомое решение, а также дисперсии и ковариации коэффициентов регрессии выражаются через элементы этой обратной матрицы:

В случае вырождения предлагается обратную матрицу A^-1 (и ее элементы) во всех приведенных выше формулах заменить на так называемую псевдообратную (“обобщенную обратную”) матрицу A^-, для которой всегда имеют место соотношения:

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

Прежде чем приступить к детальному изучению алгоритма с псевдообратной матрицей, обсудим свойства нового обобщенного решения, которые нам могут и не понравиться. При вырождении система нормальных уравнений имеет не одно определенное, а бесчисленное множество решений (так как система нормальных уравнений всегда совместна). С помощью псевдообратной матрицы выделяется решение с наименьшими значениями коэффициентов регрессии (в модель включаются все переменные и общий вклад “размазывается” по большому числу членов модели). Но малое значение коэффициента регрессии, как правило, означает малую его значимость по критерию Стьюдента. Не получится ли так, что новое обобщенное решение будет значимо в целом (по критерию Фишера), но состоять только из незначимых членов?

Произведение матрицы A на вектор X (вектор - это матрица-столбец) означает преобразование исходного вектора X в новый вектор: Y=AX. Поскольку любой вектор может быть разложен по системе координатных, достаточно рассмотреть преобразование этих координатных векторов: I₁ = (1, 0, 0, ..., 0)^T ; I₂ = (0, 1, 0, ..., 0)^T и т.д. до I_m = (0, 0, 0, ..., 1)^T . Нетрудно непосредственно убедиться, что преобразованные координатные векторы представляют собой столбцы матрицы A: AI_k = A_k (напоминаем, что матрицы перемножаются по правилу “строка на столбец").

Разным матрицам соответствуют различные линейные преобразования векторного пространства. Среди этих преобразований выделим "ортогональные" преобразования, которые не изменяют ни длин векторов, ни углов между ними (к ортогональным преобразованиям относятся всевозможные повороты и зеркальные отражения объектов). По виду матрицы A можно сразу установить, является ли преобразование с этой матрицей ортогональным: для этого необходимо, чтобы новая система координатных векторов (столбцы матрицы A) осталась ортонормированной (суммы квадратов элементов любого столбца равны единице, суммы произведений элементов любых двух столбцов равны нулю). Отсюда следует, что для ортогональных матриц имеет место соотношение A^TA = I, где I - единичная матрица. Иными словами, транспортирование ортогональной матрицы приводит к ее обращению A^T = A^-1. Ввиду коммутативности произведения взаимообратных матриц для ортогональных матриц получаем также A^TA = AA^T = I .

Любую симметричную квадратную матрицу A можно представить в виде произведения матричных сомножителей

где l _i - так называемые "собственные" или “характеристические" числа, а U - ортогональная матрица собственных векторов.

Если матрица A - невырожденная, то все ее собственные числа отличны от нуля. Нетрудно убедиться, что в этом случае матрица

является обратной к матрице A.

Если матрица A - вырожденная, то одно или несколько ее собственных чисел (последних) равны нулю. Рассмотрим систему нормальных уравнений с такой матрицей:

Умножим это матричное равенство слева на U^T, учтем, что U^TU = I, и временно переобозначим U^TX^TY = Z , U^T b = d . Получим такую систему уравнений:

Здесь (для примера) принято, что последнее собственное число равно нулю (l _m = 0); ввиду совместности системы нормальных уравнений автоматически получится z_m = 0.

d_m -произвольное.

Если положить произвольное значение d_m нулю, то вектор d будет иметь минимально возможную длину. Поскольку вектор коэффициентов регрессии b = U d получается ортогональным преобразованием из вектора d, его длина будет той же самой, т.е. минимально возможной.

Здесь по диагонали стоят числа, обратные ненулевым собственным числам матрицы A; остальные элементы матрицы - нули. Тогда

которая является псевдообратной (обобщенной обратной) к матрице A (непосредственно убеждаемся, что AA^-A = A и A^-AA^- = A^-). Если все собственные числа отличны от нуля, то A^- = A^-1.

Представим, что все собственные числа отличны от нуля, но одно или несколько чисел близки к нулю. Пусть l _m ^»  0. Ввиду совместимости системы нормальных уравнений автоматически получится, что z_m ^»  0 и формально можно вычислить компоненту . Здесь и числитель и знаменатель близки к нулю, т.е. это выражение напоминает раскрытие неопределенности типа “нуль на нуль" . Но из-за неизбежных вычислительных погрешностей результат такого отношения является неустойчивым и значение d_m будет получено с большой ошибкой, которая будет добавлена к каждому коэффициенту регрессии (к вектору b добавляется вектор d_mU_m). В этом случае (в случае плохой обусловленности матрицы A) разумнее малое собственное число l _m ^»  0 сразу положить нулю и перейти к обобщенному устойчивому решению с псевдообратной матрицей A^-. Вопрос, какое собственное число l _m ^»  0 считать малым, а какое - нет, буден рассмотрен в дополнении к этой лекции.

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

Заменим исходную систему переменных x₁ , x₂ , ... , x_m на некие их линейные комбинации (компоненты):

где .

Оказывается, набор коэффициентов u_ij для каждой компоненты образует собственный вектор корреляционной матрицы (матрицы системы нормальных уравнений), а соответствующее ему собственное число l _i равно дисперсии i-ой компоненты . Отбрасываем все компоненты с дисперсиями, меньшими единицы (l _i < 1), если сумма оставшихся собственных чисел составляет не менее 75% от числа исходных переменных. При переходе к такой системе главных компонент происходит существенное понижение размерности задачи.

Составляем уравнение регрессии относительно главных компонент:

где k £ m - число главных компонент. Ввиду некоррелированности компонент система нормальных уравнений распадается на отдельные уравнения:

Получаем вектор d - вектор коэффициентов регрессии в разложении результативного признака по главным компонентам. Осталось подставить вместо компонент их выражения через исходные переменные и перегруппировать члены, что эквивалентно матричному преобразованию b=U d .

Таким образом, временным переходом к системе главных компонент обходят все затруднения, связанные с вырождениями. Другие приложения главных компонент составили содержание нового статистического анализа - так называемого, "факторного анализа".

y = b₀ + b₁t + b₂t² + e ,

где t — последовательные годы (т.е. числа порядка 2000) за 20-летний период.

Здесь коэффициенты нормальной системы имеют порядок от а₁₁ = 1 до а_33 »  2000⁴ = 16× 10¹², причем недопустима никакая потеря младших разрядов. Дж. Форсайт показал, что расчет этого примера на ЭВМ даже с двойной точностью не приводит к правильной аппроксимации зависимости квадратичной моделью.

Конечно, затруднения именно этого примера легко преодолеваются - достаточно в качестве t задать последовательные номера (без сотен и тысяч лет); еще лучше предварительно центрировать переменную t ; наконец здесь достаточно просто можно произвести ортогонализацию членов квадратной модели. Однако предупреждение Дж. Форсайта остается в силе - существуют ситуации, когда сам переход к системе нормальных уравнений может привести к потере точности. Мы настолько привыкли составлять и решать систему нормальных уравнений для минимизации суммы квадратов невязок, что сама мысль о том, что этот путь не является единственно возможным (и даже не наилучшим), кажется неожиданной. Несколько десятилетий назад, когда все расчеты выполнялись вручную, другого пути реализации МНК практически не существовало, но сейчас, при доступности персональных компьютеров, при наличии таких пакетов прикладных программ, как Матлаб или Статграф, просто желательно рассматривать решение системы нормальных уравнений всего лишь как первое приближение, которое подлежит дальнейшему уточнению. И нередко специалист с изумлением убеждается, что стандартные МНК-оценки параметров были получены с большими смещениями (см. также пример в лекции 8).

В методе SVD Дж. Форсайта используется возможность представления любой прямоугольной матрицы размером n ´  m в виде произведения (сингулярного разложения):

где U - ортогональная матрица размера m´  m; V - ортогональная матрица размера n´  n; D₆ - ортогональная матрица размера n´  m, в которой отличны от нуля не более m элементов по главной диагонали.

Нетрудно убедиться, что матрица U состоит из собственных векторов матрицы X^TX, матрица V состоит из собственных векторов матрицы XX^T, квадраты чисел s _i² = l _i есть собственные числа матрицы X^TX (и ненулевые собственные числа матрицы XX^T).

умножим равенство слева на V^T и введем обозначения: W=V^TY; O=V^TE; U^T b=d . Получим:

Здесь k - количество ненулевых чисел s _i (k£ m). Ввиду ортогональности матрицы V суммы квадратов элементов векторов E и O одинаковы, причем можно первые k элементов вектора O сделать равными нулю, полагая . При таких значениях параметров достигается минимум суммы квадратов невязок .

При наличии мультиколинеарности несколько последних чисел s _i равны нулю (или близки к нулю), поэтому последние компоненты вектора d могут принимать любые значения. Положим их нулю, тем самым получим решение с наименьшими значениями коэффициентов регрессии. Ортогональное преобразование b=U d с матрицей U не изменяет величины суммы квадратов коэффициентов регрессии .

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

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

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

Здесь m - число подключенных переменных; c_ii - диагональные элементы обратной матрицы.

Вычисляем дисперсионное отношение для наибольшего коэффициента детерминации

Это значение оказалось меньше табличного F_0.05(2,17)=3.59, поэтому все внутренние корреляционные связи при n=20 признаются незначимыми. Однако, если бы n было бы не 20, а 200, то та же самая мультиколинеарная связь была бы уже значимой, так как в этом случае дисперсионное отношение

уже значительно превышает верхнюю границу F_0.01(2,197)=4.60 .

Ну, и что из этого следует? Хоть что-нибудь изменилось в наших представлениях? Мы и так не сомневались в том, что переменная x₄ коррелирует с переменой x₃ , для этого и решалась система нормальных уравнений.

Критерий Феррара-Глобера оценивает значимость внутренних связей между членами модели, а не степень вырожденности корреляционной матрицы, а ведь нас интересует только это. Не ту сторону явления освещает критерий Феррара-Глобера, не разрешает он поставленной проблемы. Настоящий критерий мультиколинеарности должен бы оценивать допустимость (или недопустимость) существующей тесноты внутренней связи, но не с точки зрения ее статистической значимости или незначимости; статистические представления здесь вообще ни при чем.

Теперь рассмотрим второй важный аспект поставленной проблемы. Как уже упоминалось выше, в процессе подключения дополнительных членов иногда наблюдается явление “ прохлопывания” , когда после подключения в модель очередного члена часть коэффициентов (или даже все коэффициенты) регрессии резко меняют свои значения (они могут изменить даже знаки). Эффект этот был замечен не только при обработке статистических данных, но также при аппроксимации функций при достижении некоторого критического числа членов модели. Неизвестно, связан ли эффект “ прохлопывания” с вырождением матрицы, но он является отличительным признаком наличия мультиколинеарности. Что-то подобное нам встретилось (см. предыдущую лабораторную работу) при составлении сопряженой модели, когда в качестве результативного признака была выбрана переменная x₄ . В табл. 9.1 продублирована корреляционная матрица, где переменная y переобозначена через x₀ , чтобы она не отвлекала внимание своим внешним видом.

Корреляционная матрица

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

• подключена первая переменная x₃ :

;

• дополнительно подключена переменная x₂ :

• дополнительно подключена переменная x₀ :

При подключении переменной x₀ предыдущие b -коэффициенты резко изменились; переменная x₀ оказалась значимой, ее вклад в несколько раз превышает вклад обоих ранее включенных членов. Интересно, что мультиколинеарные связи совершенно скрыли это влияние x₀ , оно проявилось только при фиксировании на постоянном уровне переменных x₂ , x₃ (сравните и ). Здесь примечателен еще один момент: b -коэффициенты по модулю оказались большими единицы. Некоторые специалисты до сих пор считают, что этого не может быть.

Рассмотрим разложение переменной x₃ по x₀ и x₂ (разложение многомерного вектора x₃ по базису x₀ , x₂); все переменные предполагаются стандартизированными:

и решим ее: b ₀ = -1.3376; b ₂ = 1.0834 .

Коэффициент детерминации (вклад двух переменных x₀ , x₂)

значительно превышает вклад только одного показателя x₀ ( =0.2116), т.е. двухфакторная модель - вполне доброкачественная. В процессе решения системы не встретилось никаких вырождений; определитель системы det = 1 - 0.81² = 0.3439 был отличен от нуля (и даже не близок к нулю). Единственная особенность этого (и предыдущего) примера в том, что здесь внутренняя связь между аргументами была теснее любой внешней связи с результативным признаком ( больше и ). Это и надо называть мультиколинеарностью?

Дадим геометрическую интерпретацию полученным результатам. После стандартизации все “многомерные вектора” имеют одинаковую “длину”. Могут ли проекции единичного вектора x₃ при разложении по единичному базису (x₀ , x₂) быть больше единицы (по модулю)? Конечно, нет, если базис ортогональный . Однако в косоугольной системе координат это возможно (b ₀ , b ₂ - косоугольные проекции). Анализируя структуру коэффициента корреляции, замечаем, что его можно интерпретировать, как косинус угла между многомерными векторами. Находим углы между нашими векторами:

На рис. 9.1 изображена плоскость, проходящая через единичные вектора x₀ , x₂ (угол между которыми j ₀₂ = 35,9° ). Вектор x₃ не лежит в этой плоскости, его проекция на координатную плоскость располагается где-то в заштрихованной области (в секторе между 117.4° и 90°  + 35.9° = 125.9° ). Его длина равна . На рис. 9.1 изображен также вектор с косоугольными проекциями (b ₀ , b ₂)=(-1.3376; 1.0834), найденными в аналитическом расчете. Этот вектор полностью соответствует проекции вектора x₃ на плоскость векторов x₀ , x₂ . Таким образом, косоугольные проекции (b ₀ , b ₂) вектора x₃ действительно оказались по модулю большими единицы (третьей проекции вектора на рис. 9.1 не видно - это случайная ошибка).

Итак, при разложении единичного вектора по единичному, но не ортогональному, базису косоугольные проекции могут быть большими единицы (по модулю). Иными словами, из-за внутренних (вполне допустимых) мультикорреляционных связей между членами модели в уравнении регрессии могут появиться b -коэффициенты большие единицы, причем это происходит на этапе подключения какой-то очередной переменной. С точки зрения математики здесь все в порядке. Модель с такими значениями b -коэффициентами правильно воспроизводит значения результативного признака. Интересно, что эффект “ прохлопывания” оказался не связанным с вырождениями. Пока мы не располагаем никакими критериями, которые бы предупреждали нас о возможности появления рассмотренного эффекта (а может быть, и не нужны такие критерии?).

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

Диагональные элементы s_ii - есть дисперсии невязок при разложении каждой не включенной еще в модель переменной по системе уже учтенных показателей. Представим себе, что одна из оставшихся вне модели переменных x_k практически точно описывается учтенными показателями, тогда все ее невязки являются всего лишь малыми ошибками и, следовательно, ее дисперсия и все ее частные ковариации с другими оставшимися вне модели переменными будут также близки к нулю: . В формуле для вычисления коэффициента частной корреляции и числитель, и знаменатель будут близки к нулю, но значение коэффициента частной корреляции может быть каким угодно, в том числе и близким к единице. В этом случае наш алгоритм подключит в модель переменную x_k с разрешающим элементом s_{kk »} 0, близким к нулю, в результате чего элементы обратной матрицы резко возрастут; очень часто при этом проявляется эффект “ прохлопывания” , при котором b -коэффициенты изменяются совершенно произвольным образом. Предлагаем дополнить алгоритм подключения следующим элементарным контролем: если какой либо диагональный элемент s_kk окажется меньше, скажем, 5% исходной величины (т.е. меньше 0.05), то такую переменную x_i следует исключить из дальнейшего анализа. Это правило надо распространить и на результативный признак: - если s₀₀ < 0. 05 , то процесс дальнейшего подключения надо прервать, так как новым переменным уже нечего объяснять.

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

Поскольку определитель квадратной матрицы A можно разложить по элементам любой строки D=|a_ij| = a_k1 A_k1 + a_k2 A_k2 + ... + a_kn A_kn , то частные производные определителя по любому элементу матрицы будут равны

где A_km - алгебраические дополнения, которые выражаются через элементы обратной матрицы c_mk .

Таким образом, относительная погрешность определителя при изменении элемента a_km на величину D a_km будет равна c_mk^.D a_km .

Продемонстрируем это на примере обращения плохо обусловленной матрицы Гильберта, элементы которой составляются по принципу . Известно, что точные значения элементов обратной матрицы Гильберта должны быть целочисленными. Так, обращая в рациональных числах матрицу

, получим .

.

Совершенно ничтожные ошибки округления в коэффициентах a₁₃ , a₂₂ , a₃₁ (абсолютная погрешность D  a₁₃ = D  a₂₂ = D  a₃₁ = 1/3000; относительная погрешность менее 0.1%) привели к заметным погрешностям в элементах обратной матрицы (сравни c₃₃ = 195.77 c точным значением 180).

Оценим последствия малых изменений элементов матрицы a₁₃ , a₂₂ , a₃₁ по предлагаемой выше методике. Вычислим ожидаемую относительную погрешность определителя: (c₃₁ + c₂₂ + c₁₃)× D  a₁₃ = -(33.28 + 210.19 + 33.28) / 3000 = -0.09 (т.е. 9%).

Результирующая погрешность возросла в 270 раз по сравнению с погрешностью в исходных коэффициентах (сравните 9% и 100/3000=1/30% )!
И этой цифре можно верить. В данном примере мы точно знали погрешность каждого элемента матрицы, поэтому мы оценили существующую погрешность определителя (а не ту, которая может появиться). Мы установили, что определитель из-за ошибок округления уменьшился на 9% , следовательно, решение системы и элементы обратной матрицы увеличились в 1.09 раз. Откорректируем полученную обратную матрицу, уменьшив все ее элементы в 1.09 раз:

.

Итак, мы имеем возможность оценивать последствия малых погрешностей в исходных коэффициентах (т.е. оценивать степень неустойчивости решения). Применительно к нашей задаче заметим, что все исходные коэффициенты корреляции (коэффициенты нормальной системы) вычислялись с четырьмя знаками после запятой; один знак считается резервным, поэтому разумно предположить, что погрешности недиагональных элементов исходной матрицы не превосходят по модулю единицы третьего знака, т.е. 0.001. Следовательно, относительную погрешность решения не будет превосходить (суммируются модули всех недиагональных элементов обратной матрицы). Оцениваем устойчивость решения при подключении в модель переменных x₂ , x₃ , x₄ . Относительная погрешность решения (см. табл. 6.7)

оказалась того же порядка, что и исходные погрешности в коэффициентах системы. В процессе подключения все разрешающие элементы были отличны от нуля (в табл. 6.7 разрешающие элементы взяты в рамочку), никакого вырождения здесь не встретилось, поэтому исходные погрешности не увеличились. А вот в одном из вариантов метода подключения-исключения нам пришлось обращать расширенную матрицу, причем на одном из этапов мы произвели преобразование Жордана-Гаусса с разрешающим элементом s₀₀=0.0648, близким к критической границе 0.05. Элементы обратной матрицы при этом заметно возросли (см. табл. 7.2). Оценим возможные последствия. Относительная погрешность всех дальнейших результатов здесь может достигать величины

Исходные погрешности увеличились в 80 раз! Погрешность в 8% уже заметная (хотя еще терпимая). При меньших же величинах разрешающего элемента (меньших 0.05) относительная погрешность будет уже порядка (10¸ 15)% , что уже считается недопустимым.