Assignment 1

STAT462 2024-S2

Braking Distance

This analysis aims to predict the braking distance of a car traveling at various speeds, utilizing a historical dataset from 1930. The dataset includes speed measurements in miles per hour and braking distances in feet.

For better legibility and relevance to a New Zealand audience, the units were converted to kilometers per hour (km/h) for speed and meters (m) for distance.

Two predictive models were employed: a simple linear regression model and a k-Nearest Neighbors (k-NN) model.

\(Q_{1}\) Data Prepration

In the data preparation steps, the dataset is loaded, the units were converted to kilometers per hour (km/h) for speed and meters (m) for distance and the dataset divided into training and test datasets with a ratio of 8:2.

#Load the dataset
df <- read.csv("braking.csv")

Unit conversion: To maintain the original dataframe, a new dataframe was created, and the previous dataset along with the converted unit versions of speed and braking distance were added to it.

#Unit conversion
# A function to convert mph to km/h
convert_mph_to_kph <- function(mph) {
  kph <- mph * 1.60934
  return(kph)
}

# A function to convert feet to meters
convert_feet_to_meters <- function(feet) {
  meters <- feet * 0.3048
  return(meters)
}

# call two columns of the original dataset (df) and convert the units and save the file into a new dataframe (df_si)
df_si <- df %>%
  mutate(sp_kph=convert_mph_to_kph(df$speed))
df_si <- df_si %>%
  mutate(dist_m=convert_feet_to_meters(df$dist))

The dataset was divided into a training set (80%) and a test set (20%) to facilitate model development and evaluation.

# making an index: for sampling we use a random 8:2 ratio of dataset which refers to different rows of the datset. 
nums <- NROW(df_si)    #for a dataset with any number of rows
set.seed(1)            # To generate a repetitive sample dataset
train_ind <- sample(1:nums, size = nums*0.8)
df_train <- df_si[train_ind, ]
df_test <- df_si[-train_ind, ]

A visualisation was also applied to get an idea of the distribution shape of the datasets.

df_train$group <- "Train"
df_test$group <- "Test"
combined_df <- rbind(df_train, df_test)
graph_breaking_1 <- ggplot(combined_df, aes(x=sp_kph, y=dist_m, color=group)) +
                      geom_point(alpha=0.3) +
                      labs(title="Braking distance vs. Speed",
                      y= "Braking distance [m]",
                      x= "Speed [kph]",
                      color="Data Source") +
                    scale_color_manual(values=c("red", "blue")) +
                    theme_minimal()


graph_breaking_1

#Saving the graph
ggsave("stopping distance [m] vs speed [kph].png", plot = graph_breaking_1 , width = 8, height = 6, dpi = 300)

\(Q_{2}\) Simple Linear Regression

Calculate \(\hat{b_0}\) and \(\hat{b_1}\) for \(Y=b_0 + b_1X\) where:

\(\hat{b_1} = \frac{cov(\underline{x},\underline{y})}{var(\underline{x})}[=corr(\underline{x},\underline{y})\frac{sd(\underline{y})}{sd(\underline{x})}]\)

\(\hat{b_0} = mean(\underline{y}) - \hat{b_1} \cdot mean (\underline{x})\)

#Calculate mean of x_mean and y_mean, Note x and y are sampled so they are x_bar and y_bar
x <- df_train$sp_kph
y <- df_train$dist_m

x_mean=mean(x)
y_mean=mean(y)


# A function to calculate b0 and b1 from any data frame
calculate_b0_b1 <- function(x,y) {
b1.hat <- cov(x,y)/var(x)
b0.hat <- y_mean-b1.hat*x_mean
return(list(b0.hat = b0.hat, b1.hat = b1.hat))
}

# calling function to calculate b0_hat and b1_hat
results <- calculate_b0_b1(x, y)

# Extract b0_hat and b1_hat from the results
b0.hat <- results$b0.hat
b1.hat <- results$b1.hat

cat("Slope (b1):", b1.hat, "\n")

## Slope (b1): 0.7643392

cat("Intercept (b0):", b0.hat, "\n")

## Intercept (b0): -5.748954

\(𝑓(𝑥)\)=-5.7489537 + 0.7643392\(x\)

\(Q_{2a}\) Simple Linear Regression Model (5km/h)

The change in speed and its impact on the model can be found by \(b_1 \times x\).

increase_speed <- 5
expected_increase_dist <- b1.hat * increase_speed
expected_increase_dist

## [1] 3.821696

Expected increase in braking distance for a 5 km/h increase in speed is 3.821696 meters

\(Q_{2b}\) Variation

To determine how much of the variation in the data can be explained by the linear regression model, the coefficient of determination needs to be calculated including \(R^{2}\).This value represents the proportion of the variance in the dependent variable that is predictable from the independent variable.

# Predict the distances using the model
df_train$predicted_dist <- b0.hat + b1.hat * x

TSS <- sum((y_mean-y)^2) #Total Sum of Squares
ESS_slr <- sum((df_train$predicted_dist - y_mean)^2)  #Explained sum of squares
Rsq_slr=ESS_slr/TSS
cat("R^2=", Rsq_slr)

## R^2= 0.5991287

An \(R^{2}\) value of 0.60 means that 60% of the variance in the dependent variable (the braking distance) can be explained by the independent variable (speed) in the above linear regression model. An \(R^{2}\) value of 0.60 indicates a moderate fit of the model to the data. It suggests that while the model does capture some of the relationship between speed and braking distance, there are other factors influencing braking distance that are not accounted for in the model. The remaining 40% of the variation in braking distance is due to other factors not included in the model,random noise, or inherent variability in the data. Additional variables or more complex models may better capture the relationship and explain more of the variance in the data.

\(Q_{2c}\) Speed as a significant predictor

residuals=y-df_train$predicted_dist
calculate_stats <- function(speed, residuals, x, x_mean, confidence_level) {
  n <- length(residuals)
  
  # Calculate the standard error of b1.hat
  se_b1.hat <- sqrt(sum(residuals^2) / (n - 2) / sum((x - x_mean)^2))
  
  # Calculate the standard error of the predicted y
  se_predicted_y <- sqrt(sum(residuals^2) / (n - 2)) * sqrt(1 + (1 / n) + ((speed - x_mean)^2 / sum((x - x_mean)^2)))
  
  # Calculate the critical t-value for the given confidence level
  alpha <- 1 - confidence_level
  t_critical <- qt(1 - alpha / 2, df = n - 2)
  
  return(list(se_b1.hat = se_b1.hat, se_predicted_y = se_predicted_y, t_critical = t_critical))
}

speed <- 30
# call function to obtain the statistical values for a specific confidence level
stats <- calculate_stats(speed, residuals, x, x_mean, confidence_level = 0.95)
print(stats)

## $se_b1.hat
## [1] 0.1014231
## 
## $se_predicted_y
## [1] 5.247293
## 
## $t_critical
## [1] 2.024394

# Calculate the t-statistic
t_stat <- b1.hat / stats$se_b1.hat

# Determine significance
is_significant <- abs(t_stat) > stats$t_critical


cat("t-statistic:", t_stat, "\n")

## t-statistic: 7.536147

cat("Critical t-value:", stats$t_critical, "\n")

## Critical t-value: 2.024394

cat("Is speed a significant predictor for distance at the 95% confidence level?", is_significant, "\n")

## Is speed a significant predictor for distance at the 95% confidence level? TRUE

Speed is a significant predictor for distance at the 95% confidence level because it does not include 0 in the interval of [0.56166 0.9670]

\(Q_{2d}\) Simple Linear Regression Model (30km/h)

#call stat function f to find the statistical information for a specific confidence interval. 
speed <- 30
stats <- calculate_stats(speed, residuals, x, x_mean, confidence_level = 0.80)
# calculate the predication interval for a specific x* and corresponding y*
predicted_dist_star <- b0.hat+b1.hat*speed
lower_bound_y.hat <- predicted_dist_star-stats$t_critical*stats$se_predicted_y
upper_bound_y.hat <- predicted_dist_star+stats$t_critical*stats$se_predicted_y

lower_bound_y.hat

## [1] 10.33754

upper_bound_y.hat

## [1] 24.0249

Predicted braking distance for 30km/h: 17.181222 meters 80% Prediction Interval: 10.3375444 to 24.0248997 meters

\(Q_{3}\) Fit a k_NN

Predict the braking distance for a car going at 30km/h using this model.

# provide speed, get the corresponded y estimation
kNN <- function(sp,k,df_train){
  neighborhood_vals <- df_train %>% 
    # generate distance between new x and each existing x
    mutate(dist = abs(df_train$sp_kph-sp)) %>% 
    # sort the distance 
    arrange(dist) %>% 
    # subset the first k rows
    slice(1:k) %>% 
    # medv is our target response
    pull(dist_m) 
  
  # take the average response as the estimate
  return (sum(neighborhood_vals)/k)
}
#Example
#print(kNN(20,21,df_train))

In order to find the best k for the kNN model, as a performance metric, MSE is calculated to find the best k.

#pred_test_knn <- sapply(df_test$sp_kph, function(speed) kNN(speed, k, df_train))

#MSE function
mse_fun <- function(y,predicted_y){
  mse_val <- mean((y - predicted_y)^2)
  return(mse_val)
}
#####################finding MSE for the range of k=1:20 (half of the data points#####
# initialise variables
k_values <- 1:40
mse_vector <- numeric(length(k_values))

# Actual y values
y <-  df_test$dist_m

#loop over k values 
for (k in k_values){
  #  Find the predicted value for df_test
  pred_test_knn <- sapply(df_test$sp_kph, FUN=kNN,k,df_train)

# find the mse_val
  mse_val <- mse_fun(y,pred_test_knn)
  mse_vector[k] <- mse_val
#  print(k)
#  print(mse_val)
#   print(pred_test_knn)
#  print(y)
}
#The best k calculated from 2 (which exclude 1 due to bizarre answer) to 20 which is half of the datapoints in df_train
best_k <- which.min(mse_vector[2:20])
print(best_k)

## [1] 7

cat("Best value for k is: ", best_k)

## Best value for k is:  7

cat("Predicted braking distance for a car going at 30km/h using kNN model (k=7) is: ",kNN(30,best_k,df_train), "m")

## Predicted braking distance for a car going at 30km/h using kNN model (k=7) is:  17.76549 m

#

#print(kNN(30,21,df_train))

Based on the MSE values calculated during cross-validation, the best k value to minimize the MSE is 7. When this k value is applied in the trained k-NN model, the predicted braking distance for a car traveling at 30 kph is approximately 17.76 meters. However, this prediction does not align optimally with the actual averaged braking distance of 15.25 meters for a car traveling at 30 kph in the original dataset.

Interestingly, using k=20 provides a prediction that closely matches the actual averaged braking distance for 30 kph. This discrepancy suggests that while k=7 minimizes the MSE across the dataset, a larger k might yield more accurate predictions for certain specific speeds, such as 30 kph.

As a generalisable model that works beyond specific data is prefered. Hence, k=7 was chosen for the kNN model.

\(Q_4\) Visualisation and Performance

# Actual y values
y <-  df_test$dist_m

  #  Find the predicted value for df_test using both models
pred_test_knn <- sapply(df_test$sp_kph, FUN=kNN,best_k,df_train)
pred_test_slr <- b0.hat + b1.hat * df_test$sp_kph


# Combining results for graph
df_test_results <- df_test
df_test_results$pred_knn <- pred_test_knn
df_test_results$pred_slr <- pred_test_slr

model_plot <- ggplot(df_test_results, aes(x = sp_kph)) +
  geom_point(aes(y = dist_m, color = "Actual"), size = 3) +  # Plot actual values
  geom_line(aes(y = pred_knn, color = "k-NN"), size = 1) +
  geom_point(aes(y = pred_knn, color = "k-NN"), size = 3) + 
  geom_line(aes(y = pred_slr, color = "Linear Regression"), size = 1) +
  geom_point(aes(y = pred_slr, color = "Linear Regression"), size = 3) + # Plot Linear Regression predictions
  labs(title = "Comparison of k-NN and Linear Regression Predictions",
       x = "Speed [kph]",
       y = "Distance [m]") +
  scale_color_manual(name = "Model Fit", values = c("Actual" = "black", "k-NN" = "red", "Linear Regression" = "blue")) +
  theme_minimal()

## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

model_plot

#Saving the graph
ggsave("stopping distance [m] vs speed [kph], Models.png", plot = model_plot , width = 8, height = 6, dpi = 300)

# find the mse_val
  mse_knn <- mse_fun(y,pred_test_knn)
  mse_slr <- mse_fun(y,pred_test_slr)
  best_model <- ifelse(mse_knn < mse_slr, "k-NN", "Linear Regression")

  
best_model  

## [1] "Linear Regression"

The best model is Linear Regression. However, both models’ behave approximately the same around the middle of data points.

Filipino Household Income

This project is about predicting household income from the number of underage children living in the household. The first two steps including loading, cleaning up and preparing training and test datasets (8:1).

\(Q_{1}\) Data Clean-up

For Question 1, the dataset is uploaded. A new dataframe with the two relevant columns was made, and these columns were relabeled to “income” and “children”. This was chosen over deleting the columns from the read dataset so all original data was still accessible.

# Load the dataset
income_csv <- read.csv("income.csv")

# New dataset with relevant columns only and columns relabeled to children and income
income_df <- income_csv[c("Total.Household.Income", "Members.with.age.5...17.years.old")]
colnames(income_df) <- c("income", "children")

\(Q_{2}\) Training & Test Datasets

The data was then divided into training and test datasets with a ratio of 8:2. set.seed was used to ensure that the same training and test data was used and make the results reproducible.

# Calculate number or rows
n_total <- nrow(income_df)

# Split data set into train and test, set seed so data is reproducible
set.seed(1)
train_ind <- sample(1:n_total, size = n_total*0.8)
income_train <- income_df[train_ind, ]
income_test <- income_df[-train_ind, ]

#set n for training sample
n <- nrow(income_train)

\(Q_{3}\)

\(Q_{3a}\) Linear Regression Model

The lm function in R was used to solve for \(y=b_0 + b_1x\) and find \(\hat{b_0}\) and \(\hat{b_1}\), where y = income and x= children. \(\hat{b_0}\) and \(\hat{b_1}\) were found by indexing the resulting coefficients from the result of the lm on the training sample.

# Calculate b0 and b1 using lm, then indexing to lm coefficients

linear_model = lm(income ~ children, data = income_train)

fhi_b0.hat = linear_model$coefficients[1]
fhi_b1.hat = linear_model$coefficients[2]

# Print results
cat("Slope (b1):", fhi_b1.hat, "\n")

## Slope (b1): -10125.03

cat("Intercept (b0):", fhi_b0.hat, "\n")

## Intercept (b0): 261025.6

# Format for formula display
b1_round <- format(fhi_b1.hat, scientific = FALSE)
b0_round <- format(fhi_b0.hat, scientific = FALSE)

b0_round <- round(fhi_b0.hat, 2)
b1_round <- round(fhi_b1.hat, 2)

# Check the sign of the coefficients
b0_display <- ifelse(b0_round >= 0, paste0(b0_round), paste0("-", abs(b0_round)))
b1_display <- ifelse(b1_round >= 0, paste0(b1_round), paste0("-", abs(b1_round)))

Following the results of the lm function the affine-linear model for the dataset is as follows:

income =261025.56 -10125.03 \(\cdot\) children

\(Q_{3b}\) Predicted Mean Income of a Household with n Children

To predict the mean income of a household using the model from Question 3a for n children, a sequence of numbers from 0 to 8 was created (n_children). A function was then used to calculate the prediction interval. To enable this function to be re-used for Question 4, \(\hat{b_0}\) and \(\hat{b_1}\) were calculated again within the function. The parameters were an x variable, y variable, alpha (\(\alpha\) was set to 0.10 to calculate 90% prediction intervals) and an explanatory variable (which was shortened to x_star). The function output is a dataframe with four columns: Number of Children, Lower Predicted Income, Predicted Income and Upper Predicted Income.

# Prediction Function
prediction_intervals <- function(x, y, alpha=0.10, x_star){
  result <- data.frame('Number of Children' = numeric(), 'Lower Predicted Income' = numeric() , 'Predicted Income' = numeric(), 'Upper Predicted Income' = numeric())
  linear_model = lm( y ~ x)
  b0 = linear_model$coefficients[1]
  b1 = linear_model$coefficients[2]
  TSS <- sum((y - mean(y))^2)
  ESS <- sum((b0 + b1 * x - mean(y))^2)
  n <- length(x)
  t_critical <- qt(1-alpha/2, n-1)
  y_hat <- b1 * x_star + b0
  tau <- t_critical * sqrt((TSS - ESS) / (n - 2)) * sqrt(1 + 1 / n + (x_star - mean(y))^2 / sum((y - mean(y))^2)) 
  y_hat_lower <- y_hat - tau
  y_hat_upper <- y_hat + tau
  
  result <- data.frame('Number of Children' = numeric(), 'Lower Predicted Income' = numeric() , 'Predicted Income' = numeric(), 'Upper Predicted Income' = numeric())
  result <- data.frame('Number of Children' = x_star, 'Lower Predicted Income' = y_hat_lower, 'Predicted Income' = y_hat, 'Upper Predicted Income' = y_hat_upper) 
  
  return(result)
}

# Create n number of children
n_children <- seq(0, 8)

# Run Function
y_star <- prediction_intervals(x=income_train$children , y=income_train$income , alpha=0.10, x_star= n_children)

# Return table
knitr::kable(y_star, caption = '$Q_{3.b}$ Prediction Intervals for Mean Income with `n` Children', col.names = c('Number of Children', 'Lower Predicted','Predicted Income','Upper Predicted'), align = "cccc")

\(Q_{3.b}\) Prediction Intervals for Mean Income with n Children

Number of Children	Lower Predicted	Predicted Income	Upper Predicted
0	-207716.7	261025.6	729767.9
1	-217841.8	250900.5	719642.8
2	-227966.8	240775.5	709517.8
3	-238091.8	230650.5	699392.8
4	-248216.9	220525.4	689267.7
5	-258341.9	210400.4	679142.7
6	-268466.9	200275.4	669017.7
7	-278591.9	190150.4	658892.6
8	-288717.0	180025.3	648767.6

\(Q_{3c}\) Percentage of Test Data Points Within 90% Prediction Intervals

To check predictions of the test data against the train, the prediction interval function was run again. Then a Boolean statement comparing the test income to the predicted incomes was made. The result was TRUE if: upper predicted <= test income >= lower predicted income. The number of TRUE statements was then divided by the total number of rows in the test data to gain a percentage.

# Running the prediction interval function
pred_test <- prediction_intervals(x=income_train$children , y=income_train$income , alpha=0.10, x_star=income_test$children)

# Checking for accuracy
accuracy <- sum(income_test$income <= pred_test$Upper.Predicted.Income & income_test$income >= pred_test$Lower.Predicted.Income)
true_datapoint <- (accuracy / (length(income_test$income)) * 100)

true_datapoint

## [1] 95.31833

Percent of datapoints within the 90% prediction intervals for the test data: 95.3183295%.

\(Q_{4}\) Linear Regression on log(income)

Steps 3a to 3c were run again using the log(income).

# Running the prediction interval function on n-children
log_prediction <- prediction_intervals(x=income_train$children , y=log(income_train$income) , alpha=0.10, x_star= n_children)

# Create table
knitr::kable(log_prediction, caption = 'Prediction Intervals for the Log of Mean Income with `n` Children', col.names = c('Number of Children', 'Lower Predicted','Predicted Income','Upper Predicted'), align = "cccc")

Prediction Intervals for the Log of Mean Income with n Children

Number of Children	Lower Predicted	Predicted Income	Upper Predicted
0	10.81236	12.08374	13.35512
1	10.81554	12.08617	13.35681
2	10.81865	12.08861	13.35857
3	10.82170	12.09104	13.36039
4	10.82469	12.09348	13.36227
5	10.82761	12.09592	13.36422
6	10.83047	12.09835	13.36623
7	10.83326	12.10079	13.36831
8	10.83599	12.10322	13.37045

# Running the predicition interval function
pred_test_log <- prediction_intervals(x=income_train$children , y=log(income_train$income) , alpha=0.10, x_star=income_test$children)

# Checking for accuracy
accuracy_log <- sum(log(income_test$income) <= pred_test$Upper.Predicted.Income & income_test$income >= pred_test$Lower.Predicted.Income)
true_datapoint_log <- (accuracy_log / (length(income_test$income)) * 100)

true_datapoint_log

## [1] 100

Percent of datapoints within the 90% prediction intervals for the test data: 100%

\(Q_{5}\) Summary of Analysing “Filipino Household Income”

The Prediction Intervals were very large for the model and included negative intervals. This could be argued to be nonsensical when dealing with mean income. A simple scatter plot for the train data shows that there appears to be a significant right or positive skew in the data indicating that income is not a Gaussian variable. This skew would impact the prediction intervals. To counteract this the Log(income) can be used because the mean incomes do not contain negative values. How the Log(income) impacts the distribution can be shown through density plots as it creates a more normal or Gaussain distribution for the income.

income_plot <- ggplot(income_train, aes(x = children, y = income)) +
  geom_point(size = 1) +
  labs(title = "Scatter Plot showing Income Train Data",
       x = "Number of Children",
       y = "Mean Income")
   
income_plot + scale_y_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE))

ggsave("Scatter Plot showing Income Train Data.png", plot = income_plot , width = 8, height = 6, dpi = 300)

mean_income <- mean(income_train$income)

density_plot <- ggplot(income_train, aes(x = income,)) +
  geom_density(size = 1) +
  geom_vline(xintercept = mean_income, color = "red", size = 1) +
  labs(title = "Density Plot of Income showing Mean Income",
       x = "Income",
       y = "Density") +
  scale_x_continuous(labels = function(x) format(x, big.mark = ",", scientific = FALSE)) +
  theme(axis.text.y = element_blank())  

density_plot_log <- ggplot(income_train, aes(x = log(income))) +
  geom_density(size = 1) +
  geom_vline(xintercept = log(mean_income), color = "red", size = 1) +
  labs(title = "Density Plot of Log(Income) showing Log(Mean Income)",
       x = "Log(Income)",
       y = "Density")+
  theme(axis.text.y = element_blank())   


print(density_plot)

ggsave("Density Plot of Income showing Mean Income.png", plot = density_plot , width = 8, height = 6, dpi = 300)
print(density_plot_log)

ggsave("Density Plot of Log(Income) showing Log(Mean Income).png", plot = density_plot_log, width = 8, height = 6, dpi = 300)

Predicting Possum Age

In this study, we analyse a dataset containing measurements of possums trapped across seven different sites in Australia, specifically within the regions of Queensland, New South Wales, and Victoria. The primary objective of this analysis is to develop a predictive model that can estimate the age of a possum based on various body measurements. By leveraging these data, we aim to uncover relationships between physical characteristics and age, providing insights that could be valuable for ecological research and wildlife management.

\(Q_{1}\) Upload the Possums dataset

Load the dataset and create a visualistion for age and total length of possums.

# Load the dataset
possum_csv <- read.csv("possums.csv")

# Define age groups
#compute correlation between age and total length
correlation <- cor(possum_csv$age, possum_csv$totlngth, use = "complete.obs")
correlation

## [1] 0.2602803

# calculate the age group vs verged total length 
###########################
#NOTE: this is an extra step which is used for visualisation at this section, other than that for the rest of questions, the age variable is considered numerical variable and the original data points are used for the further analysis and modelling) 
##########################
possums <- possum_csv
possums$age_group <- cut(possums$age, 
                         breaks = c(0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10), # Adjust these breaks as per your data
                         labels = c("0-1", "1-2", "2-3", "3-4", "4-5", "5-6","6-7" ,"7-8", "8-9", "9-10"),
                         include.lowest = TRUE)

# Calculate the mean total length and mean age for each group
age_totlgth_avg <- possums %>%
  group_by(age_group) %>%
  summarise(mean_totlgth = mean(totlngth, na.rm = TRUE),
            mean_age = mean(age, na.rm = TRUE))

age_length<- ggplot()+
  geom_point(data=possums, aes(x=totlngth, y=age, color = "Original Data"), alpha = 0.6, size = 2) +   # Original data points
  geom_point(data = age_totlgth_avg, aes(y = mean_age, x = mean_totlgth, 
             color = "Averaged Data"), size = 3) +                # Averaged data points
stat_function(fun = function(x) 6e-10 * exp(0.2589 * x), color = "black", size = 1, aes(linetype = "Exponential Fit")) +  # Add the exponential curve
  labs(title = "Age vs total length",             # Exponential fit to averaged data
         y = "Age", x = "Total Length (totlngth)") +
  ylim(0,10)+
  theme_minimal()
age_length

# Save plot
ggsave("Age vs Total Possum Length.png", plot = age_length  , width = 8, height = 6, dpi = 300)

By examining the data points and searching for trends, or specific patterns, it appears that there is a slight overall increase in total length with age if the original data is examined and age is considered as a numerical trend. This observation is supported by the computed correlation value of 0.26. If the age variable is treated as a categorical variable, and the averaged data points is calculated at each age group (1 year), an exponential trend can be fitted to the dataset. The trend between age and some other numerical variables such as belly and chest have the same patterns with different exponential coefficients.

This visualization encouraged us to explore the dataset and the relationships between its variables in a more systematic way through the following steps:

1) Check for Missing Values: In the possums file, there are two NA cells for age and 1 for foot length. The corresponding rows will be removed from dataset.

2)Review the Statistical Summary: statistical summary of all variables were reviewed. The mean, upper and lower limit of each variable revealed that the collected measurements were not recorded with the same unit. For example, chest values are in cm whereas skullw (skull width) and hdlngth (head length) are in mm. Hence a unit conversion will be also performed.

3) visualising age against a couple of variables suggested a weak correlation between the age and body part measurement. To quantify our observation, we generated a “Correlation Matrix” to possibly evaluate the correlation between age and other variables and illustrate our findings. Result of this step and corresponding interpretation will be presented after one-hot encoding and final step of data prepration.

Following code is written to perform a primary of data clean up, check for missing values, obtain a summary statistics, and convert units.

# ##################### Data Prepration #################
#make a new dataframe to remain the original data intact
possum_df <- possum_csv
#remove the unnecessary fields (case & Pop)
possum_df <- possum_df[, !names(possum_df) %in% c("case", "Pop")]
#removing rows with NA field
possum_df <- na.omit(possum_df)

# 1. Check for missing values
summary(possum_df)

##       site           sex                 age           hdlngth      
##  Min.   :1.000   Length:101         Min.   :1.000   Min.   : 82.50  
##  1st Qu.:1.000   Class :character   1st Qu.:2.000   1st Qu.: 90.70  
##  Median :4.000   Mode  :character   Median :3.000   Median : 92.90  
##  Mean   :3.673                      Mean   :3.822   Mean   : 92.73  
##  3rd Qu.:6.000                      3rd Qu.:5.000   3rd Qu.: 94.80  
##  Max.   :7.000                      Max.   :9.000   Max.   :103.10  
##      skullw         totlngth         taill          footlgth       earconch    
##  Min.   :50.00   Min.   :75.00   Min.   :32.00   Min.   :60.3   Min.   :41.30  
##  1st Qu.:55.00   1st Qu.:84.50   1st Qu.:36.00   1st Qu.:64.5   1st Qu.:44.80  
##  Median :56.40   Median :88.00   Median :37.00   Median :67.9   Median :46.80  
##  Mean   :56.96   Mean   :87.27   Mean   :37.05   Mean   :68.4   Mean   :48.13  
##  3rd Qu.:58.10   3rd Qu.:90.00   3rd Qu.:38.00   3rd Qu.:72.5   3rd Qu.:52.00  
##  Max.   :68.60   Max.   :96.50   Max.   :43.00   Max.   :77.9   Max.   :56.20  
##       eye            chest           belly      
##  Min.   :12.80   Min.   :22.00   Min.   :25.00  
##  1st Qu.:14.40   1st Qu.:25.50   1st Qu.:31.00  
##  Median :14.90   Median :27.00   Median :32.50  
##  Mean   :15.05   Mean   :27.06   Mean   :32.64  
##  3rd Qu.:15.70   3rd Qu.:28.00   3rd Qu.:34.00  
##  Max.   :17.80   Max.   :32.00   Max.   :40.00

# 2. Summary statistics
summary_stats <- possum_df %>% select_if(is.numeric) %>% summary()
print(summary_stats)

##       site            age           hdlngth           skullw     
##  Min.   :1.000   Min.   :1.000   Min.   : 82.50   Min.   :50.00  
##  1st Qu.:1.000   1st Qu.:2.000   1st Qu.: 90.70   1st Qu.:55.00  
##  Median :4.000   Median :3.000   Median : 92.90   Median :56.40  
##  Mean   :3.673   Mean   :3.822   Mean   : 92.73   Mean   :56.96  
##  3rd Qu.:6.000   3rd Qu.:5.000   3rd Qu.: 94.80   3rd Qu.:58.10  
##  Max.   :7.000   Max.   :9.000   Max.   :103.10   Max.   :68.60  
##     totlngth         taill          footlgth       earconch          eye       
##  Min.   :75.00   Min.   :32.00   Min.   :60.3   Min.   :41.30   Min.   :12.80  
##  1st Qu.:84.50   1st Qu.:36.00   1st Qu.:64.5   1st Qu.:44.80   1st Qu.:14.40  
##  Median :88.00   Median :37.00   Median :67.9   Median :46.80   Median :14.90  
##  Mean   :87.27   Mean   :37.05   Mean   :68.4   Mean   :48.13   Mean   :15.05  
##  3rd Qu.:90.00   3rd Qu.:38.00   3rd Qu.:72.5   3rd Qu.:52.00   3rd Qu.:15.70  
##  Max.   :96.50   Max.   :43.00   Max.   :77.9   Max.   :56.20   Max.   :17.80  
##      chest           belly      
##  Min.   :22.00   Min.   :25.00  
##  1st Qu.:25.50   1st Qu.:31.00  
##  Median :27.00   Median :32.50  
##  Mean   :27.06   Mean   :32.64  
##  3rd Qu.:28.00   3rd Qu.:34.00  
##  Max.   :32.00   Max.   :40.00

#3. Unit conversion
possum_df$hdlngth <- possum_df$hdlngth/10
possum_df$skullw <- possum_df$skullw/10
possum_df$footlgth <- possum_df$footlgth/10
possum_df$earconch <- possum_df$earconch/10
possum_df$eye <- possum_df$eye/10

As presented in the statistical summary of the possums dataset, the central point of the age data tends toward the minimum values. In other words, the age distribution of possums is right-tailed, meaning it is positively skewed. In this distribution, most possums are relatively young, with a smaller number of older possums creating a long tail on the right side. Compared to other numerical variables, the age data has a relatively higher mean than median, further indicating positive skewness, which will be observed in the density map of the age in the following graph. This may suggest applying a logarithmic transformation in the MLR, which will be investigated later in \(Q_{7}\).

\(Q_{2b,c}\) One-Hot Encoding

In one-hot encoding the categorical fields were changed to factors so R would not read them as numerical classes. Then using model.matrix new columns for each site (and sex) were created. If the possum was present at that site, site$n$ = 1. If the possum was not at that site, site$n$ = 0. This was also done for sex, where if sexf = 1 the possum was female. Unnecessary columns were then removed.

A column for the Site 7 was also removed because if all other sites were 0, then this would indicate that the possum was located at site 7. As such this field was redundant.

## Check the columns are stored as Factors for Categorical Data
possum_df$site <- as.factor(possum_df$site)
possum_df$sex <- as.factor(possum_df$sex)

## Perform one-hot encoding on the 'site' column
sites_code <- model.matrix(~ site - 1, data = possum_df)

## Combine
possum_df <- cbind(possum_df, sites_code)

## Perform one-hot encoding on the 'sex' column
sex_codes <- model.matrix(~ sex - 1, data = possum_df)

## Combine
possum_df <- cbind(possum_df, sex_codes)

## Remove columns
possum_df <- possum_df[, !names(possum_df) %in% c("site7", "sexm", "site", "sex")]

As mentioned above, a “Correlation Matrix” is computed to evaluate the correlation between age and variables and illustrate our findings.

# Data containing numeric columns
numeric_cols <- possum_df %>% select_if(is.numeric)  
# Visualizations
ggpairs(numeric_cols)

pairwise_plot <- ggpairs(
  numeric_cols,     
  title = "Pairwise Plot Matrix of Numeric Variables", 
  upper = list(continuous = wrap("cor", size = 4)),     
  lower = list(continuous = "smooth"),                 
  diag = list(continuous = "densityDiag")             
)
print(pairwise_plot)

Observations

*The density plots of the numerical variables printed in diagonal cells of the pairwise plot show that not a few variables follow the normal distribution such as head length, eye and belly.

*There are some variables with distinct bimodal peaks (modes), such as ear conch and foot length, and others with a merged bump, such as tail and total length. This may indicate that the data is composed of two different subgroups or populations. Since the possum dataset includes both genders, the impact of sexual dimorphism—where males and females have different physical characteristics—could be reflected in variables like ear conch, as discussed in \(Q_{7}\).

Bimodal distributions could also result from geographic locations, which is plausible given that only a few data points were collected from certain sites. Additionally, different genetic variations might be another reason for the bimodality in the data, although this cannot be explored further here, as no such information is provided in this dataset.

# Correlation matrix for all variables
cor_matrix <- cor(numeric_cols, use = "complete.obs")
corrplot(cor_matrix, method = "circle")

As seen on the circle correlation graph, there is a strong correlation between some body parts measurements such as However, between age and other numeric variables, there is a maximum correlation of 0.35 with belly. The chest, head length and total length have slightly less correlation. The least correlation is between age and earconch (0.053), implying there are independent variables. Some variables have negative correlation values, indicating an inverse relationship between them such as relationship with site and foot length.**

\(Q_{2d}\) Split dataset into Training, Validation & Test Datasets (8:1:1)

To split dataset a function is defined as the impact of generated seed number on the predicted model will be evaluated later.

n <- length(possum_df$age)
#Define a function for changing the seed number
idx_gen <- function(seedNum){
set.seed(seedNum)
train_ind <- sample(1:n, size = n*0.8)
not_train_ind <- setdiff(1:n, train_ind)
valid_ind <- sample(not_train_ind, size = (n*0.2 / 2))
return(list(train_ind=train_ind,valid_ind=valid_ind))
}
indices <- idx_gen(seedNum=1)
train_ind <- indices$train_ind
valid_ind <- indices$valid_ind
df_train <- possum_df[train_ind, ]
df_valid <- possum_df[valid_ind, ]
df_test <- possum_df[-c(train_ind, valid_ind), ]

\(Q_{3}\) The Number of Possible Linear Regression models

With 16 columns of data, 2^16 represents the total number of combinations. There are 65,536 possible multiple linear regression models.

# Number of predictors (columns excluding 'age')
n_predictors <- ncol(df_train) - 1

# Number of possible models
n_models <- 2^n_predictors
n_models

## [1] 65536

\(Q_{4}\) Step-wise Forward Feature Selection

To perform stepwise forward feature selection, the approach used in lab 4 was followed. Where the results for the best MSE result and associated feature were recorded for the training data and for the validation. These were set up to first be empty vectors based on the number of variables plus a NULL set.

Step #1: Initialisation

# upload the name of the features into a vector
features <- colnames(df_train)[colnames(df_train) != "age"]

# number of features in the model, start with 0, the null model (#counter)
p <- 0
# vector to store the best features in order, 1 is assigned to the null model.
best_features <- vector(length=length(features)+1)
#vector to store the MSEs of the best models in each complexity on training set
MSE_train_best <- vector(length=length(features)+1)

# vector to store the MSEs of the best models in each complexity on validation set
MSE_valid_best <- vector(length = length(features) +1)

The MSE was first calculated for the NULL model, where it is predicted that none of the available variables will impact \(Y\) (age). This is then used as a baseline for comparison of later models where variables are included.

Step #2: Null model

null_fun <- function(df_train){
# for the null model, we set feature as 1
predictor <- c("1")

# concatenate the response with the predictor to generate a formula
formula <- paste("age ~", predictor, sep = " ")

# use lm() to fit the model
fit_models <- lm(formula = formula, data = df_train)

# find the MSE
MSEs_train <- mean((fit_models$residuals)^2)

# store the minimum MSE, only one value at this stage though
MSE_train_best[1] <- min(MSEs_train)


# find the MSE on validation set
MSE_valid_best[1] <- mean((predict(fit_models, newdata = df_valid) - df_valid$age)^2)
return(list(MSE_train_best=MSE_train_best[1],MSE_valid_best=MSE_valid_best[1]))
}
null_fun(df_train)

## $MSE_train_best
## [1] 3.902344
## 
## $MSE_valid_best
## [1] 1.597656

#cat("Null Model, MSE training:",MSE_train_best, "\n")
#cat("Null Model, MSE valid:",MSE_valid_best, "\n")

The MSE values for the training and validation datasets in the first run of the null model (using set.seed(1)) were 3.9 and 1.59, respectively. This gap between the training and validation MSE could indicate a potential issue with data splitting or model simplicity. By changing the random number generator seed, the discrepancy between the MSE values of the training and validation datasets significantly decreased, as reported below.

The recent reported values suggests that, on average, the squared difference between the actual ages and the mean age of the training and validation datasets are 3.699375 and 2.380625, respectively. These values of the null model provide a baseline for model performance. Any other predictive model that we obtain later should aim to have a lower MSE than these baseline values. If the multiple linear regression model has an MSE significantly lower than these null model MSEs, it indicates that our model is doing a better job at predicting age based on the features than just predicting the mean age. However, the discrepancy between MSEs computed from other models with additional feature are not significantly lower than the MSEs of null model. This could be due to using additional insignificant features, over fitting, or improper use of model.

Step #3: Testing multiple seeds

By changing the set.seed value, we can check if the discrepancy between training and validation MSE is specific to the initial data split or if it’s a general issue. Consistency across different seeds would suggest that the model and data split are more robust, while persistent discrepancies might indicate the need for further investigation or model adjustments.

#
seedNum <- 2
indices <- idx_gen(seedNum)
train_ind <- indices$train_ind
valid_ind <- indices$valid_ind
df_train <- possum_df[train_ind, ]
df_valid <- possum_df[valid_ind, ]
df_test <- possum_df[-c(train_ind, valid_ind), ]

MSE_null <- null_fun(df_train)

cat("Null Model, seed #:", seedNum, "\n")

## Null Model, seed #: 2

cat("MSE training:", MSE_null$MSE_train_best, "\n")

## MSE training: 3.431094

cat("MSE valid:", MSE_null$MSE_valid_best, "\n")

## MSE valid: 3.903906

#####################################################
seedNum <- 100
indices <- idx_gen(seedNum)
train_ind <- indices$train_ind
valid_ind <- indices$valid_ind
df_train <- possum_df[train_ind, ]
df_valid <- possum_df[valid_ind, ]
df_test <- possum_df[-c(train_ind, valid_ind), ]

MSE_null <- null_fun(df_train)

cat("Null Model, seed #:", seedNum, "\n")

## Null Model, seed #: 100

cat("MSE training:", MSE_null$MSE_train_best, "\n")

## MSE training: 3.699375

cat("MSE valid:", MSE_null$MSE_valid_best, "\n")

## MSE valid: 2.380625

#####################################################
seedNum <- 500
indices <- idx_gen(seedNum)
train_ind <- indices$train_ind
valid_ind <- indices$valid_ind
df_train <- possum_df[train_ind, ]
df_valid <- possum_df[valid_ind, ]
df_test <- possum_df[-c(train_ind, valid_ind), ]

MSE_null <- null_fun(df_train)

cat("Null Model, seed #:", seedNum, "\n")

## Null Model, seed #: 500

cat("MSE training:", MSE_null$MSE_train_best, "\n")

## MSE training: 3.494375

cat("MSE valid:", MSE_null$MSE_valid_best, "\n")

## MSE valid: 3.495625

By applying three seed number, it is obvious that the observed discrepancy was due to the initial data split that is very often when the number of datapoints is low. Seed number is chosen 2 for the rest of the model computation process.

seedNum <- 2
indices <- idx_gen(seedNum)
train_ind <- indices$train_ind
valid_ind <- indices$valid_ind
df_train <- possum_df[train_ind, ]
df_valid <- possum_df[valid_ind, ]
df_test <- possum_df[-c(train_ind, valid_ind), ]

MSE_null <- null_fun(df_train)

Step #4: One feature

After the Null mode, a linear regression model with one variable is chosen. Again the ‘best’ variable is chosen by calculating MSE. There are sixteen possible formulas, with sixteen variables present in the dataset. lappy() and similar functions were used to loop through and find the best MSE result and and therefore most suitable feature.

# number of features in the model
p <- 1

# use combn() to generate all combinations of features, take p at a time, return the result as a list
predictor <- combn(features, p, simplify = FALSE)

# generate formula respectively
formula <- sapply(predictor, function(x) paste("age", "~", paste0(x, collapse = "+"), sep = " "))

# fit model respectively
fit_models <- lapply(formula, function(x) lm(x, data = df_train))

# get MSE respectively
MSEs_train <- sapply(1:length(fit_models), function(x) mean(fit_models[[x]]$residuals^2))
#MSEs_train

# the feature results the smallest MSE will be stored
best_features[1 + p] <-predictor[[which.min(MSEs_train)]]

# stored the smallest MSE
MSE_train_best[1 + p] <- min(MSEs_train)  

# find out the best performed feature and store its MSE on validation set
MSE_valid_best[1 + p] <- mean((predict(fit_models[[which.min(MSEs_train)]], newdata = df_valid) - df_valid$age)^2)


cat("MSE training:",MSE_train_best[1 + p], "for", best_features[1 + p], "\n")

## MSE training: 2.776551 for belly

cat("MSE valid:",MSE_valid_best[1 + p],"for", best_features[1 + p], "\n")

## MSE valid: 4.677446 for belly

After the ‘best’ feature is chosen, this process is then repeated. A new variable is added each time and trained and validated on how it interacts with the previously chosen ‘best’ variable. Again this is measure by the MSE. This can result in a variable that independently would not create a good model but combined with other variables creates the best result. This process can be completed from the NULL model to the final available model using a similar method as above.

Step #5: Two Features and More

# Initialisation
features <- colnames(df_train)[colnames(df_train) != "age"]
best_features <- vector(length = length(features) +1)
MSE_train_best <- vector(length = length(features) +1)
MSE_valid_best <- vector(length = length(features) +1)
predictor <- c("1")

# null model
# concatenate the response with the predictor to generate the formula
formula <- paste("age ~", predictor, sep = " ")
# use lm() to find the model
fit_models <- lm(formula = formula, data = df_train)
# find the MSE
MSEs_train <- mean((fit_models$residuals)^2)
# store the minimum MSE, only one value at stage though
MSE_train_best[1] <- min(MSEs_train)
# find the MSE on validation set
MSE_valid_best[1] <- mean((predict(fit_models, newdata = df_valid) - df_valid$age)^2)
best_features[1] <- predictor

# with features
for (i in 1 :length(features)) {
  p = i
  if(p == 1){
    predictor <- combn(features, p, simplify = FALSE)
  } else{
    predictor <- Filter(function(x) all(best_features[2:p]  %in% x), combn(features, p, simplify = FALSE))
  }
  
  formula <- sapply(predictor, function(x) paste("age", "~", paste0(x, collapse = "+"), sep = " "))
  fit_models <- lapply(formula, function(x) lm(x, data = df_train))
  MSEs_train <- sapply(1:length(fit_models), function(x) mean(fit_models[[x]]$residuals^2))
  best_features[p+1] <- setdiff(predictor[[which.min(MSEs_train)]],best_features[2:p] )
  MSE_train_best[1 + p] <- min(MSEs_train)  
  MSE_valid_best[1 + p] <- mean((predict(fit_models[[which.min(MSEs_train)]], newdata = df_valid) - df_valid$age)^2)
}

Once we have the MSE on validation set from the best model in each complexity, we can pick the single best model among them.

\(Q_{5}\) Summary of best feature selection

df_MSE <- data.frame(feature=ordered(best_features, levels = unique(best_features)), MSE_training = MSE_train_best, MSE_valid = MSE_valid_best)

knitr::kable(df_MSE, col.names = c('Feature', 'MSE Training','MSE Validation'), align = "ccc")

Feature	MSE Training	MSE Validation
1	3.431094	3.903906
belly	2.776551	4.677446
hdlngth	2.664768	4.734429
site4	2.443909	3.777170
site3	2.396543	4.283023
footlgth	2.320720	4.172309
site2	2.282702	4.235671
site1	2.234615	3.473651
totlngth	2.181925	3.235183
eye	2.141525	3.317267
taill	2.120344	3.530388
skullw	2.107239	3.405970
site5	2.101705	3.268302
chest	2.099995	3.367440
earconch	2.098247	3.413367
site6	2.097648	3.441652
sexf	2.097403	3.442616

The MSEs of the model when applied to the training and validation datasets are visualised below:

MSE_results <- ggplot(data=df_MSE, aes(x=feature)) +
  geom_point(aes(y=MSE_training, color="MSE_training")) +
  geom_line(aes(y=MSE_training, color="MSE_training", group=1)) +
  geom_point(aes(y=MSE_valid, color="MSE_validation")) + 
  geom_line(aes(y=MSE_valid, color="MSE_validation", group=1)) +
  geom_point(aes(y=min(MSE_valid), x=feature[which.min(MSE_valid)]), size = 5, color = "black", shape = 13) + 
  labs(x="Features Added in Order", y="MSE")+
  theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1))

MSE_results

ggsave("Multiple Linear Regression Results.png", plot = MSE_results, width = 8, height = 6, dpi = 300)

# Comparing Models with different selected features
model1 <- lm(age ~ belly + hdlngth + site4 + site3+ footlgth + site2 + site1 +totlngth+eye, data = df_train)
summary(model1)

## 
## Call:
## lm(formula = age ~ belly + hdlngth + site4 + site3 + footlgth + 
##     site2 + site1 + totlngth + eye, data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.8597 -1.0284 -0.1118  0.7144  4.6743 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)   
## (Intercept) -20.52607    6.52939  -3.144  0.00245 **
## belly         0.21515    0.07790   2.762  0.00733 **
## hdlngth       2.35976    0.80638   2.926  0.00462 **
## site4        -2.50019    0.86042  -2.906  0.00490 **
## site3        -1.45717    0.88686  -1.643  0.10485   
## footlgth     -2.29175    0.94921  -2.414  0.01838 * 
## site2         1.95011    0.89591   2.177  0.03288 * 
## site1         1.07144    0.77393   1.384  0.17063   
## totlngth      0.08791    0.06898   1.274  0.20676   
## eye           2.05725    1.79022   1.149  0.25440   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.564 on 70 degrees of freedom
## Multiple R-squared:  0.3758, Adjusted R-squared:  0.2956 
## F-statistic: 4.684 on 9 and 70 DF,  p-value: 7.123e-05

model2 <- lm(age ~ belly + hdlngth + site4 + footlgth + site2 , data = df_train)
summary(model2)

## 
## Call:
## lm(formula = age ~ belly + hdlngth + site4 + footlgth + site2, 
##     data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.6201 -1.0987 -0.1218  0.8993  5.2125 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -19.62351    5.42767  -3.615 0.000544 ***
## belly         0.21205    0.07617   2.784 0.006816 ** 
## hdlngth       2.34746    0.76180   3.081 0.002892 ** 
## site4        -2.33717    0.81080  -2.883 0.005161 ** 
## footlgth     -0.76131    0.50565  -1.506 0.136420    
## site2         0.61466    0.59818   1.028 0.307504    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.598 on 74 degrees of freedom
## Multiple R-squared:  0.3112, Adjusted R-squared:  0.2647 
## F-statistic: 6.687 on 5 and 74 DF,  p-value: 3.439e-05

model3 <- lm(age ~ belly + hdlngth + site4  , data = df_train)
sum3 <- summary(model3)

Exploring Coefficient Significance and Model Simplification

The next step is to identify the coefficients that are least significant. This helps to eliminate predictors that may not have a meaningful impact on the dependent variable. Examining the confidence intervals of each coefficient confirms that none of the MLR (model1) coefficients include zero, so none can be identified as insignificant based on this criterion. However, examining the p-values revealed that some coefficients have p-values greater than 0.05, indicating they can be recognized as insignificant predictors. These include “site3,” “site1,” “total length,” and “eye.”

This suggests removing these insignificant predictors from Model1, which will help simplify the model as well. This change was applied to model2. However, still there are two more insignificant features of foot length and site2 which were removed in model3 which is generated only based on three predictors of “belly , hdlngth , site4”. It is noteworthy that the R-squared decreased from 0.2956 to 0.2596 ( for the trained dataset) which is pretty normal as it is a trade-off between model complexity and performance. The simpler model with lower adjusted R-squared is preferable.

# MSE and adjusted R^2 for model on the training and validation datasets
model3_training <- lm(age ~ belly + hdlngth + site4, data = df_train)
sum3_train <- summary(model3_training)

model3_valid <- lm(age ~ belly + hdlngth + site4, data = df_valid)
sum3_valid <- summary(model3_valid)

# Calculate metrics
mse_train <- mean((df_train$age - predict(model3_training, newdata = df_train))^2)
rmse_train <- sqrt(mse_train)

mse_valid <- mean((df_valid$age - predict(model3_valid, newdata = df_valid))^2)
rmse_valid <- sqrt(mse_valid)

# Create a data frame with the results
results <- data.frame(
  Metric = c("MSE", "RMSE"),
  Train = c( mse_train, rmse_train),
  Validation = c( mse_valid, rmse_valid)
)

# Print the results table
print(results)

##   Metric    Train Validation
## 1    MSE 2.443909   2.454851
## 2   RMSE 1.563301   1.566796

Based on the lowest MSE values calculated the best predictors for the model applied to the training and validation datasets as well as considering the significance of them, the final MLR model is generated with three features of belly, head length, and site4. Other selected features do not significantly impact the model and tend to increase its complexity without providing substantial benefits.

\(Q_{6}\) Final MLR model and Assessing Model Performance

The final formula of the MLR without any feature engineer is:

\[Y=b_{0}+b_{1}X^{(1)}+b_{2}X^{(2)}+b_{3}X^{(3)}\]

\[Age=-19.08029+0.20054$\cdot$ belly +1.77337\cdot hdlngth+-2.08071\cdot site4}\]

#final model=model3 based on chosen features in the previous step
final_model <- lm(age ~ belly + hdlngth + site4  , data = df_train)

sum <- summary(final_model)
print(sum)

## 
## Call:
## lm(formula = age ~ belly + hdlngth + site4, data = df_train)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9479 -1.0789 -0.1548  0.9750  4.8992 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -19.08029    5.32367  -3.584 0.000595 ***
## belly         0.20054    0.07607   2.636 0.010155 *  
## hdlngth       1.77337    0.67280   2.636 0.010171 *  
## site4        -2.08071    0.79394  -2.621 0.010592 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.604 on 76 degrees of freedom
## Multiple R-squared:  0.2877, Adjusted R-squared:  0.2596 
## F-statistic: 10.23 on 3 and 76 DF,  p-value: 9.774e-06

According to statistical summary of model, the adjusted R-squared shows that only 25% of the predicted age values are close to the actual values which can be indicative of poor performance of MLR model for this dataset. It could be due to non-leaner relationship between response and predictor variables which has been picked up in the first step of exploration of dataset (i.e., first figure showing a possible exponential relationship between age and total length dataset). Another possible reason could be lack of data for the ages older than 7 (2 and 1 datapoints were recorded for age 8 and 9, respectively.

This can be improved by feature engineering such as interaction terms between features (total length includes head length,etc), head size(skull width x head length), classifying datasets into male and female, removing the outliers, applying log into the age as its right-tailed response, etc.

As shown on the graph the predicted data contains more values for younger possums as expected earlier in this data analysis. This can mainly be due to positive skewness of age dataset as discussed earlier in \(Q_{1}\). Indeed the relationship between the age and the predictors is non-linear. Hence, either the feature engineering can improve the model or using another model.

# Residuals
train <- df_train
valid <- df_valid
# Prediction
train$predicted <- predict(final_model, df_train)
valid$predicted <- predict(final_model, df_valid)

train$residuals <- df_train$age - train$predicted
valid$residuals <- df_valid$age - valid$predicted
max_res_train <- max(train$residuals)
min_res_train <- min(train$residuals)
max_res_valid <- max(valid$residuals)
min_res_valid <- min(valid$residuals)

cat("residuals of the training dataset are between:",max_res_train, "and", min_res_train, "\n")

## residuals of the training dataset are between: 4.899174 and -2.947948

cat("residuals of the validation dataset are between:",max_res_valid, "and", min_res_valid, "\n")

## residuals of the validation dataset are between: 2.892377 and -3.606866

The upper and lower limits of residuals of models on the training and validation datasets indicate about 5 years discrepancy between the actual age and predicated age. This 5 years within 9 years of data collection is a very large deviation from the actual values. As shown in the plot of predicated age and the actual age (above) this deviation is more toward the right side of the graph, meaning that the older ages are underestimated as the younger age. As shown in the residual plots vs the observation order, which can be age for this dataset, the largest deviation (max residuals) is for predicted value for age 9.

# Add an indicator for training and validation data
predicted_train1 <- predict(final_model, newdata = train)
train_residuals <- sum$residuals

# Create a dataframe for residuals from the training data
residuals_df <- data.frame(Fitted = predicted_train1, Residuals = train_residuals, Age = train$age)
residuals_df$Dataset <- "Training"
# Create a dataframe for residuals from the validation data
valid_df <- data.frame(Fitted = valid$predicted, Residuals = valid$residuals, Age = valid$age)
valid_df <- data.frame(Fitted = valid$predicted, Residuals = valid$residuals, Age = df_valid$age)
valid_df$Dataset <- "Validation"

# Combine the two dataframes
combined_residuals_df <- rbind(residuals_df, valid_df)

# Plot 1: Residuals vs. Fitted Values
plot1 <- ggplot(combined_residuals_df, aes(x = Fitted, y = Residuals, color = Dataset)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black") +
  labs(title = "Residuals vs. Fitted Values", x = "Fitted Values", y = "Residuals") +
  scale_x_continuous(breaks = seq(floor(min(combined_residuals_df$Fitted)), ceiling(max(combined_residuals_df$Fitted)), by = 1)) +
  scale_y_continuous(breaks = seq(floor(min(combined_residuals_df$Residuals)), ceiling(max(combined_residuals_df$Residuals)), by = 2)) +
  theme_minimal()

# Plot 2: Residuals vs. Age
plot2 <- ggplot(combined_residuals_df, aes(x = Age, y = Residuals, color = Dataset)) +
  geom_point() +
  geom_hline(yintercept = 0, color = "black") +
  labs(title = "Residuals vs. Order of data", x = "Observation order (age)", y = "Residuals") +
  scale_x_continuous(breaks = seq(floor(min(combined_residuals_df$Age)), ceiling(max(combined_residuals_df$Age)), by = 1)) +
  scale_y_continuous(breaks = seq(floor(min(combined_residuals_df$Residuals)), ceiling(max(combined_residuals_df$Residuals)), by = 2)) +
  theme_minimal()

# Combine the plots into a 2 by 1 layout
combined_plot <- plot1 / plot2

# Display the combined plot
print(combined_plot)

Model Performance Metrics

To assess the performance of the MLR model, it was applied to a test dataset. The performance metrics are considered in this assignment are: accuracy, adjusted R value, Mean Squared Error (MSE) and Root Mean Squared Error (RMSE).

#Apply model to the test dataset
y_hat_test <- predict(final_model, newdata = df_test)

Accuracy

To computing the accuracy of the model, we considered the proportion of correctly predicted age values out of the total number of predication. The following formula is used for accuracy computaion:

\[Accuracy= \frac{Total Number of Predictions}{Number of Correct Predictions}\]

To make this metric compatible with our model, we defined a safe margin based on the standard error of the residuals obtained from the final model (on the training dataset), which was about 1 year. Due to a short age range, half of this value was considered asthe safe margin (6 months). Then, the residuals of the model on the test dataset were calculated. Finally, using the accuracy formula mentioned above, the accuracy of the model on the test dataset was computed.

############### Calculate the accuracy of the model on the TEST dataset
# Find a safe margin from the residuals of the training dataset 

# Step 1: Calculate the residuals from training dataset
residuals_train <- df_train$age - sum$residuals


# Step 2: Calculate the standard deviation of residuals
residuals_sd <- sd(residuals_train)

# Print the standard deviation of residuals
cat("Standard Deviation of Residuals: ", residuals_sd, "\n")

## Standard Deviation of Residuals:  0.9998404

# residuals of test
residuals_test <- df_test$age - y_hat_test

# Define a safe margin (e.g., ±2 years)
safe_margin <- 0.5*residuals_sd

# Calculate accuracy based on this safe margin
accurate_predictions <- sum(abs(residuals_test) <= safe_margin)
total_predictions <- length(residuals_test)

accuracy <- accurate_predictions / total_predictions

cat("Accuracy within", signif(safe_margin*12,1), "months is", signif(accuracy * 100,2), "%\n")

## Accuracy within 6 months is 18 %

# Calculate residuals for the test dataset
test <- df_test
test$predicted <- y_hat_test
test$residuals <- test$age - test$predicted

# collect all predicted, and age and residual in the same data frame
test_residuals_df <- data.frame(Fitted = test$predicted, 
                                Residuals = test$residuals,
                                Age = test$age)

print(max(test$residuals))

## [1] 4.716367

print(min(test$residuals))

## [1] -2.216724

# Plot 1: Residuals vs. Fitted Values
plot1 <- ggplot(test_residuals_df, aes(x = Fitted, y = Residuals)) +
  geom_point(color = "blue") +
  geom_hline(yintercept = 0, color = "red") +
  labs(title = "Residuals vs. Fitted Values (Test Data)", x = "Fitted Values", y = "Residuals") +
  scale_x_continuous(breaks = seq(floor(min(test_residuals_df$Fitted)), ceiling(max(test_residuals_df$Fitted)), by = 1)) +
  scale_y_continuous(breaks = seq(floor(min(test_residuals_df$Residuals)), ceiling(max(test_residuals_df$Residuals)), by = 1)) +
  theme_minimal()

# Plot 2: Residuals vs. Age
plot2 <- ggplot(test_residuals_df, aes(x = Age, y = Residuals)) +
  geom_point(color = "blue") +
  geom_hline(yintercept = 0, color = "red") +
  labs(title = "Residuals vs. Age (Test Data)", x = "Age", y = "Residuals") +
  scale_x_continuous(breaks = seq(floor(min(test_residuals_df$Age)), ceiling(max(test_residuals_df$Age)), by = 1)) +
  scale_y_continuous(breaks = seq(floor(min(test_residuals_df$Residuals)), ceiling(max(test_residuals_df$Residuals)), by = 1)) +
  theme_minimal()

# Combine the plots into a 2 by 1 layout
combined_plot <- plot1 / plot2

# Display the combined plot
print(combined_plot)

#Large residual values shown in these graph illustrated the low accuracy
#of the model. Only a few datapoints(Approx.4) can be considered close to the
#actual values.
#Other performance metrics assessed in this study are MSE and RMSE. MSE
#penalizes larger errors more heavily because the errors are squared,
#making it sensitive to outlines.