Летняя школа статистики 2019 :)

Линейная регрессия линейна по параметрам

\[ Y = f(X) + \epsilon \]

\[ f(X) = \hat{\beta}_0 + \sum_{j = 1}^p{X_j\hat{\beta}_j}, \]

где \( \hat{\beta}_0 \), \( \hat{\beta}_j \) – оценки параметров; \( X_j \) – регрессоры:

• непрерывные (количественные) переменные;
  
• базисные функции (\( \mathrm{log{X}} \), \( \sqrt{X} \), \( X^2 \));
  
• полиномиальные представления: \( X_2 = X_1^2, X_3 = X_1^3 \);
  
• фиктивные переменные (dummy);
  
• взаимодействия между переменными: \( X_3 = X_1 \cdot X_2 \).

\[ \scriptsize{ RSS(\beta) = \sum_{i = 1}^n{ \Big( y_i - f(x_i) \Big) ^2} = \sum_{i = 1}^n{\Big( y_i - \beta_0 - \sum_{j = 1}^p{x_{ij}\beta_j} \Big) ^2} \rightarrow \mathrm{min} } \]

МНК-оценки: \( \hat{\beta} = (X^TX)^{-1}X^T\mathrm{y} \)

Допущения:

• остатки случайны и соответствуют условиям Гаусса-Маркова;
  
• переменные \( \mathrm{x}_1,...,\mathrm{x}_j,...,\mathrm{x}_p \) некоррелированы.
  

Теорема Гаусса-Маркова: МНК-оценки обладают наименьшей дисперсией в классе линейных несмещённых оценок:

\[ Var(\hat{\beta}) = (X^TX)^{-1}\hat{\sigma}^2, \, \mathrm{где} \, \hat{\sigma}^2 = \hat{Var}(\epsilon) \]

Однако, если пожертвовать несмещённостью, можно уменьшить дисперсию оценок параметров (LASSO, ридж-регрессия)

Пример 1 (вишнёвые деревья): trees

• \( n = 31 \), \( p = 3 \);
  
• обучающая выборка: 85%; ядро для set.seed(): 12345;
  
• Volume – объём древесины, куб. футы;
  
• Girth – диаметр ствола, дюймы;
  
• Height – высота дерева, футы.

Дисперсия оценок и устойчивость модели

График остатков (Residuals vs Fitted)

Модель \( \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_{1} \cdot \mathrm{Girth} \)

Последствия гетероскедастичности

Модель \( \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_{1} \cdot \mathrm{Girth} \), тест Бройша-Пагана

    studentized Breusch-Pagan test

data:  model.1
BP = 3.5648, df = 1, p-value = 0.05902

На уровне значимости 0.01 принимается альтернативная гипотеза (\( H_1 \): гетероскедастичность есть), следовательно:

• нельзя доверять оценкам ошибок параметров;
  
• нельзя доверять проверке гипотез.
  

Модель \( \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_{1} \cdot \mathrm{Girth} \)

# ИСХОДНАЯ МОДЕЛЬ:
coeftest(model.1)

t test of coefficients:

            Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -33.6454     3.5138 -9.5753 1.139e-09 ***
Girth         4.7814     0.2640 18.1111 1.680e-15 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# ПЕРЕОЦЕНКА ОШИБОК ПАРАМЕТРОВ: 
coeftest(model.1, vcov = vcovHC(model.1, type = 'HC3'))

t test of coefficients:

             Estimate Std. Error t value  Pr(>|t|)    
(Intercept) -33.64542    3.49103 -9.6377 1.006e-09 ***
Girth         4.78138    0.28985 16.4958 1.347e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Влиятельные наблюдения и выбросы
Сгенерированные данные

20: неаномально по \( \mathrm{x} \), аномально по \( \mathrm{y} \)
41: аномально по \( \mathrm{x} \) и по \( \mathrm{y} \)

Влиятельные наблюдения и выбросы
Сгенерированные данные

20 не смещает прямую, но увеличивает RSS
41 смещает прямую и увеличивает RSS

Влиятельные наблюдения и выбросы
trees, модель \( \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{Girth} \)

plot(model.1, 4); plot(model.1, 5)

На данных примера 1 (вишнёвые деревья)

\[ \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{Girth} + \hat{\beta}_2 \cdot \mathrm{Height} \]

(1) Существует ли зависимость между \( Y \) и \( X \)-ами?

\( H_0: \beta_{1} = \beta_{2} = 0 \)
\( H_1: \beta_{1} \neq \beta_{2} \neq 0 \)

Проверка по F-критерию (распределение Фишера):

<...>
F-statistic: 216 on 2 and 23 DF, p-value: 1.24e-15

\[ \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{Girth} + \hat{\beta}_2 \cdot \mathrm{Height} \]

(2) Насколько сильна связь \( Y \) с \( X \)-ами?

<...>
Residual standard error: 3.542 on 23 degrees of freedom
Multiple R-squared: 0.9494 Adjusted R-squared: 0.945

\[ RSS = 3.542; {RSS \over \bar{y}} = {3.542 \over 28.365} = 12.5\%; R^2_{МНОЖ} = 94.9\% \]

\[ \mathrm{Volume} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{Girth} + \hat{\beta}_2 \cdot \mathrm{Height} \]

(3) Значимо ли влияние каждого из объясняющих факторов на \( Y \)?

\( H_0: \beta_j = 0 \) против \( H_1: \beta_j \neq 0 \) для \( j = 1,...,p \)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)   -56.90       8.78   -6.48   0.0000
Girth           4.41       0.27   16.48   0.0000
Height          0.37       0.13    2.83   0.0095

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

\[ \mathrm{Volume} = -56.903 + 4.406 \cdot \mathrm{Girth} + 0.375 \cdot \mathrm{Height} \]

(4) Каков размер эффекта каждого \( X \) на \( Y \)?

Доверительные интервалы для коэффициентов:

\[ \hat{\beta} - SE(\hat{\beta}) \cdot t(\alpha, \nu) < \beta < \hat{\beta} + SE(\hat{\beta}) \cdot t(\alpha, \nu) \]

             2.5 % 97.5 %
(Intercept) -75.06 -38.75
Girth         3.85   4.96
Height        0.10   0.65

\[ \hat{\beta} - SE(\hat{\beta}) \cdot t(\alpha, \nu) < \beta < \hat{\beta} + SE(\hat{\beta}) \cdot t(\alpha, \nu) \]

Проверка на мультиколлинеарность: VIF-коэффициенты

   Girth   Height 
1.325152 1.325152 

Проверка на гетероскедастичность: тест Бройша-Пагана

    studentized Breusch-Pagan test

data:  model.all
BP = 1.5337, df = 2, p-value = 0.4645

(5) Насколько точно можно спрогнозировать \( Y \)?

Доверительный интервал для среднего (\( f(X) \)):

\[ \hat{y_i} \pm t(\alpha, \nu) \cdot se_i; \, se_i = \hat{\sigma} \sqrt{(X^f)^T(X^TX)^{-1}X^f} \]

Доверительный интервал для наблюдения
(\( f(X) + \epsilon \)):

\[ \hat{y_i} \pm t(\alpha, \nu) \cdot sf_i; \, sf_i = \hat{\sigma} \sqrt{1 + (X^f)^T(X^TX)^{-1}X^f} \]

\[ \hat{\beta} - SE(\hat{\beta}) \cdot t(\alpha, \nu) < \beta < \hat{\beta} + SE(\hat{\beta}) \cdot t(\alpha, \nu) \]

Доверительный интервал для среднего:

        fit      lwr      upr
9  21.98803 19.46414 24.51191
18 33.93068 30.68077 37.18060

Доверительный интервал для наблюдения:

        fit      lwr      upr
9  21.98803 14.23854 29.73751
18 33.93068 25.91529 41.94607

(6) Является ли связь линейной?

(7) Есть ли эффект взаимодействия \( X \)-ов?

             Estimate Std. Error t value Pr(>|t|)
(Intercept)     83.41      27.42    3.04   0.0060
Girth           -7.59       2.30   -3.31   0.0032
Height          -1.47       0.36   -4.05   0.0005
Girth:Height     0.16       0.03    5.24   0.0000

Обратите внимание на знаки!

                Estimate Std. Error t value Pr(>|t|)
(Intercept)       -27.07       3.15   -8.58   0.0000
I(Girth)            1.28       0.87    1.47   0.1551
I(Girth):Height     0.04       0.01    4.13   0.0004

    studentized Breusch-Pagan test

data:  model.inter
BP = 0.82231, df = 2, p-value = 0.6629

Нелинейность чуть меньше, нет гетероскедастичности

Пример 2 (зарплаты, Москва, 2012): wages.ru

• salary – среднемесячная зарплата после вычета налогов за последние 12 месяцев (рублей);
  
• male – пол: 1 – мужчина, 0 – женщина;
• educ – уровень образования:
  • 1 – 0-6 классов,
  • 2 – незаконченное среднее (7-8 классов),
  • 3 - незаконченное среднее плюс что-то еще,
  • 4 – законченное среднее,
  • 5 – законченное среднее специальное,
  • 6 – законченное высшее образование и выше;
• forlang - иност. язык: 1 – владеет, 0 – нет;
• exper – официальный стаж c 1.01.2002 (лет).

Гипотеза: мужчины получают больше женщин

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{male} \)

Call:
lm(formula = salary ~ male, data = df.wage.train)

Residuals:
   Min     1Q Median     3Q    Max 
-36289 -18289  -8289   9787 231711 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)    30213       4050   7.459 1.28e-11 ***
male           18076       5706   3.168  0.00193 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 32150 on 125 degrees of freedom
Multiple R-squared:  0.07433,   Adjusted R-squared:  0.06692 
F-statistic: 10.04 on 1 and 125 DF,  p-value: 0.001929

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{male} \)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)    30213    4050.45    7.46   0.0000
male           18076    5705.78    3.17   0.0019

\( \mathrm{male} = 0 \): \( \hat{\mathrm{salary}} = 29679 + 15056 \cdot 0 = 29679 \)

\( \mathrm{male} = 1 \): \( \hat{\mathrm{salary}} = 29679 + 15056 \cdot 1 = 44735 \)

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{exper} + \hat{\beta}_2 \cdot \mathrm{male} \)

Call:
lm(formula = salary ~ exper + male, data = df.wage.train)

Residuals:
   Min     1Q Median     3Q    Max 
-38045 -15772  -5494   7592 229955 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)    15123      10064   1.503  0.13549   
exper           1728       1057   1.636  0.10443   
male           17637       5674   3.108  0.00233 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31940 on 124 degrees of freedom
Multiple R-squared:  0.09388,   Adjusted R-squared:  0.07926 
F-statistic: 6.423 on 2 and 124 DF,  p-value: 0.002216

Коэффициент при exper незначим.

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{exper} + \hat{\beta}_2 \cdot \mathrm{male} + \hat{\beta}_3 \cdot \mathrm{male} \cdot \mathrm{exper} \)

Call:
lm(formula = salary ~ exper * male, data = df.wage.train)

Residuals:
   Min     1Q Median     3Q    Max 
-38270 -15478  -6387   7845 229730 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)    16857      13368   1.261    0.210
exper           1530       1460   1.048    0.297
male           13905      19679   0.707    0.481
exper:male       421       2125   0.198    0.843

Residual standard error: 32060 on 123 degrees of freedom
Multiple R-squared:  0.09417,   Adjusted R-squared:  0.07207 
F-statistic: 4.262 on 3 and 123 DF,  p-value: 0.006686

Коэффициент при exper:male незначим, эффект взаимодействия нужно исключить.

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{male} \cdot \mathrm{exper} \)

Call:
lm(formula = salary ~ male * exper - male - exper, data = df.wage.train)

Residuals:
   Min     1Q Median     3Q    Max 
-38284 -17603  -8253   9716 229716 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  30252.5     3875.9   7.805 2.05e-12 ***
male:exper    2003.2      583.8   3.431 0.000816 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 31940 on 125 degrees of freedom
Multiple R-squared:  0.08607,   Adjusted R-squared:  0.07876 
F-statistic: 11.77 on 1 and 125 DF,  p-value: 0.0008156

Все параметры модели значимы.

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{male} \cdot \mathrm{exper} \)

            Estimate Std. Error t value Pr(>|t|)
(Intercept)    30253    3875.93    7.81   0.0000
male:exper      2003     583.85    3.43   0.0008

\( \mathrm{male} = 0 \): \( \hat{\mathrm{salary}} = 30253.00 \)

\( \mathrm{male} = 1 \): \( \hat{\mathrm{salary}} = 30253.00 + 2003.00 \cdot \mathrm{exper} \)

\( \hat{\mathrm{salary}} = \hat{\beta}_0 + \hat{\beta}_1 \cdot \mathrm{male} \cdot \mathrm{exper} \)