suppressMessages(library(MASS))
data("Boston")
summary(Boston)
## crim zn indus chas
## Min. : 0.00632 Min. : 0.00 Min. : 0.46 Min. :0.00000
## 1st Qu.: 0.08204 1st Qu.: 0.00 1st Qu.: 5.19 1st Qu.:0.00000
## Median : 0.25651 Median : 0.00 Median : 9.69 Median :0.00000
## Mean : 3.61352 Mean : 11.36 Mean :11.14 Mean :0.06917
## 3rd Qu.: 3.67708 3rd Qu.: 12.50 3rd Qu.:18.10 3rd Qu.:0.00000
## Max. :88.97620 Max. :100.00 Max. :27.74 Max. :1.00000
## nox rm age dis
## Min. :0.3850 Min. :3.561 Min. : 2.90 Min. : 1.130
## 1st Qu.:0.4490 1st Qu.:5.886 1st Qu.: 45.02 1st Qu.: 2.100
## Median :0.5380 Median :6.208 Median : 77.50 Median : 3.207
## Mean :0.5547 Mean :6.285 Mean : 68.57 Mean : 3.795
## 3rd Qu.:0.6240 3rd Qu.:6.623 3rd Qu.: 94.08 3rd Qu.: 5.188
## Max. :0.8710 Max. :8.780 Max. :100.00 Max. :12.127
## rad tax ptratio black
## Min. : 1.000 Min. :187.0 Min. :12.60 Min. : 0.32
## 1st Qu.: 4.000 1st Qu.:279.0 1st Qu.:17.40 1st Qu.:375.38
## Median : 5.000 Median :330.0 Median :19.05 Median :391.44
## Mean : 9.549 Mean :408.2 Mean :18.46 Mean :356.67
## 3rd Qu.:24.000 3rd Qu.:666.0 3rd Qu.:20.20 3rd Qu.:396.23
## Max. :24.000 Max. :711.0 Max. :22.00 Max. :396.90
## lstat medv
## Min. : 1.73 Min. : 5.00
## 1st Qu.: 6.95 1st Qu.:17.02
## Median :11.36 Median :21.20
## Mean :12.65 Mean :22.53
## 3rd Qu.:16.95 3rd Qu.:25.00
## Max. :37.97 Max. :50.00
colSums(is.na(Boston))
## crim zn indus chas nox rm age dis rad
## 0 0 0 0 0 0 0 0 0
## tax ptratio black lstat medv
## 0 0 0 0 0
n = nrow(Boston)
# creating training/test sets. 70/30 split:
index = sample(1:n, n*0.70, replace = FALSE)
train = Boston[index,]
test = Boston[-index,]
# fitting MLR model for comparison.
lm.fit = lm(medv ~ .,data=train)
pred.lm = predict(lm.fit, newdata=test)
# calculating MSE of MLR:
MSE.lm = sum((test$medv - pred.lm)^2)/nrow(test)
MSE.lm
## [1] 19.52297
# scaling data based on training set: using min-max scaling
train.colmins = apply(train,2,min)
train.colmaxs = apply(train,2,max)
scaled.train = as.data.frame(scale(train,center = train.colmins, scale = train.colmaxs - train.colmins))
scaled.test = as.data.frame(scale(test,center = train.colmins, scale = train.colmaxs - train.colmins))
library(neuralnet)
##
## Attaching package: 'neuralnet'
## The following object is masked from 'package:dplyr':
##
## compute
n = names(scaled.train)
# first create formula of "outcome~predictors"
f = as.formula(paste("medv~", paste(n[!n %in% "medv"], collapse = "+")))
# fit the ANN model on training set
# 2 hidden payers = 5 & 3 neorons in each
nn53 = neuralnet(f, data=scaled.train, hidden=c(5,3), linear.output=T)
plot(nn53)
# 1 hidden layer = 5 neurons in it
nn5 = neuralnet(f, data=scaled.train, hidden=5, linear.output=T)
plot(nn5)
# prediction using ANN model
pred.nn53 = compute(nn53, scaled.test[,1:13])
# filtering prediction values only and scaling back to original scale
pred.nn53_ = pred.nn53$net.result * (max(Boston$medv) - min(Boston$medv)) + min(Boston$medv)
pred.nn5 = compute(nn5, scaled.test[,1:13])
pred.nn5_ = pred.nn5$net.result * (max(Boston$medv) - min(Boston$medv)) + min(Boston$medv)
# scaling test set outcome variable to original scale.
test.r = (scaled.test$medv) * (max(Boston$medv)-min(Boston$medv)) + min(Boston$medv)
# MSE for both models.
MSE.nn53 = sum((test.r - pred.nn53_)^2) / nrow(scaled.test)
MSE.nn5 = sum((test.r - pred.nn5_)^2) / nrow(scaled.test)
# MSE comparison:
MSE.df = data.frame(MLR = round(MSE.lm,1), NN53 = round(MSE.nn53,1), NN5 = round(MSE.nn5,1))
rownames(MSE.df) = "Test-MSE"
MSE.df
## MLR NN53 NN5
## Test-MSE 19.5 7.3 9.4
# as you can see, ANN model with mode hidden layer and neuron has better error rate and prediction power. but that comes at the cost of overfitting when new data is applied due to high flexibility.
# one data frame to have y, y-hat, type of model:
plot.data = data.frame(medv=c(test$medv, test$medv, test$medv),
pred = c(pred.lm, pred.nn5_, pred.nn53_),
model = c(rep('lm', nrow(test)), rep('ann5', nrow(test)), rep('ann53', nrow(test)))
)
# plot the y and y-hat values (original and predicted) by model type.
ggplot(plot.data, aes(x=medv, y=pred, color=model)) +
geom_point() +
geom_abline(slope=1, intercept=0) +
theme_bw()
