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


Линейная регрессия - в сущности, самая простая вещь в эконометрике и машинном обучении. Все просто: есть независимые переменные, есть зависимая переменная, есть точные формулы расчета коэффициентов регрессии, их значимости, значимости всего уравнения в целом - эти вычисления легко реализуются в 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 может еще много чего интересного. Но об этом поговорим в другой раз.

2 Комментарии

Андрей написал(а)…
Здравствуйте. Спасибо за интересную и очень полезную статью.
У меня такой вопрос. По какому методу оценивает гетероскедастичность и предположение о нормальности распределения остатков пакет glvma? Вопрос связан с тем, что я применил его к своей линейной модели. В результате модель можно применять. Однако, Breusch-Pagan test (ncvTest) показывал гетероскедастичность для этой модели, либо после преобразований ненормальное распределение остатков модели (shapiro.test(resid(m)) при отсутствии гетероскедастичности. Кстати все было похоже на экспоненциальное уравнение. Не подскажите, как их выводить в R с расчетом показателей уравнения?
Артем написал(а)…
Добрый день, Андрей. Интересный вопрос. Давайте по частям:
1) По поводу гетероскедастичности - в gvlma и тесте Бройша-Пагана реализованы просто разные подходы к определению этих явлений. Упрощенно говоря, в пакете гомо/гетероскедастичность определяется только на основании значений самих ошибок (исследуется связь с номером наблюдения), а в тесте Бройша-Погана исследуется взаимосвязь с исходными предикторами. По-моему мнению, если применять оба этих теста в связке, то об отсутствии гетероскедастичности можно говорить, только если оба этих теста дали такой результат
2) О нормальности ошибок - в gvlma расчитываются, опять же, упрощенно, скорректированные значения асимметрии и эксцесса распределения и они сравниваются с нормальным распределением. Возможно, я чего-то не понял, но если результаты применения теста Шапиро-Уилка и пакета gvlma, то желательно на график распределения ошибок посмотреть глазами непосредственно. Также, в литературе (у Кобзаря, например), указано, что в случае распределения, похожего на симметричное или нормальное, самым мощным тестом является тест Дэвида-Хартли-Пирсона - есть, наверное, смысл еще посмотреть на результаты применения иных тестов
3) Опять же, я мог не понять вопроса по экспоненциальной регрессии - но для этого вам необходимо просто взять логарифм зависимой переменной, или, наоборот, потенцировать независимые переменные (но могут получаться очень большие числа). Алгоритм проверки потом такой же.
Если что-то не ясно, задавайте вопросы еще.
Более новые Старые