Сложно о простом, или снова о предпосылках Гаусса-Маркова

Автор: Артем Черёмухин, 
Нижегородский государственный инженерно-экономический университет 

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

Однако не все предпочитают помнить о том, что условия эффективности оценок метода наименьших квадратов (МНК) достаточно жестки, и после построения модели (сколь бы хорошей она на первый взгляд не получилась), желательно вернуться на шаг назад выполнить диагностику лежащих в основе этого метода допущений. Тема диагностики линейных моделей неоднократно затрагивалась в этом блоге (см., например, здесь и здесь). В этой статье мы поговорим о применении двух пакетов для R, предназначенных для определения выполнимости условий Гаусса-Маркова - gvlma (от "global validation of linear model assumptions"), и преобразования моделей с целью добиться выполнения этих условий в случае их несоблюдения - trafo (от "transforming linear regression models"). Оба эти пакета очень хорошо могут работать и в паре. За обоими пакетами стоит очень интересная математика. Для интересующихся - ссылки на статьи по определению выполнимости условий Гаусса-Маркова (Peña & Slate 2012) и по трансформации уравнений регрессии (Medina et al. 2019).

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

library(datasets)
library(tidyverse)
library(gvlma)
library(trafo)

Собственно, применение пакета gvlma достаточно простое – нам необходимо взять набор данных, построить линейную модель, и диагностировать ее. В качестве первого примера воспользуемся набором данных Pyromycid, в котором можно проанализировать связь между концентрации вещества (conc) и показателем скорости реакции клетки на эту концентрацию (rate). Код будет таким:

model_1 <- lm(rate ~ conc, data = Puromycin)
# Анализируем выполнимость условий Гаусса-Маркова:
model_1_glvma <- gvlma(model_1)
# Смотрим на результаты:
model_1_glvma

## Call:
## lm(formula = rate ~ conc, data = Puromycin)
## 
## Coefficients:
## (Intercept)         conc  
##       93.92       105.40  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_1) 
## 
##                       Value  p-value                   Decision
## Global Stat        11.53374 0.021177 Assumptions NOT satisfied!
## Skewness            0.03486 0.851898    Assumptions acceptable.
## Kurtosis            0.46425 0.495644    Assumptions acceptable.
## Link Function      10.57938 0.001144 Assumptions NOT satisfied!
## Heteroscedasticity  0.45525 0.499853    Assumptions acceptable.

Что мы в итоге видим? У нас есть модель, и несколько предпосылок выполняются в отношении ее стандартизированных остатков, а именно предположение о нормальности распределения (Skewness, Kurtosis) и однородности дисперсии остатков и отсутствия в них корреляции (Heteroscedasticity). Результаты теста на соответствие формы связи (Link Function) позволяют судить о выполнимости предпосылки о линейности модели. В нашем случае большие значения статистики Link Function говорят о том, что данная предпосылка не выполняется. В целом модель следует считать незначимой (большие значения параметра Global Stat).

Посмотрим на стандартную сводку с результатами подгонки модели:

summary(model_1)

## Call:
## lm(formula = rate ~ conc, data = Puromycin)
## 
## Residuals:
##    Min      1Q  Median      3Q     Max 
## -49.861 -15.247  -2.861  15.686  48.054 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)    93.92       8.00   11.74 1.09e-10 ***
## conc          105.40      16.92    6.23 3.53e-06 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 28.82 on 21 degrees of freedom
## Multiple R-squared:  0.6489,	Adjusted R-squared:  0.6322 
## F-statistic: 38.81 on 1 and 21 DF,  p-value: 3.526e-06

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

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

Учтем характер этой зависимости и построим новую модель:

model_1_1 <- lm(rate ~ log(conc), data = Puromycin)
model_1_1_glvma <- gvlma(model_1_1)
model_1_1_glvma

## Call:
## lm(formula = rate ~ log(conc), data = Puromycin)
## 
## Coefficients:
## (Intercept)    log(conc)  
##       190.1         33.2  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_1_1) 
## 
##                      Value p-value                Decision
## Global Stat        1.80241  0.7720 Assumptions acceptable.
## Skewness           0.05150  0.8205 Assumptions acceptable.
## Kurtosis           0.75004  0.3865 Assumptions acceptable.
## Link Function      0.09868  0.7534 Assumptions acceptable.
## Heteroscedasticity 0.90219  0.3422 Assumptions acceptable.

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

Рассмотрим другой случай. Возьмем данные по преступлениям в США (USArrests) и проанализируем зависимость количества случаев насилия (Rape) от доли городского населения (UrbanPop):

model_2 <- lm(Rape ~ UrbanPop, data = USArrests)
summary(model_2)

## Call:
## lm(formula = Rape ~ UrbanPop, data = USArrests)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -18.644  -5.476  -1.216   5.885  27.937 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  3.78707    5.71128   0.663    0.510   
## UrbanPop     0.26617    0.08513   3.127    0.003 **
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 8.626 on 48 degrees of freedom
## Multiple R-squared:  0.1692,	Adjusted R-squared:  0.1519 
## F-statistic: 9.776 on 1 and 48 DF,  p-value: 0.003001

model_2_glvma <- gvlma(model_2)
model_2_glvma

## Call:
## lm(formula = Rape ~ UrbanPop, data = USArrests)
## 
## Coefficients:
## (Intercept)     UrbanPop  
##      3.7871       0.2662  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_2) 
## 
##                      Value  p-value                   Decision
## Global Stat        13.4206 0.009393 Assumptions NOT satisfied!
## Skewness            5.3467 0.020762 Assumptions NOT satisfied!
## Kurtosis            3.8634 0.049351 Assumptions NOT satisfied!
## Link Function       0.1553 0.693511    Assumptions acceptable.
## Heteroscedasticity  4.0553 0.044034 Assumptions NOT satisfied!

Тут все по-иному: предпосылка линейности выполняется, а остальные предпосылки - нет. Уравнение регрессии в целом значимо, хотя его свободный член статистически не отличается от нуля. Самый логичный путь – это убрать из модели свободный член, и посмотреть на нее снова. Но тут есть два "но": 1) пакет gvlma не работает с моделями с одной независимой переменной и без свободного члена; 2) для более подходящей аппроксимации данных что-то нужно знать о породившем их процессе.

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

model_21 <- lm(log(Rape) ~ UrbanPop, data = USArrests)
summary(model_21)

## Call:
## lm(formula = log(Rape) ~ UrbanPop, data = USArrests)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13433 -0.24555  0.01242  0.29998  1.07483 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.068451   0.272542   7.589 9.24e-10 ***
## UrbanPop    0.013588   0.004062   3.345  0.00161 ** 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.4116 on 48 degrees of freedom
## Multiple R-squared:  0.189,	Adjusted R-squared:  0.1721 
## F-statistic: 11.19 on 1 and 48 DF,  p-value: 0.001605

model_21_glvma <- gvlma(model_21)
model_21_glvma

## Call:
## lm(formula = log(Rape) ~ UrbanPop, data = USArrests)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.13433 -0.24555  0.01242  0.29998  1.07483 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 2.068451   0.272542   7.589 9.24e-10 ***
## UrbanPop    0.013588   0.004062   3.345  0.00161 ** 
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.4116 on 48 degrees of freedom
## Multiple R-squared:  0.189,	Adjusted R-squared:  0.1721 
## F-statistic: 11.19 on 1 and 48 DF,  p-value: 0.001605

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

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

Airquality <- na.omit(airquality)
model_30 <- lm(Ozone ~ Temp + I(Temp^2), data = Airquality)
summary(model_30)

## Call:
## lm(formula = Ozone ~ Temp + I(Temp^2), data = Airquality)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -37.270 -12.462  -3.072   9.439 123.618 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 292.95885  123.92019   2.364 0.019861 *  
## Temp         -9.22680    3.25493  -2.835 0.005476 ** 
## I(Temp^2)     0.07602    0.02116   3.593 0.000494 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 22.71 on 108 degrees of freedom
## Multiple R-squared:  0.5426,	Adjusted R-squared:  0.5342 
## F-statistic: 64.07 on 2 and 108 DF,  p-value: < 2.2e-16

model_30_glvma <- gvlma(model_30)  
model_30_glvma

## Call:
## lm(formula = Ozone ~ Temp + I(Temp^2), data = Airquality)
## 
## Coefficients:
## (Intercept)         Temp    I(Temp^2)  
##   292.95885     -9.22680      0.07602  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_30) 
## 
##                        Value p-value                   Decision
## Global Stat        471.74977  0.0000 Assumptions NOT satisfied!
## Skewness            95.92542  0.0000 Assumptions NOT satisfied!
## Kurtosis           373.46965  0.0000 Assumptions NOT satisfied!
## Link Function        2.29185  0.1301    Assumptions acceptable.
## Heteroscedasticity   0.06285  0.8020    Assumptions acceptable.

И снова все хорошо, только не выполняется условие нормального распределения остатков и модель не проходит глобальный тест. Для устранения этих проблем можно воспользоваться пакетом trafo, в котором реализовано множество методов преобразований входных данных. Применим его:

# Выясняем, какие преобразования лучше всего работают в нашем случае:
assumptions(model_30, method = "ml", plotit = FALSE)

## The default lambdarange for the Log shift opt transformation is calculated dependent on the data range. The lower value is set to 0 and the upper value to 167
## 
## The default lambdarange for the Square-root shift transformation is calculated dependent on the data range. The lower value is set to 0 and the upper value to 167
## 
## Test normality assumption 
##               Skewness Kurtosis Shapiro_W Shapiro_p
## logshiftopt     0.2316   3.6886    0.9877    0.4113
## neglog         -0.2375   3.9129    0.9861    0.3098
## manly           0.3768   3.5127    0.9845    0.2262
## modulus         0.3059   4.0153    0.9816    0.1307
## yeojohnson      0.3059   4.0153    0.9816    0.1307
## gpower          0.3176   4.0989    0.9801    0.0969
## boxcox          0.3251   4.1555    0.9789    0.0763
## bickeldoksum    0.3250   4.1555    0.9789    0.0763
## glog           -0.4391   4.4473    0.9787    0.0726
## dual            0.3271   4.3290    0.9774    0.0567
## log            -0.5369   4.8554    0.9732    0.0245
## sqrtshift       0.9970   5.8095    0.9425    0.0001
## untransformed   2.2771  11.9861    0.8231    0.0000
## reciprocal     -7.5395  72.2639    0.4059    0.0000
## 
## Test homoscedasticity assumption 
##               BreuschPagan_V BreuschPagan_p
## dual                  0.6225         0.7325
## boxcox                0.7916         0.6731
## bickeldoksum          0.7916         0.6731
## gpower                0.8659         0.6486
## modulus               0.9824         0.6119
## yeojohnson            0.9824         0.6119
## logshiftopt           1.5724         0.4556
## sqrtshift             1.9392         0.3792
## untransformed         2.7731         0.2499
## manly                 3.5720         0.1676
## neglog                3.7980         0.1497
## glog                  5.5998         0.0608
## log                   6.2208         0.0446
## reciprocal           10.7378         0.0047

Параллельно будут выводится графики следующего вида (при нажатии на Enter):

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

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

# Преобразуем модель и сохраним преобразованную переменную в столбец х:
x <- logshiftopt(object = model_30, method = "div.ks", plotit = FALSE)
Airquality1 <- as.data.frame(x, row.names = NULL, std = FALSE)

# Новой преобразованной переменной автоматически присваивается исходное имя 
# с добавлением в конце буквы t:
model_30_1 <- lm(Ozonet ~ Temp + I(Temp^2), data = Airquality1)
summary(model_30_1)

## Call:
## lm(formula = Ozonet ~ Temp + I(Temp^2), data = Airquality1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.92571 -0.24386 -0.03017  0.23484  1.30209 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)   
## (Intercept)  5.5009744  2.0453446   2.690  0.00829 **
## Temp        -0.0942196  0.0537237  -1.754  0.08230 . 
## I(Temp^2)    0.0009172  0.0003493   2.626  0.00989 **
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.3749 on 108 degrees of freedom
## Multiple R-squared:  0.5982,	Adjusted R-squared:  0.5908 
## F-statistic:  80.4 on 2 and 108 DF,  p-value: < 2.2e-16

model_30_1_glvma <- gvlma(model_30_1)
model_30_1_glvma

## Call:
## lm(formula = Ozonet ~ Temp + I(Temp^2), data = Airquality1)
## 
## Coefficients:
## (Intercept)         Temp    I(Temp^2)  
##   5.5009744   -0.0942196    0.0009172  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_30_1) 
## 
##                     Value p-value                   Decision
## Global Stat        10.773 0.02924 Assumptions NOT satisfied!
## Skewness            3.424 0.06424    Assumptions acceptable.
## Kurtosis            4.523 0.03344 Assumptions NOT satisfied!
## Link Function       1.825 0.17672    Assumptions acceptable.
## Heteroscedasticity  1.001 0.31714    Assumptions acceptable.

По выведенному результату видим, что коэффициент при Temp не значим на 5%-ном уровне значимости. Попробуем убрать эту переменную из модели:

model_30_2 <- lm(Ozonet ~ I(Temp^2), data = Airquality1)
summary(model_30_2)

## Call:
## lm(formula = Ozonet ~ I(Temp^2), data = Airquality1)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.97586 -0.27015 -0.01469  0.23484  1.25610 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 1.924e+00  1.554e-01   12.38   <2e-16 ***
## I(Temp^2)   3.062e-04  2.461e-05   12.44   <2e-16 ***
## ---
## Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
## 
## Residual standard error: 0.3784 on 109 degrees of freedom
## Multiple R-squared:  0.5868,	Adjusted R-squared:  0.583 
## F-statistic: 154.8 on 1 and 109 DF,  p-value: < 2.2e-16

model_30_2_glvma <- gvlma(model_30_2)
model_30_2_glvma

## Call:
## lm(formula = Ozonet ~ I(Temp^2), data = Airquality1)
## 
## Coefficients:
## (Intercept)    I(Temp^2)  
##   1.9240568    0.0003062  
## 
## 
## ASSESSMENT OF THE LINEAR MODEL ASSUMPTIONS
## USING THE GLOBAL TEST ON 4 DEGREES-OF-FREEDOM:
## Level of Significance =  0.05 
## 
## Call:
##  gvlma(x = model_30_2) 
## 
##                    Value p-value                Decision
## Global Stat        7.387  0.1168 Assumptions acceptable.
## Skewness           1.518  0.2179 Assumptions acceptable.
## Kurtosis           2.302  0.1292 Assumptions acceptable.
## Link Function      2.511  0.1130 Assumptions acceptable.
## Heteroscedasticity 1.057  0.3040 Assumptions acceptable.

С преобразованными переменными все получилось. Если вам необходимо увидеть сами результаты оптимизации, то установите plotit = TRUE:

# Выводим график определения оптимального значения параметра преобразования:
x <- logshiftopt(object = model_30, method = "div.ks", plotit = TRUE)

# Выводим наиболее вероятные оптимальные значения параметра преобразования:
x$lambdavector[x$measvector == x$measoptim]

# [1] 11.31849 11.32500 11.35000

С учетом документации по пакету trafo и характера примененного преобразования мы теперь можем явно написать уравнение полученной зависимости: \[ \ln(y + 11.35) = 1.924 + 3.063 \cdot 10^{-4} \cdot x^2.\]

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