Занятие 6. Корреляционный анализ. Линейная регрессия.
Алла Тамбовцева
Файлы для работы можно скачать здесь
Прежде чем обсуждать, как строить регрессионные модели и как их тестировать, вспомним о том, какие основные типы регрессий есть и в каких случаях они используются.
Если зависимая переменная количественная, то используется обычная линейная регрессия (ordinary least squares regression). Если зависимая переменная качественная, то используется логистическая регрессия (logistic regression) или пробит-регрессия (probit regression). При этом нужно понимать, что у логистической регрессии много разновидностей. В самом простом случае, когда нас интересует наличие/отсутствие признака (бинарный отклик), используется обычная логит-регрессия, для более сложных случаев, где зависимая переменная принимает больше двух значений, существует мультиномиальная логистическая регрессия и условная логистическая регрессия.
Есть еще один тип регрессий, который не так часто обсуждается, но может быть полезен, - тобит-регрессия. Используется она в случае, когда отклик является количественным, но при этом принимает значения на ограниченном промежутке. В связи с этим другое название модели – censored regression. Один из частых примеров ее использования – когда зависимая переменная представляет собой доли чего-либо или показатель, который принимает значение от 0 до 1. Имеет смысл использовать тобит-регрессию и тогда, когда распределение зависимой переменной сильно скошено, причем таким образом, что переменная часто принимает нулевые или сильно близкие к нулю значения.
Так как по анализу качественных данных будет отдельный курс, мы сконцентрируемся на количественных данных и обсудим линейную регрессию.
Однако построению моделей всегда предшествует этап предварительного анализа данных, в который входит, в частности, анализ корреляций между переменными.
Корреляционный анализ
Возьмем базу данных по социально-экономическим показателям разных стран в 2008 году, с которой раньше работали на курсе “Анализ политологических данных в пакете STATA”.
install.packages('foreign')library(foreign) # не забываем обращаться к библиотеке
df <- read.dta("C:/Users/AllaT/Desktop/Rclass/soc_ec.dta")Посмотрим на переменные:
summary(df)## country gdp gdpgr cpi
## Length:208 Min. : 321 Min. :-4 Min. :1.000
## Class :character 1st Qu.: 2210 1st Qu.: 1 1st Qu.:2.500
## Mode :character Median : 7388 Median : 2 Median :3.400
## Mean :12639 Mean : 3 Mean :4.022
## 3rd Qu.:16865 3rd Qu.: 5 3rd Qu.:5.200
## Max. :78599 Max. :12 Max. :9.300
## NA's :42 NA's :39 NA's :28
## fhpress gastil va ps
## Min. : 9.00 Min. :1.000 Min. :-2.240000 Min. :-3.280000
## 1st Qu.:23.00 1st Qu.:1.000 1st Qu.:-0.790000 1st Qu.:-0.595000
## Median :45.50 Median :2.000 Median : 0.050000 Median : 0.130000
## Mean :46.77 Mean :1.756 Mean :-0.003689 Mean :-0.002802
## 3rd Qu.:66.00 3rd Qu.:2.000 3rd Qu.: 0.890000 3rd Qu.: 0.830000
## Max. :98.00 Max. :3.000 Max. : 1.530000 Max. : 1.520000
## NA's :14 NA's :15 NA's :2 NA's :1
## gove rq rl
## Min. :-2.510000 Min. :-2.770000 Min. :-2.690000
## 1st Qu.:-0.737500 1st Qu.:-0.680000 1st Qu.:-0.770000
## Median :-0.140000 Median :-0.145000 Median :-0.200000
## Mean : 0.004418 Mean :-0.007549 Mean :-0.008164
## 3rd Qu.: 0.735000 3rd Qu.: 0.780000 3rd Qu.: 0.820000
## Max. : 2.530000 Max. : 2.000000 Max. : 1.960000
## NA's :2 NA's :4 NA's :1
## cc
## Min. :-1.900000
## 1st Qu.:-0.772500
## Median :-0.255000
## Mean :-0.005539
## 3rd Qu.: 0.682500
## Max. : 2.340000
## NA's :4Комментарий к переменным: gdp - ВВП, gdpgr - рост ВВП, cpi - индекс восприятия коррупции, fhpress - индекс свободы прессы (Freedom House), gastil - индекс Гастила, показатели Worldwide Governance indicators: va - подотчетность (Voice & Accountability), ps - политическая стабильность (Political Stability & Lack of Violence), gove - эффективность управления (Government Effectiveness), качество управления (Regulatory Quality), rl - верховенство закона (Rule of Law), cc - борьба с коррупцией (control of corruption).
Для удобства “закрепим” переменные, чтобы к ним было проще обращаться:
attach(df)Коэффициент корреляции
Для начала посмотрим на корреляцию каких-нибудь двух переменных:
cor(gdp, gdpgr) # ВВП и рост ВВП## [1] NAПочему коэффицент корреляции не рассчитывается? Потому что в базе есть пропущенные значения! Можем прописать следующую опцию:
cor(gdp, gdpgr, use = "complete.obs") # использовать все, кроме NA (complete observations)## [1] -0.318576Но лучше удалим все пропущенные значения сейчас, чтобы не было неожиданностей.
df <- na.omit(df)
attach(df)Как известно, существуют разные коэффициенты корреляции. Самые распространенные - линейный коэффициент корреляции Пирсона, коэффициент ранговой корреляции Спирмена и коэффициент ранговой корреляции Кендалла. По умолчанию считается коэффициент Пирсона, остальные можно получить, прописав дополнительный аргумент:
cor(va, ps, method = "spearman") # voice and accountability, political stability## [1] 0.6867763cor(va, ps, method = "kendal")## [1] 0.5099086Kорреляционная и ковариационная матрицы
Если хотим посмотреть на корреляцию “всего со всем”, можем указать столбцы в базе (переменные) и получить корреляционную матрицу:
cor(df[7:11])## va ps gove rq rl
## va 1.0000000 0.6347854 0.7630947 0.7864847 0.7835651
## ps 0.6347854 1.0000000 0.6713709 0.5874213 0.7419493
## gove 0.7630947 0.6713709 1.0000000 0.9279800 0.9488619
## rq 0.7864847 0.5874213 0.9279800 1.0000000 0.8846102
## rl 0.7835651 0.7419493 0.9488619 0.8846102 1.0000000Аналогично можем получить ковариационную матрицу:
cov(df[7:11])## va ps gove rq rl
## va 0.9393136 0.5686409 0.7167090 0.7246989 0.7334991
## ps 0.5686409 0.8543026 0.6013503 0.5161994 0.6623679
## gove 0.7167090 0.6013503 0.9391138 0.8549874 0.8881398
## rq 0.7246989 0.5161994 0.8549874 0.9039074 0.8123311
## rl 0.7334991 0.6623679 0.8881398 0.8123311 0.9329071Однако, чтобы получить значимость коэффициентов корреляции, потребуется пакет Hmisc.
install.packages('Hmisc')library(Hmisc)
# Внимание: функция rcorr требует на входе матрицу, а не data.fame!
rcorr(as.matrix(df[7:11]), type = 'spearman')## va ps gove rq rl
## va 1.00 0.69 0.76 0.76 0.78
## ps 0.69 1.00 0.68 0.58 0.77
## gove 0.76 0.68 1.00 0.91 0.93
## rq 0.76 0.58 0.91 1.00 0.85
## rl 0.78 0.77 0.93 0.85 1.00
##
## n= 158
##
##
## P
## va ps gove rq rl
## va 0 0 0 0
## ps 0 0 0 0
## gove 0 0 0 0
## rq 0 0 0 0
## rl 0 0 0 0Визуализация связи переменных
Часто связь между переменными хочется визуализировать. Строить диаграммы рассеяния (scatterplots) мы уже умеем.
plot(va, ps, xlab = 'Voice and accountability', ylab = 'Political stability')
Теперь добавим линию регрессии. Пусть подотчетность (va) будет независимой переменной, а политическая стабильность (ps) – зависимой.
plot(va, ps, xlab = 'Voice and accountability', ylab = 'Political stability')
# lm - от linear model,
# y ~ x, зависимая ~ независимая
abline(lm(ps ~ va), col = 'red')
Если корреляционный анализ - основная часть работы, то могут потребоваться и более продвинутые графики. Для этого нам потребуются пакеты car и gclus.
install.packages("car")
install.packages("gclus")Сначала построим матрицу диаграмм рассеяния (graph matrix), у которой на диагонали будут гистограммы для соответствующих переменных.
library(car)
# diagonal - что стоит на диагонали
# smooth = FALSE - обычная регрессия,
# smooth = TRUE - loess регрессия
scatterplotMatrix(df[7:11], diagonal = "histogram", smooth = FALSE)
При желании можем поставить на диагонали что-то еще, например, график плотности распределения. Заодно поменяем тип регрессии на loess (lowess).
LOWESS (locally weighted scatterplot smoothing) - взвешенная нелинейная регрессия, при которой подгонка модели происходит не на всей выборке, а на подвыборках. В результате мы получаем взвешенные оценки коэффициентов; визуально (для парной регрессии) - получаем уже не прямую, а кривую, которая показывает связь переменных при разных значениях независимой переменной. Если значения независимой переменной зависят от времени (или сама переменная - показатель времени), можно сказать, что кривая (lowess curve) показывает тренд.
# smooth= TRUE - lowess regression
scatterplotMatrix(df[7:11], diagonal = "density", smooth = TRUE)
А теперь делаем почти идеальный график и двигаемся дальше:
# создадим вектор с названиями переменных
labs = c("Voice & Accountability", "Political Stability",
"Government Effectiveness", "Regulatory Quality", "Rule of Law")
# reg.line = FALSE - нет регрессионной прямой, оставим только smooth (lowess curve)
# var.labels - названия переменных на диагонали
# cex.labels - размер шрифта для labels
# main - название графика
scatterplotMatrix(df[7:11], diagonal = "histogram", smooth = TRUE,
reg.line = FALSE,
var.labels = labs, cex.labels = 0.75,
main = "Correlations of Worldwide Governance Indicators")
Еще один вариант красивого графика для корреляций - разноцветный график, созданный с помощью пакета gclus.
library(gclus)
# получим вектор коэффициентов (по модулю)
coeffs <- abs(cor(df[7:11]))
# зададим цвета (автоматическое разбиение по вектору коэффициентов)
colors <- dmat.color(coeffs)
# отсортируем так, чтобы графики, где связь переменных наибольшая,
# были ближе к диагонали
order <- order.single(coeffs)
# строим сам график
# gap - расстояние между графиками в матрице
cpairs(df[7:11], order, panel.colors = colors, gap = .5,
main = "Correlations of Worldwide Governance Indicators" )
Теперь мы наконец можем перейти к самим регрессиям.
Линейная регрессия (OLS regression)
Для начала просто построим регрессию. Запишем уравнение регрессии, используя символику R:
gdpgr ~ ps + rl + gove + cc
Эта запись означает, что рост ВВП (gdpgr) - зависимая переменная, остальные (политическая стабильность, верховенство закона, эффективность государственного управления, борьба с коррупцией - показатели WGI) - независимые. Оценим модель и назовем ее m1:
# lm - от linear model
m1 <- lm(formula = gdpgr ~ ps + rl + gove + cc)Если просто выведем на экран m1, увидим только коэффициенты:
m1##
## Call:
## lm(formula = gdpgr ~ ps + rl + gove + cc)
##
## Coefficients:
## (Intercept) ps rl gove cc
## 2.8749 0.6791 -1.9981 1.3700 -1.3499Для полноценной выдачи нужна функция summary():
summary(m1)##
## Call:
## lm(formula = gdpgr ~ ps + rl + gove + cc)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.638 -1.883 -0.232 1.731 10.081
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.8749 0.2351 12.226 <2e-16 ***
## ps 0.6791 0.3787 1.793 0.0749 .
## rl -1.9981 1.0488 -1.905 0.0586 .
## gove 1.3700 0.8017 1.709 0.0895 .
## cc -1.3499 0.8530 -1.583 0.1156
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.905 on 153 degrees of freedom
## Multiple R-squared: 0.238, Adjusted R-squared: 0.2181
## F-statistic: 11.95 on 4 and 153 DF, p-value: 1.788e-08Как видим, показатели у модели плохие. И это не случайно. Чуть позже мы увидим, почему (в конце концов, на плохой модели удобно показывать, какие проблемы существуют, и к чему это приводит).
Посмотрим, какие показатели можно вытащить из выдачи:
attributes(m1)## $names
## [1] "coefficients" "residuals" "effects" "rank"
## [5] "fitted.values" "assign" "qr" "df.residual"
## [9] "xlevels" "call" "terms" "model"
##
## $class
## [1] "lm"Если понадобится, можем вызывать их по отдельности:
m1$coefficients## (Intercept) ps rl gove cc
## 2.8749158 0.6790695 -1.9981140 1.3700189 -1.3498723Как уже выяснили, модель плохая. Если вспоминать теорию по линейной регрессии, проблемы возникают, когда нарушаются условия Гаусса-Маркова. Основные проблемы, наличие которых надо обязательно проверять, следующие:
мультиколлинеарность
гетероскедастичность
наличие влиятельных наблюдений.
Обо всем по порядку.
Mультиколлинеарность
Мультиколлинеарность возникает, когда две и более независимых переменных сильно связаны друг с другом (разумеется, имеется в виду линейная связь).
Последствия мультиколлинеарности:
Если мультиколлинеарность строгая (коэффициент корреляции переменных по модулю равен 1), то невозможно оценить коэффициенты и рассчитать ошибки.
Если мультиколлинеарность нестрогая – неверные оценки коэффициентов и большие стандартные ошибки.
Неустойчивость оценок: возможны неверные выводы из-за больших стандартных ошибок (получаем заниженные t–статистики). Из-за этого при удалении нескольких или даже одного наблюдения из выборки можно получить совсем другие оценки коэффициентов, даже с противоположными знаками.
Индикаторы мультиколлинеарности:
Высокий R2 и незначимые коэффициенты.
Сильная парная корреляция предикторов.
Сильные частные корреляции предикторов.
Высокий VIF – variance inflation factor.
Как получить два первых индикатора – понятно: нужно посмотреть в выдачу или построить корреляционную матрицу для предикторов (и scatterplot заодно). С остальными будем разбираться.
Для расчета частного коэффициента корреляции нам понадобится пакет ggm.
install.packages("ggm")Смысл частного коэффициента корреляции:
Есть несколько переменных x1, x2, x3, x4, x5, все как-то связаны между собой. Хотим оценить связь между x1 и x2 “в чистом виде”, то есть получить коэффициент корреляции, “очищенный” от влияния других x. “Контролируем” эффект остальных переменных и считаем частный коэффициент корреляции.
library(ggm)
# смотрим на ЧКК переменных rl и gove,
# "контролируем" на эффект остальных
# var(df) - обязательно должны указать ковариационную матрицу
pcor(c("rl","gove","cc","ps"), var(df))## [1] 0.5341819Конечно, одних корреляций мало. Поэтому нужен более строгий показатель - VIF (variance inflation factor), показатель “вздутия” дисперсии. Он рассчитывается для каждого предиктора и показывает, насколько увеличивается дисперсия оценки коэффициента при предикторе из-за наличия связи с другими переменными. Средний VIF>5 (в крайнем случае VIF>10) говорит о том, что присутствует проблема мульколлинеарности.
library(car)
vif(m1)## ps rl gove cc
## 2.279098 19.087674 11.225831 13.280766В отличие от STATA, R не показывает средний VIF. Но его можно легко посчитать:
mean(vif(m1))## [1] 11.46834Способы борьбы с мультиколлинеарностью:
Избавление от предикторов, которые сильно коррелируют друг с другом (НО: осторожно, если исключаем много и не по делу, может возникнуть проблема пропущенных переменных).
Увеличение числа наблюдений (НО: возможно столкнуться с проблемой: увеличение числа наблюдений ведет к увеличению их изменчивости).
Метод главных компонент (если присутствуют переменные, которые логически связаны): объединить несколько связанных переменных в одну.
Ridge regression.
Гетероскедастичность
Гетероскедастичность - нарушение постоянства дисперсии ошибок.
Последствия гетероскедастичности:
Неверная проверка гипотез (и неверные выводы).
Получение менее эффективных оценок (эффективные оценки – с минимальной дисперсией).
Методы выявления гетероскедастичности:
Визуально (например, scatterplot независимая переменная vs остатки).
Анализ остатков: насколько их распределение близко к нормальному (проверка нормальности + гистограммы).
Тест Бреуша–Пагана:
H0: нет гетероскедастичности (гомоскедастичность), H1: есть гетероскедастичность.
- Тест Уайта:
H0: нет гетероскедастичности (гомоскедастичность), H1: есть гетероскедастичность.
Построим графики (к пункту 1):
plot(ps, m1$residuals)
plot(rl, m1$residuals)
plot(gove, m1$residuals)
plot(cc, m1$residuals)
По графикам сложно сказать, есть ли гетероскедастичность или нет. Явной неоднородности дисперсии ошибок не наблюдается, скорее всего, проблемы гетероскедастичности нет.
Посмотрим на более точный тест - тест Бреуша-Пагана (с тестом Уайта в R проблемы).
ncvTest(m1)## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 2.08712 Df = 1 p = 0.1485459Тест подтверждает наши догадки - нет оснований отвергнуть гипотезу о гомоскедастичности.
Тем не менее, надо знать, как с гетероскедастичностью бороться.
Способы борьбы с гетероскедастичностью:
Обобщенный метод наименьших квадратов (ОМНК)
Реализуемый метод наименьших квадратов (РОМНК)
“Взвешенный” метод наименьших квадратов – использование “весов”, которые влияют на оценки коэффициентов.
Использование стандартных ошибок, устойчивых к гетероскедастичности (вместо обычных стандартных ошибок по умолчанию).
Первые три метода используются редко, потому что они требуют довольно серьезных теоретических допущений, которые на практике выполняются нечасто.
Последний метод - стандартный. Однако тут надо помнить следующую вещь. Существуют 2 вида ошибок:
hc0 – если выборка большая
hc3 – если выборка маленькая (n<300).
Для использования таких ошибок необходим пакет sandwich и пакет lmtest.
install.packages("sandwich")
install.packages("lmtest")library(sandwich)
library(lmtest)
# по умолчанию - используются ошибки hc3
# выдает ковариационную матрицу для модели с такими ошибками
vcovHC(m1)## (Intercept) ps rl gove cc
## (Intercept) 0.056733392 0.01479564 0.008791928 -0.05536674 0.02857882
## ps 0.014795636 0.14907935 -0.178482912 0.03232969 0.05425082
## rl 0.008791928 -0.17848291 1.410536963 -0.51760961 -0.75130388
## gove -0.055366743 0.03232969 -0.517609613 0.80171807 -0.27662578
## cc 0.028578823 0.05425082 -0.751303882 -0.27662578 0.98711106# ковариационная матрица для обычной модели
vcov(m1)## (Intercept) ps rl gove cc
## (Intercept) 0.055292506 0.002315476 0.03452547 -0.02311236 -0.008790379
## ps 0.002315476 0.143428259 -0.14190269 0.04576574 -0.003563882
## rl 0.034525468 -0.141902686 1.10001374 -0.44913828 -0.523195308
## gove -0.023112361 0.045765740 -0.44913828 0.64266362 -0.198653290
## cc -0.008790379 -0.003563882 -0.52319531 -0.19865329 0.727567864Видим, что матрицы отличаются друг друга. А как включить эти ошибки в модель и переоценить ее?
# относится к пакету lmtest
coeftest(m1, vcov. = vcovHC)##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.87492 0.23819 12.0700 < 2e-16 ***
## ps 0.67907 0.38611 1.7588 0.08062 .
## rl -1.99811 1.18766 -1.6824 0.09453 .
## gove 1.37002 0.89539 1.5301 0.12806
## cc -1.34987 0.99353 -1.3587 0.17626
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Если нужны hc0 ошибки, указываем это и проделываем аналогичные операции.
vcovHC(m1, type = "HC0")## (Intercept) ps rl gove cc
## (Intercept) 0.052810703 0.01448119 0.006606439 -0.04831565 0.02432063
## ps 0.014481191 0.13601225 -0.161336973 0.02738409 0.04977251
## rl 0.006606439 -0.16133697 1.273368503 -0.46072579 -0.68068450
## gove -0.048315645 0.02738409 -0.460725792 0.71756670 -0.25213483
## cc 0.024320635 0.04977251 -0.680684496 -0.25213483 0.89564930coeftest(m1, vcov. = vcovHC(m1, type = "HC0"))##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.87492 0.22981 12.5102 < 2e-16 ***
## ps 0.67907 0.36880 1.8413 0.06751 .
## rl -1.99811 1.12844 -1.7707 0.07860 .
## gove 1.37002 0.84709 1.6173 0.10787
## cc -1.34987 0.94639 -1.4263 0.15581
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Здесь мы видим, что значимость коэффициентов не изменилась, изменились только стандартные ошибки коэффициентов, причем несильно. Но, конечно, так бывает далеко не всегда.
Выбросы и влиятельные наблюдения
Для начала просто посмотрим, есть ли выбросы:
outlierTest(m1)##
## No Studentized residuals with Bonferonni p < 0.05
## Largest |rstudent|:
## rstudent unadjusted p-value Bonferonni p
## 14 3.670067 0.00033504 0.052936Видно, что у нас есть один выброс – наблюдение номер 14. Однако одного отклонения от модельных значений недостаточно. Для того, чтобы наблюдение могло значительно исказить оценки коэффициентов регрессии, оно должно обладать большим потенциалом влияния. Одним из показателей влиятельности является мера Кука. Ее можно довольно легко посчитать:
cooks.distance(m1)## 1 2 3 4 5
## 1.190125e-03 2.676924e-03 2.575390e-02 1.349610e-03 3.593773e-03
## 6 7 8 9 10
## 3.913364e-03 1.313518e-03 1.226237e-02 1.982968e-03 2.892188e-02
## 11 12 13 14 15
## 1.892047e-04 1.669836e-03 1.867609e-03 8.562711e-02 3.364360e-06
## 16 17 18 19 20
## 3.420997e-03 9.063500e-03 1.111828e-04 3.853166e-03 4.935803e-04
## 21 22 23 24 25
## 4.832962e-03 2.925947e-03 2.223188e-03 6.885735e-04 6.073095e-03
## 26 27 28 29 30
## 4.760291e-02 8.703351e-04 1.112919e-02 6.687213e-03 3.939422e-02
## 31 32 33 34 35
## 1.322002e-03 1.271999e-04 3.822145e-04 1.827712e-02 1.749945e-03
## 36 37 38 39 40
## 6.964955e-04 3.368539e-03 8.165473e-04 1.614148e-04 7.206572e-05
## 41 42 43 44 45
## 1.687532e-04 4.605230e-03 4.435805e-03 5.408185e-03 4.805269e-02
## 46 47 48 49 50
## 1.326894e-02 2.324888e-02 2.422580e-04 2.671982e-04 1.966304e-02
## 51 52 53 54 55
## 2.471872e-04 1.262112e-04 6.808301e-04 1.366429e-04 6.194654e-05
## 56 57 58 59 60
## 8.371859e-03 2.143355e-02 1.462540e-03 9.691804e-04 1.058949e-02
## 61 62 63 64 65
## 3.074107e-03 2.936137e-04 1.393129e-03 1.596047e-02 4.573917e-04
## 66 67 68 69 70
## 3.379143e-05 1.228811e-02 3.553848e-03 8.736940e-03 3.232236e-02
## 71 72 73 74 75
## 2.522978e-03 1.348480e-04 1.287171e-02 6.890706e-03 1.444127e-04
## 76 77 78 79 80
## 2.079038e-03 5.038442e-04 2.531663e-02 1.570791e-02 2.527251e-07
## 81 82 83 84 85
## 3.810971e-03 2.191387e-04 9.196616e-04 3.197416e-03 6.958842e-04
## 86 87 88 89 90
## 7.646462e-04 7.425513e-03 1.867223e-04 6.705437e-05 5.789158e-04
## 91 92 93 94 95
## 6.571413e-03 9.439846e-03 7.241238e-03 4.776405e-03 1.074121e-03
## 96 97 98 99 100
## 2.204763e-04 1.085104e-03 6.092543e-04 6.892500e-03 7.130862e-03
## 101 102 103 104 105
## 1.473704e-03 1.050703e-03 4.178520e-04 9.220305e-04 2.107374e-03
## 106 107 108 109 110
## 3.945369e-03 2.725412e-04 1.692067e-04 2.011121e-02 4.489812e-03
## 111 112 113 114 115
## 1.944579e-03 9.969010e-04 7.994095e-03 6.969255e-03 1.494113e-02
## 116 117 118 119 120
## 4.444750e-05 4.991321e-05 5.920451e-02 6.215471e-05 2.200377e-05
## 121 122 123 124 125
## 4.166798e-03 1.454531e-03 3.291936e-03 5.688769e-02 5.114098e-03
## 126 127 128 129 130
## 1.658287e-08 7.396276e-04 3.822929e-03 5.732211e-04 7.329356e-05
## 131 132 133 134 135
## 3.462456e-02 3.904125e-03 2.145459e-03 3.084581e-03 3.395352e-04
## 136 137 138 139 140
## 1.250739e-03 3.740153e-03 5.929908e-04 8.079504e-05 6.650757e-04
## 141 142 143 144 145
## 1.542579e-02 2.708065e-02 9.854119e-03 7.700348e-04 3.170239e-04
## 146 147 148 149 150
## 4.574658e-05 4.939014e-03 2.634564e-03 5.009881e-04 1.005538e-05
## 151 152 153 154 155
## 6.490462e-02 1.809722e-05 2.196858e-03 8.003688e-03 9.255831e-03
## 156 157 158
## 1.373626e-04 4.768805e-03 2.834198e-04Однако в таком виде оценивать влиятельность наблюдений не очень удобно. Конечно, можно найти максимальное значение меры Кука, посмотреть на наблюдения с наибольшими значениями, но есть еще более удобный, стандартный способ - график “leverage vs residuals”.
# id.method = noteworthy - отметить наблюдения с наибольшими отклонениями
influencePlot(m1, id.method = "noteworthy", main = "Influence Plot",
sub = "Circle size is proportial to Cook's Distance")
## StudRes Hat CookD
## 14 3.670067 0.03323388 0.2926211
## 30 -1.339572 0.09937221 0.1984798Так как влиятельность наблюдения определяется отклонением от предсказанных значений и потенциалом влияния (leverage), влиятельные наблюдения надо искать в правом верхнем углу данного графика (большие остатки + высокий потенциал влияния).
В нашем случае влиятельных наблюдений нет. Есть выброс (14) и наблюдение с большим потенциалом влияния (30).
Кроме этого привычного графика в R есть целый набор графиков, позволяющих увидеть влиятельность наблюдений.
# после исполнения команды нужно перейти в окно консоли и следовать инструкциям
plot(m1, which = c(1:5))
Первый и третий графики позволяют не только выявить наблюдения, которые сильно отклоняются от предсказанных значений, но и визуально оценить, есть ли гетероскедастичность или нет. Последний график - тот же, что мы строили перед этим (leverage vs residuals), только на нем величина меры Кука отмечена не с помощью кругов разных размеров, а с помощью пунктирной линии. Здесь ее нет, потому что нет влиятельных наблюдений; вообще линия отчерчивает область в правом верхнем углу.
Можно вызвать любой график из пяти:
plot(m1, which = 4)
Основной метод борьбы с влиятельными наблюдениями - их исключение. Но при этом, конечно, исключать наблюдения надо обоснованно (если речь идет о межстрановых исследованиях, и стран не так много, понятно, что раскидываться странами будет опрометчиво). Выбросы тоже желательно удалять, но опять же с осторожностью - в идеале, для каждого нетипичного/ влиятельного наблюдения нужно проводить небольшое case study, чтобы понять, чем обусловленна нетипичность и стоит ли наблюдение удалять.
Выгрузка таблиц в LaTex
R позволяет выгружать в LaTex не только таблицы с результатами регрессии, но и почти любые таблицы, например, таблицу с описательными статистиками.
Для этого необходим пакет stargazer.
install.packages("stargazer")library(stargazer)
# команда ниже выдает код для LaTex
stargazer(df[7:11]) ##
## % Table created by stargazer v.5.1 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu
## % Date and time: Пт, фев 12, 2016 - 18:05:02
## \begin{table}[!htbp] \centering
## \caption{}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lccccc}
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## Statistic & \multicolumn{1}{c}{N} & \multicolumn{1}{c}{Mean} & \multicolumn{1}{c}{St. Dev.} & \multicolumn{1}{c}{Min} & \multicolumn{1}{c}{Max} \\
## \hline \\[-1.8ex]
## va & 158 & $-$0.049 & 0.969 & $-$2.200 & 1.530 \\
## ps & 158 & $-$0.098 & 0.924 & $-$2.610 & 1.520 \\
## gove & 158 & $-$0.046 & 0.969 & $-$1.890 & 2.530 \\
## rq & 158 & $-$0.036 & 0.951 & $-$2.130 & 1.920 \\
## rl & 158 & $-$0.096 & 0.966 & $-$1.680 & 1.960 \\
## \hline \\[-1.8ex]
## \end{tabular}
## \end{table}Полученный код можно скопировать в tex-файл и редактировать дальше, если есть необходимость.
Пусть теперь у нас есть несколько регрессий (не претендую на их осмысленность, просто для примера выгрузки таблицы).
m1 <- lm(gdpgr ~ ps + rl)
m2 <- lm(gdpgr ~ ps + rl + gove + cc)
m3 <- lm(gdpgr ~ ps + rl + gove + cc + gastil)
# align=TRUE - объдиняем все в одну таблицу
stargazer(m1, m2, m3, title = "Regression Results", align = TRUE)##
## % Table created by stargazer v.5.1 by Marek Hlavac, Harvard University. E-mail: hlavac at fas.harvard.edu
## % Date and time: Пт, фев 12, 2016 - 18:05:02
## % Requires LaTeX packages: dcolumn
## \begin{table}[!htbp] \centering
## \caption{Regression Results}
## \label{}
## \begin{tabular}{@{\extracolsep{5pt}}lD{.}{.}{-3} D{.}{.}{-3} D{.}{.}{-3} }
## \\[-1.8ex]\hline
## \hline \\[-1.8ex]
## & \multicolumn{3}{c}{\textit{Dependent variable:}} \\
## \cline{2-4}
## \\[-1.8ex] & \multicolumn{3}{c}{gdpgr} \\
## \\[-1.8ex] & \multicolumn{1}{c}{(1)} & \multicolumn{1}{c}{(2)} & \multicolumn{1}{c}{(3)}\\
## \hline \\[-1.8ex]
## ps & 0.596 & 0.679^{*} & 0.754^{**} \\
## & (0.377) & (0.379) & (0.381) \\
## & & & \\
## rl & -1.961^{***} & -1.998^{*} & -1.915^{*} \\
## & (0.361) & (1.049) & (1.046) \\
## & & & \\
## gove & & 1.370^{*} & 1.435^{*} \\
## & & (0.802) & (0.800) \\
## & & & \\
## cc & & -1.350 & -1.251 \\
## & & (0.853) & (0.852) \\
## & & & \\
## gastil & & & 0.574 \\
## & & & (0.392) \\
## & & & \\
## Constant & 2.902^{***} & 2.875^{***} & 1.890^{***} \\
## & (0.234) & (0.235) & (0.712) \\
## & & & \\
## \hline \\[-1.8ex]
## Observations & \multicolumn{1}{c}{158} & \multicolumn{1}{c}{158} & \multicolumn{1}{c}{158} \\
## R$^{2}$ & \multicolumn{1}{c}{0.217} & \multicolumn{1}{c}{0.238} & \multicolumn{1}{c}{0.249} \\
## Adjusted R$^{2}$ & \multicolumn{1}{c}{0.207} & \multicolumn{1}{c}{0.218} & \multicolumn{1}{c}{0.224} \\
## Residual Std. Error & \multicolumn{1}{c}{2.926 (df = 155)} & \multicolumn{1}{c}{2.905 (df = 153)} & \multicolumn{1}{c}{2.894 (df = 152)} \\
## F Statistic & \multicolumn{1}{c}{21.485$^{***}$ (df = 2; 155)} & \multicolumn{1}{c}{11.948$^{***}$ (df = 4; 153)} & \multicolumn{1}{c}{10.059$^{***}$ (df = 5; 152)} \\
## \hline
## \hline \\[-1.8ex]
## \textit{Note:} & \multicolumn{3}{r}{$^{*}$p$<$0.1; $^{**}$p$<$0.05; $^{***}$p$<$0.01} \\
## \end{tabular}
## \end{table}Здесь можно посмотреть, как будут выглядеть две таблицы.
Внимание:
в случае с таблицей регрессий в преамбуле в tex-файле необходимо прописать usepackage{dcolumn}, иначе файл не скомпилируется. Стоит обращать внимание на то, что прописано в качестве комментария в начале LaTex-овского кода, который появляется в консоли: там указаны пакеты, которые необходимо указать в преамбуле к документу (tex).