Кілька слів про «лінійної» регресії

    Іноді так буває: завдання можна вирішити чи не арифметично, а на розум перш за все приходять всякі інтеграли Лебега і функції Бесселя. Ось починаєш навчати нейронну мережу, потім додаєш ще парочку прихованих шарів, експериментуєш з кількістю нейронів, функціями активації, потім згадуєш про SVM і Random Forest і починаєш все спочатку. І все ж, незважаючи на прямо таки достаток цікавих статистичних методів навчання, лінійна регресія залишається одним з популярних інструментів. І для цього є свої передумови, не останнє місці серед яких займає интуитивность в інтерпретації моделі.
 
 
 Трохи формул
У простому випадку лінійну модель можна представити так:
 
y i = a 0 + a 1 x i + ε i
 
де a 0 — математичне сподівання залежної змінної y i , коли змінна x i дорівнює нулю; a 1 — очікувана зміна залежної змінної y i при зміні x i на одиницю (цей коефіцієнт підбирають таким чином, щоб величина ½Σ (y i i ) 2 була мінімальна — це так звана «функція невязки» ); ε i — випадкова помилка.
При цьому коефіцієнти a 1 і a 0 можна виразити через матан коефіцієнт кореляції Пірсона , стандартні відхилення і середні значення змінних x і y:
 
â 1 = cor (y, x) σ y / σ x
 
â 0 = ȳ — â 1
 
 Діагностика і помилки моделі
Щоб модель була коректною, необхідно виконання умов Гаусса-Маркова , тобто помилки повинні бути гомоскедастічни з нульовим математичним очікуванням. Графік залишків e i = y i — ŷ i допомагає визначити, наскільки адекватна побудована модель (e i можна вважати оцінкою ε i ).
Подивимося на графік залишків у випадку простої лінійної залежності y 1 ~ x:
 Прихований текст
set.seed(1)
n <- 100
x <- runif(n)
y1 <- x + rnorm(n, sd=.1)
fit1 <- lm(y1 ~ x)
par(mfrow=c(1, 2))
plot(x, y1, pch=21, col="black", bg="lightblue", cex=.9)
abline(fit1)
plot(x, resid(fit1), pch=21, col="black", bg="lightblue", cex=.9)
abline(h=0)

 
 
Залишки більш-менш рівномірно розподілені відносно горизонтальної осі, що говорить про «відсутність систематичної зв'язку між значеннями випадкового члена в будь-яких двох спостереженнях». А тепер досліджуємо такий же графік, але побудований для лінійної моделі, яка насправді не є лінійною:
 Прихований текст
y2 <- log(x) + rnorm(n, sd=.1)
fit2 <- lm(y2 ~ x)
plot(x, y2, pch=21, col="black", bg="lightblue", cex=.9)
abline(fit2)
plot(x, resid(fit2), pch=21, col="black", bg="lightblue", cex=.9)
abline(h=0)

 
 
За графіком y 2 ~ x начебто можна припустити лінійну залежність, але у залишків є патерн, а значить, чиста лінійна регресія тут не пройде . А ось що насправді означає гетероскедастичності :
 Прихований текст
y3 <- x + rnorm(n, sd=.001*x)
fit3 <- lm(y3 ~ x)
plot(x, y3, pch=21, col="black", bg="lightblue", cex=.9)
abline(fit3)
plot(x, resid(fit3), pch=21, col="black", bg="lightblue", cex=.9)
abline(h=0)

 
 
Лінійна модель з такими «раздувающимися» залишками некоректна. Ще іноді буває корисно побудувати графік квантилів залишків проти квантилів, які можна було б очікувати за умови, що залишки нормально розподілені:
 Прихований текст
qqnorm(resid(fit1))
qqline(resid(fit1))
qqnorm(resid(fit2))
qqline(resid(fit2))

 
 
На другому графіку чітко видно, що припущення про нормальність залишків можна відкинути (що знову таки говорить про некоректність моделі). А ще бувають такі ситуації:
 Прихований текст
x4 <- c(9, x)
y4 <- c(3, x + rnorm(n, sd=.1))
fit4 <- lm(y4 ~ x4)
par(mfrow=c(1, 1))
plot(x4, y4, pch=21, col="black", bg="lightblue", cex=.9)
abline(fit4)

 
 
Це так званий «викид» , який може сильно спотворити результати і привести до помилкових висновків. В R є кошти для його виявлення — за допомогою стандартизованность заходи dfbetas і hat values:
 
> round(dfbetas(fit4), 3)
    (Intercept)      x4
1        15.987 -26.342
2        -0.131   0.062
3        -0.049   0.017
4         0.083   0.000
5         0.023   0.037
6        -0.245   0.131
7         0.055   0.084
8         0.027   0.055
.....

 
> round(hatvalues(fit4), 3)
1     2     3     4     5     6     7     8     9     10...
0.810 0.012 0.011 0.010 0.013 0.014 0.013 0.014 0.010 0.010...

Як видно, перший член вектора x4 надає помітно більший вплив на параметри регресійної моделі, ніж решта, будучи, таким чином, викидом.
 
 Вибір моделі при множинної регресії
Природно, що при множинної регресії виникає питання: чи варто враховувати всі змінні? З одного боку, здавалося б, що варто, тому будь-яка змінна потенційно несе корисну інформацію. Крім того, збільшуючи кількість змінних, ми збільшуємо і R 2 (до речі, саме з цієї причини цю міру не можна вважати надійною при оцінці якості моделі). З іншого боку, коштувати пам'ятати про такі речі, як AIC і BIC , які вводять штрафи за складність моделі. Абсолютне значення інформаційного критерію само по собі не має сенсу, тому треба порівнювати ці значення у кількох моделей: в нашому випадку — з різною кількістю змінних. Модель з мінімальним значенням інформаційного критерію буде найкращою (хоча тут є про що посперечатися).
Розглянемо датасета UScrime з бібліотеки MASS:
 
library(MASS)
data(UScrime)
stepAIC(lm(y~., data=UScrime))

Модель з найменшим значенням AIC має такі параметри:
 
Call:
lm(formula = y ~ M + Ed + Po1 + M.F + U1 + U2 + Ineq + Prob, 
    data = UScrime)

Coefficients:
(Intercept)            M           Ed          Po1          M.F           U1           U2         Ineq         Prob  
  -6426.101        9.332       18.012       10.265        2.234       -6.087       18.735        6.133    -3796.032  

Таким чином, оптимальна модель з урахуванням AIC буде такою:
 
fit_aic <- lm(y ~ M + Ed + Po1 + M.F + U1 + U2 + Ineq + Prob, data=UScrime)
summary(fit_aic)

 
...
Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -6426.101   1194.611  -5.379 4.04e-06 ***
M               9.332      3.350   2.786  0.00828 ** 
Ed             18.012      5.275   3.414  0.00153 ** 
Po1            10.265      1.552   6.613 8.26e-08 ***
M.F             2.234      1.360   1.642  0.10874    
U1             -6.087      3.339  -1.823  0.07622 .  
U2             18.735      7.248   2.585  0.01371 *  
Ineq            6.133      1.396   4.394 8.63e-05 ***
Prob        -3796.032   1490.646  -2.547  0.01505 *  
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Якщо уважно придивитися, то виявиться, що у змінних MF і U1 досить високе значення p-value, що як би натякає нам, що ці змінні не так вже й важливі. Але p-value — досить неоднозначна міра при оцінки важливості тієї чи іншої змінної для статистичної моделі. Наочно цей факт демонструє приклад:
 
data <- read.table("http://www4.stat.ncsu.edu/~stefanski/NSF_Supported/Hidden_Images/orly_owl_files/orly_owl_Lin_9p_5_flat.txt")
fit <- lm(V1~. -1, data=data)
summary(fit)$coef

 
Estimate Std. Error   t value     Pr(>|t|)
V2  1.1912939  0.1401286  8.501431 3.325404e-17
V3  0.9354776  0.1271192  7.359057 2.568432e-13
V4  0.9311644  0.1240912  7.503873 8.816818e-14
V5  1.1644978  0.1385375  8.405652 7.370156e-17
V6  1.0613459  0.1317248  8.057300 1.242584e-15
V7  1.0092041  0.1287784  7.836752 7.021785e-15
V8  0.9307010  0.1219609  7.631143 3.391212e-14
V9  0.8624487  0.1198499  7.196073 8.362082e-13
V10 0.9763194  0.0879140 11.105393 6.027585e-28

p-values ​​у кожної змінної — практично нуль, і можна припустити, що всі змінні важливі для цієї лінійної моделі. Але насправді, якщо придивитися до залишків, виходить якось так:
 Прихований текст
plot(predict(fit), resid(fit), pch=".")

 
 
І все ж, альтернативний підхід грунтується на дисперсійному аналізі , в якому значення p-value грають ключову роль. Порівняємо модель без змінної MF з моделлю, побудованою з урахуванням тільки AIС ,:
 
fit_aic0 <- update(fit_aic,  ~ . - M.F)
anova(fit_aic0, fit_aic)

 
Analysis of Variance Table

Model 1: y ~ M + Ed + Po1 + U1 + U2 + Ineq + Prob
Model 2: y ~ M + Ed + Po1 + M.F + U1 + U2 + Ineq + Prob
  Res.Df     RSS Df Sum of Sq      F Pr(>F)
1     39 1556227                           
2     38 1453068  1    103159 2.6978 0.1087

Враховуючи P-значення, рівне 0.1087, при рівні значущості α = 0.05 ми можемо зробити висновок, що немає статистично значимого свідчення на користь альтернативної гіпотези, тобто на користь моделі з додатковою змінною MF
 
 Посилання:
 
     
  1. Чудовий курс Regression Models
  2.  
  3. Residual Plots and Data Sets
  4.  
  

    
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.