Routine Prelim Analysis for Boston

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. Census

2. 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 FORWARD

3. 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

# 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. Routine Exploratory Analysis for Numeric Variables

Orig <- Boston2  # Set your tibble here

  ### Check for N/A. Remove if there's any.
  summary(Orig)

##       resp            crim                zn             indus      
##  Min.   : 5.00   Min.   : 0.00632   Min.   :  0.00   Min.   : 0.46  
##  1st Qu.:17.02   1st Qu.: 0.08205   1st Qu.:  0.00   1st Qu.: 5.19  
##  Median :21.20   Median : 0.25651   Median :  0.00   Median : 9.69  
##  Mean   :22.53   Mean   : 3.61352   Mean   : 11.36   Mean   :11.14  
##  3rd Qu.:25.00   3rd Qu.: 3.67708   3rd Qu.: 12.50   3rd Qu.:18.10  
##  Max.   :50.00   Max.   :88.97620   Max.   :100.00   Max.   :27.74  
##       chas              nox               rm             age        
##  Min.   :0.00000   Min.   :0.3850   Min.   :3.561   Min.   :  2.90  
##  1st Qu.:0.00000   1st Qu.:0.4490   1st Qu.:5.886   1st Qu.: 45.02  
##  Median :0.00000   Median :0.5380   Median :6.208   Median : 77.50  
##  Mean   :0.06917   Mean   :0.5547   Mean   :6.285   Mean   : 68.57  
##  3rd Qu.:0.00000   3rd Qu.:0.6240   3rd Qu.:6.623   3rd Qu.: 94.08  
##  Max.   :1.00000   Max.   :0.8710   Max.   :8.780   Max.   :100.00  
##       dis              rad              tax           ptratio     
##  Min.   : 1.130   Min.   : 1.000   Min.   :187.0   Min.   :12.60  
##  1st Qu.: 2.100   1st Qu.: 4.000   1st Qu.:279.0   1st Qu.:17.40  
##  Median : 3.207   Median : 5.000   Median :330.0   Median :19.05  
##  Mean   : 3.795   Mean   : 9.549   Mean   :408.2   Mean   :18.46  
##  3rd Qu.: 5.188   3rd Qu.:24.000   3rd Qu.:666.0   3rd Qu.:20.20  
##  Max.   :12.127   Max.   :24.000   Max.   :711.0   Max.   :22.00  
##      black            lstat      
##  Min.   :  0.32   Min.   : 1.73  
##  1st Qu.:375.38   1st Qu.: 6.95  
##  Median :391.44   Median :11.36  
##  Mean   :356.67   Mean   :12.65  
##  3rd Qu.:396.23   3rd Qu.:16.95  
##  Max.   :396.90   Max.   :37.97

  sum(is.na(Orig))

## [1] 0

  # If there is na in the data, run below
  dim(Orig)

## [1] 506  14

  Orig <- Orig %>% na.omit()
  dim(Orig)

## [1] 506  14

  ###  Correlation Check
  # Pick columns that are numeric
  library(corrplot)    # install.packages("corrplot")

## corrplot 0.95 loaded

    Orig_num <- Orig %>% select_if(is.numeric)
    cor(Orig_num)

##               resp        crim          zn       indus         chas         nox
## resp     1.0000000 -0.38830461  0.36044534 -0.48372516  0.175260177 -0.42732077
## crim    -0.3883046  1.00000000 -0.20046922  0.40658341 -0.055891582  0.42097171
## zn       0.3604453 -0.20046922  1.00000000 -0.53382819 -0.042696719 -0.51660371
## indus   -0.4837252  0.40658341 -0.53382819  1.00000000  0.062938027  0.76365145
## chas     0.1752602 -0.05589158 -0.04269672  0.06293803  1.000000000  0.09120281
## nox     -0.4273208  0.42097171 -0.51660371  0.76365145  0.091202807  1.00000000
## rm       0.6953599 -0.21924670  0.31199059 -0.39167585  0.091251225 -0.30218819
## age     -0.3769546  0.35273425 -0.56953734  0.64477851  0.086517774  0.73147010
## dis      0.2499287 -0.37967009  0.66440822 -0.70802699 -0.099175780 -0.76923011
## rad     -0.3816262  0.62550515 -0.31194783  0.59512927 -0.007368241  0.61144056
## tax     -0.4685359  0.58276431 -0.31456332  0.72076018 -0.035586518  0.66802320
## ptratio -0.5077867  0.28994558 -0.39167855  0.38324756 -0.121515174  0.18893268
## black    0.3334608 -0.38506394  0.17552032 -0.35697654  0.048788485 -0.38005064
## lstat   -0.7376627  0.45562148 -0.41299457  0.60379972 -0.053929298  0.59087892
##                  rm         age         dis          rad         tax    ptratio
## resp     0.69535995 -0.37695457  0.24992873 -0.381626231 -0.46853593 -0.5077867
## crim    -0.21924670  0.35273425 -0.37967009  0.625505145  0.58276431  0.2899456
## zn       0.31199059 -0.56953734  0.66440822 -0.311947826 -0.31456332 -0.3916785
## indus   -0.39167585  0.64477851 -0.70802699  0.595129275  0.72076018  0.3832476
## chas     0.09125123  0.08651777 -0.09917578 -0.007368241 -0.03558652 -0.1215152
## nox     -0.30218819  0.73147010 -0.76923011  0.611440563  0.66802320  0.1889327
## rm       1.00000000 -0.24026493  0.20524621 -0.209846668 -0.29204783 -0.3555015
## age     -0.24026493  1.00000000 -0.74788054  0.456022452  0.50645559  0.2615150
## dis      0.20524621 -0.74788054  1.00000000 -0.494587930 -0.53443158 -0.2324705
## rad     -0.20984667  0.45602245 -0.49458793  1.000000000  0.91022819  0.4647412
## tax     -0.29204783  0.50645559 -0.53443158  0.910228189  1.00000000  0.4608530
## ptratio -0.35550149  0.26151501 -0.23247054  0.464741179  0.46085304  1.0000000
## black    0.12806864 -0.27353398  0.29151167 -0.444412816 -0.44180801 -0.1773833
## lstat   -0.61380827  0.60233853 -0.49699583  0.488676335  0.54399341  0.3740443
##               black      lstat
## resp     0.33346082 -0.7376627
## crim    -0.38506394  0.4556215
## zn       0.17552032 -0.4129946
## indus   -0.35697654  0.6037997
## chas     0.04878848 -0.0539293
## nox     -0.38005064  0.5908789
## rm       0.12806864 -0.6138083
## age     -0.27353398  0.6023385
## dis      0.29151167 -0.4969958
## rad     -0.44441282  0.4886763
## tax     -0.44180801  0.5439934
## ptratio -0.17738330  0.3740443
## black    1.00000000 -0.3660869
## lstat   -0.36608690  1.0000000

    corrplot::corrplot(cor(Orig_num))

    corrplot::corrplot(cor(Orig_num), method="number")

  ### Visualization
  # pairs(Orig) equivalent
  library(GGally)      # install.packages("GGally")
    GGally::ggpairs(Orig[, c(2:7, 1)], aes(alpha=1))

    GGally::ggpairs(Orig[, c(8:14, 1)], aes(alpha=1))

  ### Chi-sq test of association for each column
  # Pick columns that has more than 1 unique value
  list_cols <- Orig %>% summarise_all(function(x) length(unique(x))) %>%
        gather() %>% filter(value>1) %>% pull(key)

  ### Apply chisq.test to those columns
  ChiSq.pval <- Orig %>% select(list_cols) %>%
        summarise_all(funs(chisq.test(., Orig$resp)$p.value))

## Warning: Using an external vector in selections was deprecated in tidyselect 1.1.0.
## ℹ Please use `all_of()` or `any_of()` instead.
##   # Was:
##   data %>% select(list_cols)
## 
##   # Now:
##   data %>% select(all_of(list_cols))
## 
## See <https://tidyselect.r-lib.org/reference/faq-external-vector.html>.
## This warning is displayed once per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

## Warning: There were 14 warnings in `summarise()`.
## The first warning was:
## ℹ In argument: `resp = chisq.test(resp, Orig$resp)$p.value`.
## Caused by warning in `chisq.test()`:
## ! Chi-squared approximation may be incorrect
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 13 remaining warnings.

  ChiSq.pval <- Orig %>% summarise_all(funs(chisq.test(., Orig$resp)$p.value))

## Warning: `funs()` was deprecated in dplyr 0.8.0.
## ℹ Please use a list of either functions or lambdas:
## 
## # Simple named list: list(mean = mean, median = median)
## 
## # Auto named with `tibble::lst()`: tibble::lst(mean, median)
## 
## # Using lambdas list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## There were 14 warnings in `summarise()`. The first warning was: ℹ In argument:
## `resp = chisq.test(resp, Orig$resp)$p.value`. Caused by warning in
## `chisq.test()`: ! Chi-squared approximation may be incorrect ℹ Run
## `dplyr::last_dplyr_warnings()` to see the 13 remaining warnings.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

  ChiSq.pval

## # A tibble: 1 × 14
##    resp  crim            zn  indus  chas    nox    rm   age   dis      rad   tax
##   <dbl> <dbl>         <dbl>  <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>    <dbl> <dbl>
## 1     0 0.174 0.00000000365 0.0461 0.269 0.0478 0.604 0.963 0.814 0.000286 0.113
## # ℹ 3 more variables: ptratio <dbl>, black <dbl>, lstat <dbl>

  barplot(t(t(as.matrix(ChiSq.pval))), las=2, cex.names=1,
      main="p-values from Chi-sq test of association")
  abline(h=.05, col='red')

  which(ChiSq.pval < .05)             # Col num of variables w <.05 p-value

## [1]  1  3  4  6 10

  list_cols[which(ChiSq.pval < .05)]  # Top list of 'important' variables

## [1] "resp"  "zn"    "indus" "nox"   "rad"

  which(ChiSq.pval > .6)              # Col num of variables w >.6 p-value

## [1]  7  8  9 13

  list_cols[which(ChiSq.pval > .6)]   # Top list of 'unimportatnt' variables

## [1] "rm"    "age"   "dis"   "black"

### End of Routine Exploratory Analysis.
### Try repeating with Boston3.

CODE ONLY