1. Data loading and prep

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


##----------------------------------------
### Data Prep

library(tidyverse)            # install.packages("tidyverse")
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.2.1     ✔ readr     2.2.0
## ✔ forcats   1.0.1     ✔ stringr   1.6.0
## ✔ ggplot2   4.0.2     ✔ tibble    3.3.1
## ✔ lubridate 1.9.5     ✔ tidyr     1.3.2
## ✔ purrr     1.2.1     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ✖ dplyr::select() masks MASS::select()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Turn Boston into tibble
Boston <- as_tibble(Boston)
Boston

# Rename "medv" column as "resp"
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
Boston3 <- Boston2 %>% mutate( chas=as.factor(chas) )
Boston3


2. Three variable model with surface plot

library(plotly)  # install.packages('plotly', repos='https://cran.case.edu/')

# 1. Prepare data with a 2-level categorical factor
DF = Boston3

# 2. Fit the model (2 continuous predictors + 1 categorical predictor)
Fit01 <- lm(resp ~ zn + rm + chas, data = DF)
summary(Fit01)
## 
## Call:
## lm(formula = resp ~ zn + rm + chas, data = DF)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -24.705  -3.287   0.239   2.778  39.400 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -30.46972    2.65423 -11.480  < 2e-16 ***
## zn            0.06657    0.01285   5.182 3.19e-07 ***
## rm            8.26353    0.42787  19.313  < 2e-16 ***
## chas1         4.52123    1.12556   4.017 6.80e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.38 on 502 degrees of freedom
## Multiple R-squared:  0.5217, Adjusted R-squared:  0.5189 
## F-statistic: 182.5 on 3 and 502 DF,  p-value: < 2.2e-16
# 3. Setup the coordinate grid for the continuous variables
x_range <- seq(min(DF$zn), max(DF$zn), length.out = 30)
y_range <- seq(min(DF$rm), max(DF$rm), length.out = 30)
grid <- expand.grid(zn = x_range, rm = y_range)

# 4. Generate distinct prediction matrices for both levels
grid_A <- transform(grid, chas = factor("0", levels=levels(DF$chas)))
grid_B <- transform(grid, chas = factor("1", levels=levels(DF$chas)))

z_A <- matrix(predict(Fit01, newdata = grid_A), nrow = 30, ncol = 30)
z_B <- matrix(predict(Fit01, newdata = grid_B), nrow = 30, ncol = 30)

# 5. Define matching colors (V-Engine = Orange, Straight = Blue)
colors_pts   <- c("#E69F00", "#56B4E9")
color_A_pln  <- list(c(0, 1), c("#E69F00", "#E69F00"))
color_B_pln  <- list(c(0, 1), c("#56B4E9", "#56B4E9"))

# 6. Build the 3D plot
plot_ly(DF, x = ~zn, y = ~rm, z = ~resp) %>%
  # Scatter points
  add_markers(color = ~chas, colors=colors_pts, marker=list(size=6), opacity=0.9) %>%
  # V-Engine plane
  add_surface(x = x_range, y = y_range, z = z_A, name = "chas 0 plane",
              colorscale = color_A_pln, opacity = 0.4, showscale = FALSE) %>%
  # Straight-Engine plane
  add_surface(x = x_range, y = y_range, z = z_B, name = "chas 1 plane",
              colorscale = color_B_pln, opacity = 0.4, showscale = FALSE) %>%
  # Labels and style
  layout(scene = list(
    xaxis = list(title = "zn"),
    yaxis = list(title = "rn"),
    zaxis = list(title = "resp")
  ))

\[\begin{align*} Y = -30.5 + 0.067 (zn) + 8.26 (rm) + 4.52 (chas) \end{align*}\]

# Scatter Plots for each variables
plot(DF$zn, DF$resp)

plot(DF$rm, DF$resp)

plot(DF$chas, DF$resp)



COCE ONLY

##--------------------------------------
## 1. Data loading and prep

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


##----------------------------------------
### Data Prep

library(tidyverse)            # install.packages("tidyverse")

# Turn Boston into tibble
Boston <- as_tibble(Boston)
Boston

# Rename "medv" column as "resp"
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
Boston3 <- Boston2 %>% mutate( chas=as.factor(chas) )
Boston3


##--------------------------------------
## 2. Three variable model with surface plot

library(plotly)  # install.packages('plotly', repos='https://cran.case.edu/')

# a. Prepare data with a 2-level categorical factor
DF = Boston3

# b. Fit the model (2 continuous predictors + 1 categorical predictor)
Fit01 <- lm(resp ~ zn + rm + chas, data = DF)
summary(Fit01)

# c. Setup the coordinate grid for the continuous variables
x_range <- seq(min(DF$zn), max(DF$zn), length.out = 30)
y_range <- seq(min(DF$rm), max(DF$rm), length.out = 30)
grid <- expand.grid(zn = x_range, rm = y_range)

# d. Generate distinct prediction matrices for both levels
grid_A <- transform(grid, chas = factor("0", levels=levels(DF$chas)))
grid_B <- transform(grid, chas = factor("1", levels=levels(DF$chas)))

z_A <- matrix(predict(Fit01, newdata = grid_A), nrow = 30, ncol = 30)
z_B <- matrix(predict(Fit01, newdata = grid_B), nrow = 30, ncol = 30)

# e. Define matching colors (V-Engine = Orange, Straight = Blue)
colors_pts   <- c("#E69F00", "#56B4E9")
color_A_pln  <- list(c(0, 1), c("#E69F00", "#E69F00"))
color_B_pln  <- list(c(0, 1), c("#56B4E9", "#56B4E9"))

# f. Build the 3D plot
plot_ly(DF, x = ~zn, y = ~rm, z = ~resp) %>%
  # Scatter points
  add_markers(color = ~chas, colors=colors_pts, marker=list(size=6), opacity=0.9) %>%
  # V-Engine plane
  add_surface(x = x_range, y = y_range, z = z_A, name = "chas 0 plane",
              colorscale = color_A_pln, opacity = 0.4, showscale = FALSE) %>%
  # Straight-Engine plane
  add_surface(x = x_range, y = y_range, z = z_B, name = "chas 1 plane",
              colorscale = color_B_pln, opacity = 0.4, showscale = FALSE) %>%
  # Labels and style
  layout(scene = list(
    xaxis = list(title = "zn"),
    yaxis = list(title = "rn"),
    zaxis = list(title = "resp")
  ))


# Y = -30.5 + 0.067 (zn) + 8.26 (rm) + 4.52 (chas)


##--------------------------------------
# Scatter Plots for each variables
plot(DF$zn, DF$resp)
plot(DF$rm, DF$resp)
plot(DF$chas, DF$resp)