Корреляция и регрессия

Из R.Cody, J.Smith. Applied statistics and the SAS programming Language, Pretence Hall, 1991


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

Корреляция.

Мы записали значения "ПОЛ" (GENDER), "РОСТ" (HEIGHT) "ВЕС"(WEIGHT) и "ВОЗРАСТ"( AGE) семи человек и хотим подсчитать корреляцию между "РОСТОМ" (HEIGHT) и "ВЕСОМ"(WEIGHT) и "ВЕСОМ"(WEIGHT) и "ВОЗРАСТОМ"( AGE). Программа, которая может это выполнить будет выглядеть следующим образом:

DATA CORR_EG;
   INPUT GENDER $ HEIGHT WEIGHT AGE;
DATALINES;
M 68 155 23
F 61  99 20
F 63 115 21
M 70 205 45
M 69 170  .
F 65 125 30
M 72 220 48
;
PROC CORR DATA=CORR_EG;
   TITLE 'Example of a Correlation Matrix';
   VAR HEIGHT WEIGHT AGE;
RUN;

Результаты работы этой программы показаны ниже:

Correlation Analysis

            3 'VAR' Variables:  HEIGHT   WEIGHT   AGE


                        Simple Statistics

 Variable             N          Mean       Std Dev           Sum

 HEIGHT               7     66.857143      3.976119    468.000000
 WEIGHT               7    155.571429     45.796132   1089.000000
 AGE                  6     31.166667     12.416387    187.000000

            Simple Statistics

 Variable       Minimum          Maximum

 HEIGHT       61.000000        72.000000
 WEIGHT       99.000000       220.000000
 AGE          20.000000        48.000000


  Pearson Correlation Coefficients / Prob > |R| under Ho: Rho=0
  / Number of Observations

                     HEIGHT            WEIGHT               AGE

   HEIGHT           1.00000           0.97165           0.86614
                     0.0               0.0003            0.0257
                          7                 7                 6

   WEIGHT           0.97165           1.00000           0.92496
                     0.0003            0.0               0.0082
                          7                 7                 6

   AGE              0.86614           0.92496           1.00000
                     0.0257            0.0082            0.0
                          6                 6                 6

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

В данной распечатке мы видим, что корреляция между "РОСТОМ" (HEIGHT) и "ВЕСОМ"(WEIGHT) составляет 0,97165, уровень достоверности составляет 0,0003; корреляция между "РОСТОМ" (HEIGHT) и "ВОЗРАСТОМ"( AGE) составляет 0,86614 (p = 0,0257); корреляция между "ВЕСОМ"(WEIGHT) и "ВОЗРАСТОМ"( AGE) составляет 0,92496 (p = 0,0082). Давайте на мгновение сконцентрируемся на корреляции между "РОСТОМ" (HEIGHT) и "ВЕСОМ"(WEIGHT). Небольшое p-значение, подсчитанное для этого коэффициента корреляции, указывает, что крайне мало вероятно, что мы получили бы подобную большую величину корреляционного коэффициента за счет действия случайных, если бы выборка из семи пациентов была бы взята из популяции, где истинный коэффициент корреляции равен нулю. Помните, что на самом деле это всего лишь пример. Корреляции подобных размеров крайне редко встречаются в реальных данных.

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

PROC CORR options;
   VAR  список переменных;
RUN;
 

Термин "список переменных" должен быть замещен списком, состоящим из названий переменных, которые отделены друг от друга пробелами. Если не указаны никакие опции будет рассчитан Пирсоновский корреляционный коэффициент, а так же будет подсчитана простейшая описательная статистика. Как мы обсудим позднее имеется также несколько непараметрических корреляционных коэффициентов. Они вызываются при помощи ключевых слов SPEARMAN, которое приведет к расчету рангового коэффициента корреляции Спирмана, KENDALL, что приведет к расчету коэффициента Тау-б Кендала и HOEFFDING, которое приведет к расчету D -статистики Хоффдинга. Если Вы используете одну из этих опций, например Спирман (SPEARMAN), Вы должны использовать также опцию Пирсон (PEARSON), если Вы хотите, чтобы были рассчитаны и Пирсоновские коэффициенты корреляции. Если Вы не хотите, чтобы обсчитывались и распечатывались простые статистические показатели, включите в список Ваших опций ключевое слово NOSIMPLE.

Например, если бы Вы хотели, чтобы были рассчитаны как коэффициент корреляции Спирмана, так и коэффициент корреляции Пирсона, но не хотите, чтобы были рассчитаны простые статистические показатели, Вам необходимо было бы написать следующий код:

PROC CORR DATA=CORR_EG PEARSON SPEARMAN NOSIMPLE;
   TITLE 'Example of a Correlation Matrix';
   VAR HEIGHT WEIGHT AGE;
RUN;

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

Если все, что Вы хотите - это корреляции между одной или несколькими переменными и другим набором переменных, то для этого необходимо будет использовать команду WITH. Процедура CORR тогда рассчитает корреляционный коэффициент между каждой переменной, которая перечислена в списке, идущим следом за командой WITH и с каждой переменной, которая указана в команде VAR. Например, у нас есть переменные IQ и GPA (средняя оценка за экзамен) в наборе данных, который имеет название RESULTS. Мы, кроме того, подсчитали оценки студентов по десяти тестам (TEST 1- TEST 10). Если все, что мы хотим - это посмотреть на корреляции между IQ и GPA и этими десятью тестами, то тогда необходимо будет использовать следующий синтаксис:

 PROC CORR DATA=RESULTS;
   VAR IQ GPA;
   WITH TEST1-TEST10;
RUN;
 

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

Значимость коэффициента корреляции .

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

Значимость коэффициента корреляции является производной от величины коэффициента корреляции и размера выборки. Когда имеется большое количество точек, даже небольшой коэффициент корреляции может быть значимым. Например, когда имеется десять точек корреляционный коэффициент 0,63 или больше является достоверным (на уровне 0,05); со ста точками корреляционный коэффициент 0,195 уже будет достоверным. Обратите внимание, что отрицательная корреляция имеет абсолютно аналогичную силу связи, такую же как и положительная корреляция (несмотря на то, что взаимоотношения между двумя показателями противоположны). Корреляция - 0,40 имеет такую же силу , как и корреляция +0,40.

Чрезвычайно важно помнить, что корреляция указывает только на силу взаимоотношений - она не предполагает, что существуют причинно-следственные связи. Например, мы абсолютно с одинаковой вероятностью найдем высокую положительную корреляцию между количеством больниц в каждом из 50 штатов и количеством домашних животных в каждом штате. Означает ли это, что домашние животные приводят к тому, что люди заболевают и поэтому требуется строить новые больницы? Сомнительно. Наиболее вероятное объяснение заключается в том, что обе переменные количество и домашних животных и количество больниц связаны с размерами населения.

Быть ЗНАЧИМЫМ не равняется быть ВАЖНЫМ или СИЛЬНЫМ. Иными словами определение коэффициента корреляции не дает нам очень большой информации. После того, как мы узнаем, что наш коэффициент корреляции значимо отличается от нуля, мы должны дальше искать возможность интерпретировать важность этой корреляции. Давайте на секундочку мы отступим от основной темы обсуждения и зададим себе вопрос: "Что мы имеем в виду, когда говорим о значимости коэффициента корреляции?" Предположим у нас имеется популяция, в которой были измерены значения x и y и корреляция между этими точками является нулевой. Мы можем представить себе, что у нас будет график этой популяции, показанный ниже:

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

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

Как интерпретировать коэффициент корреляции.

Один из самых хороших способов интерпретации коэффициента корреляции (r) - это посмотреть на квадрат коэффициента (r - квадрат); r - квадрат может быть интерпретирован как пропорция вариабельности одной переменной, объясняемая вариабельностью другой переменной. В качестве примера корреляция между ростом и весом, которую мы рассчитывали раньше, составляет 0,97. Соответственно r - квадрат равняется 0,94. Мы можем сказать, что 94% вариабельности веса может быть объяснена вариабельностью роста (и наоборот). Кроме того (1 - 0,94) или 6% дисперсии веса является следствием действия других фактором, кроме вариабельности роста. Помня про эту интерпретацию, если мы посмотрим на коэффициент корреляции 0,4, мы обнаружим, что только 0,16 или 16% дисперсии данной переменной объясняется вариабельностью других переменных.

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

Важное допущение, про которое мы должны помнить, когда мы обсуждаем коэффициент корреляции, заключается в то, что пара точек х и у независимы от всех остальных точек. Иными словами каждая пара точек была измерена у разных людей. Например у нас имеется 20 человек и мы измерили две переменные, скажем частоту сердечных сокращений и артериальное давление у каждого из этих 20 человек. Затем мы можем подсчитать коэффициент корреляции, используя 20 пар точек. А теперь предположим, что мы сделали 10 измерений артериального давления и частоты сердечных сокращений для каждого пациента. Мы не можем подсчитать достоверный коэффициент корреляции, используя все 200 единиц данных (10 измерений для каждого из 20 человек), поскольку 10 измерений для каждого пациента не являются независимыми. Наилучшим методом в этом случае будет подсчитать среднее артериальное давление и среднюю частоту сердечных сокращений для каждого пациента, а затем использовать средние значения для того, чтобы подсчитать коэффициент корреляции. После того, как мы это сделаем, при интерпретации результатов мы должны будем помнить, что мы связывали друг с другом средние значения для каждого пациента.

Частные коэффициенты корреляции.

Исследователь может захотеть определить силу взаимосвязи между двумя переменными, когда удаляется влияние других переменных. Один из способов добиться этого заключается в подсчете частного коэффициента корреляции. Для того, чтобы удалить эффект одной или более переменной из корреляции, нам необходимо использовать команду PARTIAL для того, чтобы перечислить те переменные, чей эффект мы хотели бы удалить. Используя, описанный выше, набор данных CORR_EG мы можем подсчитать частные коэффициенты корреляции между ростом (HEIGHT) и весом (WEIGHT), удалив при этом эффект возраста ( AGE) следующим образом:

PROC CORR DATA=CORR_EG NOSIMPLE;
   TITLE 'Example of a Partial Correlation';
   VAR HEIGHT WEIGHT;
   PARTIAL AGE;
RUN;

Как Вы можете видеть в распечатке, представленной ниже частный коэффициент корреляции между ростом и весом теперь немножко ниже (0,91934), но все еще остается значимым (p=0,0272).

                       Correlation Analysis

              1 'PARTIAL' Variables:  AGE
              2 'VAR'     Variables:  HEIGHT   WEIGHT


          Pearson Partial Correlation Coefficients
          / Prob > |R| under Ho: Partial Rho=0 / N = 6

                              HEIGHT            WEIGHT

            HEIGHT           1.00000           0.91934
                              0.0               0.0272

            WEIGHT           0.91934           1.00000
                              0.0272            0.0

Линейная регрессия

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

WEIGHT |
       |
   225 +
       |                                                        A
       |
       |
       |                                              A
   200 +
       |
       |
       |
       |
   175 +
       |                                         A
       |
       |
       |                                    A
   150 +
       |
       |
       |
       |
   125 +                     A
       |
       |           A
       |
       |
   100 + A
       |
       +-+----+----+----+----+----+----+----+----+----+----+----+-
        61   62   63   64   65   66   67   68   69   70   71   72

                                  HEIGHT

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

Процедура REG (сокращение от слова регрессия) имеет следующую общую форму:

PROC REG опции;
   MODEL зависимая переменная (ые) = независимые переменные / опции;
RUN;

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

PROC REG DATA=CORR_EG;
   TITLE 'Regression Line for Height-Weight Data';
   MODEL WEIGHT=HEIGHT;
RUN;

Результат этой процедуры показан ниже:

Model: MODEL1
Dependent Variable: WEIGHT


                       Analysis of Variance

                         Sum of         Mean
Source          DF      Squares       Square      F Value
Prob>F

Model            1  11880.32724  11880.32724       84.451
0.0003
Error            5    703.38705    140.67741
C Total          6  12583.71429

    Root MSE      11.86075     R-square       0.9441
    Dep Mean     155.57143     Adj R-sq       0.9329
    C.V.           7.62399

                       Parameter Estimates

                       Parameter      Standard    T for H0:
      Variable  DF      Estimate         Error   Parameter=0

      INTERCEP   1   -592.644578   81.54217462        -7.268
      HEIGHT     1     11.191265    1.21780334         9.190


      Variable  DF    Prob > |T|

      INTERCEP   1        0.0008
      HEIGHT     1        0.0003

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

y = a + bx

где а = постоянный член уравнения, b = ее наклон.

Мы можем написать уравнение для наилучшей формы прямой линии определенной следующим образом:

ВЕС = - 592,64 + 11,19 x РОСТ

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

ВЕС = - 592,64 + 11,19 x 70 = 190,66 фунтов

Справа от оценок параметра находится столбец, обозначенный Standart Error, T for H0 : Parameter = 0, и Prob > |T|. Т- значение и связанные с ним вероятности (Вероятность (Prob) > |T|)), тестируют гипотезу о том, что параметр на самом деле является нулем. Иными словами, если бы истинный наклон и постоянный член уравнения были бы нулем, какова бы была вероятность получения, просто за счет действия случайных факторов, значения равные или больше тем, что мы обнаружили? Стандартная ошибка (Standart Error) может быть представлена примерно таким же образом как и стандартная ошибка среднего. Она отражает точность, с которой мы знаем истинный показатель наклона прямой и постоянный член уравнения.

В нашем случае наклон кривой составляет 11,19, а стандартная ошибка этого наклона 1,22. Поэтому мы можем сформировать 95% доверительный интервал для наклона прямой, взяв две (примерно), стандартных ошибки выше и ниже среднего значения. 95% доверительный интервал нашего значения наклона кривой в нашем случае составляет от 8,75 до 13,63. На самом деле, поскольку количество точек в нашем примере достаточно мало (n = 7), мы бы должны были бы пойти взять таблицу t-оценок и обнаружить количество стандартных ошибок выше и ниже среднего для 95% доверительного интервала. (Что справедливо, когда n меньше 30). Взяв эту самую таблицу t-оценок, мы смотрим на количество степеней свободы (df), которое равно n - 2 и уровень значимости, который был бы равен 0,05. Значение t для пяти степеней свободы и р = 0,05 - это 2,57. Таким образом наш 95% доверительный интервал составляет 11,19 плюс или минус 2,57 x 1,22 = 3,14.

Анализируя следующую часть распечатки, мы видим следующие значения для Root MSE, R-square, Dep Mean, Adj R-sq., и C.V. Root MSE - это квадратный корень дисперсии ошибки. Иными словами это стандартное отклонение остатков. R-square (R- квадрат) является ни чем иным, как квадратом множественного коэффициента корреляции. Поскольку у нас имеется всего лишь одна независимая переменная (рост) R- квадрат является квадратом Пирсоновского коэффициента корреляции между ростом и весом. Как упоминалось в предыдущем разделе квадрат корреляционного коэффициента информирует нас о том, какой процент вариации в зависимой переменной может быть объяснен вариабельностью независимой переменной. Когда имеется более одной переменной R- квадрат будет отражать вариабельность зависимой переменной, которая может быть объяснена линейной комбинацией всех независимых переменных. Dep Mean - это средняя зависимой переменной: в данном случае веса. C.V. - коэффициент вариации. Ну и наконец, Adj R-sq.- это квадрат коэффициента корреляции откорректированный на количество независимых переменных, внесенных в уравнение. Эта корректировка обычно уменьшает значение R- квадрат. Различия обычно не очень большие и они становятся более значимыми, когда мы работаем с большим количеством независимых переменных.

Разделение общей суммы квадратов.

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

Общая сумма квадратов - эта сумма, возведенных в квадрат отклонений веса каждого человека от среднего значения. Эта общая сумма квадратов может быть разделена на две части. Сумма квадратов вследствие регрессии или модели и сумма квадратов ошибки (иногда называемой остаточной). Одна часть, которая в распечатке называется Sum of Squares (ERROR) отражает отклонение каждого веса от предсказанного веса. Другая часть отражает отклонение между предсказанными значениями и средними. Эта часть называется суммой квадратов в результате действия модели. Столбец, помеченный средний квадрат (Mean Square) - это сумма квадратов, деленная на количество степеней свободы. В нашем случае имеется семь точек данных. Общее количество степеней свободы равняется n - 1 или 6. Модель имеет одну степень свободы. Количество степеней свободы для ошибки равняется общему количеству степеней свободы минус количество степеней свободы для модели, что дает нам 5. Мы можем рассматривать средний квадрат как соответствующую дисперсию (квадрат стандартного отклонения) для каждого из этих двух частей вариабельности. Наши интуиция подсказывает нам, что если отклонения вокруг регрессионной линии являются небольшими (небольшой средний квадрат ошибки) по сравнению с отклонениями между предсказанными значениями и средним (средние квадраты модели) - у нас имеется хорошая регрессионная линия. Для того, чтобы сравнивать средние квадраты формируется отношение:

F= Средний квадрат для модели/ Средний квадрат для ошибки.

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

Составление графиков точек на регрессионной линии.

Для того, чтобы нарисовать график наших данных по росу и весу мы можем использовать процедуру PLOT следующим образом:

PROC PLOT DATA=CORR_EG;
   PLOT WEIGHT*HEIGHT;
RUN;

Однако, в этом случае нам бы пришлось нарисовать регрессионную линию вручную, используя данные о постоянном члене уравнения и наклоне кривой, полученной при исполнении процедуры REG. Было бы желательно, чтобы программа SAS нарисовала регрессионную линию (на самом деле предсказанные регрессией значения) за нас. Процедура REG позволяет использовать команду PLOT. Форма этой команды следующая:

PLOT переменная y * переменная x = символ / опции ;

Переменная у и переменная х - это названия переменных, которые должны быть отложены по осям у и х, соответственно. Символ может быть либо какой-то отдельной буквой в одиночных кавычках, или названием переменной, чьи значения будут использоваться в качестве символа для рисования (так же как в процедуре PLOT). Имеется однако некие особые "имена переменных ", которые могут использоваться с командой PLOT. В особенности нас интересуют имена PREDICTED. (предсказанный) и RESIDUAL. (остатки). Обратите внимание, что после каждого имени идет точка! Эти слова используются для того, чтобы нарисовать предсказанные и остаточные значения. Довольно частым ключевым словом является слово OVERLAY, которое используется для того, чтобы нарисовать более одного графика в одной и той же системе координат.

Предположим, что мы хотели бы увидеть график взаимосвязи веса и роста в примере процедуры REG, описанном выше. В добавление мы бы хотели, чтобы эти точки были указаны на регрессионной линии. Мы хотели бы, чтобы значки, которыми мы обозначаем точки веса и роста, были бы звездочками (*), а символы, показывающие предсказанные значения на регрессионной линии, обозначались бы буквой 'P'. Процедура, которая сопровождается командой PLOT, показана в следующем отрывке кода, а за ним идет сгенерированный рисунок:

 PROC REG DATA=CORR_EG;
   MODEL WEIGHT = HEIGHT;
   PLOT PREDICTED.*HEIGHT = 'P' WEIGHT*HEIGHT='*' / OVERLAY;
RUN;

           +-----+------+------+------+------+------+------+-----†
  PRED     |                                                     |
P          |                                                     |
r          |                                                     |
e      225 +                                                     +
d          |                                               *     |
i          |                                               P     |
c          |                                        *            |
t      200 +                                                     +
e          |                                        P            |
d          |                                                     |
           |                                     P               |
V      175 +                                                     +
a          |                                 P   *               |
l          |                                                     |
u          |                                 *                   |
e      150 +                                                     +
           |                                                     |
o          |                       P                             |
f          |                                                     |
       125 +                       *                             +
W          |                                                     |
E          |                ?                                    |
I          |                                                     |
G      100 +         *                                           +
H          |                                                     |
T          |         P                                           |
           |                                                     |
        75 +                                                     +
           |                                                     |
           |                                                     |
           +-----+------+------+------+------+------+------+-----+
                60     62     64     66     68     70     72

                                   HEIGHT

 

Рисунки, содержащие остатки и доверительные интервалы.

Прежде, чем мы закончим этот раздел, мы хотели бы упомянуть, что программа SAS имеет возможности нарисовать различные расчетные значения, кроме предсказанных значений, зависимой переменной. Наиболее полезные статистики, за исключением предсказанных значений - это следующие:

 

RESIDUAL.

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

 

L95.

U95.

Нижний или верхний 95% доверительный интервал для индивидуального значения. (то есть 95% значений независимой переменной будут находиться между этими двумя значениями).

 

L95M.

U95M.

Нижний или верхний 95% доверительный интервал для среднего зависимой переменной данного значения независимой переменной.

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

(а) График оригинальных данных, предсказанного значения и верхнего и нижнего 95% доверительного интервала для среднего веса.

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

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

PROC REG DATA=CORR_EG;
   MODEL WEIGHT = HEIGHT;
   PLOT PREDICTED.*HEIGHT = 'P'
        U95M.*HEIGHT='-' L95M.*HEIGHT='-'
        WEIGHT*HEIGHT='*' / OVERLAY;
   PLOT RESIDUAL.*HEIGHT='o';
RUN;

Два графика приведены ниже:

           +-----+------+------+------+------+------+------+-----†
P PRED     |                                                     |
r      300 +                                                     +
e          |                                                     |
d          |                                                     |
i          |                                                     |
c          |                                               -     |
t          |                                               ?     |
e      200 +                                     -  ?      -     +
d          |                                 -   P  ?            |
           |                                 P   ?               |
V          |                       -         ?                   |
a          |                -      ?                             |
l          |         -      ?      -                             |
u      100 +         *      -                                    +
e          |         P                                           |
           |         -                                           |
o          |                                                     |
f          |                                                     |
           |                                                     |
W        0 +                                                     +
E          +-----+------+------+------+------+------+------+-----+
I               60     62     64     66     68     70     72
G
                                   HEIGHT

На графике звездочки (*) - это оригинальные данные. Р представляет собой предсказанные значения, а тире (-) верхний и нижний 95% доверительный интервал от среднего веса для данного роста. Вопросительные знаки (?) представляют собой множественные наблюдения (по всей вероятности, когда исходные данные и предсказанные значения очень близки друг к другу). Следующий график будет обсуждаться в следующем разделе.

Добавление квадратичного значения к уравнению регрессии.

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

           +-----+------+------+------+------+------+------+-----†
  RESIDUAL |                                                     |
           |                                                     |
        20 +                                                     +
           |                                                     |
           |                                        o            |
           |                                                     |
R       10 +         o                                           +
e          |                                               o     |
s          |                                                     |
i          |                o                                    |
d        0 +                                                     +
u          |                                                     |
a          |                                                     |
l          |                                                     |
       -10 +                       o             o               +
           |                                 o                   |
           |                                                     |
           |                                                     |
       -20 +                                                     +
           |                                                     |
           +-----+------+------+------+------+------+------+-----+
                60     62     64     66     68     70     72

                                   HEIGHT

или

HEIGHT2 = HEIGHT * * 2;

Затем Вы можете написать команды MODEL и PLOT:

PROC REG DATA=CORR_EG;
  MODEL WEIGHT=HEIGHT HEIGHT2
  PLOT RESIDUAL. * HEIGHT = 'o';
RUN;

Когда Вы запустите данную модель, Вы получите, что r-квадрат равняется 0,9743, что лучше, по сравнению с r-квадрат равным 0,9441, полученным в рамках линейной модели. Изучение остаточного графика показывает, что распределение остатков выглядит значительно более случайным, чем в более раннем варианте. Здесь правда необходимо сделать несколько предупреждений. Во-первых, помните, что это пример с очень небольшим количеством данных и оригинальный коэффициент корреляции уже был достаточно высоким. Во-вторых, хотя возможно ввести даже кубическую зависимость и так далее, необходимо помнить, что результаты должны быть легко интерпретируемыми.

           +-----+------+------+------+------+------+------+-----†
  RESIDUAL |                                                     |
        20 +                                                     +
           |                                                     |
           |                                                     |
           |                                        o            |
           |                                                     |
R          |                                                     |
e       10 +                                                     +
s          |                                                     |
i          |                                                     |
d          |                                                     |
u          |                o                                    |
a          |                                                     |
l        0 +                                                     +
           |         o             o                             |
           |                                                     |
           |                                     o         o     |
           |                                 o                   |
           |                                                     |
       -10 +                                                     +
           +-----+------+------+------+------+------+------+-----+
                60     62     64     66     68     70     72

                                   HEIGHT

Трансформирование данных .

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

Пациент

Доза лекарства

ЧСС

1

2

60

2

2

58

3

4

63

4

4

62

5

8

67

6

8

65

7

16

70

8

16

70

9

32

74

10

32

73

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

DATA HEART;
   INPUT DOSE HR;
DATALINES;        
 2        60
 2        58
 4        63
 4        62
 8        67
 8        65
16        70
16        70
32        74
32        73
;
PROC PLOT DATA=HEART;
   PLOT HR*DOSE='o';
RUN;

PROC REG DATA=HEART;
   MODEL HR = DOSE;
RUN;

Результирующий график представлен ниже:

              Plot of HR*DOSE.  Symbol used is 'o'.
HR |
   |
74 +                                                            o
73 +                                                            o
72 +
71 +
70 +                            o
69 +
68 +
67 +            o
66 +
65 +            o
64 +
63 +    o
62 +    o
61 +
60 +o
59 +
58 +o
   |
   ++---+---+---+---+---+---+---+---+---+---+---+---+---+---+---+
    2   4   6   8  10  12  14  16  18  20  22  24  26  28  30  32

                                DOSE

NOTE: 1 obs hidden.
Model: MODEL1
Dependent Variable: HR

                      Analysis of Variance
                         Sum of         Mean
Source      DF      Squares       Square      F Value  Prob>F

Model       1    233.48441    233.48441       49.006   0.0001
Error       8     38.11559      4.76445
C Total     9    271.60000

    Root MSE       2.18276     R-square       0.8597
    Dep Mean      66.20000     Adj R-sq       0.8421
    C.V.           3.29722

                       Parameter Estimates

                      Parameter      Standard    T for H0:
     Variable  DF      Estimate         Error   Parameter=0

     INTERCEP   1     60.708333    1.04491764        58.099
     DOSE       1      0.442876    0.06326447         7.000

     Variable  DF    Prob > |T|

Либо за счет клинической оценки результатов, либо внимательно изучая график, мы можем придти к выводу, что взаимоотношения между двумя показателями не являются линейными. Мы видим примерно равное увеличение в частоте сердечных сокращений каждый раз, когда доза удваивается. Поэтому мы можем нарисовать логарифм дозы против частоты сердечных сокращений тогда, когда мы ожидаем линейное взаимоотношение. В программе SAS имеется некоторое количество встроенных функций, таких как логарифмы и тригонометрические функции. Мы можем писать математические уравнения для того, чтобы определить новые переменные, размещая эти команды между командами INPUT и DATALINES. В программах SAS мы представляем сложение и вычитание, умножение и деление знаками +, - , * , и /, соответственно. Возведение в степень записывается как **. Для того, чтобы создать новую переменную, которая равна логарифму дозы, мы пишем следующее:

 DATA HEART;
   INPUT DOSE HR;
   LDOSE = LOG(DOSE);
DATALINES;

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

Мы можем теперь начертить зависимость между логарифмом дозы и частотой сердечных сокращений и подсчитать новое регрессионное уравнение:

PROC PLOT DATA=HEART;
   PLOT HR*LDOSE='o';
RUN;

PROC REG DATA=HEART;
   MODEL HR=LDOSE;
RUN;

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

                Plot of HR*LDOSE.  Symbol used is 'o'.
    HR |
       |
    74 +                                                  o
    73 +                                                  o
    72 +
    71 +
    70 +                                      o
    69 +
    68 +
    67 +                          o
    66 +
    65 +                          o
    64 +
    63 +              o
    62 +              o
    61 +
    60 +  o
    59 +
    58 +  o
       |
       +--+-----------+-----------+-----------+-----------+--
       0.6931      1.3863      2.0794      2.7726      3.4657

                                LDOSE

NOTE: 1 obs hidden.
Model: MODEL1
Dependent Variable: HR


                      Analysis of Variance

                         Sum of         Mean
Source          DF      Squares       Square      F Value Prob>F

Model            1    266.45000    266.45000      413.903 0.0001
Error            8      5.15000      0.64375
C Total          9    271.60000

    Root MSE       0.80234     R-square       0.9810
    Dep Mean      66.20000     Adj R-sq       0.9787
    C.V.           1.21199

                       Parameter Estimates

                      Parameter      Standard    T for H0:
     Variable  DF      Estimate         Error   Parameter=0

     INTERCEP   1     55.250000    0.59503151        92.852
     LDOSE      1      5.265837    0.25883212        20.345


     Variable  DF    Prob > |T|

     INTERCEP   1        0.0001
     LDOSE      1        0.0001

Обратите внимание на то, что теперь точки находятся значительно ближе к регрессионной линии. MEAN SQUARE ERROR меньше и r-квадрат больше, что подтверждает наше заключение о том, что дозировка зависит от частоты сердечных сокращений скорее в виде логарифмической кривой, а не просто линейной зависимости.

Подсчет индивидуальных кривых наклона.

В некоторых случая нам оказывается необходимо подсчитать кривые наклона для каждого пациента, входящего в исследование ( в том случае, когда у нас было несколько измерений), или графики наклона для подгруппы из всей популяции. Примером могло бы быть исследование, в котором каждый пациент относится либо к группе А либо к группе В. Мы затем три раза измеряем его или ее деятельность в рамках определенного теста, а далее мы бы хотели изучить скорость обучения (наклон кривой изменения оценок в зависимости от времени) и сравнить ее в двух группах. До 6 версии SAS эта проблема разрешалась достаточно сложно. Благодаря ответу компании SAS на запросы пользователей, мы сейчас можем формировать новый набор данных, содержащий наклоны кривых, полученных при помощи процедуры REG.

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

 DATA TEST;
   INPUT ID GROUP $ TIME SCORE;
DATALINES;
1 A 1 2
1 A 2 5
1 A 3 7
2 A 1 4
2 A 2 6
2 A 3 9
3 B 1 8
3 B 2 6
3 B 3 2
4 B 1 8
4 B 2 7
4 B 3 3
;
PROC SORT DATA=TEST;
   BY ID GROUP;
RUN;

PROC REG OUTEST=SLOPES DATA=TEST;
   BY ID GROUP;
   MODEL SCORE = TIME / NOPRINT;
RUN;

PROC PRINT DATA=SLOPES;
   TITLE 'Listing of the Data Set SLOPES';
RUN;

PROC TTEST DATA=SLOPES;
   TITLE 'Comparing Slopes Between Groups';
   CLASS GROUP;
   VAR TIME;
RUN;


Listing of the Data Set SLOPES                

OBS ID GROUP _MODEL_ _TYPE_ _DEPVAR_  _RMSE_ INTERCEP TIME SCORE

 1   1   A   MODEL1  PARMS   SCORE   0.40825  -0.3333  2.5   -1
 2   2   A   MODEL1  PARMS   SCORE   0.40825   1.3333  2.5   -1
 3   3   B   MODEL1  PARMS   SCORE   0.81650  11.3333 -3.0   -1
 4   4   B   MODEL1  PARMS   SCORE   1.22474  11.0000 -2.5   -1


Comparing Slopes Between Groups
               
                         TTEST PROCEDURE
Variable: TIME

GROUP       N             Mean          Std Dev        Std Error
----------------------------------------------------------------
A           2       2.50000000      0.00000E+00      0.00000E+00
B           2      -2.75000000      3.53553E-01      2.50000E-01

Variances        T       DF    Prob>|T|
---------------------------------------
Unequal    21.0000      1.0      0.0303
Equal      21.0000      2.0      0.0023

NOTE: All values are the same for one CLASS level.
 

Набор данных SLOPES, который бы создан процедурой REG, и соответствующие результаты t-теста, показаны далее. (Здесь мы себя немножко опережаем. Можно почитать по поводу сравнения средних при помощи процедуры t-тест в других разделах). Переменная время (TIME) представляет собой наклон кривой оценок в зависимости от времени (то есть коэффициенты времени). Причина, по которой переменная GROUP включена после команды ID, заключается в том, что нам необходимо, чтобы эта переменная была бы включена в, формирующийся в результате деятельности процедуры REG, файл. Можно также включить и другие переменные по которым пациент имеет только одно значение, например пол. Поскольку наклон кривой в данном примере является расчетом всех наблюдений для данного пациента, любая другая переменная, которая могла бы иметь другие значения для индивидуума, не может быть использована. Поэтому зависимая переменная в каждой модели имеет значение - 1.