# code chunks: Boston dataset and Neural Networks 
# Date: 28.11.2025 



# Load and prepare data 
library("neuralnet")
library("MASS")
library("dplyr")
library("mlbench")
library(magrittr)
library(patchwork)
library(ggplot2)
library("reshape2")

set.seed("3456")

data = Boston 
head(data)

# Min max normalization - Rescaling 
summary(data)
max_data <- apply(data, 2, max)
min_data <- apply(data, 2, min)
data_scaled <- as.data.frame(scale(data, center = min_data, scale = max_data - min_data))


# Split into a simple train and test dataset 
split  <-  sample(1:nrow(data),round(0.75*nrow(data)))

train_data <- as.data.frame(data_scaled[split,])
test_data <- as.data.frame(data_scaled[-split, ])

nrow(test_data)
nrow(train_data)
summary(train_data)


# Set up the neural net 
n = names(data)
f = as.formula(paste("medv ~ ", paste(n[!n %in% "medv"], collapse = "+" 
)))

net = neuralnet(f, data = train_data, hidden = c(5,3), linear.output = T)
plot(net)
summary(net)



# predict with neural net on test data
predict_net <-  predict(net, test_data[,1:13])
predict_net_unscaled <-  as.data.frame(predict_net * (max(data$medv)- min(data$medv)) + min (data$medv))

# compute RMSEs
rmse.net <- caret::RMSE(predict_net, test_data$medv)
print(paste("RMSE neural net:", rmse.net))
test_unscaled <- as.data.frame((test_data$medv) * (max(data$medv)- min(data$medv)) + min (data$medv))
# Compare distribution of actual and predicted testing dataset 
summary(predict_net_unscaled)
summary(test_unscaled)


# Set up a regression model for comparison 
Regression_Model <- glm(medv ~ ., data = train_data)
summary(Regression_Model)

# predict with regression model 
predict_lm <- predict(Regression_Model, test_data[,1:13])
predict_lm_unscaled <-  as.data.frame(predict_lm * (max(data$medv)- min(data$medv)) + min (data$medv))
rmse.lm <- caret::RMSE(predict_lm, test_data$medv)


# Compare
print(paste("Root Mean Squared Error Linear Model:", rmse.lm))
print(paste("Root Mean Squared Error Neural Net:", rmse.net))







# Plot the distribution of predicted values and compare to actual values -----
# of test data 


# Combine the three data frames into one
df_combined <- data.frame(
  predict_lm_unscaled = predict_lm_unscaled[[1]],
  predict_net_unscaled = predict_net_unscaled[[1]],
  test_unscaled = test_unscaled[[1]]
)

# Convert to long format for ggplot2
df_long <- melt(df_combined, variable.name = "model", value.name = "prediction")

# Boxplot
ggplot(df_long, aes(x = model, y = prediction, fill = model)) +
  geom_boxplot() +
  theme_minimal() +
  labs(title = "Distribution of Values", x = "Model", y = "Predicted Value") +
  theme(legend.position = "none")


# Scatter 
ggplot(df_combined, aes(x = test_unscaled)) +
  geom_point(aes(y = predict_lm_unscaled, color = 'Linear Model'), size = 2, alpha = 0.7) +
  geom_point(aes(y = predict_net_unscaled, color = 'Neural Net 1'), size = 2, alpha = 0.7) +
  geom_abline(slope = 1, intercept = 0, linewidth = 1) +
  scale_color_manual(values = c('Linear Model' = 'blue', 
                                'Neural Net 1' = 'red')) +
  labs(title = 'real vs predicted values',
       x = 'Actual medv',
       y = 'Predicted medv',
       color = 'Model') +
  theme_minimal()