Регрессионный анализ
Как выглядят данные
## state metro_res white hs_grad poverty female_house
## 1 Alabama 55.4 71.3 79.9 14.6 14.2
## 2 Alaska 65.6 70.8 90.6 8.3 10.8
## 3 Arizona 88.2 87.7 83.8 13.3 11.1
## 4 Arkansas 52.5 81.0 80.9 18.0 12.1
## 5 California 94.4 77.5 81.1 12.8 12.6
## 6 Colorado 84.5 90.2 88.7 9.4 9.6Описание переменных:
state - название штата
metro_res - процент людей, проживающих в столичной области этого штата
white - процент белого населения
hs_grad - процент людей с высшим образованием
proverty - процент людей за чертой бедности
female_house - процент домохозяек
Будем исследовать влияние всех факторов на уровень бедности.
Строим скаттерплоты
Прологорифмируем переменные, распределение которых отличается от нормального
df$hs_grad <- log(df$hs_grad)
df$poverty <- log(df$poverty)
df$female_house <- log(df$female_house)
Обратим внимание на выбросы и на достаточно большую отрицательную корреляцию между female_house - hs_grade и female_house - white
Строим модель множественной линейной регрессии
fit <- lm(poverty~ hs_grad + female_house + white + metro_res, df)##
## Call:
## lm(formula = poverty ~ hs_grad + female_house + white + metro_res,
## data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.35393 -0.13438 -0.00669 0.11311 0.38447
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 25.296875 4.382403 5.772 6.34e-07 ***
## hs_grad -4.818029 0.850596 -5.664 9.20e-07 ***
## female_house -0.264650 0.272951 -0.970 0.3373
## white -0.006363 0.002931 -2.171 0.0351 *
## metro_res -0.003906 0.001809 -2.159 0.0361 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 0.175 on 46 degrees of freedom
## Multiple R-squared: 0.6195,
## Adjusted R-squared: 0.5864
## F-statistic: 18.72 on 4 and 46 DF, p-value: 3.4e-09Незначимым получился коэфициент при female_house 
Outliers
Найдем выбросы по Куку
cook <- cooks.distance(fit)
out_cook <- as.numeric(names(cook[cook > 4/(nrow(df)- ncol(df) - 2)]))
out_cook## [1] 9 12 30Найдем выбросы по Махаланобису
dist_mah <- mahalanobis(df[,-1],colMeans(df[,-1]),cov(df[,-1]))
dist_mah## [1] 5.192060 4.132341 4.033433 4.478637 5.711998 2.255827 3.537904
## [8] 4.227607 18.060438 2.001490 2.591124 33.111897 3.478634 1.447457
## [15] 1.550477 2.599909 1.306799 2.357585 4.436872 4.623875 6.620225
## [22] 3.234320 2.365627 4.219125 11.134768 1.954967 7.756225 2.333047
## [29] 4.034742 8.508540 4.274844 5.543738 4.589793 3.856887 6.911343
## [36] 1.925615 2.370502 1.946499 1.085312 6.745484 4.518224 2.989707
## [43] 2.137484 10.362756 3.722855 6.097001 1.792793 2.671126 9.300806
## [50] 1.177789 2.681490Видим, что сильно отличаются 9 и 12 элементы. Они же у нас были в outliers по Куку.
Уберем выбросы и построим скаттерплоты
Визуально видно, что выбросов стало меньше.
Redundancy
red <- redun(~metro_res + white + hs_grad + female_house, df2, type = "ordinary", r2 = 0.0000001)
pcors <- pcor(df2[-1])$estimate[4,-4]
semi.pcors <- spcor(df2[-1])$estimate[4,-4]
red.pcors <- data.frame(cbind(red$rsq1, pcors,semi.pcors))
colnames(red.pcors) <- c("R_square","partial correlation", "semi-partial correlation")
rownames(red.pcors) <- c("metro_res", "white", "hs_grade", "female_house")
print(red.pcors)## R_square partial correlation semi-partial correlation
## metro_res 0.4050353 -0.3927009 -0.25061973
## white 0.8058021 -0.1066424 -0.06295022
## hs_grade 0.6372115 -0.6582902 -0.51326542
## female_house 0.8911231 -0.1304459 -0.07722199Тут можно убрать признаки female_house и white. Уберем только female_house
Снова построим линейную регрессию
fit2 <- lm(poverty~hs_grad+ white + metro_res, df2)##
## Call:
## lm(formula = poverty ~ hs_grad + white + metro_res, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.24815 -0.12485 -0.02112 0.11280 0.33805
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.865814 2.461111 8.885 2.25e-11 ***
## hs_grad -4.291186 0.574672 -7.467 2.38e-09 ***
## white 0.000139 0.002734 0.051 0.95969
## metro_res -0.005265 0.001568 -3.357 0.00163 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 0.1563 on 44 degrees of freedom
## Multiple R-squared: 0.6496,
## Adjusted R-squared: 0.6257
## F-statistic: 27.18 on 3 and 44 DF, p-value: 4.207e-10Здесь оказался незначимым коэффициент metro_res. Уберем его из модели
fit2 <- lm(poverty~hs_grad + metro_res, df2)##
## Call:
## lm(formula = poverty ~ hs_grad + metro_res, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2480 -0.1249 -0.0210 0.1128 0.3382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.823186 2.288032 9.538 2.24e-12 ***
## hs_grad -4.278756 0.514247 -8.320 1.18e-10 ***
## metro_res -0.005280 0.001522 -3.470 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 0.1546 on 45 degrees of freedom
## Multiple R-squared: 0.6495,
## Adjusted R-squared: 0.634
## F-statistic: 41.7 on 2 and 45 DF, p-value: 5.683e-11Получили модель \(poverty = -4.2\cdot hs\_grad -0.005\cdot metro\_res + 21.8\)
Теперь выполним пошаговую регрессию.
fit_stepwise <- lm(poverty~hs_grad + female_house + white + metro_res ,df2)
forward(fit_stepwise, 0.05, full=TRUE)## Forward selection, alpha-to-enter: 0.05
##
## Full model: poverty ~ hs_grad + female_house + white + metro_res
## <environment: 0x77f41e8>
##
## --= Step 1 =--
## Single term additions
##
## Model:
## poverty ~ 1
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 3.0684 -130.00
## hs_grad 1 1.70528 1.3631 -166.95 57.5459 1.214e-09 ***
## female_house 1 0.49623 2.5722 -136.47 8.8745 0.004605 **
## white 1 0.19320 2.8752 -131.12 3.0910 0.085380 .
## metro_res 1 0.33864 2.7298 -133.62 5.7065 0.021061 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## --= Step 2 =--
## Single term additions
##
## Model:
## poverty ~ hs_grad
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1.3631 -166.95
## female_house 1 0.069675 1.2935 -167.47 2.4240 0.126495
## white 1 0.012365 1.3508 -165.38 0.4119 0.524243
## metro_res 1 0.287759 1.0754 -176.33 12.0415 0.001159 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## --= Step 3 =--
## Single term additions
##
## Model:
## poverty ~ hs_grad + metro_res
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1.0754 -176.33
## female_house 1 0.0062015 1.0692 -174.61 0.2552 0.6159
## white 1 0.0000631 1.0753 -174.33 0.0026 0.9597##
## Call:
## lm(formula = poverty ~ hs_grad + metro_res, data = df2)
##
## Coefficients:
## (Intercept) hs_grad metro_res
## 21.82319 -4.27876 -0.00528backward(fit_stepwise, 0.05, full=TRUE)## Backward elimination, alpha-to-remove: 0.05
##
## Full model: poverty ~ hs_grad + female_house + white + metro_res
## <environment: 0x6787d60>
##
## --= Step 1 =--
## Single term deletions
##
## Model:
## poverty ~ hs_grad + female_house + white + metro_res
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1.0570 -173.16
## hs_grad 1 0.80835 1.8654 -147.89 32.8840 8.887e-07 ***
## female_house 1 0.01830 1.0753 -174.33 0.7444 0.393054
## white 1 0.01216 1.0692 -174.61 0.4946 0.485652
## metro_res 1 0.19273 1.2497 -167.12 7.8403 0.007623 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## --= Step 2 =--
## Single term deletions
##
## Model:
## poverty ~ hs_grad + female_house + metro_res
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1.0692 -174.61
## hs_grad 1 0.91336 1.9825 -146.97 37.5879 2.167e-07 ***
## female_house 1 0.00620 1.0754 -176.33 0.2552 0.615950
## metro_res 1 0.22429 1.2935 -167.47 9.2301 0.003995 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## --= Step 3 =--
## Single term deletions
##
## Model:
## poverty ~ hs_grad + metro_res
## Df Sum of Sq RSS AIC F value Pr(>F)
## <none> 1.0754 -176.33
## hs_grad 1 1.65440 2.7298 -133.62 69.230 1.183e-10 ***
## metro_res 1 0.28776 1.3631 -166.95 12.041 0.001159 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1##
## Call:
## lm(formula = poverty ~ hs_grad + metro_res, data = df2)
##
## Coefficients:
## (Intercept) hs_grad metro_res
## 21.82319 -4.27876 -0.00528Пошаговая регрессия вперед и назад “включила” в модель два фактора, оказывающих значимое влияние на зависимую переменную: hs_grad и metro_res.
summary(fit2)##
## Call:
## lm(formula = poverty ~ hs_grad + metro_res, data = df2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.2480 -0.1249 -0.0210 0.1128 0.3382
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 21.823186 2.288032 9.538 2.24e-12 ***
## hs_grad -4.278756 0.514247 -8.320 1.18e-10 ***
## metro_res -0.005280 0.001522 -3.470 0.00116 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## s: 0.1546 on 45 degrees of freedom
## Multiple R-squared: 0.6495,
## Adjusted R-squared: 0.634
## F-statistic: 41.7 on 2 and 45 DF, p-value: 5.683e-11plot(fit2)
