Aziz Al-Najjar, Assignment 4

https://rpubs.com/Azizkhaled/1167315

Q1: Environment Canada Daily Weather

#libraries
options(warn = -1)
library(tidyverse)
library(lubridate)
library(RColorBrewer)
library(tidymodels) # a framework for modelling in tidyverse
library(GGally) # pairs plot
library(glue) # new favourite library
library(fda) # for a nice plot about splines
library(mgcv) # library we will use for splines
library(dplyr)

Step 1: Obtain Environment Canada’s daily temperature daily for any weather station. (I chose SAANICHTON CDA BC) Spanning 1970 to 2022.

Reason for choosing SAANICHTON: I tried multiple stations, SAANICHTON was among the few stations that has data between 1970 and 2022 with few Max.Temp missing values, similar to the Ottawa CDA data.

Data acquisition:

Similar to the tutorial provided, changing the stationID to 87 for SAANICHTON

base_url = c("https://climate.weather.gc.ca/climate_data/bulk_data_e.html?format=csv&stationID=87&Prov=BC&Year=", "&Month=1&Day=1&timeframe=2&submit=Download+Data")

destfile = list()
for(year in 1970:2022){
  download.file(paste0(
                      base_url[1],year,
                      base_url[2]),
                      destfile = paste0("Data/SAANICHTON",year,"-daily.csv"))
}

Creating the data set:

# Find all csv files in the "Data" folder.
files2load = list.files(path="Data",pattern="daily")
data_daily=NULL
for(lp in 1:length(files2load)){
  data_daily = rbind(data_daily,read_csv(paste0("Data/",files2load[lp])))
}

Step 2: Data Cleaning and Exploratory Data plots

Adding the day.of.year variable

data_daily <- data_daily |> 
  rename_all(make.names) |>
  mutate(Date.Time = as.Date(Date.Time, format = "%Y-%m-%d"),
         day.of.year = yday(Date.Time)) # Add day.of.year variable

Lets generate some plots to understand the data better

data_daily <- data_daily |> 
  dplyr::select(Date.Time, Month, Year, `day.of.year`, `Max.Temp...C.`)
# data_daily <- na.omit(data_daily)
data_daily |> 
  ggpairs(title = "Daily Data",
          axisLabels = "show")

We can see the plots above are not very pretty, but you still can see the relationship between day.of.year and Max.Temp. It hard to see the other relationship. Let’s try plotting them separately for better interpetation

Max temp vs date plot:

ggplot(data_daily, aes(x = Date.Time, y = Max.Temp...C., colour=factor(Year))) +
  geom_line(alpha = 0.7) + 
  geom_smooth(se = TRUE, color = "red") + 
  ggtitle("Time Series of Max Temperature with Trend") +
  xlab("Date") + ylab("Max Temperature (C)") +
  theme(legend.position = "none") 

We can see that the max temperature was fairly consistant accross the years. However, This plot does not provide much information due to huge data variation

Focusing on one month to see the trend over the years.

# Plot means of the Septemper max temperatures
data_daily |> 
  filter(Month == "09") |> 
  group_by(Year) |> 
  summarise(MeanOfSepMax = mean(Max.Temp...C.)) |>
  ungroup() |>
  ggplot() +
  geom_line(aes(x = Year, y = MeanOfSepMax, color = "red")) +
  ggtitle("Means of the Sep Max Temperatures") +
  labs(x = "Year", y = "degrees C") +
  theme(legend.position = "none") 

Despite the gaps indicating missing data, the remaining data seems to be consistent and shows a reasonable level of variation from year to year, suggesting no abrupt changes in temperature trends.

Plotting September data ignoring the missing days in some Septs

# Plot means of the June max and min temperatures
# Plot means of the Septemper max temperatures
data_daily |> 
  filter(Month == "09") |> 
  group_by(Year) |> 
  summarise(MeanOfSepMax = mean(Max.Temp...C., na.rm = TRUE)) |>
  ungroup() |>
  ggplot() +
  geom_line(aes(x = Year, y = MeanOfSepMax, color = "red")) +
  ggtitle("Means of the Sep Max Temperatures") +
  labs(x = "Year", y = "degrees C") +
  theme(legend.position = "none") 

This plot shows that there is no missing gaps or obvious outliers when ignoring the missing days. Might indicate that ignoring the missing values might be a safe approach since we lack the knowledge on how to import missing values.

Max Temp vs day.of.year

# Scatter Plot with Day of Year
ggplot(data_daily, aes(x = day.of.year, y = Max.Temp...C., , color=factor(Year))) +
  geom_point(alpha = 0.8) + 
  ggtitle("Max Temperature vs. Day of Year") +
 theme(legend.position = "none") 

From this plot, we can see that there is a clear seasonal pattern where temperatures peak in the middle of the year and are lower at the beginning and end, which corresponds to summer and winter seasons, respectively. The overlapping points suggest variability in maximum temperatures for the same days across different years. This indicates that Spline models might be a good fit here.

Step 3: Find best model to fit to Max.Temp…C. with Year and day.of.year as predictors.

data_daily <- data_daily |> 
  dplyr::select(Year, `day.of.year`, `Max.Temp...C.`)
data_daily <- na.omit(data_daily)

#going with option 1
# Split the data into training and testing sets
set.seed(123)  # Set a seed for reproducibility
data_split <- initial_split(data_daily, prop = 0.5)
train_data <- training(data_split)
test_data <- testing(data_split)

Linear models- Simple Linear :

lm_model = linear_reg() |>
set_engine("lm")

Simple_basic_recipe <- recipe(Max.Temp...C. ~ Year + `day.of.year` , data = train_data) 

# Create and fit the workflow with the training data
Simple_fitted_model <- workflow() |>
  add_recipe(Simple_basic_recipe) |>
  add_model(lm_model) |>
  fit(data = train_data)
summary(Simple_fitted_model|> extract_fit_engine())

Call:
stats::lm(formula = ..y ~ ., data = data)

Residuals:
     Min       1Q   Median       3Q      Max 
-21.9199  -4.6154  -0.6071   4.7348  22.1548 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -2.378e+01  8.527e+00  -2.789  0.00529 ** 
Year         1.805e-02  4.273e-03   4.225 2.41e-05 ***
day.of.year  9.230e-03  6.096e-04  15.142  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 6.275 on 9479 degrees of freedom
Multiple R-squared:  0.0255,    Adjusted R-squared:  0.02529 
F-statistic:   124 on 2 and 9479 DF,  p-value: < 2.2e-16

# Simple_fitted_model |> extract_fit_engine() |> plot(scheme = 2, residuals = TRUE)


# Generate predictions on the testing set
predictions <- predict(Simple_fitted_model, new_data = test_data) |>
  bind_cols(test_data)

# Plot the results
ggplot() +
  geom_point(aes(x = `day.of.year`, y = `Max.Temp...C.`, color = "Actual"), data = test_data) +
  geom_line(aes(x = `day.of.year`, y = .pred, color = "Predicted"), data = predictions) +
  labs(title = "Max Temperature Predictions by Day of Year (Simple_Linear)",
       x = "Day of Year",
       y = "Max Temperature (C)",
       color = "Legend") +
  theme(legend.position = "bottom")

The simple linear model roughly fitted a straight line over the average Max temp as in ‘Max Temperature Predictions by Day of Year (Simple_Linear)’ plot, which shows that the model is not capturing the variability in the data well. Not the best modelling choice, as we can also see in the low R-squared value of 0.025

Linear Models: step_cut

# Update the recipe with the training data

step_cut_recipe <- recipe(Max.Temp...C. ~ Year + `day.of.year` , data = train_data) |> 
  step_cut(`day.of.year`, breaks = seq(1, 365, length.out = 12+1))

# Create and fit the workflow with the training data
step_cut_fitted_model <- workflow() |>
  add_recipe(step_cut_recipe) |>
  add_model(lm_model) |>
  fit(data = train_data)

summary(step_cut_fitted_model|> extract_fit_engine())

Call:
stats::lm(formula = ..y ~ ., data = data)

Residuals:
    Min      1Q  Median      3Q     Max 
-17.663  -1.965  -0.142   1.810  16.549 

Coefficients:
                         Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -28.307360   4.227640  -6.696 2.27e-11 ***
Year                     0.017649   0.002118   8.333  < 2e-16 ***
day.of.year(31.3,61.7]   1.366141   0.155394   8.791  < 2e-16 ***
day.of.year(61.7,92]     3.475982   0.153416  22.657  < 2e-16 ***
day.of.year(92,122]      6.209084   0.155705  39.877  < 2e-16 ***
day.of.year(122,153]     9.580320   0.155928  61.441  < 2e-16 ***
day.of.year(153,183]    12.589525   0.153755  81.881  < 2e-16 ***
day.of.year(183,213]    15.239769   0.154984  98.332  < 2e-16 ***
day.of.year(213,244]    15.293723   0.155497  98.354  < 2e-16 ***
day.of.year(244,274]    12.119983   0.155445  77.970  < 2e-16 ***
day.of.year(274,304]     6.862771   0.155239  44.208  < 2e-16 ***
day.of.year(304,335]     2.293135   0.153797  14.910  < 2e-16 ***
day.of.year(335,365]    -0.053149   0.154632  -0.344    0.731    
day.of.year(365,366]    -1.512043   1.180135  -1.281    0.200    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.109 on 9468 degrees of freedom
Multiple R-squared:  0.761, Adjusted R-squared:  0.7607 
F-statistic:  2319 on 13 and 9468 DF,  p-value: < 2.2e-16

step_cut_fitted_model |> extract_fit_engine() |> plot(scheme = 2, residuals = TRUE)

Warning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameter

Warning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameter

Warning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameterWarning: "scheme" is not a graphical parameterWarning: "residuals" is not a graphical parameter

# Generate predictions on the testing set
predictions <- predict(step_cut_fitted_model, new_data = test_data) |>
  bind_cols(test_data)

# Plot the results
ggplot() +
  geom_point(aes(x = `day.of.year`, y = `Max.Temp...C.`, color = "Actual"), data = test_data) +
  geom_line(aes(x = `day.of.year`, y = .pred, color = "Predicted"), data = predictions) +
  labs(title = "Max Temperature Predictions by Day of Year with 12 step cut",
       x = "Day of Year",
       y = "Max Temperature (C)",
       color = "Legend") +
  theme(legend.position = "bottom")

With the step cut, the model effectively captures the seasonal variation, resulting in a pattern that aligns with the cyclical nature of temperature changes over the year. As we increase the number of step cuts, the model will fit better to the data. 12 cuts here roughly represents the months.

Linear Models: step_poly

# Update the recipe with the training data

poly_recipe <- recipe(Max.Temp...C. ~ Year + `day.of.year` , data = train_data) |> 
  step_poly(`day.of.year`, degree = 10)

# Create and fit the workflow with the training data
poly_model <- workflow() |>
  add_recipe(poly_recipe) |>
  add_model(lm_model) |>
  fit(data = train_data)

summary(poly_model|> extract_fit_engine())
poly_model |> extract_fit_engine() |> plot(scheme = 2, residuals = TRUE)


# Generate predictions on the testing set
predictions <- predict(poly_model, new_data = test_data) |>
  bind_cols(test_data)

# Plot the results
ggplot() +
  geom_point(aes(x = `day.of.year`, y = `Max.Temp...C.`, color = "Actual"), data = test_data) +
  geom_line(aes(x = `day.of.year`, y = .pred, color = "Predicted"), data = predictions) +
  labs(title = "Max Temperature Predictions by Day of Year with step_poly",
       x = "Day of Year",
       y = "Max Temperature (C)",
       color = "Legend") +
  theme(legend.position = "bottom")

Aiming for a smoother fit, the 10-degree poly_cut model delivers whats promised. The summary statistics show an improvement with a higher R-squared value of 0.7817, indicating that the polynomial model explains more of the variance in the data than the previous linear model with step cuts. The lower residual standard error also suggests better predictive accuracy. Higher degrees might fit the training data better, but might result in an overfitting

GAM models:

gam_model = gen_additive_mod( mode = "regression" ) |>
  set_engine("mgcv")
gam_fit = workflow() |>
  add_variables(outcomes = Max.Temp...C., predictors = c(Year, `day.of.year`)) |>
    add_model(gam_model,
      formula = Max.Temp...C. ~ s(Year) +
        s(`day.of.year`)) |>
fit(data = train_data)
gam_fit |> extract_fit_engine() |> summary()

Family: gaussian 
Link function: identity 

Formula:
Max.Temp...C. ~ s(Year) + s(day.of.year)

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 13.92762    0.03026   460.3   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                 edf Ref.df       F p-value    
s(Year)        8.913  8.998   26.12  <2e-16 ***
s(day.of.year) 8.771  8.985 3823.74  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

R-sq.(adj) =  0.785   Deviance explained = 78.5%
GCV = 8.6984  Scale est. = 8.6813    n = 9482

# gam_fit |> extract_fit_engine() |>
#   plot(residuals = TRUE, cex = 5, shade = TRUE)

gam_fit |> extract_fit_engine() |> plot(scheme = 2, residuals = TRUE)

# Generate predictions on the testing set
predictions <- predict(gam_fit, new_data = test_data) |>
  bind_cols(test_data)

# Plot the results
ggplot() +
  geom_point(aes(x = `day.of.year`, y = `Max.Temp...C.`, color = "Actual"), data = test_data) +
  geom_line(aes(x = `day.of.year`, y = .pred, color = "Predicted"), data = predictions) +
  labs(title = "Max Temperature Predictions by Day of Year with gam_fit",
       x = "Day of Year",
       y = "Max Temperature (C)",
       color = "Legend") +
  theme(legend.position = "bottom")

The GAM is fitting a non-linear relationship between the day of the year and temperature, allowing for a smooth curve that captures the seasonal trend. Compared to the linear model, the GAM provides a more flexible fit, accommodating the cyclical nature of temperature changes, which is evidenced by the smooth prediction curve closely tracking the actual data points.

The summary of the GAM shows that both the smooth term for the year (s(Year)) and the day of the year (s(day.of.year)) are statistically significant, with the model explaining 78.5% of the deviance (a good level of model fit for this type of data). The model’s adjusted R-squared value is 0.785, indicating that it accounts for a substantial proportion of the variability in the maximum temperatures, which is an improvement over simpler linear models and justifies the reason behind selecting it as the best fit to our data.

Q2: Option wine: An R attempt (Solution is submetted in jupyter notebook)

library(tidyverse) # readable code that is easier to troubleshoot
library(tidymodels) # a framework for modelling in tidyverse
library(vip) # variable importance
library(ISLR2) # Boston dataset

Load the dataset and perform data exploration

# Load the wine quality dataset
wine_data <- read_csv("Assignment_4_wine.csv")

Rows: 6497 Columns: 13── Column specification ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr  (1): type
dbl (12): fixed.acidity, volatile.acidity, citric.acid, residual.sugar, chlorides, free.sulfur.dioxide, total.sulfur.dioxide, density, pH, sulphates, alcohol, quality
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

quality_classes <- unique(wine_data$quality)
num_classes <- length(quality_classes)
print(num_classes)

[1] 7

quality_counts <- table(wine_data$quality)
print(quality_counts)

   3    4    5    6    7    8    9 
  30  216 2138 2836 1079  193    5 

ggplot(wine_data, aes(x = quality)) +
  geom_histogram(binwidth = 1, boundary = 0.5, fill = "blue", color = "black") +
  labs(title = "Histogram of Wine Quality", x = "Quality", y = "Count") +
  theme_minimal()

We can see how the quality data is uniformly distributed, which shows great class imbalance. We need to consider this when splitting the data between training and test. Using stratify will ensure that both testing and training sets has all classes

Data split and model initialization

# Split the data into training and testing sets
set.seed(42)
wine_data$quality <- as.factor(wine_data$quality)
data_split <- initial_split(wine_data, prop = 0.8, strata = quality)
train_data <- training(data_split)
test_data <- testing(data_split)
data_folds = vfold_cv(train_data, v = 5, strata = quality)

# Define the model and recipe
forest_model <- rand_forest(trees = tune(), mtry = tune(), mode = "classification") |> 
  set_engine("ranger", importance = "impurity") |> 
  set_mode("classification")

wine_recipe <- recipe(quality ~ ., data = train_data) 

# Create a workflow
model_flow = workflow() |>
  add_recipe( wine_recipe ) |>
  add_model( forest_model )

Model Fine-tuning with grid_search:

Considering the massive data imbalance, we include the precision and recall in the analysis.

my_metrics <- metric_set(accuracy, precision, recall, f_meas)
grid_ctrl = control_grid(save_pred = FALSE,
  allow_par = TRUE, parallel_over = "everything",
  save_workflow = TRUE)
T1 = Sys.time()
# Perform the grid search with the specified metrics
grid_tuned <- model_flow |> tune_grid(
  resamples = data_folds,
  grid = grid_regular(
    mtry(range = c(2, ncol(train_data) - 1)),
    trees(range = c(1, 600)),
    levels = 10
  ),
  metrics = my_metrics, 
  control = grid_ctrl
)
(elapsed = Sys.time()-T1)
# Print the elapsed time and grid search results

print(elapsed)

Time difference of 21.82304 mins

print(grid_tuned)

# Tuning results
# 10-fold cross-validation using stratification 

There were issues with some computations:

  - Warning(s) x1: While computing multiclass `precision()`, some levels had no predicted events (i.e. `true_positive + false_positive = 0`). Precision is und...

Run `show_notes(.Last.tune.result)` for more information.

best_results = grid_tuned |> select_best(metric =  "f_meas")
final_forest_fit = model_flow |>
  finalize_workflow(best_results) |>
  last_fit(split = data_split)
best_results

After fine-tuning, we see that the best f1 score performance occured with mtry = 2 and trees = 100

Plotting the variable importance:

final_workflow <- finalize_workflow(model_flow, best_results)
# Fit the final model to the entire training dataset
final_fit <- fit(final_workflow, data = train_data)

# Extract the fitted model object
Fitobject <- final_fit |> extract_fit_engine()

# Define the prediction function for the permutation importance
pred_fun <- function(model, newdata) {
  predict(model, newdata = newdata, type = "prob")$.pred_class
}

# Calculate permutation-based variable importance
set.seed(123) # For reproducibility
vip_perm <- vi(object = Fitobject, 
               feature_names = names(select(train_data, -quality)), 
               target = "quality",
               metric = "f_meas",
               train = train_data,
               pred_wrapper = pred_fun,
               nsim = 1000,
               parallel = TRUE,
               verbose = TRUE)

# Plot the variable importance
vip(vip_perm)

We can see that the alcohol level has larger importance compared to the other features. However, I am not sure why some features are not in the plot. I assume that ‘type’ feature has very little importance in predicting the quality, but to confirm this. I performed the same approach on python using scikit-learn. The results confirm the plots above.

Q3: The analysis and plot is in the ipynb file

results: we can confirm that following are the important variables: X11 , X1 , X3 , X47 ||| - X6 - X78 - X38 - X36 - X29 - X54