4.1 Initial fit

Let us begin by fitting a simple MARS model with default parameters (degree as 1) and observe the summary.

# let us use MARS (Multivariate Adaptive Regression Spline)
library(earth)    # for MARS
library(pdp)      # for partial dependence plots
library(vip)      # for variable importance plots

# a simple fit
boston_mars <- earth(cmedv ~ ., data = boston,
                     degree = 1  # tuning parameter
)
# MARS summary
summary(boston_mars) 

## Call: earth(formula=cmedv~., data=boston, degree=1)
## 
##                 coefficients
## (Intercept)        -2.743460
## h(-70.99-lon)      23.712668
## h(lon- -70.99)     38.145330
## h(lat-42.2835)     43.008986
## h(4.42228-crim)    -0.572711
## h(crim-4.42228)    -0.409119
## h(crim-18.811)      0.330218
## h(0.488-nox)      -24.939405
## h(nox-0.488)      -19.484251
## h(rm-6.431)         5.998860
## h(rm-7.313)         8.569941
## h(rm-7.929)       -18.441106
## h(dis-1.5106)      34.885916
## h(2.4358-dis)      36.199897
## h(dis-2.4358)     -35.834647
## h(300-tax)          0.029591
## h(14.7-ptratio)     2.136646
## h(ptratio-14.7)    -0.458695
## h(386.75-b)        -0.005791
## h(b-386.75)        -0.116258
## h(6.12-lstat)       2.105861
## h(lstat-6.12)      -0.609421
## h(lstat-23.79)      0.505551
## 
## Selected 23 of 26 terms, and 10 of 15 predictors
## Termination condition: Reached nk 31
## Importance: rm, lstat, ptratio, dis, crim, nox, tax, lon, b, lat, ...
## Number of terms at each degree of interaction: 1 22 (additive model)
## GCV 11.46153    RSS 4813.864    GRSq 0.8643274    RSq 0.8869394

OBSERVATIONS:

1. Rsq of the best model as per default parameters is 0.887. This indicates that there has been an improvement in comparison to the linear model.

2. Unlike the case of Linear Regression, p-values or confidence intervals are not provided in the summary of MARS for the coefficients.

3. MARS has forward pass where variable + knot combination that gives the greatest improvement to the current model are added and pruning pass where basis functions are removed one at a time. Thus, in the summary we can observe something like this : (crim-4.42228) along with its coefficient.

4.2 Tuning the parameter: degree

Let us tune the parameter degree, to obtain its optimal values for modelling. For doing the same, I shall use the caret package. The parameter is optimized keeping Rsquared as the metric (the higher the better) using 5-fold Cross Validation repeated thrice.

# tuning parameters
library(caret)
getModelInfo("earth")$earth$parameters 

##   parameter   class          label
## 1    nprune numeric         #Terms
## 2    degree numeric Product Degree

# Tune a MARS model. Setting nprune max as 100
set.seed(103)  # for reproducibility
boston_mars_tune <- train(
  x = subset(boston, select = -cmedv),
  y = boston$cmedv,
  method = "earth",
  metric = "Rsquared",
  trControl = trainControl(method = "repeatedcv", 
                           number = 5, repeats = 3),
  tuneGrid = expand.grid(degree = 1:5, nprune = 100) # tuning only degree and not nprune
)

# Print model tuning summary
print(boston_mars_tune)

## Multivariate Adaptive Regression Spline 
## 
## 506 samples
##  15 predictor
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 3 times) 
## Summary of sample sizes: 406, 404, 404, 405, 405, 405, ... 
## Resampling results across tuning parameters:
## 
##   degree  RMSE      Rsquared   MAE     
##   1       3.660745  0.8398051  2.514046
##   2       3.657539  0.8497680  2.318028
##   3       3.755125  0.8414312  2.370777
##   4       3.764644  0.8402414  2.381942
##   5       3.764644  0.8402414  2.381942
## 
## Tuning parameter 'nprune' was held constant at a value of 100
## Rsquared was used to select the optimal model using the largest value.
## The final values used for the model were nprune = 100 and degree = 2.

# Plot model tuning summary
ggplot(boston_mars_tune) + theme_light()

The above plot indicates that degree 2 provides the highest Rsquared was obtained for degree equal to 2. Corresponding RMSE and MAE values are low as well.

4.3 The best model

We know that best degree parameter is 2. Let us find out the best model using pmethod as cv. The selection of number of terms, and predictors will be determined by “cv” or Cross Validation (10 folds done thrice) rather than by default standard GCV statistic.

# the best parameter: degree = 2
set.seed(103)
fit1 <- earth(cmedv~., data=boston, ncross=3, nfold=10, degree = 2, pmethod="cv",keepxy=TRUE)

summary(fit1)

## Call: earth(formula=cmedv~., data=boston, pmethod="cv", keepxy=TRUE, degree=2,
##             nfold=10, ncross=3)
## 
##                               coefficients
## (Intercept)                       28.96957
## h(42.2248-lat)                   -18.68634
## h(lat-42.2248)                   -15.76396
## h(crim-4.42228)                   -0.11324
## h(6.431-rm)                       -0.45963
## h(rm-6.431)                       10.30938
## h(age-18.5)                       -0.07599
## h(1.3567-dis)                    941.03891
## h(dis-1.3567)                     -0.50888
## h(lstat-6.12)                     -0.79467
## h(lstat-22.11)                     0.62181
## h(1.3567-dis) * b                 -2.34986
## h(-70.965-lon) * h(age-18.5)       0.31074
## h(lon- -70.965) * h(age-18.5)      1.32895
## h(4.42228-crim) * h(tax-224)      -0.00158
## h(4.42228-crim) * h(224-tax)       0.01795
## h(0.624-nox) * h(rm-6.431)        -5.26121
## h(nox-0.624) * h(rm-6.431)       -88.36153
## h(0.713-nox) * h(lstat-6.12)       2.05404
## h(nox-0.713) * h(lstat-6.12)       1.02961
## h(6.431-rm) * h(dis-1.8209)       -0.49396
## h(rm-6.431) * h(ptratio-18.6)     -3.05364
## h(rm-6.431) * h(18.6-ptratio)      0.62770
## h(307-tax) * h(6.12-lstat)         0.02077
## h(tax-307) * h(6.12-lstat)         0.02198
## 
## Selected 25 of 29 terms, and 11 of 15 predictors (pmethod="cv")
## Termination condition: Reached nk 31
## Importance: rm, lstat, nox, tax, dis, ptratio, b, crim, lon, age, lat, ...
## Number of terms at each degree of interaction: 1 10 14
## GRSq 0.9096501  RSq 0.929844  mean.oof.RSq 0.8352213 (sd 0.0827)
## 
## pmethod="backward" would have selected:
##     22 terms 11 preds,  GRSq 0.9113226  RSq 0.928802  mean.oof.RSq 0.8326279

par(mfrow = c(1, 2))

#model selection
q1 <- plot(fit1, which=1,
col.mean.infold.rsq="blue", col.infold.rsq="lightblue",
col.grsq=0, col.rsq=0, col.vline=0, col.oof.vline=0)

# model selection
q2 <- plotres(fit1, which=1, info = TRUE)

INTERPRETATION:

1. The plot on the left indicates training and testing performance (Rsq on Y axis) obtained from the 10 fold cross validation performed thrice. The performance on training data (blue curve) increases as we increase model complexity; on independent data the performance (pink curve) peaks and then decreases.

2. The plot on the right shows the best model selection (25 of 29 terms, 11 of 15 predictors using pmethod=“cv”) using green line, indicating optimal terms as 25.

5.1 Residual vs Fitted values

# residuals vs fitted
plot(fit1, which = 3) # which = 3 for residuals vs fitted

The Residuals vs Fitted graph shows the residual for each value of the predicted response. The red line is a lowess fit. In this instance, the red line is almost constant and lying around 0. Thus, the mean of residuals is almost 0.

Visually, the variance of the residuals is constant and does not spread along the fitted values. Thus, this model is fulfilling the homoscedasticity test. However, in earth, constant variance of the residuals isn’t as important as it is in linear models.

However, this plot does indicate some cases of outliers. Those residuals are marked with numbers (372, 427, 407) indicating extreme values, or leverage / influential points.

5.2 Cumulative Distribution of residuals

plotres(fit1, which = 2, info = TRUE) # which = 2 for cumulative distribution

The Cumulative Distribution graph shows the cumulative distribution of the absolute values of residuals. What we would ideally like to see is a graph that starts at 0 and shoots up quickly to 1.

We see that 95% of the absolute values of residuals are less than about 4.61 (look at the vertical gray line for 95%). So in the training data, 95% of the time the predicted value is within 4.61 units of the observed value.

5.3 Residual QQ plot

plotres(fit1, which = 4, info = TRUE) # which = 4 for residual qqplot

The QQ (quantile-quantile) plot compares the distribution of the residuals to a normal distribution. If the residuals are distributed normally they will lie on the line.

Here, the distribution is approximately normal (look at the actual vs normal at the bottom) with heavy tailed data. Normality of the residuals often isn’t too important for earth models, but the graph is useful for discovering outlying residuals and other anomalies.

As seen earlier, observation points 427, 407, 372 do seem to be outliers/leverage/influential points. Ideally, they should be removed from the data before modelling is implemented or further investigation should be done after those observation points. For simplicity, we ignore those outliers and continue ahead.

The best model passes the residual analyses tests. We can proceed ahead to communicate our results to the stakeholders.

6.1 Cross Validation performance of the best model

A simple 10 fold cross validation has been performed on the data and following metrics were calculated to assess the performance of the best model.

cv <- modelr::crossv_kfold(boston, 10)
models <- purrr::map(cv$train, ~ fit1)

#RMSE
error_RMSE <- purrr::map2_dbl(models, cv$test, modelr::rmse)

#MAE
error_MAE <- purrr::map2_dbl(models, cv$test, modelr::mae)

#MAPE
error_MAPE <- purrr::map2_dbl(models, cv$test, modelr::mape)

# summary of errors
tibble("RMSE" = mean(error_RMSE),
        "MAE" = mean(error_MAE),
       "MAPE" = mean(error_MAPE))

## # A tibble: 1 x 3
##    RMSE   MAE   MAPE
##   <dbl> <dbl>  <dbl>
## 1  2.41  1.79 0.0958

The above are the metrics for the best model.

1. RMSE interpretation: Typical errors are around, on average, 2.41 USD (in 1000’s) away from actual value.
2. MAPE interpretation: Typical errors are around 9.6 % of the actual value.

Ideally one does not perform hyperparameter tuning, model selection and cross validation on the same dataset. However, due to low number of observations for this dataset, I have performed the above three on the same dataset. This might lead to biased results.

6.2 Important Variables

The following plot shows us the relative importance as per the best model fit by earth.

# Variable importance plot
vip(
  fit1, 
  num_features = 15
)

To properly understand the association of these variables with cmedv, let us look at the partial dependence plots.

6.3 Association of the top three variables (rm, lstat and nox) on cmedv

The following plots show association of predictor variables one at a time with cmedv.

# Partial dependence of cmedv on rm
p1 <- fit1 %>%
  pdp::partial(pred.var = "rm") %>%
  autoplot(color = "red2", size = 1) +
  geom_point(data = boston, aes(x = rm, y = cmedv), alpha = 0.1) +
  theme_light()

# Partial dependence of cmedv on lstat
p2 <- fit1 %>%
  pdp::partial(pred.var = "lstat") %>%
  autoplot(color = "red2", size = 1) +
  geom_point(data = boston, aes(x = lstat, y = cmedv), alpha = 0.1) +
  theme_light()

# Partial dependence of cmedv on nox
p3 <- fit1 %>%
  pdp::partial(pred.var = "nox") %>%
  autoplot(color = "red2", size = 1) +
  geom_point(data = boston, aes(x = nox, y = cmedv), alpha = 0.1) +
  theme_light()

# 3 figures arranged in 2 rows and 2 columns
annotate_figure(ggarrange(p1, p2, p3,  ncol = 2, nrow = 2),
                top = text_grob("Partial dependence plots"))

1. cmedv increases as rm increases. Both these variables seem to be directly proportional. This is logically coherent with the fact that as the average number of rooms increase, the housing price will increase.

2. Opposite is the case with lstat. They seem to be indirectly proportional. The more the percentage of lower status of population in that area, lower will be the housing price.

3. nox and cmedv are indirectly proportional. This is logically coherent with the fact that housing prices will be high (and buyers would be ready to pay more) where emissions (Nitric oxides concentration due to pollution) are lower.

Let us look at the combined effect of two variables (like rm and lstat on cmedv).

# Partial dependence of cmedv on rm and lstat
fit1 %>%
  pdp::partial(pred.var = c("rm", "lstat"), chull = TRUE) %>%
  autoplot() +
  theme_light()

This 2-D plot gives us an idea about how the association of two variables at a time on cmedv. On the X- Axis we have rm and on the Y-Axis we have lstat. The variation is house prices cmedv is indicated with the help of color scale (yellow indicating high prices and dark violet indicating low prices).

The following 3-D plot is a good visualization tool to understand the association but is difficult to interpet.

# Partial dependence of cmedv on both rm and nox
pd <- pdp::partial(
  fit1, 
  pred.var = c("rm", "nox"),
  chull = TRUE
)
# Interactive 3-D plot
plotly::plot_ly(
  x = ~rm, 
  y = ~nox, 
  z = ~yhat, 
  data = pd,
  type = "mesh3d"
)