knitr::opts_chunk$set(echo = TRUE)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.4.3
library(splines)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tidyr)
# Load data bawaan
data("airquality")
air <- na.omit(airquality)  # hapus NA
# Library untuk spline
library(splines)
# Fungsi untuk evaluasi model
eval_model <- function(model, data, response){
  y_true <- data[[response]]
  y_pred <- predict(model, data)
  mse <- mean((y_true - y_pred)^2)        # hitung MSE manual
  aic <- AIC(model)                       # AIC
  adj_r2 <- summary(model)$adj.r.squared  # Adjusted R2
  return(c(MSE = mse, AIC = aic, AdjR2 = adj_r2))
}

### ----------------------------
### 1. Ozone ~ Temp
### ----------------------------
# Linear
mod_lin_temp <- lm(Ozone ~ Temp, data = air)
res_lin_temp <- eval_model(mod_lin_temp, air, "Ozone")
# Tangga (piecewise)
mod_step_temp <- lm(Ozone ~ cut(Temp, breaks=3), data = air)
res_step_temp <- eval_model(mod_step_temp, air, "Ozone")
# Spline biasa
mod_spline_temp <- lm(Ozone ~ bs(Temp, df=3), data = air)
res_spline_temp <- eval_model(mod_spline_temp, air, "Ozone")
# Natural spline
mod_ns_temp <- lm(Ozone ~ ns(Temp, df=3), data = air)
res_ns_temp <- eval_model(mod_ns_temp, air, "Ozone")
par(mar = c(4, 4, 2, 1))   # bottom, left, top, right
plot(mod_ns_temp, which = 1, main = "Residuals Ozone ~ Temp (Natural Spline)")

—————————-

2. Ozone ~ Wind

—————————-

# Linear
mod_lin_wind <- lm(Ozone ~ Wind, data = air)
res_lin_wind <- eval_model(mod_lin_wind, air, "Ozone")
# Tangga (piecewise)
mod_step_wind <- lm(Ozone ~ cut(Wind, breaks=3), data = air)
res_step_wind <- eval_model(mod_step_wind, air, "Ozone")
# Spline biasa
mod_spline_wind <- lm(Ozone ~ bs(Wind, df=3), data = air)
res_spline_wind <- eval_model(mod_spline_wind, air, "Ozone")
# Natural spline
mod_ns_wind <- lm(Ozone ~ ns(Wind, df=3), data = air)
res_ns_wind <- eval_model(mod_ns_wind, air, "Ozone")
par(mar = c(4, 4, 2, 1))   # bottom, left, top, right
plot(mod_ns_temp, which = 1, main = "Residuals Ozone ~ Wind (Natural Spline)")

### —————————- ### Gabungkan hasil ### —————————-

hasil <- rbind(
  Temp_Linear  = res_lin_temp,
  Temp_Step    = res_step_temp,
  Temp_Spline  = res_spline_temp,
  Temp_NSpline = res_ns_temp,
  Wind_Linear  = res_lin_wind,
  Wind_Step    = res_step_wind,
  Wind_Spline  = res_spline_wind,
  Wind_NSpline = res_ns_wind
)

print(round(hasil, 3))
##                  MSE      AIC AdjR2
## Temp_Linear  561.869 1023.775 0.483
## Temp_Step    541.609 1021.699 0.497
## Temp_Spline  496.085 1013.953 0.535
## Temp_NSpline 490.279 1012.647 0.541
## Wind_Linear  685.655 1045.876 0.369
## Wind_Step    706.858 1051.257 0.344
## Wind_Spline  531.415 1021.590 0.502
## Wind_NSpline 525.450 1020.337 0.508
ggplot(air, aes(x = Temp, y = Ozone)) +
  geom_point(color = "seashell4") +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "peru") +
  geom_smooth(method = "lm", formula = y ~ ns(x, df=3), se = FALSE, color = "rosybrown") +
  geom_smooth(method = "lm", formula = y ~ bs(x, df=3), se = FALSE, color = "palegreen4") +
  labs(
    title = "Hubungan Ozone ~ Temp",
    subtitle = "Biru = Linear, Merah = Natural Spline, Hijau = Spline",
    x = "Temperature", 
    y = "Ozone"
  )

# Ozone ~ Wind
ggplot(air, aes(x = Wind, y = Ozone)) +
  geom_point(color = "black") +
  geom_smooth(method = "lm", formula = y ~ x, se = FALSE, color = "skyblue") +
  geom_smooth(method = "lm", formula = y ~ ns(x, df=3), se = FALSE, color = "pink") +
  geom_smooth(method = "lm", formula = y ~ bs(x, df=3), se = FALSE, color = "palegreen") +
  labs(
    title = "Hubungan Ozone ~ Wind",
    subtitle = "Biru = Linear, Merah = Natural Spline, Hijau = Spline",
    x = "Wind", 
    y = "Ozone"
  )