Task1
df <- read.csv(file = "Angell.csv", header = TRUE, as.is = FALSE) #Скачиваем csv, сохраняем в папку, где лежит r markdown и считываем данные.
head(df) #Head() Показывает первые шесть строк данных.3. Описание данных
Рассматриваемые данные содержат различную информацию о 43 городах в Америке около 1950г.
moral - Moral integration - Composite of crime rate and welfare expenditures
hetero - Ethnic Heterogenity - From percentages of nonwhite and foreign-born white residents
mobility - Geographic Mobility - From percentages of residents moving into and out of the city
region - A factor with levels: E Northeast; MW Midwest; S Southeast; W West.
#library(dplyr)
library(lattice)
panel.hist.splom <- function(x, ...) {
yrng <- current.panel.limits()$ylim
h <- hist(x, plot = FALSE)
breaks <- h$breaks
nB <- length(breaks)
y <- h$counts
y <- yrng[1] + 0.95 * diff(yrng) * y / max(y)
panel.rect(breaks[-nB], yrng[1], breaks[-1], y, col="cyan", ...)
}
splom(~subset(df, select=c(moral:mobility)),
data = df, groups = region,
diag.panel = panel.hist.splom, pch = 21, auto.key = list(space="top", columns=4))
dfl <- transform(df, hetero=log(hetero))
splom(~subset(dfl, select=c(moral:mobility)),
data = dfl, groups = region,
diag.panel = panel.hist.splom, pch = 21, auto.key = list(space="right", columns=1, title = "Region"))
dflo <- dfl #сохраняем датасет
dflo[38,] <- NA #удаляем аутлаер относительно линейной зависимости между hetero и mobility
splom(~dflo[,2:4],
data = dflo, subset = region == 'S')
dfpart <- na.omit(subset(dflo, region == "S"))4.1 Линейная регрессия
Модель: \[y = b_0 + b_1 x_1 + \dots + b_p x_p +\varepsilon\]
summary(lm(mobility ~ hetero, data = dfpart)) #Функция, строящая линейную регрессию. Предсказываем mobility по hetero. В полученных результатах смотрим сначала на самую нижнюю строчку, в которой находится p-value, полученное в результате проверки гипотезы о значимости регрессии. Далее смотрим на таблицу Coefficients. В первом столбце показаны коэффициенты свободного члена (b_0) и коэффициенты перед другими предсказывающими признаками (b_1, b_2, etc). Самый правый столбец в данной таблице есть p-value, полученное в результате проверги гипотезы о том, что коэффициент перед признаком равняется нулю.##
## Call:
## lm(formula = mobility ~ hetero, data = dfpart)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.8946 -3.7644 0.0698 1.0746 8.5405
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 90.669 10.051 9.021 2.05e-06 ***
## hetero -14.660 2.559 -5.729 0.000132 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.111 on 11 degrees of freedom
## Multiple R-squared: 0.749, Adjusted R-squared: 0.7261
## F-statistic: 32.82 on 1 and 11 DF, p-value: 0.0001325summary(mdl <- lm(mobility ~ hetero + moral, data = dfpart)) #Здесь все то же самое, но уже два предсказывающих признака. В данном случае, p-value проверки гипотезы о том, коэффициент перед moral равняется нулю, есть 0.02. Одна из причин, из-за которой результаты регрессии могут быть неправильными - сильно зависимые "независимые" признаки. В нашем случае moral и hetero не коррелируют, поэтому данный вариант мы отбрасываем. Далее, стоит посмотреть на аутлаеры.##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6684 -1.5507 -0.9228 1.3104 5.4488
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 105.3438 9.6566 10.909 7.12e-07 ***
## hetero -16.0823 2.1032 -7.646 1.75e-05 ***
## moral -1.1792 0.4346 -2.713 0.0218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.272 on 10 degrees of freedom
## Multiple R-squared: 0.8554, Adjusted R-squared: 0.8265
## F-statistic: 29.58 on 2 and 10 DF, p-value: 6.318e-05# Стандартизованные коэффициенты регрессии (так называются коэффициенты регрессии, если признаки сначала стандартизовать, а потом строить регрессию). Их правильнее анализировать, так как их значение не зависит от масштаба.
library(lm.beta)## Warning: package 'lm.beta' was built under R version 4.0.3summary(lm.beta(mdl))##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6684 -1.5507 -0.9228 1.3104 5.4488
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 105.3438 0.0000 9.6566 10.909 7.12e-07 ***
## hetero -16.0823 -0.9494 2.1032 -7.646 1.75e-05 ***
## moral -1.1792 -0.3369 0.4346 -2.713 0.0218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.272 on 10 degrees of freedom
## Multiple R-squared: 0.8554, Adjusted R-squared: 0.8265
## F-statistic: 29.58 on 2 and 10 DF, p-value: 6.318e-05## Аутлаеры в регрессионной модели
library(MASS)
r <- stdres(mdl)
jackknife.res <- studres(mdl) #deleted residuals
plot(jackknife.res ~ r)
abline(coef = c(0,1))
d1 <- cooks.distance(mdl)
o <- cbind(dfpart, d1, r, jackknife.res)
osorted <- o[order(-d1),]
head(osorted)o[d1 > 4*mean(d1),] #Один из критериев, по которому можем искать аутлаеры. Смотрим на расстояния Кука у индивидов, которые больше четырех средних расстояний Кука по всем индивидам.Тульса является городом, бизнессмен которого является “отцом” шоссе 66. По этой причине город заслужил прозвище “место рождения шоссе 66”, что оказало очень сильное влияние на становление города.
#Перестраиваем модель с убранным аутлаером.
summary(mdl2 <- lm(mobility ~ hetero + moral, data = dfpart[-c(11,12),])) ##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart[-c(11,
## 12), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.4024 -1.2432 -0.2792 1.2871 5.3981
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 108.744 12.252 8.876 2.05e-05 ***
## hetero -17.096 2.388 -7.158 9.63e-05 ***
## moral -1.128 0.502 -2.248 0.0548 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 2.825 on 8 degrees of freedom
## Multiple R-squared: 0.8734, Adjusted R-squared: 0.8418
## F-statistic: 27.6 on 2 and 8 DF, p-value: 0.0002566summary(lm.beta(lm(mobility ~ hetero, data = dfpart[-c(11,12),])))##
## Call:
## lm(formula = mobility ~ hetero, data = dfpart[-c(11, 12), ])
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.0958 -2.9358 0.7436 1.9711 5.0351
##
## Coefficients:
## Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) 87.8446 0.0000 9.6085 9.142 7.51e-06 ***
## hetero -14.1600 -0.8908 2.4079 -5.881 0.000235 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.402 on 9 degrees of freedom
## Multiple R-squared: 0.7935, Adjusted R-squared: 0.7705
## F-statistic: 34.58 on 1 and 9 DF, p-value: 0.0002347Практически ничего не изменилось (в целом, возможно, это не был outlier)
stepAIC(mdl) #backward## Start: AIC=33.41
## mobility ~ hetero + moral
##
## Df Sum of Sq RSS AIC
## <none> 107.05 33.408
## - moral 1 78.82 185.86 38.581
## - hetero 1 625.87 732.92 56.417##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart)
##
## Coefficients:
## (Intercept) hetero moral
## 105.344 -16.082 -1.179По результатам видно, что от убирания каждой из переменных значение AIC возрастает, т.е. каждая нужна. Но от удаления из модели moral возрастает несильно, а вот hetero точно нельзя убирать.
#Теперь посмотрим, как зависимые переменные всё портят. VIF (variance inflation factor) - 1/(1-Rj^2), где Rj - коэффициент детерминации (множественный коэффициент корреляции) j-го предиктора по остальным предикторам.
summary(mdl)##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart)
##
## Residuals:
## Min 1Q Median 3Q Max
## -4.6684 -1.5507 -0.9228 1.3104 5.4488
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 105.3438 9.6566 10.909 7.12e-07 ***
## hetero -16.0823 2.1032 -7.646 1.75e-05 ***
## moral -1.1792 0.4346 -2.713 0.0218 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.272 on 10 degrees of freedom
## Multiple R-squared: 0.8554, Adjusted R-squared: 0.8265
## F-statistic: 29.58 on 2 and 10 DF, p-value: 6.318e-05library(car)## Warning: package 'car' was built under R version 4.0.3## Loading required package: carDatavif(mdl)## hetero moral
## 1.066209 1.066209dfpart.art <- dfpart
set.seed(100)
dfpart.art$hetero2 <- dfpart$hetero + rnorm(nrow(dfpart), 0, 0.1)
summary(mdl.art <- lm(mobility ~ hetero + hetero2 + moral, data = dfpart.art))##
## Call:
## lm(formula = mobility ~ hetero + hetero2 + moral, data = dfpart.art)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.3471 -1.0757 -0.6504 1.6103 5.2251
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 100.9791 12.8386 7.865 2.53e-05 ***
## hetero -30.7553 27.0887 -1.135 0.286
## hetero2 15.8551 29.1759 0.543 0.600
## moral -1.2082 0.4539 -2.662 0.026 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.394 on 9 degrees of freedom
## Multiple R-squared: 0.86, Adjusted R-squared: 0.8133
## F-statistic: 18.43 on 3 and 9 DF, p-value: 0.0003498vif(mdl.art)## hetero hetero2 moral
## 164.399659 163.570295 1.081132stepAIC(mdl.art)## Start: AIC=34.99
## mobility ~ hetero + hetero2 + moral
##
## Df Sum of Sq RSS AIC
## - hetero2 1 3.401 107.05 33.408
## - hetero 1 14.845 118.49 34.728
## <none> 103.64 34.988
## - moral 1 81.597 185.24 40.537
##
## Step: AIC=33.41
## mobility ~ hetero + moral
##
## Df Sum of Sq RSS AIC
## <none> 107.05 33.408
## - moral 1 78.82 185.86 38.581
## - hetero 1 625.87 732.92 56.417##
## Call:
## lm(formula = mobility ~ hetero + moral, data = dfpart.art)
##
## Coefficients:
## (Intercept) hetero moral
## 105.344 -16.082 -1.179