Regression Analysis for Boston
Nao Mimoto – Class Web Page
1. Preliminary
The dataset “Boston” comes with package MASS. When you load the package, it automatically loads the dataset as well.
library(MASS) # install.packages("MASS")
Boston # see the data
head(Boston) # another way to see the data
?Boston # see explanation for variables
# 1 crim per capita crime rate by town.
# 2 zn proportion of residential land zoned for lots over 25,000 sq.ft.
# 3 indus proportion of non-retail business acres per town.
# 4 chas Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
# 5 nox nitrogen oxides concentration (parts per 10 million).
# 6 rm average number of rooms per dwelling.
# 7 age proportion of owner-occupied units built prior to 1940.
# 8 dis weighted mean of distances to five Boston employment centres.
# 9 rad index of accessibility to radial highways.
# 10 tax full-value property-tax rate per $10,000.
# 11 ptratio pupil-teacher ratio by town.
# 12 black 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.
# 13 lstat lower status of the population (percent).
# 14 medv median value of owner-occupied homes in \$1000s.
# 13 lstat = lower status of the population (percent)
# Proportion of population that is lower status
# = 1/2 (proportion of adults without,
# some high school education and proportion of male
# workers classified as laborers). The logarithmic
# specification implies that socioeconomic status
# distinctions mean more in the upper brackets of
# society than in the lower classes. Source: 1970 U. S. Census2. Turn Boston into a tibble
Boston comes as data.frame. It is a good habit to turn all the data you work with into a tibble for consistency and conveinience. A tibble is a modern version of data.frame that can be worked on using pakages in tidyverse.
class(Boston) # show the class of Boston
# - it's data.frame not tibble
# it's better to turn it into tibble
library(tidyverse) # install.packages("tidyverse")
Boston <- as_tibble(Boston)
Boston
class(Boston) # now it should be a tibble
is_tibble(Boston) # is Boston tibble now?
is.data.frame(Boston) # it's still data.frame too
Boston # only first 10 rows will print by default
# print(Boston, n=100) # if you want to see more rows
# view(Boston) # spreadsheet-like view NOT RECOMMENDED
# CLOSE TO MOVE FORWARD3. Set Response variable
I like to name response variable ‘resp’ so that I don’t have to adjust the code every time. It also makes clear what is Y (resp), and what is X (everything else).
library(tidyverse)
### set medv as the response variable
# Rename "medv" column as "resp" to streamline analysis.
Boston2 <- Boston %>% rename(resp=medv)
Boston2## # A tibble: 506 × 14
## crim zn indus chas nox rm age dis rad tax ptratio black
## <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl>
## 1 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3 397.
## 2 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8 397.
## 3 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8 393.
## 4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7 395.
## 5 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7 397.
## 6 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7 394.
## 7 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2 396.
## 8 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2 397.
## 9 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2 387.
## 10 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2 387.
## # ℹ 496 more rows
## # ℹ 2 more variables: lstat <dbl>, resp <dbl>## # A tibble: 506 × 14
## resp crim zn indus chas nox rm age dis rad tax ptratio
## <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 24 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3
## 2 21.6 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8
## 3 34.7 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8
## 4 33.4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7
## 5 36.2 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7
## 6 28.7 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7
## 7 22.9 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2
## 8 27.1 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2
## 9 16.5 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2
## 10 18.9 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2
## # ℹ 496 more rows
## # ℹ 2 more variables: black <dbl>, lstat <dbl># turn "chas" column into 0/1 factor (not used here)
Boston3 <- Boston2 %>% mutate( chas=as.factor(chas) )
Boston3## # A tibble: 506 × 14
## resp crim zn indus chas nox rm age dis rad tax ptratio
## <dbl> <dbl> <dbl> <dbl> <fct> <dbl> <dbl> <dbl> <dbl> <int> <dbl> <dbl>
## 1 24 0.00632 18 2.31 0 0.538 6.58 65.2 4.09 1 296 15.3
## 2 21.6 0.0273 0 7.07 0 0.469 6.42 78.9 4.97 2 242 17.8
## 3 34.7 0.0273 0 7.07 0 0.469 7.18 61.1 4.97 2 242 17.8
## 4 33.4 0.0324 0 2.18 0 0.458 7.00 45.8 6.06 3 222 18.7
## 5 36.2 0.0690 0 2.18 0 0.458 7.15 54.2 6.06 3 222 18.7
## 6 28.7 0.0298 0 2.18 0 0.458 6.43 58.7 6.06 3 222 18.7
## 7 22.9 0.0883 12.5 7.87 0 0.524 6.01 66.6 5.56 5 311 15.2
## 8 27.1 0.145 12.5 7.87 0 0.524 6.17 96.1 5.95 5 311 15.2
## 9 16.5 0.211 12.5 7.87 0 0.524 5.63 100 6.08 5 311 15.2
## 10 18.9 0.170 12.5 7.87 0 0.524 6.00 85.9 6.59 5 311 15.2
## # ℹ 496 more rows
## # ℹ 2 more variables: black <dbl>, lstat <dbl>4. Full Model: Regression with all variables
##
## Call:
## lm(formula = resp ~ ., data = Boston2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.595 -2.730 -0.518 1.777 26.199
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.646e+01 5.103e+00 7.144 3.28e-12 ***
## crim -1.080e-01 3.286e-02 -3.287 0.001087 **
## zn 4.642e-02 1.373e-02 3.382 0.000778 ***
## indus 2.056e-02 6.150e-02 0.334 0.738288
## chas 2.687e+00 8.616e-01 3.118 0.001925 **
## nox -1.777e+01 3.820e+00 -4.651 4.25e-06 ***
## rm 3.810e+00 4.179e-01 9.116 < 2e-16 ***
## age 6.922e-04 1.321e-02 0.052 0.958229
## dis -1.476e+00 1.995e-01 -7.398 6.01e-13 ***
## rad 3.060e-01 6.635e-02 4.613 5.07e-06 ***
## tax -1.233e-02 3.760e-03 -3.280 0.001112 **
## ptratio -9.527e-01 1.308e-01 -7.283 1.31e-12 ***
## black 9.312e-03 2.686e-03 3.467 0.000573 ***
## lstat -5.248e-01 5.072e-02 -10.347 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7338
## F-statistic: 108.1 on 13 and 492 DF, p-value: < 2.2e-16
# How to looking at all the things in the output
str(Reg01) # or you can click on environment tab in Rstudio5. Removing variables by hand
From the Full Model, remove age
Remove indus
##
## Call:
## lm(formula = resp ~ . - age - indus, data = Boston2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5984 -2.7386 -0.5046 1.7273 26.2373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.341145 5.067492 7.171 2.73e-12 ***
## crim -0.108413 0.032779 -3.307 0.001010 **
## zn 0.045845 0.013523 3.390 0.000754 ***
## chas 2.718716 0.854240 3.183 0.001551 **
## nox -17.376023 3.535243 -4.915 1.21e-06 ***
## rm 3.801579 0.406316 9.356 < 2e-16 ***
## dis -1.492711 0.185731 -8.037 6.84e-15 ***
## rad 0.299608 0.063402 4.726 3.00e-06 ***
## tax -0.011778 0.003372 -3.493 0.000521 ***
## ptratio -0.946525 0.129066 -7.334 9.24e-13 ***
## black 0.009291 0.002674 3.475 0.000557 ***
## lstat -0.522553 0.047424 -11.019 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7348
## F-statistic: 128.2 on 11 and 494 DF, p-value: < 2.2e-166. Adding variables by hand
# Using variables that had low p-value in Chi-square test
Reg04 <- lm(resp ~ zn + indus + nox + rad, data=Boston2)
summary(Reg04)##
## Call:
## lm(formula = resp ~ zn + indus + nox + rad, data = Boston2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.993 -4.736 -1.601 3.028 33.309
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 29.95772 2.35422 12.725 < 2e-16 ***
## zn 0.05360 0.01830 2.929 0.00355 **
## indus -0.39030 0.08522 -4.580 5.88e-06 ***
## nox -4.22098 5.06033 -0.834 0.40460
## rad -0.14095 0.05302 -2.659 0.00810 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.924 on 501 degrees of freedom
## Multiple R-squared: 0.2635, Adjusted R-squared: 0.2576
## F-statistic: 44.81 on 4 and 501 DF, p-value: < 2.2e-16Remove nox from the model
# Using variables that had low p-value in Chi-square test
Reg05 <- lm(resp ~ zn + indus + rad, data=Boston2)
summary(Reg05)##
## Call:
## lm(formula = resp ~ zn + indus + rad, data = Boston2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -13.046 -4.723 -1.503 2.970 33.332
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 28.13613 0.87912 32.005 < 2e-16 ***
## zn 0.05684 0.01788 3.180 0.00157 **
## indus -0.42850 0.07185 -5.964 4.65e-09 ***
## rad -0.15468 0.05038 -3.070 0.00226 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 7.922 on 502 degrees of freedom
## Multiple R-squared: 0.2625, Adjusted R-squared: 0.258
## F-statistic: 59.54 on 3 and 502 DF, p-value: < 2.2e-16
7. Exaustive Search
Go through all possible combination of variables, and find the best model using AIC, BIC, or Adjusted R-squared.
library(leaps) # install.packages('leaps', repos='https://cran.case.edu/')
## nvmax sets max number of variables in a subset
allFits = regsubsets(resp ~ ., data=Boston2,
nvmax=15,
method = "exhaustive")
modelList = summary(allFits)
modelList## Subset selection object
## Call: regsubsets.formula(resp ~ ., data = Boston2, nvmax = 15, method = "exhaustive")
## 13 Variables (and intercept)
## Forced in Forced out
## crim FALSE FALSE
## zn FALSE FALSE
## indus FALSE FALSE
## chas FALSE FALSE
## nox FALSE FALSE
## rm FALSE FALSE
## age FALSE FALSE
## dis FALSE FALSE
## rad FALSE FALSE
## tax FALSE FALSE
## ptratio FALSE FALSE
## black FALSE FALSE
## lstat FALSE FALSE
## 1 subsets of each size up to 13
## Selection Algorithm: exhaustive
## crim zn indus chas nox rm age dis rad tax ptratio black lstat
## 1 ( 1 ) " " " " " " " " " " " " " " " " " " " " " " " " "*"
## 2 ( 1 ) " " " " " " " " " " "*" " " " " " " " " " " " " "*"
## 3 ( 1 ) " " " " " " " " " " "*" " " " " " " " " "*" " " "*"
## 4 ( 1 ) " " " " " " " " " " "*" " " "*" " " " " "*" " " "*"
## 5 ( 1 ) " " " " " " " " "*" "*" " " "*" " " " " "*" " " "*"
## 6 ( 1 ) " " " " " " "*" "*" "*" " " "*" " " " " "*" " " "*"
## 7 ( 1 ) " " " " " " "*" "*" "*" " " "*" " " " " "*" "*" "*"
## 8 ( 1 ) " " "*" " " "*" "*" "*" " " "*" " " " " "*" "*" "*"
## 9 ( 1 ) "*" " " " " "*" "*" "*" " " "*" "*" " " "*" "*" "*"
## 10 ( 1 ) "*" "*" " " " " "*" "*" " " "*" "*" "*" "*" "*" "*"
## 11 ( 1 ) "*" "*" " " "*" "*" "*" " " "*" "*" "*" "*" "*" "*"
## 12 ( 1 ) "*" "*" "*" "*" "*" "*" " " "*" "*" "*" "*" "*" "*"
## 13 ( 1 ) "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*" "*"Best by min AIC
## [1] 11## (Intercept) crim zn indus chas nox
## TRUE TRUE TRUE FALSE TRUE TRUE
## rm age dis rad tax ptratio
## TRUE FALSE TRUE TRUE TRUE TRUE
## black lstat
## TRUE TRUE ### Run the best model using lm()
ixIsItIn = modelList$which[ixA, -1] # -1 is to remove intercept
ixIsItIn## crim zn indus chas nox rm age dis rad tax
## TRUE TRUE FALSE TRUE TRUE TRUE FALSE TRUE TRUE TRUE
## ptratio black lstat
## TRUE TRUE TRUE## [1] "crim" "zn" "chas" "nox" "rm" "dis" "rad"
## [8] "tax" "ptratio" "black" "lstat"## [1] "crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat"## [1] "resp ~ crim + zn + chas + nox + rm + dis + rad + tax + ptratio + black + lstat"##
## Call:
## lm(formula = as.formula(formula_string), data = Boston2)
##
## Coefficients:
## (Intercept) crim zn chas nox rm
## 36.341145 -0.108413 0.045845 2.718716 -17.376023 3.801579
## dis rad tax ptratio black lstat
## -1.492711 0.299608 -0.011778 -0.946525 0.009291 -0.522553## [1] 3023.726## [1] 3078.671## [1] 0.7348058Best by max Adj R-sq
## [1] 11## (Intercept) crim zn indus chas nox
## TRUE TRUE TRUE FALSE TRUE TRUE
## rm age dis rad tax ptratio
## TRUE FALSE TRUE TRUE TRUE TRUE
## black lstat
## TRUE TRUE ### Run the best model using lm()
ixIsItIn = modelList$which[ixB, -1] # -1 is to remove intercept
listOfVars = paste(names(ixIsItIn)[ixIsItIn], collapse=" + ")
formula_string = paste("resp ~", listOfVars)
Fit11 = lm(as.formula(formula_string), data=Boston2)
summary(Fit11)##
## Call:
## lm(formula = as.formula(formula_string), data = Boston2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.5984 -2.7386 -0.5046 1.7273 26.2373
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 36.341145 5.067492 7.171 2.73e-12 ***
## crim -0.108413 0.032779 -3.307 0.001010 **
## zn 0.045845 0.013523 3.390 0.000754 ***
## chas 2.718716 0.854240 3.183 0.001551 **
## nox -17.376023 3.535243 -4.915 1.21e-06 ***
## rm 3.801579 0.406316 9.356 < 2e-16 ***
## dis -1.492711 0.185731 -8.037 6.84e-15 ***
## rad 0.299608 0.063402 4.726 3.00e-06 ***
## tax -0.011778 0.003372 -3.493 0.000521 ***
## ptratio -0.946525 0.129066 -7.334 9.24e-13 ***
## black 0.009291 0.002674 3.475 0.000557 ***
## lstat -0.522553 0.047424 -11.019 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared: 0.7406, Adjusted R-squared: 0.7348
## F-statistic: 128.2 on 11 and 494 DF, p-value: < 2.2e-16## [1] 3023.726## [1] 3078.671## [1] 0.7348058COCE ONLY
##--------------------------------------
## 1. Preliminary
library(MASS) # install.packages("MASS")
Boston # see the data
head(Boston) # another way to see the data
?Boston # see explanation for variables
# 1 crim per capita crime rate by town.
# 2 zn proportion of residential land zoned for lots over 25,000 sq.ft.
# 3 indus proportion of non-retail business acres per town.
# 4 chas Charles River dummy variable (= 1 if tract bounds river; 0 otherwise).
# 5 nox nitrogen oxides concentration (parts per 10 million).
# 6 rm average number of rooms per dwelling.
# 7 age proportion of owner-occupied units built prior to 1940.
# 8 dis weighted mean of distances to five Boston employment centres.
# 9 rad index of accessibility to radial highways.
# 10 tax full-value property-tax rate per $10,000.
# 11 ptratio pupil-teacher ratio by town.
# 12 black 1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town.
# 13 lstat lower status of the population (percent).
# 14 medv median value of owner-occupied homes in \$1000s.
# 13 lstat = lower status of the population (percent)
# Proportion of population that is lower status
# = 1/2 (proportion of adults without,
# some high school education and proportion of male
# workers classified as laborers). The logarithmic
# specification implies that socioeconomic status
# distinctions mean more in the upper brackets of
# society than in the lower classes. Source: 1970 U. S. Census
##--------------------------------------
## 2. Turn Boston into a tibble
class(Boston) # show the class of Boston
# it's data.frame not tibble
# it's better to turn it into tibble
library(tidyverse) # install.packages("tidyverse")
Boston <- as_tibble(Boston)
Boston
class(Boston) # now it should be a tibble
is_tibble(Boston) # is Boston tibble now?
is.data.frame(Boston) # it's still data.frame too
Boston # only first 10 rows will print by default
# print(Boston, n=100) # if you want to see more rows
# view(Boston) # spreadsheet-like view NOT RECOMMENDED
# CLOSE TO MOVE FORWARD
##--------------------------------------
## 3. Set Response variable
library(tidyverse)
### set medv as the response variable
# Rename "medv" column as "resp" to streamline analysis.
Boston2 <- Boston %>% rename(resp=medv)
Boston2
# move "resp" columnm to the 1st column
Boston2 <- Boston2 %>% relocate(resp)
Boston2
# turn "chas" column into 0/1 factor (not used here)
Boston3 <- Boston2 %>% mutate( chas=as.factor(chas) )
Boston3
##--------------------------------------
## 4. Full Model: Regression with all variables
Reg01 <- lm(resp ~ . , data=Boston2)
summary(Reg01)
plot(Reg01)
# How to looking at all the things in the output
str(Reg01) # or you can click on environment tab in Rstudio
##--------------------------------------
## 5. Removing variables by hand
# From the full model, remove age
Reg02 <- lm(resp ~ . -age, data=Boston2)
summary(Reg02)
### Remove indus
Reg03 <- lm(resp ~ . -age -indus, data=Boston2)
summary(Reg03)
##--------------------------------------
## 6. Adding variables by hand
# Using variables that had low p-value in Chi-square test
Reg04 <- lm(resp ~ zn + indus + nox + rad, data=Boston2)
summary(Reg04)
### Remove nox from the model
# Using variables that had low p-value in Chi-square test
Reg05 <- lm(resp ~ zn + indus + rad, data=Boston2)
summary(Reg05)
plot(Reg05)
##--------------------------------------
## 7. Exaustive Search
library(leaps) # install.packages('leaps', repos='https://cran.case.edu/')
## nvmax sets max number of variables in a subset
allFits = regsubsets(resp ~ ., data=Boston2,
nvmax=15,
method = "exhaustive")
modelList = summary(allFits)
modelList
### Best by min AIC
# (AIC is directly proportional to Mallow’s Cp)
ixA = which.min(modelList$cp)
print(ixA)
modelList$which[ixA, ]
### Run the best model using lm()
ixIsItIn = modelList$which[ixA, -1] # -1 is to remove intercept
ixIsItIn
names(ixIsItIn)[ixIsItIn]
listOfVars = paste(names(ixIsItIn)[ixIsItIn], collapse=" + ")
listOfVars
formula_string = paste("resp ~", listOfVars)
formula_string
Fit10 = lm(as.formula(formula_string), data=Boston2)
Fit10
AIC(Fit10)
BIC(Fit10)
summary(Fit10)$adj.r.sq
### Best by max Adj R-sq
ixB = which.max(modelList$adjr2)
print(ixB)
modelList$which[ixB, ]
### Run the best model using lm()
ixIsItIn = modelList$which[ixB, -1] # -1 is to remove intercept
listOfVars = paste(names(ixIsItIn)[ixIsItIn], collapse=" + ")
formula_string = paste("resp ~", listOfVars)
Fit11 = lm(as.formula(formula_string), data=Boston2)
summary(Fit11)
AIC(Fit11)
BIC(Fit11)
summary(Fit11)$adj.r.sq
### Best by min BIC
ixC = which.min(modelList$bic)
print(ixC)
modelList$which[ixC, ]