top of page

Regression Predictive Modelling Problem In R | Investigate The Boston House Price Dataset

Problem Definition

For this project we will investigate the Boston House Price dataset. It is already included in the library mlbench, so you just need to install the package, as follows:


install.library("mlbench") 
# load the package 
library(mlbench) 
# list the contents of the package 
library(help = "mlbench") 
# Boston Housing Data 
data(BostonHousing) 

Each record in the database describes a Boston suburb or town. The data was drawn from the Boston Standard Metropolitan Statistical Area (SMSA) in 1970. The attributes are defined as follows (taken from the UCI Machine Learning Repository)


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 town2

4. CHAS: Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)

5. NOX: nitric 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 distances to five Boston employment centers

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. B: 1000(𝐵𝑘 − 0.63) 2 , where 𝐵𝑘 is the proportion of blacks by town

13. LSTAT: % lower status of the population

14. MEDV: Median value of owner-occupied homes in $1000s We can see that the input attributes have a mixture of units


Load the Dataset

The dataset is available in the mlbench package. Let's start by offloading the required packages and loading the dataset.


# load packages 
library(mlbench) 
library(caret) 
library(corrplot) 
# attach the BostonHousing dataset 
data(BostonHousing) 

Validation Dataset

It is a good idea to use a validation hold out set. This is a sample of the data that we hold back from our analysis and modelling. We use it right at the end of our project to confirm the accuracy of our final model. It is a smoke test that we can use to see if we messed up and to give us confidence on our estimates of accuracy on unseen data.

# Split out validation dataset 
# create a list of 80% of the rows in the original dataset we can use for  training 
set.seed(7)
validationIndex <- createDataPartition(BostonHousing$medv, p=0.80,  list=FALSE) 
# select 20% of the data for validation 
validation <- BostonHousing[-validationIndex,] 
# use the remaining 80% of data to training and testing the models 
dataset <- BostonHousing[validationIndex,] 

Analyse Data

The objective of this step in the process is to better understand the problem.


Descriptive Statistics

Let's start off by confirming the dimensions of the dataset, e.g., the number of rows and columns.

# dimensions of dataset 
dim(dataset) 

We have 407 instances to work with and can confirm the data has 14 attributes including the class attribute medv.

407 14


Let's also look at the data types of each attribute.

# list types for each attribute 
sapply(dataset, class) 

We can see that one of the attributes (chas) is a factor while all of the others are numeric.


crim zn indus chas nox rm age dis rad tax

ptratio b

"numeric" "numeric" "numeric" "factor" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric" "numeric"

lstat medv

"numeric" "numeric"


Let's now take a peak at the first 20 rows of the data.

# take a peek at the first 5 rows of the data 
head(dataset, n=20)

We can confirm that the scales for the attributes are all over the place because of the differing units. We may benefit from some transforms later on.


crim zn indus chas nox rm age dis rad tax ptratio b lstat medv

2 0.02731 0.0 7.07 0 0.469 6.421 78.9 4.9671 2 242 17.8 396.90 9.14 21.6 3 0.02729 0.0 7.07 0 0.469 7.185 61.1 4.9671 2 242 17.8 392.83 4.03 34.7 4 0.03237 0.0 2.18 0 0.458 6.998 45.8 6.0622 3 222 18.7 394.63 2.94 33.4 5 0.06905 0.0 2.18 0 0.458 7.147 54.2 6.0622 3 222 18.7 396.90 5.33 36.2 6 0.02985 0.0 2.18 0 0.458 6.430 58.7 6.0622 3 222 18.7 394.12 5.21 28.7 7 0.08829 12.5 7.87 0 0.524 6.012 66.6 5.5605 5 311 15.2 395.60 12.43 22.9 8 0.14455 12.5 7.87 0 0.524 6.172 96.1 5.9505 5 311 15.2 396.90 19.15 27.1 9 0.21124 12.5 7.87 0 0.524 5.631 100.0 6.0821 5 311 15.2 386.63 29.93 16.5 13 0.09378 12.5 7.87 0 0.524 5.889 39.0 5.4509 5 311 15.2 390.50 15.71 21.7 14 0.62976 0.0 8.14 0 0.538 5.949 61.8 4.7075 4 307 21.0 396.90 8.26 20.4 15 0.63796 0.0 8.14 0 0.538 6.096 84.5 4.4619 4 307 21.0 380.02 10.26 18.2 16 0.62739 0.0 8.14 0 0.538 5.834 56.5 4.4986 4 307 21.0 395.62 8.47 19.9 17 1.05393 0.0 8.14 0 0.538 5.935 29.3 4.4986 4 307 21.0 386.85 6.58 23.1 18 0.78420 0.0 8.14 0 0.538 5.990 81.7 4.2579 4 307 21.0 386.75 14.67 17.5 19 0.80271 0.0 8.14 0 0.538 5.456 36.6 3.7965 4 307 21.0 288.99 11.69 20.2 20 0.72580 0.0 8.14 0 0.538 5.727 69.5 3.7965 4 307 21.0 390.95 11.28 18.2 23 1.23247 0.0 8.14 0 0.538 6.142 91.7 3.9769 4 307 21.0 396.90 18.72 15.2 25 0.75026 0.0 8.14 0 0.538 5.924 94.1 4.3996 4 307 21.0 394.33 16.30 15.6 26 0.84054 0.0 8.14 0 0.538 5.599 85.7 4.4546 4 307 21.0 303.42 16.51 13.9 27 0.67191 0.0 8.14 0 0.538 5.813 90.3 4.6820 4 307 21.0 376.88 14.81 16.6


Let's summarize the distribution of each attribute.

# summarize attribute distributions 
summary(dataset) 

We can note that chas is a pretty unbalanced factor. We could transform this attribute to numeric to make calculating descriptive statistics and plots easier.


crim zn indus chas nox rm age

Min. : 0.00906 Min. : 0.00 Min. : 0.46 0:376 Min. :0.3850 Min. :3.863 Min. : 2.90


1st Qu.: 0.08556 1st Qu.: 0.00 1st Qu.: 5.19 1: 31 1st Qu.:0.4530 1st Qu.:5.873 1st Qu.: 45.05


Median : 0.28955 Median : 0.00 Median : 9.90 Median :0.5380 Median :6.185 Median : 77.70


Mean : 3.58281 Mean :10.57 Mean :11.36 Mean :0.5577 Mean :6.279 Mean : 68.83


3rd Qu.: 3.50464 3rd Qu.: 0.00 3rd Qu.:18.10 3rd Qu.:0.6310 3rd Qu.:6.611 3rd Qu.: 94.55


Max. :88.97620 Max. :95.00 Max. :27.74 Max. :0.8710 Max. :8.780 Max. :100.00


dis rad tax ptratio b lstat medv

Min. : 1.130 Min. : 1.000 Min. :188.0 Min. :12.60 Min. : 0.32 Min. : 1.730 Min. : 5.00


1st Qu.: 2.031 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:374.50 1st Qu.: 6.895 1st Qu.:17.05


Median : 3.216 Median : 5.000 Median :330.0 Median :19.00 Median :391.13 Median :11.500 Median :21.20


Mean : 3.731 Mean : 9.464 Mean :405.6 Mean :18.49 Mean :357.88 Mean :12.827 Mean :22.61


3rd Qu.: 5.100 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.27 3rd Qu.:17.175 3rd Qu.:25.00


Max. :10.710 Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90 Max. :37.970 Max. :50.00


Let's go ahead and convert chas to a numeric attribute.


dataset[,4] <- as.numeric(as.character(dataset[,4]))


Now, let's now take a look at the correlation between all of the numeric attributes. cor(dataset[,1:13])


This is interesting. We can see that many of the attributes have a strong correlation (e.g.> 0:70 or < 0:70). For example:

  • nox and indus with 0.77

  • dist and indus with 0.71

  • tax and indus with 0.72

  • age and nox with 0.72

  • dist and nox with 0.76

crim zn indus chas nox rm age dis rad tax ptratio b lstat

crim 1.00000000 -0.19790631 0.40597009 -0.05713065 0.4232413 -0.21513269 0.3543819 -0.3905097 0.64240501 0.60622608 0.2892983 -0.3021185 0.47537617


zn -0.19790631 1.00000000 -0.51895069 -0.04843477 -0.5058512 0.28942883 - 0.5707027 0.6561874 -0.29952976 -0.28791668 -0.3534121 0.1692749 -0.39712686


indus 0.40597009 -0.51895069 1.00000000 0.08003629 0.7665481 -0.37673408 0.6585831 -0.7230588 0.56774365 0.68070916 0.3292061 -0.3359795 0.59212718


chas -0.05713065 -0.04843477 0.08003629 1.00000000 0.1027366 0.08252441 0.1093812 -0.1114242 -0.00901245 -0.02779018 -0.1355438 0.0472442 -0.04569239


nox 0.42324132 -0.50585121 0.76654811 0.10273656 1.0000000 -0.29885055 0.7238371 -0.7708680 0.58516760 0.65217875 0.1416616 -0.3620791 0.58196447


rm -0.21513269 0.28942883 -0.37673408 0.08252441 -0.2988506 1.00000000 - 0.2325359 0.1952159 -0.19149122 -0.26794733 -0.3200037 0.1553992 -0.62038075


age 0.35438190 -0.57070265 0.65858310 0.10938121 0.7238371 -0.23253586 1.0000000 -0.7503321 0.45235421 0.50164657 0.2564318 -0.2512574 0.59321281


This is collinearity and we may see better results with regression algorithms if the correlated attributes are removed.


Unimodal Data Visualizations

Let's look at visualizations of individual attributes. It is often useful to look at your data using multiple different visualizations in order to spark ideas. Let's look at histograms of each attribute to get a sense of the data distributions.


# histograms each attribute 
par(mfrow=c(2,7)) 
for(i in 1:13) { 
hist(dataset[,i], main=names(dataset)[i]) 
}


We can see that some attributes may have an exponential distribution, such as crim, zn, age and b. We can see that others may have a bimodal distribution such as rad and tax


Result:












Let's look at the same distributions using density plots that smooth them out a bit.


# density plot for each attribute 
par(mfrow=c(2,7)) 
for(i in 1:13) { 
plot(density(dataset[,i]), main=names(dataset)[i]) 
} 

Output Result:











This perhaps adds more evidence to our suspicion about possible exponential and bimodal distributions. It also looks like nox, rm and lsat may be skewed Gaussian distributions, which might be helpful later with transforms.


Let's look at the data with box and whisker plots of each attribute.


# boxplots for each attribute 
par(mfrow=c(2,7)) 
for(i in 1:13) { 
boxplot(dataset[,i], main=names(dataset)[i]) 
}

This helps point out the skew in many distributions so much so that data looks like outliers (e.g., beyond the whisker of the plots).


Output Result:











Multi-modal Data Visualisations

Let's look at some visualisations of the interactions between variables. The best place to start is a scatterplot matrix.



# scatterplot matrix 
pairs(dataset[,1:13]) 

We can see that some of the higher correlated attributes do show good structure in their relationship. Not linear, but nice predictable curved relationships


Output Result:

Scatterplot Matrix of Boston House Dataset Input Attributes.













# correlation plot 
correlations <- cor(dataset[,1:13]) 
corrplot(correlations, method="circle") 

The larger darker blue dots confirm the positively correlated attributes we listed early (not the diagonal). We can also see some larger darker red dots that suggest some negatively correlated attributes. For example, tax and rad. These too may be candidates for removal to better improve accuracy of models later on


Output Result:














Summary of Ideas

There is a lot of structure in this dataset. We need to think about transforms that we could use later to better expose the structure which in turn may improve modelling accuracy. So far it would be worth trying:

  • Feature selection and removing the most correlated attributes.

  • Normalizing the dataset to reduce the effect of differing scales.

  • Standardizing the dataset to reduce the effects of differing distributions.

  • Box-Cox transform to see if flattening out some of the distributions improves accuracy.

With lots of additional time I would also explore the possibility of binning (discretization) of the data. This can often improve accuracy for decision tree algorithms.


Evaluate Algorithms: Baseline

We have no idea what algorithms will do well on this problem. Gut feel suggests regression algorithms like GLM and GLMNET may do well. It is also possible that decision trees and even SVM may do well. I have no idea. Let's design our test harness. We will use 10-fold cross-validation (each fold will be about 360 instances for training and 40 for test) with 3 repeats


The dataset is not too small, and this is a good standard test harness configuration. We will evaluate algorithms using the RMSE and R2 metrics. RMSE will give a gross idea of how wrong all predictions are (0 is perfect) and R2 will give an idea of how well the model has fit the data (1 is perfect, 0 is worst).


# Prepare the test harness for evaluating algorithms. 
# Run algorithms using 10-fold cross validation 
trainControl <- trainControl(method="repeatedcv", number=10, repeats=3) metric <- "RMSE"

Let's create a baseline of performance on this problem and spot-check a number of different algorithms. We will select a suite of different algorithms capable of working on this regression problem. The 6 algorithms selected include:

  • Linear Algorithms: Linear Regression (LR), Generalized Linear Regression (GLM) and Penalized Linear Regression (GLMNET)

  • Non-Linear Algorithms: Classification and Regression Trees (CART), Support Vector Machines (SVM) with a radial basis function and k-Nearest Neighbours (KNN)


We know the data has differing units of measure so we will standardize the data for this baseline comparison. This will those algorithms that prefer data in the same scale (e.g., instance-based methods and some regression algorithms) a chance to do well.


# Estimate accuracy of machine learning algorithms

# LM

set.seed(7)

fit.lm <- train(medv~., data=dataset, method="lm", metric=metric, preProc=c("center", "scale"), trControl=trainControl)


# GLM

set.seed(7)

fit.glm <- train(medv~., data=dataset, method="glm", metric=metric, preProc=c("center", "scale"), trControl=trainControl)


# GLMNET

set.seed(7)

fit.glmnet <- train(medv~., data=dataset, method="glmnet", metric=metric, preProc=c("center", "scale"), trControl=trainControl)


# SVM

set.seed(7)

fit.svm <- train(medv~., data=dataset, method="svmRadial", metric=metric,13 preProc=c("center", "scale"), trControl=trainControl)


# CART

set.seed(7) grid <- expand.grid(.cp=c(0, 0.05, 0.1))

fit.cart <- train(medv~., data=dataset, method="rpart", metric=metric, tuneGrid=grid, preProc=c("center", "scale"), trControl=trainControl)


# KNN

set.seed(7)

fit.knn <- train(medv~., data=dataset, method="knn", metric=metric, preProc=c("center", "scale"), trControl=trainControl)


The algorithms all use default tuning parameters, except CART which is fussy on this dataset and has 3 default parameters specified. Let's compare the algorithms. We will use a simple table of results to get a quick idea of what is going on. We will also use a dot plot to show the 95% confidence level for the estimated metrics.


# Collect resample statistics from models and summarize results. 
# Compare algorithms 
results <- resamples(list(LM=fit.lm, GLM=fit.glm, GLMNET=fit.glmnet,  SVM=fit.svm, CART=fit.cart, KNN=fit.knn)) 
summary(results) 
dotplot(results)

Output Result:

Dotplot Compare Machine Learning Algorithms on the Boston House Price Dataset with Box-Cox Power Transform














We can see that this indeed decrease the RMSE and increased the 𝑅 2 on all except the CART algorithms. The RMSE of SVM dropped to an average of 3.761.


# Output of estimated accuracy of models on transformed dataset.

Models: LM, GLM, GLMNET, SVM, CART, KNN20

Number of resamples: 30


RMSE

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's

LM 3.404 3.811 4.399 4.621 5.167 7.781 0

GLM 3.404 3.811 4.399 4.621 5.167 7.781 0

GLMNET 3.312 3.802 4.429 4.611 5.123 7.976 0

SVM 2.336 2.937 3.543 3.761 4.216 8.207 0

CART 2.797 3.434 4.272 4.541 5.437 9.248 0

KNN 2.474 3.608 4.308 4.563 5.080 8.922 0


Rsquared

Min. 1st Qu. Median Mean 3rd Qu. Max. NA's

LM 0.5439 0.7177 0.7832 0.7627 0.8257 0.8861 0

GLM 0.5439 0.7177 0.7832 0.7627 0.8257 0.8861 0

GLMNET 0.5198 0.7172 0.7808 0.7634 0.8297 0.8909 0

SVM 0.5082 0.8249 0.8760 0.8452 0.8998 0.9450 0

CART 0.3614 0.6733 0.8197 0.7680 0.8613 0.9026 0

KNN 0.4065 0.7562 0.8073 0.7790 0.8594 0.9043 0


Improve Results with Tuning

We can improve the accuracy of the well performing algorithms by tuning their parameters. In this section we will look at tuning the parameters of SVM with a Radial Basis Function (RBF). with more time it might be worth exploring tuning of the parameters for CART and KNN. It might also be worth exploring other kernels for SVM besides the RBF. Let's look at the default parameters already adopted.


# Display estimated accuracy of a model. 
print(fit.svm) 

The C parameter is the cost constraint used by SVM. Learn more in the help for the ksvm function? ksvm. We can see from previous results that a C value of 1.0 is a good starting point.


# Output of estimated accuracy of a model.

Support Vector Machines with Radial Basis Function Kernel

407 samples

13 predictor


Pre-processing: centered (13), scaled (13), Box-Cox transformation (11)

Resampling: Cross-Validated (10 fold, repeated 3 times)

Summary of sample sizes: 366, 367, 366, 366, 367, 367, ...

Resampling results across tuning parameters


C RMSE Rsquared RMSE SD Rsquared SD

0.25 4.555338 0.7906921 1.533391 0.11596196

0.50 4.111564 0.8204520 1.467153 0.10573527

1.00 3.761245 0.8451964 1.323218 0.09487941


Tuning parameter 'sigma' was held constant at a value of 0.07491936

RMSE was used to select the optimal model using the smallest value.

The final values used for the model were sigma = 0.07491936 and C = 1.


Let's design a grid search around a C value of 1. We might see a small trend of decreasing RMSE with increasing C, so let’s try all integer C values between 1 and 10. Another parameter that caret lets us tune is the sigma parameter. This is a smoothing parameter. Good sigma values are often start around 0.1, so we will try numbers before and after.


# Tune the parameters of a model. 
# tune SVM sigma and C parametres 
trainControl <- trainControl(method="repeatedcv", number=10, repeats=3) 
metric <- "RMSE" 
set.seed(7) 
grid <- expand.grid(.sigma=c(0.025, 0.05, 0.1, 0.15), .C=seq(1, 10, by=1)) 
fit.svm <- train(medv~., data=dataset, method="svmRadial", metric=metric,  tuneGrid=grid, preProc=c("BoxCox"), trControl=trainControl) 
print(fit.svm) 
plot(fit.svm) 

Output Result:

Algorithm Tuning Results for SVM on the Boston House Price Dataset.














We can see that the sigma values flatten out with larger C cost constraints. It looks like we might do well with a sigma of 0.05 and a C of 10. This gives us a respectable RMSE of 2.977085.


# Output of tuning the parameters of a model.

Support Vector Machines with Radial Basis Function Kernel

407 samples

13 predictor


Pre-processing: Box-Cox transformation (11)

Resampling: Cross-Validated (10 fold, repeated 3 times)


19.7. Ensemble Methods 174

Summary of sample sizes: 366, 367, 366, 366, 367, 367, ...


Resampling results across tuning parameters:

sigma C RMSE Rsquared RMSE SD Rsquared SD

0.025 1 3.889703 0.8335201 1.4904294 0.11313650

0.025 2 3.685009 0.8470869 1.4126374 0.10919207

0.025 3 3.562851 0.8553298 1.3664097 0.10658536

0.025 4 3.453041 0.8628558 1.3167032 0.10282201

0.025 5 3.372501 0.8686287 1.2700128 0.09837303

------

RMSE was used to select the optimal model using the smallest value.

The final values used for the model were sigma = 0.05 and C = 10.


If we wanted to take this further, we could try even more fine tuning with more grid searches. We could also explore trying to tune other parameters of the underlying ksvm() function. Finally and as already mentioned, we could perform some grid searches on the other non-linear regression methods.


Ensemble Methods

We can try some ensemble methods on the problem and see if we can get a further decrease in our RMSE. In this section we will look at some boosting and bagging techniques for decision trees. Additional approaches you could look into would be blending the predictions of multiple well performing models together, called stacking. Let's take a look at the following ensemble methods:

  • Random Forest, bagging (RF).

  • Gradient Boosting Machines boosting (GBM).

  • Cubist, boosting (CUBIST).


# Estimate accuracy of ensemble methods. 
# try ensembles 
trainControl <- trainControl(method="repeatedcv", number=10, repeats=3) 
metric <- "RMSE" 

# Random Forest 
set.seed(seed) 
fit.rf <- train(medv~., data=dataset, method="rf", metric=metric,  preProc=c("BoxCox"), trControl=trainControl) 

# Stochastic Gradient Boosting set.seed(seed) 
fit.gbm <- train(medv~., data=dataset, method="gbm", metric=metric,  preProc=c("BoxCox"), trControl=trainControl, verbose=FALSE) 

# Cubist set.seed(seed) 
fit.cubist <- train(medv~., data=dataset, method="cubist", metric=metric, preProc=c("BoxCox"), trControl=trainControl) 

# Compare algorithms 
ensembleResults <- resamples(list(RF=fit.rf, GBM=fit.gbm,  CUBIST=fit.cubist)) 
summary(ensembleResults) 
dotplot(ensembleResults)

Output:

Ensemble Methods on the Boston House Price Dataset.














Get R Programming Project help, R Programming Homework Help with an affordable price. Send your requirement details at realcode4you@gmail.com and get instant help with an affordable prices.

Comentários


bottom of page