Диагностика в линейной регрессии, часть 1 “Выбросы”
Раздел 6.1 учебника Фокса (стр. 192-201): Unusual Data
Представим себе модель y=b0+bx+e. Добавим в нее член gd, где g – коэффициент, а d – двоичная переменная, принимающая значение 1 для первого случая и 0 для остальных. t-тест, проверяющий гипотезу о том, что g=0, есть стьюдентизированный остаток для первого случая. Если |t|<2, то остаток явно мал, т.е. первый случай не отличается существенно от остальных случаев в выборке. Таким же образом получаются остатки и для всех остальных случаев.
Рассматриваем знакомый пример, анализ престижа профессий.
library(car)
## Loading required package: MASS
## Loading required package: nnet
attach(Duncan)
Duncan[1:10, ]
## type income education prestige
## accountant prof 62 86 82
## pilot prof 72 76 83
## architect prof 75 92 90
## author prof 55 90 76
## chemist prof 64 86 90
## minister prof 21 84 87
## professor prof 64 93 93
## dentist prof 80 100 90
## reporter wc 67 87 52
## engineer prof 72 86 88
scatterplot(prestige ~ income, id.method = list("identify"), lab = row.names(Duncan))
scatterplot(prestige ~ education, id.method = list("identify"), lab = row.names(Duncan))
mod.duncan <- lm(prestige ~ income + education)
qqPlot(mod.duncan, simulate = T, id.method = list("identify"), lab = row.names(Duncan))
Из mod.duncan берутся именно стьюдентизированные остатки, а аргумент simulate=T очерчивает 95%-ый доверительный участок (envelope) вокруг линии оптимума методом бутстрэп (bootstrap). Видим, что только один случай находится на границе участка, остальные вписываются. Этот случай – minister (пастор).
outlierTest(mod.duncan)
##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 6 3.135 0.003177 0.143
Функция outlier.test проверяет на значимость самый большой выброс в указанной модели (mod.duncan) методом Бонферрони, в котором двусторонняя вероятность нулевой гипотезы умножается на размер выборки. Смысл процедуры Бонферрони в том, что по теории вероятности мы должны ожидать какое-то количество выбросов, например 5% из 100%; чем больше выборка, тем больше будет выбросов. Если у нас всего одно наблюдение и вероятность выброса 5%, то при 2-х наблюдениях вероятность выброса уже 10% и т.д. Вероятность по Бонферрони поэтому правильнее. Видим, что вероятность по Бонферрони незначима (14%>5%); следовательно, поскольку это самый большой выброс, можно заключить, что проблем с выбросами в данной модели нет.
Функция outlierTest выдает номер строки в файле, для которой получен самый большой остаток: №6
Duncan[6, ] # это, конечно, пастор (minister)
## type income education prestige
## minister prof 21 84 87
Leverage – “рычаг” – влияние отдельного наблюдения на наклон регрессионной линии, зависящее в т.ч. от того, насколько далеко от центра наблюдений находится данное наблюдение. Стандартной мерой “рычага” являются hat values, названные так потому, что они связывают отдельное наблюдение Y с предсказанным значением Y^ (y-hat). То, как они рассчитываются, находится за пределами нашего понимания :-) Среднее значение hat value = (k+1)/n, к-во коэффициентов в модели, деленное на размер выборки. Величины, превышающие среднее значение в два раза (в малых выборках – в три раза), считаются проблемными в том смысле, что такие наблю- дения сильно влияют на регрессионные коэффициенты.
plot(hatvalues(mod.duncan))
На графике по оси абсцисс номера случаев, по ординатам – величины hat values.
Чтобы решить, какие случаи являются проблемными, проведем горизонтальные линии на уровне две и три средние величины hat values. Для этого умножим 2 и 3 на 3/45, т.к. к-во коэффициентов в модели 3, а размер выборки – 45.
plot(hatvalues(mod.duncan)) + abline(h = c(2, 3) * 3/45, lty = 2)
## numeric(0)
Теперь установим, какие случаи можно считать проблемными
plot(hatvalues(mod.duncan)) + abline(h = c(2, 3) * 3/45, lty = 2)
## numeric(0)
identify(1:45, hatvalues(mod.duncan), row.names(Duncan))
## integer(0)
(Первый аргумент определяет шкалу X, второй – шкалу Y, третий – ярлыки.)
Получается, что хотя самым большим выбросом является minister, наибольшим рычагом обладает RR.engineer (машинист); у conductor (проводника) рычаг тоже больше, чем у пастора.
Мера, которая объединяет величины выброса и рычага, – dfbeta(ij)=b(j-i)-b(j), где b(j) – регрессионный коэффициент для X(j), а b(j-i) – он же, но подсчитанный без i-того случая. (Независимые переменные: 1,2 … j … k; а случаи – 1,2, … i … n.)
dfbs.duncan <- dfbetas(mod.duncan)
создаем объект dfbs.duncan, в котором записаны все 45*3 dfbeta из данной модели.
dfbs.duncan[1:5, ]
## (Intercept) income education
## 1 -0.0225344 6.662e-04 0.0359439
## 2 -0.0254350 5.088e-02 -0.0081183
## 3 -0.0091867 6.484e-03 0.0056193
## 4 -0.0000472 -6.018e-05 0.0001398
## 5 -0.0658168 1.700e-02 0.0867771
plot(dfbs.duncan[, c(2, 3)])
identify(dfbs.duncan[, 2], dfbs.duncan[, 3], row.names(Duncan))
## integer(0)
Cлучай пастора увеличивает эффект образования и уменьшает эффект дохода на престиж. Случай машиниста, наоборот, уменьшает эффект образования и увеличивает эффект дохода на престиж профессии.
Проблема с dfbeta – их очень много. На их основе можно получить суммарный индекс – расстояния Кука (по-прежнему для каждого случая, но уже не для каждой независимой переменной в отдельности, а в целом для данного случая).
Расстояние Кука показывает, насколько данный случай меняет коэффициенты регрессии. Считается, что величина cooks.distance, превышающая 4/(n-k-1), является проблемной.
plot(cooks.distance(mod.duncan))
abline(h = 4/(45 - 2 - 1), lty = 2)
identify(1:45, cooks.distance(mod.duncan), row.names(Duncan))
## integer(0)
Далее Фокс предлагает красивый график, сочетающий информацию о стьюдентизированных остатках, рычагах и расстояниях Кука.
plot(hatvalues(mod.duncan), rstudent(mod.duncan), type = "n")
задаются координаты, абсциссы – рычаги, ординаты – остатки; аргумент type='n' указывает, что точки наносить на график не нужно (они будут нанесены позднее)
plot(hatvalues(mod.duncan), rstudent(mod.duncan), type = "n")
abline(h = c(-2, 0, 2), lty = 2)
чертятся три горизонтальных пунктира на высоте -2, 0 и +2; т.к. по ординатам остатки, пунктиры соответствуют пределам и средней стьюдентизированных остатков.
plot(hatvalues(mod.duncan), rstudent(mod.duncan), type = "n")
abline(h = c(-2, 0, 2), lty = 2)
abline(v = c(2, 3) * 3/45, lty = 2)
чертятся вертикальные пунктиры, соответствующие пределам рычагов
cook <- sqrt(cooks.distance(mod.duncan))
создается вектор корней из расстояний Кука
plot(hatvalues(mod.duncan), rstudent(mod.duncan), type = "n")
abline(h = c(-2, 0, 2), lty = 2)
abline(v = c(2, 3) * 3/45, lty = 2)
points(hatvalues(mod.duncan), rstudent(mod.duncan), cex = 10 * cook/max(cook))
теперь наносятся точки, не нанесенные в строке plot; их размер умножается на отношения расстояния Кука для данного наблюдения к max расстояний Кука. Таким образом размер точки соответствует расстоянию Кука. (Расстояния берутся из вектора cook.)
теперь можно идентифицировать проблемные наблюдения на графике
plot(hatvalues(mod.duncan), rstudent(mod.duncan), type = "n")
abline(h = c(-2, 0, 2), lty = 2)
abline(v = c(2, 3) * 3/45, lty = 2)
points(hatvalues(mod.duncan), rstudent(mod.duncan), cex = 10 * cook/max(cook))
identify(hatvalues(mod.duncan), rstudent(mod.duncan), row.names(Duncan))
## integer(0)
Такой график для переменной X1 образуется из остатков двух регрессионных моделей: 1) зависимость Y от всех X, кроме X1; 2) зависимость X1 от всех остальных X. Эти остатки распологаются в двумерной системе координат, где (1) откладывается по ординатам, а (2) – по абсциссам. В этой системе проводится регрессионная линия, которая имеет наклон b1 (как в полной регрессии), ее остатки – остатки многомерной модели.
Такие графики строятся далее для всех остальных X по очереди, включая двоичные переменные. В сущности этот график показывает частное влияние X на Y в многомерной регрессионной модели.
На таких графиках удобно наблюдать не только отдельные выбросы, но пары (или даже большие множества) наблюдений, который “раскачивают” регрессию.
avPlots(mod.duncan, id.method = list("identify"), labels = row.names(Duncan))
на графике income|others видно, что minister и conductor соместно “опускают” регрессию, т.е. снижают эффект дохода на престиж; RR.engineer ее слегка поднимает, но с меньшим успехом.
# -------------------- Задание (без проверки и оценки). Тем, кто
# использует в своей работе линейные модели, проверить наличие выбросов,
# сильно влияющих на регрессионные коэффициенты.
detach(Duncan)
ls()
## [1] "cook" "dfbs.duncan" "mod.duncan"
rm(cook, dfbs.duncan, mod.duncan)
library(car)
attach(Ornstein)
Ornstein[1:10, ]
## assets sector nation interlocks
## 1 147670 BNK CAN 87
## 2 133000 BNK CAN 107
## 3 113230 BNK CAN 94
## 4 85418 BNK CAN 48
## 5 75477 BNK CAN 66
## 6 40742 FIN CAN 69
## 7 40140 TRN CAN 46
## 8 26866 BNK CAN 16
## 9 24500 TRN CAN 77
## 10 23700 MIN US 6
table(interlocks)
## interlocks
## 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17
## 28 19 14 11 8 14 11 6 12 7 4 12 9 8 4 3 6 3
## 18 19 20 21 22 23 25 27 28 29 30 31 32 33 34 35 36 39
## 9 2 4 3 2 3 4 4 5 5 2 2 2 3 1 1 1 1
## 40 42 43 44 46 48 51 55 66 69 77 87 94 107
## 2 1 1 1 1 1 1 1 1 1 1 1 1 1
plot(table(interlocks))
mod.ornstein <- lm(interlocks + 1 ~ assets + nation + sector)
Вместо пуассоновской регрессии на этот раз запущена линейная регрессия (так делал сам автор исследования). Единица к зав. переменной добавлена Фоксом, чтобы было легче управляться с последующими преобразованиями переменной.
summary(mod.ornstein)
##
## Call:
## lm(formula = interlocks + 1 ~ assets + nation + sector)
##
## Residuals:
## Min 1Q Median 3Q Max
## -25.00 -6.60 -1.63 4.78 40.73
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.13e+01 1.56e+00 7.22 7.4e-12 ***
## assets 8.10e-04 6.12e-05 13.23 < 2e-16 ***
## nationOTH -1.24e+00 2.70e+00 -0.46 0.6456
## nationUK -5.78e+00 2.67e+00 -2.16 0.0318 *
## nationUS -8.62e+00 1.50e+00 -5.76 2.6e-08 ***
## sectorBNK -1.78e+01 5.91e+00 -3.02 0.0028 **
## sectorCON -4.71e+00 4.73e+00 -1.00 0.3203
## sectorFIN 5.15e+00 2.65e+00 1.95 0.0527 .
## sectorHLD 8.78e-01 4.00e+00 0.22 0.8267
## sectorMAN 1.15e+00 2.06e+00 0.56 0.5785
## sectorMER 1.49e+00 2.64e+00 0.57 0.5721
## sectorMIN 4.88e+00 2.07e+00 2.36 0.0190 *
## sectorTRN 6.17e+00 2.76e+00 2.24 0.0263 *
## sectorWOD 8.23e+00 2.68e+00 3.07 0.0024 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 9.83 on 234 degrees of freedom
## Multiple R-squared: 0.646, Adjusted R-squared: 0.627
## F-statistic: 32.9 on 13 and 234 DF, p-value: <2e-16
Anova(mod.ornstein)
## Anova Table (Type II tests)
##
## Response: interlocks + 1
## Sum Sq Df F value Pr(>F)
## assets 16904 1 175.06 < 2e-16 ***
## nation 3449 3 11.91 2.8e-07 ***
## sector 3706 9 4.26 3.9e-05 ***
## Residuals 22595 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
qqPlot(mod.ornstein, sim = T)
qq.plot остатков уже рассматривался в первой части занятия но еще нагляднее непараметрическая оценка плотности распределения, которую мы ранее использовали для оценки нормальности переменных.
plot(density(rstudent(mod.ornstein)), main = "Studentized residuals variation density")
#'studetized residuals variation density' = 'Плотность распределения стьюдентизированных остатков'
Остатки смещены вправо, что можно скорректировать понижением степени Y. Фокс пишет, что методом проб и ошибок можно прийти к тому, что подходящим преобразованием для Y в данном случае является квадратный корень.
library(MASS)
boxcox(mod.ornstein)
функция boxcox (не путать с boxCoxPowers!) подыскивает наилучшее преобразование для Y вида (Yл - 1)/л. Если наилучшее л=0, это означает, что наилучшим преобразованием будет натуральный логарифм Y. На графике указана оценка и 95%-ый интервал для л; по вертикальной оси отложена функция правдоподобия (ее можно понимать, как вероятность). Видно, что наилучшая оценка находится где-то между нулем и единицей. Чтобы точнее оценить эту величину, можно уменьшить диапазон л, который по умолчанию равен (-2; 2).
boxcox(mod.ornstein, lambda = seq(0.1, 0.5, by = 0.01))
теперь видно, что оценка л около 0,3 (примерно 0,29), а 95%-ый интервал оценки – примерно (0,18; 0,41).
plot(fitted.values(mod.ornstein), rstudent(mod.ornstein))
abline(h = 0, lty = 2)
видно, что размах остатков на уровне Y=20 больше, чем при меньших значениях Y
spreadLevelPlot(mod.ornstein)
## Warning: 1 negative fitted value removed
##
## Suggested power transformation: 0.4097
в отличие от предыдущего графика, в spread.level.plot по вертикальной оси отложены абсолютные величины (“модуль”) остатков. Кроме того, в основном окне предлагают степень трансформации Y, которая решила бы проблему гетероскедастичности. Эта степень примерно такая же, что и предложенная процедурой boxcox.
y <- 3 * ((interlocks)^(1/3) - 1)
mod.ornstein.y <- lm(y ~ assets + nation + sector)
spreadLevelPlot(mod.ornstein.y)
##
## Suggested power transformation: 1.204
ncvTest(mod.ornstein)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 46.99 Df = 1 p = 7.152e-12
данный тест формально проверяет гетероскедастичность; тест значим.
ncvTest(mod.ornstein.y)
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 6.843 Df = 1 p = 0.008899
с преобразованной Y дисперсию остатков можно считать постоянной
ncvTest(mod.ornstein, ~assets + nation + sector)
## Non-constant Variance Score Test
## Variance formula: ~ assets + nation + sector
## Chisquare = 74.74 Df = 13 p = 1.066e-10
здесь проверяется зависимость дисперсии ошибки не только от уровня Y, но и от трех независимых переменных. Разница в двух тестах – 12 ст.св. и примерно 27,75 пунктов в хи-квадрате. Критическое значение хи-кв при 12 ст.св на 5%-м уровне значимости 21,03. Значение 27>21, из чего следует вывод, что дисперсия ошибки зависит не только от Y, но и хотя бы от одной из трех независимых переменных.
Фокс пишет, что отчасти проблема связана с нелинейной зависимостью Y и assets. Проверим:
scatterplot(interlocks ~ assets)
нелинейность видна в начале диапазона, где сосредоточена основная масса наблюдений
логарифмическая трансформация решает проблему.
scatterplot(log(interlocks + 1) ~ log(assets))
Основной подход, предлагаемый Фоксом – трансформация Y. Традиционный подход – weighted least squares.
ncvTest(mod.ornstein, ~assets)
## Non-constant Variance Score Test
## Variance formula: ~ assets
## Chisquare = 20.17 Df = 1 p = 7.087e-06
ncvTest(mod.ornstein, ~nation)
## Non-constant Variance Score Test
## Variance formula: ~ nation
## Chisquare = 32.49 Df = 3 p = 4.131e-07
ncvTest(mod.ornstein, ~sector)
## Non-constant Variance Score Test
## Variance formula: ~ sector
## Chisquare = 18.48 Df = 9 p = 0.03001
ncvTest(mod.ornstein, ~assets + nation)
## Non-constant Variance Score Test
## Variance formula: ~ assets + nation
## Chisquare = 44.42 Df = 4 p = 5.246e-09
scatterplot(log(abs(rstudent(mod.ornstein))) ~ log(fitted(mod.ornstein) + 2))
воспроизвели spread.level.plot
scatterplot(log(abs(rstudent(mod.ornstein))) ~ log(assets))
assets являются источником гетероскедастичности
summary(lm(interlocks ~ assets + nation + sector, weights = 1/assets))
##
## Call:
## lm(formula = interlocks ~ assets + nation + sector, weights = 1/assets)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.5727 -0.1295 -0.0172 0.1214 1.0793
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.02e+00 9.01e-01 7.79 2.2e-13 ***
## assets 1.41e-03 2.55e-04 5.54 8.1e-08 ***
## nationOTH -1.35e+00 2.09e+00 -0.65 0.5185
## nationUK -3.05e+00 1.66e+00 -1.84 0.0677 .
## nationUS -5.90e+00 9.77e-01 -6.04 6.0e-09 ***
## sectorBNK -2.95e+01 1.60e+01 -1.84 0.0667 .
## sectorCON -4.30e+00 2.45e+00 -1.75 0.0808 .
## sectorFIN 3.77e+00 5.36e+00 0.70 0.4827
## sectorHLD -1.87e+00 3.22e+00 -0.58 0.5619
## sectorMAN 1.03e+00 1.16e+00 0.89 0.3754
## sectorMER 2.21e+00 1.74e+00 1.27 0.2044
## sectorMIN 2.78e+00 1.28e+00 2.16 0.0315 *
## sectorTRN 1.30e+00 2.39e+00 0.54 0.5870
## sectorWOD 5.02e+00 1.93e+00 2.61 0.0097 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.24 on 234 degrees of freedom
## Multiple R-squared: 0.298, Adjusted R-squared: 0.259
## F-statistic: 7.66 on 13 and 234 DF, p-value: 1.36e-12
Anova(lm(interlocks ~ assets + nation + sector, weights = 1/assets))
## Anova Table (Type II tests)
##
## Response: interlocks
## Sum Sq Df F value Pr(>F)
## assets 1.77 1 30.70 8.1e-08 ***
## nation 2.18 3 12.55 1.2e-07 ***
## sector 1.25 9 2.41 0.012 *
## Residuals 13.52 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod.ornstein)
## Anova Table (Type II tests)
##
## Response: interlocks + 1
## Sum Sq Df F value Pr(>F)
## assets 16904 1 175.06 < 2e-16 ***
## nation 3449 3 11.91 2.8e-07 ***
## sector 3706 9 4.26 3.9e-05 ***
## Residuals 22595 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
обратите внимание, как изменилась сумма квадратов остатков!
еще один способ – корректировка коэффициентов методом Уайта с помощью heteroscedacticity-consistent standard errors. На основе этого метода создана функция hccm (heteroscedascticity-consistent covariance matrix) в пакете car. Имеется 5 способов корректировки: hc0, hc1, hc2, hc3, hc4. (Для понимания метода нужна матричная алгебра.)
help(hccm)
Anova(mod.ornstein, white.adjust = "hc3")
## Analysis of Deviance Table (Type II tests)
##
## Response: interlocks + 1
## Df F Pr(>F)
## assets 1 46.56 7.5e-11 ***
## nation 3 12.99 7.0e-08 ***
## sector 9 3.89 0.00013 ***
## Residuals 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Anova(mod.ornstein)
## Anova Table (Type II tests)
##
## Response: interlocks + 1
## Sum Sq Df F value Pr(>F)
## assets 16904 1 175.06 < 2e-16 ***
## nation 3449 3 11.91 2.8e-07 ***
## sector 3706 9 4.26 3.9e-05 ***
## Residuals 22595 234
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
результаты несколько отличаются: эффект assets выглядит слабее, хотя все равно является самым значимым. (Такой же эффект, что и метод взвешенных наименьших квадратов.)
detach()
ls()
## [1] "mod.ornstein" "mod.ornstein.y" "y"
rm(mod.ornstein, mod.ornstein.y, y)