IPDA Assignment 2 - Code
2024-05-27
# Import all libraries
library(caret)## Loading required package: ggplot2## Loading required package: latticelibrary(ggplot2)
library(GGally)## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2library(BMA)## Loading required package: survival##
## Attaching package: 'survival'## The following object is masked from 'package:caret':
##
## cluster## Loading required package: leaps## Loading required package: robustbase##
## Attaching package: 'robustbase'## The following object is masked from 'package:survival':
##
## heart## Loading required package: inline## Loading required package: rrcov## Scalable Robust Estimators with High Breakdown Point (version 1.7-5)library(table1)##
## Attaching package: 'table1'## The following objects are masked from 'package:base':
##
## units, units<-library(tidyverse)## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ lubridate 1.9.3 ✔ tibble 3.2.1
## ✔ purrr 1.0.2 ✔ tidyr 1.3.1## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ✖ purrr::lift() masks caret::lift()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors1 - Present the study design and aims (5%)
The study is a cross-sectional investigation of house conditions in 506 suburbs in the Boston area. The aim is to develop and validate a model for predicting housing prices.
2 - Describe the characteristics of the study sample (5%)
housing_data = read.csv("C:\\Users\\User\\Downloads\\Housing prices data.csv")
table1(~ ID + crime + zone + industry + river + nox + age + distance + radial + ptratio + lstat + price + rooms, data=housing_data)| Overall (N=506) |
|
|---|---|
| ID | |
| Mean (SD) | 254 (146) |
| Median [Min, Max] | 254 [1.00, 506] |
| crime | |
| Mean (SD) | 3.61 (8.60) |
| Median [Min, Max] | 0.257 [0.00632, 89.0] |
| zone | |
| Mean (SD) | 11.4 (23.3) |
| Median [Min, Max] | 0 [0, 100] |
| industry | |
| Mean (SD) | 11.1 (6.86) |
| Median [Min, Max] | 9.69 [0.460, 27.7] |
| river | |
| Mean (SD) | 0.0692 (0.254) |
| Median [Min, Max] | 0 [0, 1.00] |
| nox | |
| Mean (SD) | 0.555 (0.116) |
| Median [Min, Max] | 0.538 [0.385, 0.871] |
| age | |
| Mean (SD) | 68.6 (28.1) |
| Median [Min, Max] | 77.5 [2.90, 100] |
| distance | |
| Mean (SD) | 3.80 (2.11) |
| Median [Min, Max] | 3.21 [1.13, 12.1] |
| radial | |
| Mean (SD) | 9.55 (8.71) |
| Median [Min, Max] | 5.00 [1.00, 24.0] |
| ptratio | |
| Mean (SD) | 18.5 (2.16) |
| Median [Min, Max] | 19.1 [12.6, 22.0] |
| lstat | |
| Mean (SD) | 12.7 (7.14) |
| Median [Min, Max] | 11.4 [1.73, 38.0] |
| price | |
| Mean (SD) | 22.5 (9.20) |
| Median [Min, Max] | 21.2 [5.00, 50.0] |
| rooms | |
| Mean (SD) | 6.28 (0.703) |
| Median [Min, Max] | 6.21 [3.56, 8.78] |
3 - Develop a model for predicting housing prices (30%)
3.1 - Examine the distribution of housing prices and the correlation among its candidate predictors. Interpret the graph.
Price distribution graph - The distribution of housing prices is skewed right, with a rather sharp drop after the peak, indicating a ‘non-normal’ graph. From this, we can see that the majority of suburbs have a median house price sitting in the $15,000-$25,000 price band. Correlation matrix - Notably, almost all variables have a strong correlation with the median house price of suburb (corr=> |±0.7|) except for the ‘river’ variable.
ggplot(housing_data, aes(x = price)) +
geom_histogram(binwidth = 1, fill = "blue", color = "black") +
theme_minimal() +
labs(title = "Distribution of Housing Prices", x = "Price (x $1000)", y = "Frequency")
# Correlation matrix
correlation_matrix = cor(housing_data %>% select(-ID))
ggcorr(correlation_matrix, label = TRUE)
### 3.2 - Develop a model for predicting housing prices using a Bayesian
Model Averaging (BMA) approach
xvars = housing_data[, c("lstat", "ptratio", "radial", "distance", "age", "rooms", "nox", "river", "industry", "zone", "crime")]
yvar = housing_data[,c("price")]
bma_model = bicreg(x = xvars, y = yvar, strict = FALSE, OR = 20)
# Summary of BMA model
summary(bma_model)##
## Call:
## bicreg(x = xvars, y = yvar, strict = FALSE, OR = 20)
##
##
## 10 models were selected
## Best 5 models (cumulative posterior probability = 0.8713 ):
##
## p!=0 EV SD model 1 model 2 model 3
## Intercept 100.0 38.066751 5.29477 39.98405 34.99981 40.80894
## lstat 100.0 -0.555997 0.04916 -0.55346 -0.54777 -0.55201
## ptratio 100.0 -0.974062 0.15049 -1.00175 -0.85618 -1.11544
## radial 49.7 0.054116 0.06207 0.10457 . 0.11877
## distance 100.0 -1.356310 0.22229 -1.45079 -1.45011 -1.19964
## age 0.0 0.000000 0.00000 . . .
## rooms 100.0 3.961303 0.42193 3.85056 4.00998 3.99022
## nox 100.0 -19.756773 3.80180 -21.37566 -17.71878 -21.95511
## river 100.0 3.175507 0.87748 3.13944 3.17045 3.09444
## industry 6.8 -0.004933 0.02359 . . .
## zone 67.1 0.025403 0.02106 0.03658 0.04109 .
## crime 76.4 -0.077851 0.05413 -0.11854 -0.08042 -0.11061
##
## nVar 9 8 8
## r2 0.727 0.724 0.723
## BIC -601.35620 -600.89285 -600.23419
## post prob 0.275 0.218 0.157
## model 4 model 5
## Intercept 37.12121 36.92263
## lstat -0.57683 -0.56984
## ptratio -0.91872 -1.00275
## radial . .
## distance -1.38632 -1.14459
## age . .
## rooms 3.93816 4.11181
## nox -18.84096 -18.74043
## river 3.30769 3.24430
## industry . .
## zone 0.03512 .
## crime . .
##
## nVar 7 6
## r2 0.720 0.716
## BIC -599.84863 -599.17436
## post prob 0.129 0.092par(mfrow = c(1,1))
imageplot.bma(bma_model)
model = lm(price ~ lstat + ptratio + radial + distance + rooms + nox + river + zone + crime, data = housing_data)
summary(model)##
## Call:
## lm(formula = price ~ lstat + ptratio + radial + distance + rooms +
## nox + river + zone + crime, data = housing_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -15.8550 -2.9880 -0.5477 1.8770 26.4105
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 39.98405 4.94523 8.085 4.78e-15 ***
## lstat -0.55346 0.04797 -11.537 < 2e-16 ***
## ptratio -1.00175 0.13049 -7.677 8.74e-14 ***
## radial 0.10457 0.04071 2.569 0.010495 *
## distance -1.45079 0.18905 -7.674 8.92e-14 ***
## rooms 3.85056 0.41105 9.368 < 2e-16 ***
## nox -21.37566 3.50093 -6.106 2.07e-09 ***
## river 3.13944 0.86975 3.610 0.000338 ***
## zone 0.03658 0.01357 2.695 0.007279 **
## crime -0.11854 0.03330 -3.559 0.000407 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 4.847 on 496 degrees of freedom
## Multiple R-squared: 0.7273, Adjusted R-squared: 0.7223
## F-statistic: 146.9 on 9 and 496 DF, p-value: < 2.2e-163.3 - Check the model’s assumptions.
# Residual diagnostics
par(mfrow = c(2, 2))
plot(model)
# Reset plotting layout
par(mfrow = c(1, 1))3.4 - Present its mathematical equation in the format of Y = I + 1X1 + 2X2 + … + nXn
price = 39.984 - 0.553lstat - 1.002ptratio + 0.105radial - 1.451distance + 3.851rooms - 21.378nox + 3.139river + 0.037zone - 0.119crime
3.5 - Interpret the findings
From a common sense perspective this seems to be an accurate model. Detractions from the housing price (lstat, ptratio, distance, nox, crime) and additions (radial, rooms, river, zone) all make immediate sense regarding housing prices. However, there are some concerns based on the model assumptions. The results of the RvsF and S-L graph suggest non-linearity and heteroscedasicity, respectively.
4 - Validate the model’s performance internally (30%)
4.1 - Identify the optimal internal validation method used to validate the model’s performance. Explain the reason(s).
I have chosen cross-validation/k-fold validation as the dataset is not extensively large (thereby not too computationally intensive/better generalisation) and the reduced chance of result variance (compared to test/train split).
4.2 - Conduct an internal validation to assess the model’s performance.
train_control = trainControl(method = "cv", number = 10)
cv_model = train(price ~ lstat + ptratio + radial + distance + rooms + nox + river + zone + crime, data = housing_data, method = "lm", trControl = train_control)4.3 - Present and interpret the findings
From the R-squared result we know the model to have an accuracy of 72.45%, which is okay but not incredible. With an MSE(adj) of ±$3485.32 and RMSE(adj) of ±4874.87, we can make a few possible judgements on the model and how it could be improved. First, there are some large errors given by the model, either several across the board or a smaller selection of larger errors. Second, we can suggest that the model may benefit from techniques to handle outliers/reduce the influence of larger errors (e.g., regularisation/transforming the target variable).
print(cv_model)## Linear Regression
##
## 506 samples
## 9 predictor
##
## No pre-processing
## Resampling: Cross-Validated (10 fold)
## Summary of sample sizes: 457, 454, 455, 455, 456, 456, ...
## Resampling results:
##
## RMSE Rsquared MAE
## 4.864676 0.7228761 3.482862
##
## Tuning parameter 'intercept' was held constant at a value of TRUE