Файл Angell в пакете car содержит данные о 43 городах США, которые были использованы в одном из первых применений множественной линейной регрессии в социологии. (См. Angell, R. C. “The moral integration of American cities.” American Journal of Sociology, 57 (1951, part 2): 1-140.) Индекс “моральной интеграции” города (т.е. социальной интеграции, противоположности аномии) сочетает информацию об уровне преступности в городе и данные о доле расходов на социальные программы в муниципальном бюджете. Высокие значения этого индекса связаны с низким уровнем преступности и высокой долей социальных расходов. Индекс этнической разнородности высчитывался по сочетанию пропорции небелого населения города и пропорции лиц, родившихся за пределами США. Индекс географической мобильности был построен на основе пропорций лиц, выехавших из города на новое место жительства и прибывших в город из других мест. Angell утверждает, что “моральная интеграция” (moral) должна уменьшаться с ростом этнической разнородности (hetero) и географической мобильности (mobility).
require("car")
А. Используя двумерные графики, рассмотрите зависимость между moral с одной стороны и (по очереди) hetero и mobility с другой стороны. В каждом случае попробуйте определить, будет ли регрессия линейной или нелинейной.
plot(Angell$moral ~ Angell$hetero)
plot(Angell$moral ~ Angell$mobility)
Судя по графикам, в первом случае регрессия, вероятно, будет линейной. Во втором - нелинейной.
Б. Попробуйте преобразовать переменные hetero и mobility. Сравните R-квадраты в моделях с трансформированными переменными и в моделях с нетрансформированными переменными. Сравните SSE и MSE в этих же моделях. Приводит ли преобразование к улучшению моделей?
summary(lm(Angell$moral ~ Angell$hetero))
##
## Call:
## lm(formula = Angell$moral ~ Angell$hetero)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.678 -2.610 0.249 2.297 6.693
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 14.4236 0.8251 17.48 < 2e-16 ***
## Angell$hetero -0.1028 0.0221 -4.64 3.5e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.93 on 41 degrees of freedom
## Multiple R-squared: 0.345, Adjusted R-squared: 0.329
## F-statistic: 21.6 on 1 and 41 DF, p-value: 3.49e-05
anova(lm(Angell$moral ~ Angell$hetero))
## Analysis of Variance Table
##
## Response: Angell$moral
## Df Sum Sq Mean Sq F value Pr(>F)
## Angell$hetero 1 185 184.7 21.6 3.5e-05 ***
## Residuals 41 351 8.6
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(Angell$moral ~ Angell$mobility))
##
## Call:
## lm(formula = Angell$moral ~ Angell$mobility)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.104 -2.253 0.425 2.355 5.543
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 16.1437 1.4518 11.12 6e-14 ***
## Angell$mobility -0.1791 0.0496 -3.61 0.00083 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.15 on 41 degrees of freedom
## Multiple R-squared: 0.241, Adjusted R-squared: 0.223
## F-statistic: 13 on 1 and 41 DF, p-value: 0.00083
anova(lm(Angell$moral ~ Angell$mobility))
## Analysis of Variance Table
##
## Response: Angell$moral
## Df Sum Sq Mean Sq F value Pr(>F)
## Angell$mobility 1 129 129.1 13 0.00083 ***
## Residuals 41 407 9.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
В первой модели R-квадрат порядка 34%, во второй - 24%, т.е обе модели плохо описывают зависимость переменной moral от переменных hetero и mobility.
Попробуем улучшить эти модели.
summary(lm(Angell$moral ~ sqrt(Angell$hetero)))
##
## Call:
## lm(formula = Angell$moral ~ sqrt(Angell$hetero))
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.496 -2.317 0.003 2.366 6.775
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 17.946 1.508 11.90 7.0e-15 ***
## sqrt(Angell$hetero) -1.261 0.269 -4.68 3.1e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.92 on 41 degrees of freedom
## Multiple R-squared: 0.348, Adjusted R-squared: 0.333
## F-statistic: 21.9 on 1 and 41 DF, p-value: 3.09e-05
anova(lm(Angell$moral ~ sqrt(Angell$hetero)))
## Analysis of Variance Table
##
## Response: Angell$moral
## Df Sum Sq Mean Sq F value Pr(>F)
## sqrt(Angell$hetero) 1 187 186.7 21.9 3.1e-05 ***
## Residuals 41 349 8.5
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
summary(lm(Angell$moral ~ log(Angell$mobility)))
##
## Call:
## lm(formula = Angell$moral ~ log(Angell$mobility))
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.783 -2.495 0.551 1.912 5.044
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.40 4.19 6.78 3.4e-08 ***
## log(Angell$mobility) -5.29 1.28 -4.13 0.00017 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.04 on 41 degrees of freedom
## Multiple R-squared: 0.294, Adjusted R-squared: 0.277
## F-statistic: 17.1 on 1 and 41 DF, p-value: 0.000173
anova(lm(Angell$moral ~ log(Angell$mobility)))
## Analysis of Variance Table
##
## Response: Angell$moral
## Df Sum Sq Mean Sq F value Pr(>F)
## log(Angell$mobility) 1 158 157.5 17.1 0.00017 ***
## Residuals 41 378 9.2
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Увеличившиеся значения R-квадрат говорят о том, что новые переменные (sqrt(Angell$hetero) и log(Angell$mobility)) лучше объясняют moral, однако, в первом случае это увеличение несущественное, во втором случае мы наблюдаем прирост значения R-квадрат всего в 5%. Суммы квадратов остатков в новых моделях уменьшились, но незначительно. В целом, нельзя сказать, что трансформация переменных привела к значимому улучшению моделей.
В. Проверьте зависимость переменной moral от переменной region (регион США: e - Восток, mw - Средний Запад, s - Юг, w - Запад).
plot(Angell$moral ~ Angell$region)
summary(aov(Angell$moral ~ factor(Angell$region)))
## Df Sum Sq Mean Sq F value Pr(>F)
## factor(Angell$region) 3 343 114.5 23.2 8.8e-09 ***
## Residuals 39 192 4.9
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
P-значение достаточно мало (8.8e-09), для того, чтобы отвергнуть нулевую ипотезу, говорющую о независимости переменной moral от региона.
Посмотрим, как зависит переменная moral от hetero и mobiliy по регионам.
library(ggplot2)
sp1 <- ggplot(Angell, aes(x = hetero, y = moral)) + geom_point(shape = 1)
sp2 <- ggplot(Angell, aes(x = mobility, y = moral)) + geom_point(shape = 1)
sp1 + facet_grid(region ~ .)
sp2 + facet_grid(region ~ .)
Г.Проверьте зависимость moral от всех переменных, учитывая возможность мультипликативных эффектов. Подберите модель, которая представляется Вам наилучшей. Объясните, почему Вы ее выбрали. Дайте общую интепретацию модели.
Составим несколько моделей
model.31 <- lm(Angell$moral ~ Angell$mobility * Angell$region)
model.32 <- lm(Angell$moral ~ Angell$hetero * Angell$region)
model.33 <- lm(Angell$moral ~ (Angell$hetero + Angell$mobility) * Angell$region)
Посмотрим, какая модель будет лучшей.
summary(model.31)
##
## Call:
## lm(formula = Angell$moral ~ Angell$mobility * Angell$region)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.267 -1.799 0.384 1.436 3.785
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 12.6932 4.9663 2.56 0.015 *
## Angell$mobility 0.1681 0.3085 0.54 0.589
## Angell$regionMW -0.6888 5.5327 -0.12 0.902
## Angell$regionS -4.2297 5.5839 -0.76 0.454
## Angell$regionW -5.8513 7.7746 -0.75 0.457
## Angell$mobility:Angell$regionMW -0.1576 0.3215 -0.49 0.627
## Angell$mobility:Angell$regionS -0.1906 0.3178 -0.60 0.553
## Angell$mobility:Angell$regionW -0.0699 0.3466 -0.20 0.841
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.32 on 35 degrees of freedom
## Multiple R-squared: 0.649, Adjusted R-squared: 0.578
## F-statistic: 9.23 on 7 and 35 DF, p-value: 2e-06
summary(model.32)
##
## Call:
## lm(formula = Angell$moral ~ Angell$hetero * Angell$region)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.79 -1.20 0.21 1.20 3.28
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 18.2562 1.8041 10.12 6.2e-12 ***
## Angell$hetero -0.1230 0.0705 -1.74 0.08976 .
## Angell$regionMW -4.1314 2.3676 -1.74 0.08977 .
## Angell$regionS -9.0486 2.3899 -3.79 0.00058 ***
## Angell$regionW -4.0026 4.3131 -0.93 0.35976
## Angell$hetero:Angell$regionMW 0.0379 0.0963 0.39 0.69660
## Angell$hetero:Angell$regionS 0.0950 0.0758 1.25 0.21839
## Angell$hetero:Angell$regionW -0.1028 0.2412 -0.43 0.67270
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.15 on 35 degrees of freedom
## Multiple R-squared: 0.699, Adjusted R-squared: 0.638
## F-statistic: 11.6 on 7 and 35 DF, p-value: 1.62e-07
summary(model.33)
##
## Call:
## lm(formula = Angell$moral ~ (Angell$hetero + Angell$mobility) *
## Angell$region)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.596 -1.414 0.157 1.365 3.438
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 22.4982 7.3894 3.04 0.0047 **
## Angell$hetero -0.1559 0.0908 -1.72 0.0960 .
## Angell$mobility -0.2182 0.3681 -0.59 0.5577
## Angell$regionMW -4.4744 8.6334 -0.52 0.6080
## Angell$regionS -8.7397 8.5679 -1.02 0.3156
## Angell$regionW -11.1468 10.8119 -1.03 0.3105
## Angell$hetero:Angell$regionMW 0.0172 0.1266 0.14 0.8929
## Angell$hetero:Angell$regionS 0.1040 0.0975 1.07 0.2943
## Angell$hetero:Angell$regionW -0.0432 0.2598 -0.17 0.8690
## Angell$mobility:Angell$regionMW 0.1131 0.3850 0.29 0.7708
## Angell$mobility:Angell$regionS 0.1166 0.3790 0.31 0.7603
## Angell$mobility:Angell$regionW 0.2840 0.3992 0.71 0.4821
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.19 on 31 degrees of freedom
## Multiple R-squared: 0.723, Adjusted R-squared: 0.624
## F-statistic: 7.34 on 11 and 31 DF, p-value: 5.41e-06
Третья модель лучше всего объясняет переменную moral, однако, ни один из коэффициентов не является значимым. Попробуем исключить переменную region из модели.
model.34 <- lm(Angell$moral ~ Angell$hetero + Angell$mobility)
summary(model.34)
##
## Call:
## lm(formula = Angell$moral ~ Angell$hetero + Angell$mobility)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.071 -1.194 -0.206 1.738 4.195
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 19.9408 1.1927 16.72 < 2e-16 ***
## Angell$hetero -0.1086 0.0170 -6.39 1.3e-07 ***
## Angell$mobility -0.1933 0.0354 -5.46 2.7e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.24 on 40 degrees of freedom
## Multiple R-squared: 0.624, Adjusted R-squared: 0.606
## F-statistic: 33.2 on 2 and 40 DF, p-value: 3.13e-09
Эта модель объясняет переменную moral немного хуже предыдущей, однако, все переменные в ней оказываются значимыми.