1.1 Sumber Dataset

Dataset ini adalah Seoul Bike Sharing Demand Dataset yang bersumber dari UCI Machine Learning Repository. Data merekam jumlah sepeda yang disewa per jam di sistem bike-sharing kota Seoul, Korea Selatan, beserta kondisi cuaca dan informasi waktu selama 1 tahun penuh (Desember 2017 – November 2018).

1.2 Load Library & Import Data

library(tidyverse)
library(ggcorrplot) 
library(car) 
library(gridExtra) 
library(skimr) 
library(readr) 
library(tibble)
library(knitr)
library(dplyr)
library(ggplot2)
library(corrplot)
library(lubridate) 
library(DT)

# Import data — sesuaikan path file CSV   
setwd("D:/Semester 4/Eksplorasi dan Visualisasi Data/Project")
seoulbike <- read.csv2("SeoulBikeData.csv", 
               sep = ";",
               fileEncoding = "latin1",
               check.names = FALSE,
               stringsAsFactors = FALSE)

# Convert ke tibble  
seoulbike <- as_tibble(seoulbike)

print(seoulbike)

## # A tibble: 8,760 × 14
##    Date       `Rented Bike Count`  Hour `Temperature (°C)` `Humidity (%)`
##    <chr>                    <int> <int>              <dbl>          <int>
##  1 01/12/2017                 254     0               -5.2             37
##  2 01/12/2017                 204     1               -5.5             38
##  3 01/12/2017                 173     2               -6               39
##  4 01/12/2017                 107     3               -6.2             40
##  5 01/12/2017                  78     4               -6               36
##  6 01/12/2017                 100     5               -6.4             37
##  7 01/12/2017                 181     6               -6.6             35
##  8 01/12/2017                 460     7               -7.4             38
##  9 01/12/2017                 930     8               -7.6             37
## 10 01/12/2017                 490     9               -6.5             27
## # ℹ 8,750 more rows
## # ℹ 9 more variables: `Wind speed (m/s)` <dbl>, `Visibility (10m)` <int>,
## #   `Dew point temperature (°C)` <dbl>, `Solar Radiation (MJ/m2)` <dbl>,
## #   `Rainfall (mm)` <dbl>, `Snowfall (cm)` <dbl>, Seasons <chr>, Holiday <chr>,
## #   `Functioning Day` <chr>

1.3 Jumlah Observasi dan Variabel

# Jumlah baris (observasi) dan kolom (variabel)
cat("Jumlah Observasi (baris):", nrow(seoulbike), "\n")

## Jumlah Observasi (baris): 8760

cat("Jumlah Variabel (kolom) :", ncol(seoulbike), "\n" )

## Jumlah Variabel (kolom) : 14

colnames(seoulbike)

##  [1] "Date"                       "Rented Bike Count"         
##  [3] "Hour"                       "Temperature (°C)"          
##  [5] "Humidity (%)"               "Wind speed (m/s)"          
##  [7] "Visibility (10m)"           "Dew point temperature (°C)"
##  [9] "Solar Radiation (MJ/m2)"    "Rainfall (mm)"             
## [11] "Snowfall (cm)"              "Seasons"                   
## [13] "Holiday"                    "Functioning Day"

1.4 Keterangan Variabel

1.5 Tujuan Analisis

Tujuan analisis adalah untuk memahami pola dan faktor-faktor yang mempengaruhi permintaan penyewaan sepeda di Seoul, meliputi :

• pengaruh kondisi cuaca (suhu, kelembaban, hujan, salju),
  
• pengaruh waktu (jam, musim, hari libur),
  
• mengidentifikasi variabel-variabel paling relevan untuk prediksi jumlah sepeda yang disewa.

2.1 Statistik Deskriptif

# Statistik Deskriptif
skim(seoulbike)

Variable type: character

Variable type: numeric

Interpretasi Statistik Deskriptif

Proses pemeriksaan dan ringkasan data dilakukan menggunakan fungsi skim() dari library skimr. Berdasarkan output ringkasan data, berikut adalah analisis mendalam mengenai karakteristik dataset Seoul Bike Sharing Demand:

A. Dimensi Data dan Kualitas Input

• Volume Data: Dataset ini merekam data operasional selama 1 tahun penuh (Desember 2017 – November 2018) dengan total 8.760 baris pengamatan. Angka ini merepresentasikan pencatatan per jam secara penuh dalam 365 hari (\(365 \times 24 \text{ jam}\)).
  
• Kandungan Informasi: Terdapat 14 variabel yang terbagi menjadi 10 variabel numerik (metrik cuaca dan volume rental) serta 4 variabel karakter/kategorikal (informasi waktu dan status kontrol).
  
• Kebersihan Data (Data Integrity): Seluruh variabel memiliki complete rate sebesar 100% (n_missing = 0). Hal ini menunjukkan bahwa dataset bersih dari missing values, sehingga tidak memerlukan tahapan imputasi data sebelum pemodelan.

B. Analisis Dinamika Variabel Target (Rented Bike Count)

• Volume Rata-rata: Secara umum, rata-rata sepeda yang disewa di Kota Seoul adalah 704,60 unit per jam.
  
• Volatilitas Tinggi: Standar deviasi yang besar (645,00) mencerminkan fluktuasi permintaan yang sangat masif antar-jam.
  
• Distribusi data: Rentang peminjaman bergerak dari 0 hingga 3.556 sepeda. Berdasarkan visualisasi histogram, bentuk distribusi condong ke kanan (right-skewed). Nilai kuartil 3 (p75) berada di angka 1.065,25, yang mengindikasikan bahwa lonjakan permintaan di atas 2.000 unit merupakan outlier berfrekuensi rendah yang dipicu oleh jam sibuk (peak hours) atau musim tertentu.

C. Analisis Variabel Kategorikal (Faktor Kalender & Sistem)

• Distribusi Waktu: Variabel Date mencakup 365 hari unik dan variabel Seasons menangkap karakteristik 4 musim di Korea Selatan secara proporsional.
• Kondisi Sistem (Functioning Day): Terdapat 2 nilai unik pada variabel ini. Fakta bahwa nilai minimum target berada di angka 0 sepeda berkaitan langsung dengan variabel kontrol ini, di mana pada hari-hari tertentu sistem tidak beroperasi (Non-functioning Day) karena pemeliharaan, sehingga menihilkan transaksi penyewaan secara artifisial.

D. Analisis Variabel Prediktor (Kondisi Mikroklimat Seoul)

Kondisi cuaca di Seoul menunjukkan variabilitas musiman yang sangat kontras, yang secara logis menjadi faktor penentu (driving factors) dalam keputusan pengguna untuk bersepeda:

• Termal (Temperature (°C) & Dew point temperature (°C)): Suhu udara di Seoul sangat ekstrem dengan rentang mencapai 57,2°C, bergerak dari -17,8°C pada puncak musim dingin hingga 39,4°C pada puncak musim panas, dengan rata-rata tahunan berada di suhu sejuk 12,88°C.

• Kelembaban & Visibilitas: Kelembaban rata-rata berada pada angka 58,23%. Menariknya, jarak pandang (Visibility) di Seoul mayoritas sangat baik dengan nilai median (p50) hingga nilai maksimum tertahan di angka 2.000 meter (batas atas sensor).

• Presipitasi (Rainfall (mm) & Snowfall (cm)): Rata-rata curah hujan (0,15 mm) dan salju (0,08 cm) per jam bernilai sangat kecil. Berdasarkan histogram, distribusi kedua variabel ini didominasi oleh nilai 0. Hal ini membuktikan bahwa cuaca buruk (hujan/salju lebat) adalah fenomena rare events (jarang terjadi), namun ketika terjadi, intensitasnya bisa ekstrem (hujan hingga 35 mm/jam dan salju hingga 8,8 cm).

E. Missing Value

Berdasarkan hasil skim() di atas, seluruh variabel menunjukkan n_missing = 0 yang artinya dataset tidak memiliki missing value dan siap digunakan untuk analisis lebih lanjut tanpa proses imputasi.

2.2 Visualisasi Distribusi

library(ggplot2)
library(patchwork)

## Warning: package 'patchwork' was built under R version 4.4.3

theme_clean <- theme_minimal(base_size = 11) +
  theme(
    plot.title = element_text(face = "bold", size = 12, margin = margin(b = 10)),
    panel.grid.major.x = element_blank(), # Hapus grid vertikal
    panel.grid.minor = element_blank(),   # Hapus semua grid minor
    axis.text.x = element_text(color = "#333333"),
    axis.text.y = element_text(color = "#333333")
  )

plot_color <- "#5B84B1" 

p1 <- ggplot(seoulbike, aes(x = Seasons)) + 
  geom_bar(fill = plot_color, width = 0.6) +
  labs(title = "Distribusi Musim", x = NULL, y = "Frekuensi") +
  theme_clean

p2 <- ggplot(seoulbike, aes(x = Holiday)) +
  geom_bar(fill = plot_color, width = 0.5) +
  labs(title = "Status Hari Libur", x = NULL, y = NULL) + # Hapus label Y
  theme_clean

p3 <- ggplot(seoulbike, aes(x = `Functioning Day`)) +
  geom_bar(fill = plot_color, width = 0.5) +
  labs(title = "Hari Operasional", x = NULL, y = NULL) + # Hapus label Y
  theme_clean

# Menggabungkan plot dengan patchwork
p1 + p2 + p3 + 
  plot_annotation(
    title = "Analisis Distribusi Variabel Kategorikal",
    theme = theme(plot.title = element_text(face = "bold", size = 14, hjust = 0.5, margin = margin(b = 15)))
  )

ggplot(seoulbike, aes(x = `Rented Bike Count`)) +
  geom_histogram(binwidth = 50, fill = "#5B84B1", color = "white") +
  labs(
    title = "Distribusi Jumlah Sepeda yang Disewa",
    x = "Jumlah Sepeda Disewa",
    y = "Frekuensi"
  ) +
  theme_minimal()

seoulbike %>%
  select(where(is.numeric), -Hour, -`Rented Bike Count`) %>%
  pivot_longer(cols = everything(),
               names_to = "Variabel",
               values_to = "Nilai") %>%
  ggplot(aes(x = Nilai)) +
  geom_histogram(bins = 30, fill = "#5B84B1",
                 color = "white", alpha = 0.7) +
  facet_wrap(~Variabel, scales = "free") +
  labs(title = "Distribusi Variabel Numerik") +
  theme_minimal() +
  theme(axis.text.x = element_text(size = 7))

seoulbike %>%
  group_by(Hour) %>%
  summarise(Rata_rata = mean(`Rented Bike Count`)) %>%
  ggplot(aes(x = Hour, y = Rata_rata)) +
  geom_line(color = "#5B84B1", linewidth = 1.2) +
  geom_point(color = "#5B84B1", size = 3) +
  labs(
    title = "Rata-rata Penyewaan Sepeda per Jam",
    x = "Jam",
    y = "Rata-rata Jumlah Sepeda"
  ) +
  scale_x_continuous(breaks = 0:23) +
  theme_minimal()

p4 <- ggplot(seoulbike, aes(x = Seasons, y = `Rented Bike Count`)) +
  geom_boxplot(fill = "#5B84B1", alpha = 0.7) +
  labs(title = "Per Musim", x = NULL, y = "Jumlah Sepeda") +
  theme_minimal()

p5 <- ggplot(seoulbike, aes(x = Holiday, y = `Rented Bike Count`)) +
  geom_boxplot(fill = "#5B84B1", alpha = 0.7) +
  labs(title = "Hari Libur vs Bukan", x = NULL, y = "Jumlah Sepeda") +
  theme_minimal()

grid.arrange(p4, p5, ncol = 2)

2.3 Analisis Korelasi

# Hitung matriks korelasi
cor_matrix <- seoulbike %>%
  select(where(is.numeric)) %>%
  cor(use = "complete.obs")

# Heatmap
ggcorrplot(
  cor_matrix,
  method = "square",
  type = "lower",
  lab = TRUE,
  lab_size = 3,
  colors = c("#D9E4EC", "white", "#5B84B1"),
  title = "Heatmap Korelasi Antar Variabel Numerik",
  ggtheme = theme_minimal()
)

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## ℹ The deprecated feature was likely used in the ggcorrplot package.
##   Please report the issue at <https://github.com/kassambara/ggcorrplot/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Analisis korelasi digunakan untuk mengukur kekuatan dan arah hubungan linear antar variabel numerik, khususnya terhadap variabel target Rented Bike Count.

Berdasarkan heatmap dan tabel korelasi di atas, Temperature memiliki korelasi positif tertinggi terhadap jumlah penyewaan sepeda (0.54), diikuti Hour (0.41) dan Dew Point Temperature (0.38). Kondisi cuaca cerah yang direpresentasikan oleh Solar Radiation (0.26) dan Visibility (0.20) juga berkontribusi positif. Sebaliknya, Humidity (-0.20), Snowfall (-0.14), dan Rainfall (-0.12) berkorelasi negatif yang mencerminkan bahwa kondisi cuaca buruk menurunkan minat bersepeda.

Tidak ada variabel dengan korelasi sangat kuat (>0.7) terhadap target, menunjukkan bahwa permintaan penyewaan sepeda ditentukan oleh kombinasi faktor cuaca dan waktu, bukan satu faktor tunggal. Yang perlu diwaspadai adalah korelasi antara Temperature dan Dew Point Temperature sebesar 0.91 antar sesama variabel prediktor, indikasi multikolinearitas yang akan dianalisis lebih lanjut pada sub-bab berikutnya.

2.4 Identifikasi Multikoliniearitas

Identifikasi multikolinearitas dilakukan untuk mendeteksi apakah ada variabel prediktor yang saling berkorelasi terlalu tinggi satu sama lain. Jika ada, variabel tersebut membawa informasi yang redundan dan dapat mengganggu analisis selanjutnya.

Multikolinearitas diuji menggunakan Variance Inflation Factor (VIF) Variabel dengan VIF > 10 mengindikasikan adanya multikolinearitas serius.

model_vif <- lm(`Rented Bike Count` ~
                  `Temperature (°C)` +
                  `Humidity (%)` +
                  `Wind speed (m/s)` +
                  `Visibility (10m)` +
                  `Dew point temperature (°C)` +
                  `Solar Radiation (MJ/m2)` +
                  `Rainfall (mm)`+
                  `Snowfall (cm)`,
                data = seoulbike )

vif_values <- vif(model_vif)

vif_df <- data.frame(
  Variabel = names(vif_values),
  VIF = round(vif_values, 2)
)

vif_df <- vif_df[order(-vif_df$VIF), ]

kable(vif_df, row.names = FALSE,
      caption = "Nilai VIF per Variabel")

Hasil uji VIF mengonfirmasi temuan pada analisis korelasi sebelumnya. Dew Point Temperature (VIF = 115.47) dan Temperature (VIF = 86.65) memiliki nilai VIF yang jauh melampaui batas toleransi 10, menunjukkan multikolinearitas yang sangat serius di antara keduanya. Hal ini konsisten dengan korelasi 0.91 yang terdeteksi sebelumnya kedua variabel membawa informasi yang hampir identik sehingga salah satunya perlu dieliminasi pada tahap Feature Selection. Humidity (VIF = 20.32) juga melampaui batas, kemungkinan karena hubungannya yang erat dengan suhu dan titik embun. Sementara itu, lima variabel lainnya memiliki VIF mendekati 1, mengindikasikan tidak ada masalah multikolinearitas.

2.5 Deteksi Outlier

Deteksi outlier dilakukan menggunakan metode IQR (Interquartile Range), yaitu dengan menghitung rentang antara kuartil pertama (Q1) dan kuartil ketiga (Q3). Data dikategorikan sebagai outlier apabila nilainya berada di luar batas berikut:

• Batas Bawah = Q1 - 1.5 × IQR
• Batas Atas = Q3 + 1.5 × IQR

Metode IQR dipilih karena bersifat robust , tidak terpengaruh oleh nilai ekstrem itu sendiri dalam perhitungannya, sehingga lebih andal dibanding metode berbasis mean dan standar deviasi (Z-score) yang rentan terhadap distribusi tidak normal. Mengingat beberapa variabel dalam dataset ini memiliki distribusi right-skewed, metode IQR lebih sesuai untuk digunakan.

num_df <- Filter(is.numeric, seoulbike)

outlier_result <- data.frame(
  Variabel = names(num_df),
  Jumlah_Outlier = sapply(num_df, function(x) {
    Q1 <- quantile(x, 0.25)
    Q3 <- quantile(x, 0.75)
    IQR_val <- Q3 - Q1
    sum(x < Q1 - 1.5 * IQR_val | x > Q3 + 1.5 * IQR_val)
  })
)

outlier_result <- outlier_result[order(-outlier_result$Jumlah_Outlier), ]

kable(outlier_result, row.names = FALSE,
      caption = "Jumlah Outlier per Variabel (Metode IQR)")

Berdasarkan hasil deteksi outlier menggunakan metode IQR, terdapat 5 variabel yang memiliki outlier dan 5 variabel lainnya bersih dari outlier.

Solar Radiation mencatat outlier terbanyak (641 data) diikuti Rainfall (528) dan Snowfall (443). Ketiganya konsisten dengan temuan distribusi sebelumnya, ketiga variabel ini didominasi nilai 0 (malam hari, tidak hujan, tidak bersalju), sehingga ketika nilainya tinggi secara otomatis terhitung sebagai outlier oleh metode IQR. Outlier pada ketiga variabel ini bukan data error, melainkan mencerminkan kondisi alam yang memang jarang terjadi.

Wind Speed (161) dan Rented Bike Count (158) juga memiliki outlier dalam jumlah lebih kecil, mencerminkan adanya jam-jam tertentu dengan angin kencang atau lonjakan permintaan yang tidak biasa.

Sementara itu, Hour, Temperature, Humidity, Visibility, dan Dew Point Temperature tidak memiliki outlier sama sekali , menunjukkan bahwa kelima variabel ini terdistribusi secara wajar tanpa nilai ekstrem.

3.1 Day of Week, Day, Month, Year and Day Category Features

Fitur Day of Week, Day, Month, dan Year dibuat untuk memecah informasi dari variabel Date menjadi komponen waktu spesifik, sehingga setiap elemen tanggal dapat dianalisis secara terpisah. Sementera fitur Day Category dibuat untuk mengidentifikasi hari ke dalam kategori hari kerja (weekday) atau akhir pekan (weekend). Variabel-variabel ini dibuat karena dataset awal tidak menyimpan informasi tersebut secara eksplisit. Padahal, pola penyewaan sepeda cenderung berbeda antara hari kerja dan akhir pekan, maupun antara bulan satu dengan bulan lainnya.

Variabel Date dalam dataset yang digunakan secara langsung, akan menghasilkan model yang hanya membaca tanggal sebagai string tanpa memahami konteks hari di baliknya. Dengan membuat variabel variabel-variabel tersebut, informasi dapat diekstrak dan direpresentasikan menjadi variabel yang lebih informatif dan mudah diinterpretasikan dalam analisis pola penyewaan sepeda.

seoulbike <- seoulbike %>%
  mutate(   
    Date = dmy(Date),
    
    # Pemisahan Date
    "Day of Week"  = weekdays(Date),
    Day = day(Date),
    Month = month((Date), label = TRUE),
    Year  = year(Date),
    
    # Kategori Hari
    `Day Category` = ifelse(`Day of Week` %in% c("Saturday", "Sunday"), "Weekend", "Weekday")
  )

DT::datatable(seoulbike)

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

# Visualisasi Ringkasan Day of Week Features
seoulbike %>%
  group_by(`Day of Week`) %>%
  summarise(
    Avg = round(mean(`Rented Bike Count`), 2),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = reorder(`Day of Week`, -Avg),
             y = Avg)) +
  geom_col(fill = "lightblue") +
  geom_text(aes(label = Avg)) +
  labs(
    title = "Average Bike Rentals by Day of Week",
    x = "Day of Week",
    y = "Average Rented Bike Count"
  ) +
  theme_minimal()

Grafik tersebut menunjukkan bahwa pada hari Jumat penyewaan sepeda cenderung mencatat rata-rata tertinggi. Rata-rata penyewaan sepeda pada hari kerja (Senin-Jumat) relatif stabil, sementara pada akhir pekan (Sabtu-Minggu) justru mengalami penuruanan. Hal ini dapat mengindikasikan bahwa sepeda lebih banyak digunakan untuk keperluan transportasi harian daripada rekreasi akhir pekan.

# Visualisasi Ringkasan Month Features
seoulbike %>%
  group_by(Month) %>%
  summarise(
    Avg = round(mean(`Rented Bike Count`), 2),
    .groups = "drop"
  ) %>%
  ggplot(aes(x = Month,
             y = Avg)) +
  geom_col(fill = "lightblue") +
  geom_text(aes(label = Avg)) +
  labs(
    title = "Average Bike Rentals by Month",
    x = "Month",
    y = "Average Rented Bike Count"
  ) +
  theme_minimal()

Grafik tersebut menunjukkan pola musiman yang sangat jelas dalam penyewaan sepeda, dimana rata-rata penyewaan terus meningkat pada bulan-bulan musim semi hingga musim panas (Maret, April, Mei, Juni, Juli, Agustus) dan mencapai puncaknya pada bulan Juni, kemudian mulai menurun secara bertahap memasuki musim gugur hingga mencapai titik terendah pada bulan-bulan musim dingin (Desember, Januari, Februari). Hal ini dapat mengindikasikan bahwa banyak orang di Seoul kurang nyaman untuk bersepeda pada bulan Desember-Februari mengingat kondisi suhu rendah dan kemungkinan cuaca bersalju.

# Visualisasi Ringkasan Day Category Features
seoulbike %>%
  group_by(`Day Category`, Hour) %>%
  summarise(Average = mean(`Rented Bike Count`), .groups = "drop") %>%
  ggplot(aes(x = Hour, 
             y = Average, 
             color = `Day Category`, 
             group = `Day Category`)) +
  geom_line() +
  geom_point() +
  labs(title = "Average Bike Rentals by Hour on Weekends and Weekdays") +
  theme_minimal()

Grafik tersebut memperlihatkan perbedaan pola penyewaan yang sangat kontras antara Weekdays dan Weekends:

• Weekdays memperlihatkan dua puncak tajam, yaitu sekitar pukul 08.00 (jam berangkat kerja/sekolah) dan sekitar pukul 18.00 (jam pulang kerja). Di luar jam-jam tersebut, permintaan turun signifikan.

• Weekends menunjukkan pola yang jauh lebih landai dan bertahap, dengan permintaan mulai meningkat perlahan sejak pagi dan mencapai puncak di siang hingga sore hari (sekitar pukul 12.00–15.00), lalu menurun kembali.

Perbedaan pola yang mencolok ini memperkuat pentingnya variabel Day Category sebagai fitur yang mampu menangkap perbedaan perilaku pengguna dalam model prediksi.

3.2 Rush Hour Feature

Variabel Rush Hour dibuat untuk merepresentasikan kondisi penyewaan sepeda pada waktu sibuk, yaitu waktu ketika aktivitas meningkat, seperti jam berangkat kerja (sekitar pukul 07.00-09.00) atau jam pulang kerja (sekitar pukul 17.00-19.00). Pada periode tersebut, kebutuhan transportasi biasanya lebih tinggi sehingga permintaan penyewaan sepeda juga cenderung meningkat.

Jika variabel Hour dalam dataset digunakan secara langsung, model hanya melihat jam sebagai angka biasa dari 0–23 tanpa memahami pola aktivitas manusia di baliknya. Dengan membuat fitur Rush Hour, informasi tersebut dapat diringkas menjadi variabel yang lebih sederhana dan lebih mudah diinterpretasikan.

# Feature Rush Hour
seoulbike <- seoulbike %>%
  mutate("Rush Hour" = ifelse(Hour %in% c(7:9, 17:19), "Rush Hour", "Non-Rush Hour"))

DT::datatable(seoulbike)

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

# Menampilkan Ringkasan Feature Rush Hour
seoulbike %>%
  group_by(`Rush Hour`)  %>%
  summarise(
    "Total Penyewaan" = sum(`Rented Bike Count`),
    Average = mean(`Rented Bike Count`)
  )

## # A tibble: 2 × 3
##   `Rush Hour`   `Total Penyewaan` Average
##   <chr>                     <int>   <dbl>
## 1 Non-Rush Hour           3944254    600.
## 2 Rush Hour               2228060   1017.

Hasil ringkasan menunjukkan bahwa total dan rata-rata penyewaan sepeda pada periode Rush Hour lebih tinggi dibandingkan dengan periode Non-Rush Hour. Hal ini mengonfirmasi asumsi awal, bahwa pada jam sibuk (07.00–09.00 dan 17.00–19.00) permintaan sepeda lebih tinggi. Hal ini kemungkinan besar terjadi karena sepeda digunakan sebagai moda transportasi alternatif untuk menghindari kemacetan. Fitur ini berhasil mengekstrak informasi penting dari variabel Hour ke dalam bentuk yang lebih ringkas dan mudah diinterpretasikan.

3.3 Feature Precipitation

Dataset memiliki dua variabel cuaca yang berkaitan dengan kondisi presipitasi (proses jatuhnya materi air dari atmosfer), yaitu Rainfall (mm) dan Snowfall (cm). Keduanya merepresentasikan kondisi yang secara langsung dapat mempengaruhi kenyamanan dan keamanan bersepeda. Fitur Precipitation dibuat sebagai variabel biner yang merangkum keduanya, dimana:

• Hujan/Salju jika terdapat curah hujan atau salju.
• Tidak Hujan/Salju jika tidak terdapat curah hujan atau salju.

Pendekatan ini dipilih karena dari perspektif pengguna sepeda, yang paling berpengaruh adalah ada atau tidaknya presipitasi, bukan seberapa besar intensitasnya.

# Feature Precipitation
seoulbike <- seoulbike %>%
  mutate(Precipitation = ifelse(`Rainfall (mm)` > 0 | `Snowfall (cm)` > 0, "Hujan/Salju", "Tidak Hujan/Salju"))

DT::datatable(seoulbike)

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

# Menampilkan Ringkasan Feature Precipitation
seoulbike %>%
  group_by(Precipitation) %>%
  summarise(
    "Total Penyewaan" = sum(`Rented Bike Count`),
    Average           = mean(`Rented Bike Count`)
  )

## # A tibble: 2 × 3
##   Precipitation     `Total Penyewaan` Average
##   <chr>                         <int>   <dbl>
## 1 Hujan/Salju                  165637    176.
## 2 Tidak Hujan/Salju           6006677    768.

Hasil ringkasan menunjukkan perbedaan yang signifikan dalam rata-rata penyewaan sepeda antara kondisi presipitasi dan tidak. Rata-rata penyewaan saat kondisi Tidak Hujan/Salju jauh lebih tinggi dibandingkan saat Hujan/Salju. Temuan ini sangat logis karena kondisi hujan atau salju membuat penggunaan sepeda menjadi tidak nyaman bahkan berisiko, sehingga banyak pengguna beralih ke moda transportasi lain. Fitur Precipitation terbukti mampu menangkap pola penurunan permintaan saat kondisi cuaca buruk secara efektif.

3.4 Encoding Variabel Kategorik

Encoding Variabel Kategorik adalah teknik untuk mengubah variabel bertipe kategorik (teks/faktor) menjadi representasi numerik agar dapat diproses oleh algoritma machine learning maupun model statistik. Hal ini dilakukan karena sebagian besar algoritma hanya dapat memproses data dalam format angka, sehingga variabel kategorik yang dihasilkan dari tahap feature engineering maupun yang sudah ada sebelumnya perlu dikonversi terlebih dahulu.

Terdapat dua teknik encoding yang digunakan dalam project ini, yaitu:

• Label Encoding: Digunakan untuk variabel dengan urutan yang bermakna (ordinal), seperti Seasons, Day of Week, dan Month. Label encoding mengubah kategori menjadi angka bulat berurutan yang berurutan.

• Binary/One-Hot Encoding: Digunakan untuk variabel biner atau yang tidak memiliki urutan, seperti Functioning Day, Holiday, Day Category, Rush Hour, dan Precipitation. Variabel ini dikonversi menjadi 0 atau 1 untuk merepresentasikan dua kondisi yang mungkin terjadi.

# Encoding Variabel Kategorik
seoulbike <- seoulbike %>%
  mutate(
    Seasons =  as.integer(factor(Seasons, levels = c("Spring", "Summer", "Autumn", "Winter"))), # Label Encoding
    
    `Functioning Day` =  ifelse(`Functioning Day` == "Yes", 1L, 0L), # One Hot Encoding
    
    Holiday =  ifelse(Holiday == "Holiday", 1L, 0L), # One Hot Encoding
    
    `Day of Week` = as.integer(factor(`Day of Week`, levels = c("Monday", "Tuesday", "Wednesday", "Thursday", "Friday", "Saturday", "Sunday"))), # Label Encoding
    
    Month = as.integer(factor(Month, levels = c("Jan", "Feb", "Mar", "Apr", "May", "Jun", "Jul", "Aug", "Sep", "Oct", "Nov", "Dec"))), # Label Encoding
    
    `Day Category` = ifelse(`Day Category` == "Weekend", 1L, 0L), # One Hot Encoding
    
    `Rush Hour` =  ifelse(`Rush Hour` == "Rush Hour", 1L, 0L), # One Hot Encoding
    
    Precipitation =  ifelse(Precipitation == "Hujan/Salju", 1L, 0L), # One Hot Encoding
  )

DT::datatable(seoulbike)

## Warning in instance$preRenderHook(instance): It seems your data is too big for
## client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

Setelah proses encoding, seluruh variabel dalam dataset kini berbentuk numerik dan siap digunakan dalam proses pemodelan. Variabel Seasons, Day of Week, dan Month telah dikonversi menggunakan label encoding sesuai urutan yang bermakna, sementara Functioning Day, Holiday, Rush Hour, dan Precipitation telah dikonversi menjadi nilai biner (0 dan 1). Dataset kini memiliki representasi yang lebih konsisten dengan berbagai metode analisis statistik maupun machine learning.

4.1 Filter Method

Filter Method adalah pendekatan feature selection yang mengevaluasi relevansi variabel berdasarkan karakteristik statistiknya. Proses filter dilakukan dengan mengevaluasi nilai Variance Inflation Factor (VIF) untuk mendeteksi multikolinearitas, yaitu kondisi di mana dua atau lebih variabel prediktor saling berkorelasi kuat satu sama lain, sehingga informasi yang terkandung dalam satu variabel sebagian besar sudah terwakili oleh variabel lainnya.

Interpretasi nilai VIF:

• VIF = 1: Tidak ada korelasi antar variabel prediktor
  
• VIF 1–5: Korelasi sedang, umumnya masih dapat ditoleransi
  
• VIF > 10: Korelasi sangat tinggi, indikasi multikolinearitas serius yang perlu ditangani

Adanya multikolinearitas dapat menyebabkan estimasi koefisien regresi menjadi tidak stabil, sulit diinterpretasikan, dan menghasilkan standar error yang besar. Oleh karena itu, penting untuk mengevaluasi dan memilih variabel yang paling relevan.

# Seleksi variabel numerik
df_num <- seoulbike %>% 
  select(where(is.numeric))

# Membuat model
model <- lm(`Rented Bike Count` ~ ., data = df_num)
summary(model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ ., data = df_num)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1272.65  -251.34   -38.91   211.99  2016.75 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -5.968e+05  4.478e+04 -13.327  < 2e-16 ***
## Hour                          2.539e+01  6.809e-01  37.292  < 2e-16 ***
## `Temperature (°C)`            2.122e+01  3.415e+00   6.216 5.35e-10 ***
## `Humidity (%)`               -6.869e+00  9.733e-01  -7.057 1.83e-12 ***
## `Wind speed (m/s)`            1.675e+01  4.748e+00   3.527 0.000422 ***
## `Visibility (10m)`            1.431e-02  9.160e-03   1.562 0.118288    
## `Dew point temperature (°C)`  2.075e+00  3.565e+00   0.582 0.560486    
## `Solar Radiation (MJ/m2)`    -4.010e+01  7.097e+00  -5.651 1.65e-08 ***
## `Rainfall (mm)`              -3.604e+01  4.271e+00  -8.438  < 2e-16 ***
## `Snowfall (cm)`               1.098e+02  1.190e+01   9.224  < 2e-16 ***
## Seasons                      -4.973e+01  5.170e+00  -9.619  < 2e-16 ***
## Holiday                      -1.170e+02  2.004e+01  -5.838 5.48e-09 ***
## `Functioning Day`             9.080e+02  2.442e+01  37.181  < 2e-16 ***
## `Day of Week`                 5.522e+00  3.502e+00   1.577 0.114932    
## Day                          -6.229e-01  4.997e-01  -1.247 0.212596    
## Month                         3.417e+01  1.825e+00  18.724  < 2e-16 ***
## Year                          2.955e+02  2.219e+01  13.317  < 2e-16 ***
## `Day Category`               -9.988e+01  1.550e+01  -6.445 1.22e-10 ***
## `Rush Hour`                   3.600e+02  1.007e+01  35.759  < 2e-16 ***
## Precipitation                -2.975e+02  1.874e+01 -15.877  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 8740 degrees of freedom
## Multiple R-squared:  0.6141, Adjusted R-squared:  0.6133 
## F-statistic: 732.1 on 19 and 8740 DF,  p-value: < 2.2e-16

# Menghitung Variance Inflation Factor (VIF)
sort(vif(model), decreasing = TRUE)

## `Dew point temperature (°C)`           `Temperature (°C)` 
##                   117.999028                    90.570602 
##               `Humidity (%)`                `Day of Week` 
##                    21.384751                     2.666182 
##               `Day Category`                        Month 
##                     2.664091                     2.155179 
##                         Year    `Solar Radiation (MJ/m2)` 
##                     2.083753                     2.069549 
##                Precipitation                      Seasons 
##                     1.836193                     1.814968 
##           `Visibility (10m)`              `Snowfall (cm)` 
##                     1.690322                     1.470232 
##           `Wind speed (m/s)`              `Rainfall (mm)` 
##                     1.318097                     1.264062 
##                         Hour            `Functioning Day` 
##                     1.209501                     1.056706 
##                          Day                  `Rush Hour` 
##                     1.052015                     1.034631 
##                      Holiday 
##                     1.024851

Berdasarkan hasil perhitungan Variance Inflation Factor (VIF), variabel Dew point temperature (°C), Temperature (°C), dan Humidity (%) memiliki nilai VIF yang sangat tinggi, yaitu dmasing-masing sekitar 117.999028, 90.570602, 21.384751. Hal ini menunjukkan bahwa ketiga variabel ini saling berkorelasi kuat sehingga menunjukkan adanya multikolinearitas dalam model.

Variabel Dew point temperature (°C) dipilih untuk dihapus dari model dengan mempertimbangkan hubungan variabel Dew point temperature (°C) dengan variabel Temperature (°C), dan Humidity (%), di mana secara konsep variabel dew point temperature dapat dihitung melalui temperature dan humidity.

# Mengeliminasi variabel penyebab multikolinearitas
df_num_2 <- df_num %>%
  select(-`Dew point temperature (°C)`)

# Membuat model baru
model_2 <- lm(`Rented Bike Count` ~ ., data = df_num_2)

# Menghitung Variance Inflation Factor (VIF) model baru
vif(model_2)

##                      Hour        `Temperature (°C)`            `Humidity (%)` 
##                  1.207183                  2.840763                  2.823433 
##        `Wind speed (m/s)`        `Visibility (10m)` `Solar Radiation (MJ/m2)` 
##                  1.316743                  1.681253                  1.976511 
##           `Rainfall (mm)`           `Snowfall (cm)`                   Seasons 
##                  1.258169                  1.469864                  1.809970 
##                   Holiday         `Functioning Day`             `Day of Week` 
##                  1.024848                  1.056623                  2.660184 
##                       Day                     Month                      Year 
##                  1.052015                  2.153991                  2.083523 
##            `Day Category`               `Rush Hour`             Precipitation 
##                  2.664002                  1.034631                  1.817553

Setelah Dew point temperature (°C) dieliminasi, seluruh variabel memiliki nilai VIF yang rendah. Hal ini menunjukkan bahwa model kini memiliki variabel-variabel prediktor yang lebih independen satu sama lain, sehingga estimasi koefisien regresi menjadi lebih stabil dan dapat diandalkan.

Dari hasil seleksi ini, variabel yang dipertahankan untuk tahap analisis adalah seluruh variabel numerik kecuali Dew point temperature (°C), yaitu:
- Hour
- Temperature (°C)
- Humidity (%)
- Wind speed (m/s)
- Visibility (10m)
- Solar Radiation (MJ/m2)
- Rainfall (mm)
- Snowfall (cm)
- Seasons
- Holiday
- Functioning Day
- Day of Week
- Month
- Year
- Day Category
- Rush Hour
- Precipitation

4.2 Stepwise Regression

Stepwise Regression adalah metode seleksi fitur berbasis wrapper yang menggunakan kriteria statistik untuk menentukan variabel mana yang harus dimasukkan atau dikeluarkan dari model. Berbeda dengan Filter Method yang mengevaluasi variabel secara individual, Stepwise Regression mempertimbangkan kombinasi variabel dan dampaknya secara bersama-sama terhadap kualitas model.

Kriteria yang digunakan adalah AIC (Akaike Information Criterion), di mana semakin kecil nilai AIC, semakin baik model karena berhasil menyeimbangkan antara akurasi prediksi dan kompleksitas model (jumlah variabel). Metode Stepwise menggunakan tiga pendekatan untuk mengevaluasi variabel secara bertahap (step by step).

• Forward Selection: Dimulai dari model kosong (null model), lalu variabel ditambahkan satu per satu selama penambahan tersebut dapat menurunkan nilai AIC.
• Backward Elimination: Dimulai dari model penuh (full model), lalu variabel dihapus satu per satu selama penghapusan tersebut dapat menurunkan nilai AIC.
• Bidirectional Stepwise: Kombinasi keduanya, dimana variabel dapat ditambahkan maupun dihapus pada setiap langkah.

# model lengkap
full_model <- lm(`Rented Bike Count` ~ ., data = df_num)
summary(full_model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ ., data = df_num)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1272.65  -251.34   -38.91   211.99  2016.75 
## 
## Coefficients:
##                                Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                  -5.968e+05  4.478e+04 -13.327  < 2e-16 ***
## Hour                          2.539e+01  6.809e-01  37.292  < 2e-16 ***
## `Temperature (°C)`            2.122e+01  3.415e+00   6.216 5.35e-10 ***
## `Humidity (%)`               -6.869e+00  9.733e-01  -7.057 1.83e-12 ***
## `Wind speed (m/s)`            1.675e+01  4.748e+00   3.527 0.000422 ***
## `Visibility (10m)`            1.431e-02  9.160e-03   1.562 0.118288    
## `Dew point temperature (°C)`  2.075e+00  3.565e+00   0.582 0.560486    
## `Solar Radiation (MJ/m2)`    -4.010e+01  7.097e+00  -5.651 1.65e-08 ***
## `Rainfall (mm)`              -3.604e+01  4.271e+00  -8.438  < 2e-16 ***
## `Snowfall (cm)`               1.098e+02  1.190e+01   9.224  < 2e-16 ***
## Seasons                      -4.973e+01  5.170e+00  -9.619  < 2e-16 ***
## Holiday                      -1.170e+02  2.004e+01  -5.838 5.48e-09 ***
## `Functioning Day`             9.080e+02  2.442e+01  37.181  < 2e-16 ***
## `Day of Week`                 5.522e+00  3.502e+00   1.577 0.114932    
## Day                          -6.229e-01  4.997e-01  -1.247 0.212596    
## Month                         3.417e+01  1.825e+00  18.724  < 2e-16 ***
## Year                          2.955e+02  2.219e+01  13.317  < 2e-16 ***
## `Day Category`               -9.988e+01  1.550e+01  -6.445 1.22e-10 ***
## `Rush Hour`                   3.600e+02  1.007e+01  35.759  < 2e-16 ***
## Precipitation                -2.975e+02  1.874e+01 -15.877  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 8740 degrees of freedom
## Multiple R-squared:  0.6141, Adjusted R-squared:  0.6133 
## F-statistic: 732.1 on 19 and 8740 DF,  p-value: < 2.2e-16

# model null
null_model <- lm(`Rented Bike Count` ~ 1, data = df_num)
summary(null_model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ 1, data = df_num)
## 
## Residuals:
##    Min     1Q Median     3Q    Max 
## -704.6 -513.6 -200.1  360.6 2851.4 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  704.602      6.891   102.2   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 645 on 8759 degrees of freedom

# forward selection 
forward_model <- step(null_model, scope = list(lower = null_model,upper = full_model), direction = "forward", trace = 0)
summary(forward_model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ `Temperature (°C)` + Hour + 
##     `Rush Hour` + `Functioning Day` + `Humidity (%)` + Precipitation + 
##     Month + Year + `Rainfall (mm)` + Seasons + `Snowfall (cm)` + 
##     `Day Category` + Holiday + `Solar Radiation (MJ/m2)` + `Wind speed (m/s)` + 
##     `Visibility (10m)` + `Day of Week`, data = df_num)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1275.08  -250.97   -38.39   211.48  2016.50 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.952e+05  4.477e+04 -13.295  < 2e-16 ***
## `Temperature (°C)`         2.313e+01  6.036e-01  38.326  < 2e-16 ***
## Hour                       2.539e+01  6.800e-01  37.344  < 2e-16 ***
## `Rush Hour`                3.601e+02  1.007e+01  35.769  < 2e-16 ***
## `Functioning Day`          9.057e+02  2.436e+01  37.184  < 2e-16 ***
## `Humidity (%)`            -6.273e+00  3.494e-01 -17.952  < 2e-16 ***
## Precipitation             -2.995e+02  1.863e+01 -16.077  < 2e-16 ***
## Month                      3.405e+01  1.823e+00  18.681  < 2e-16 ***
## Year                       2.947e+02  2.218e+01  13.284  < 2e-16 ***
## `Rainfall (mm)`           -3.616e+01  4.261e+00  -8.486  < 2e-16 ***
## Seasons                   -4.948e+01  5.162e+00  -9.585  < 2e-16 ***
## `Snowfall (cm)`            1.084e+02  1.185e+01   9.144  < 2e-16 ***
## `Day Category`            -9.971e+01  1.550e+01  -6.435 1.30e-10 ***
## Holiday                   -1.170e+02  2.004e+01  -5.841 5.37e-09 ***
## `Solar Radiation (MJ/m2)` -4.032e+01  6.916e+00  -5.831 5.71e-09 ***
## `Wind speed (m/s)`         1.681e+01  4.744e+00   3.544 0.000396 ***
## `Visibility (10m)`         1.649e-02  9.022e-03   1.828 0.067601 .  
## `Day of Week`              5.407e+00  3.498e+00   1.546 0.122258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 8742 degrees of freedom
## Multiple R-squared:  0.614,  Adjusted R-squared:  0.6133 
## F-statistic: 818.1 on 17 and 8742 DF,  p-value: < 2.2e-16

# backward elimination 
backward_model <- step(full_model, direction = "backward", trace = 0) 
summary(backward_model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ Hour + `Temperature (°C)` + 
##     `Humidity (%)` + `Wind speed (m/s)` + `Visibility (10m)` + 
##     `Solar Radiation (MJ/m2)` + `Rainfall (mm)` + `Snowfall (cm)` + 
##     Seasons + Holiday + `Functioning Day` + `Day of Week` + Month + 
##     Year + `Day Category` + `Rush Hour` + Precipitation, data = df_num)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1275.08  -250.97   -38.39   211.48  2016.50 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.952e+05  4.477e+04 -13.295  < 2e-16 ***
## Hour                       2.539e+01  6.800e-01  37.344  < 2e-16 ***
## `Temperature (°C)`         2.313e+01  6.036e-01  38.326  < 2e-16 ***
## `Humidity (%)`            -6.273e+00  3.494e-01 -17.952  < 2e-16 ***
## `Wind speed (m/s)`         1.681e+01  4.744e+00   3.544 0.000396 ***
## `Visibility (10m)`         1.649e-02  9.022e-03   1.828 0.067601 .  
## `Solar Radiation (MJ/m2)` -4.032e+01  6.916e+00  -5.831 5.71e-09 ***
## `Rainfall (mm)`           -3.616e+01  4.261e+00  -8.486  < 2e-16 ***
## `Snowfall (cm)`            1.084e+02  1.185e+01   9.144  < 2e-16 ***
## Seasons                   -4.948e+01  5.162e+00  -9.585  < 2e-16 ***
## Holiday                   -1.170e+02  2.004e+01  -5.841 5.37e-09 ***
## `Functioning Day`          9.057e+02  2.436e+01  37.184  < 2e-16 ***
## `Day of Week`              5.407e+00  3.498e+00   1.546 0.122258    
## Month                      3.405e+01  1.823e+00  18.681  < 2e-16 ***
## Year                       2.947e+02  2.218e+01  13.284  < 2e-16 ***
## `Day Category`            -9.971e+01  1.550e+01  -6.435 1.30e-10 ***
## `Rush Hour`                3.601e+02  1.007e+01  35.769  < 2e-16 ***
## Precipitation             -2.995e+02  1.863e+01 -16.077  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 8742 degrees of freedom
## Multiple R-squared:  0.614,  Adjusted R-squared:  0.6133 
## F-statistic: 818.1 on 17 and 8742 DF,  p-value: < 2.2e-16

# stepwise (bidirectional) 
stepwise_model <- step(null_model, scope = list(lower = null_model, upper = full_model), direction = "both", trace = 0) 
summary(stepwise_model)

## 
## Call:
## lm(formula = `Rented Bike Count` ~ `Temperature (°C)` + Hour + 
##     `Rush Hour` + `Functioning Day` + `Humidity (%)` + Precipitation + 
##     Month + Year + `Rainfall (mm)` + Seasons + `Snowfall (cm)` + 
##     `Day Category` + Holiday + `Solar Radiation (MJ/m2)` + `Wind speed (m/s)` + 
##     `Visibility (10m)` + `Day of Week`, data = df_num)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1275.08  -250.97   -38.39   211.48  2016.50 
## 
## Coefficients:
##                             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)               -5.952e+05  4.477e+04 -13.295  < 2e-16 ***
## `Temperature (°C)`         2.313e+01  6.036e-01  38.326  < 2e-16 ***
## Hour                       2.539e+01  6.800e-01  37.344  < 2e-16 ***
## `Rush Hour`                3.601e+02  1.007e+01  35.769  < 2e-16 ***
## `Functioning Day`          9.057e+02  2.436e+01  37.184  < 2e-16 ***
## `Humidity (%)`            -6.273e+00  3.494e-01 -17.952  < 2e-16 ***
## Precipitation             -2.995e+02  1.863e+01 -16.077  < 2e-16 ***
## Month                      3.405e+01  1.823e+00  18.681  < 2e-16 ***
## Year                       2.947e+02  2.218e+01  13.284  < 2e-16 ***
## `Rainfall (mm)`           -3.616e+01  4.261e+00  -8.486  < 2e-16 ***
## Seasons                   -4.948e+01  5.162e+00  -9.585  < 2e-16 ***
## `Snowfall (cm)`            1.084e+02  1.185e+01   9.144  < 2e-16 ***
## `Day Category`            -9.971e+01  1.550e+01  -6.435 1.30e-10 ***
## Holiday                   -1.170e+02  2.004e+01  -5.841 5.37e-09 ***
## `Solar Radiation (MJ/m2)` -4.032e+01  6.916e+00  -5.831 5.71e-09 ***
## `Wind speed (m/s)`         1.681e+01  4.744e+00   3.544 0.000396 ***
## `Visibility (10m)`         1.649e-02  9.022e-03   1.828 0.067601 .  
## `Day of Week`              5.407e+00  3.498e+00   1.546 0.122258    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 401.1 on 8742 degrees of freedom
## Multiple R-squared:  0.614,  Adjusted R-squared:  0.6133 
## F-statistic: 818.1 on 17 and 8742 DF,  p-value: < 2.2e-16

# AIC 
AIC(forward_model, backward_model, stepwise_model) # semakin kecil lebih baik

##                df      AIC
## forward_model  19 129898.4
## backward_model 19 129898.4
## stepwise_model 19 129898.4

# melihat variabel terpilih
formula(forward_model) 

## `Rented Bike Count` ~ `Temperature (°C)` + Hour + `Rush Hour` + 
##     `Functioning Day` + `Humidity (%)` + Precipitation + Month + 
##     Year + `Rainfall (mm)` + Seasons + `Snowfall (cm)` + `Day Category` + 
##     Holiday + `Solar Radiation (MJ/m2)` + `Wind speed (m/s)` + 
##     `Visibility (10m)` + `Day of Week`

formula(backward_model) 

## `Rented Bike Count` ~ Hour + `Temperature (°C)` + `Humidity (%)` + 
##     `Wind speed (m/s)` + `Visibility (10m)` + `Solar Radiation (MJ/m2)` + 
##     `Rainfall (mm)` + `Snowfall (cm)` + Seasons + Holiday + `Functioning Day` + 
##     `Day of Week` + Month + Year + `Day Category` + `Rush Hour` + 
##     Precipitation

formula(stepwise_model)

## `Rented Bike Count` ~ `Temperature (°C)` + Hour + `Rush Hour` + 
##     `Functioning Day` + `Humidity (%)` + Precipitation + Month + 
##     Year + `Rainfall (mm)` + Seasons + `Snowfall (cm)` + `Day Category` + 
##     Holiday + `Solar Radiation (MJ/m2)` + `Wind speed (m/s)` + 
##     `Visibility (10m)` + `Day of Week`

Interpretasi Hasil Stepwise Regression:

Perbandingan Nilai AIC:

Nilai AIC dari Forward, Backward, dan Bidirectional Stepwise menghasilkan nilai yang sama, menunjukkan bahwa ketiga metode berhasil menemukan kombinasi variabel yang optimal dan konsisten. Kesamaan nilai AIC ini mengindikasikan bahwa solusi yang ditemukan bersifat stabil dan tidak bergantung pada arah seleksi yang digunakan.

Variabel yang Terpilih:

Berdasarkan output formula dari ketiga model stepwise, variabel-variabel yang terpilih, yaitu:
- Hour
- Temperature (°C)
- Humidity (%)
- Wind speed (m/s)
- Visibility (10m)
- Solar Radiation (MJ/m2)
- Rainfall (mm)
- Snowfall (cm)
- Seasons
- Holiday
- Functioning Day
- Day of Week
- Month
- Year
- Day Category
- Rush Hour
- Precipitation

Variabel-variavel ini merupakan prediktor yang paling berkontribusi signifikan terhadap prediksi jumlah penyewaan sepeda. Variabel ini dipilih karena secara konsisten menurunkan nilai AIC, yang menandakan bahwa informasi yang dibawa dapat meningkatkan kemampuan prediktif model.

Variabel Dew point temperature (°C) dieliminasi oleh model stepwise karena informasi yang dikandungnya sudah sepenuhnya terwakili oleh variabel Temperature (°C) dan Humidity (%) yang tetap dipertahankan dalam model. Ketika kedua variabel tersebut sudah ada dalam model, kehadiran Dew point temperature (°C) tidak lagi memberikan informasi tambahan yang baru, sehingga penambahan variabel ini tidak menurunkan nilai AIC, melainkan hanya menambah kompleksitas model tanpa peningkatan kemampuan prediktif yang berarti. Hasil ini juga sejalan dengan temuan pada tahap Filter Method (VIF) sebelumnya, di mana Dew point temperature (°C) teridentifikasi memiliki multikolinearitas tinggi dengan Temperature (°C) dan Humidity (%). Dengan demikian, eliminasi variabel ini oleh kedua metode secara konsisten memperkuat kesimpulan bahwa Dew point temperature (°C) merupakan variabel yang redundan dalam model.

Sementara itu, variabel Day (tanggal dalam bulan, bernilai 1–31) dieliminasi karena tidak memiliki hubungan yang bermakna terhadap jumlah penyewaan sepeda. Berbeda dengan Month yang mencerminkan pola musiman atau Day of Week yang mencerminkan pola aktivitas mingguan, tanggal dalam satu bulan tidak merepresentasikan pola perilaku pengguna yang konsisten, tidak ada alasan mengapa tanggal 1 atau tanggal 15 menghasilkan jumlah penyewaan yang lebih tinggi atau lebih rendah dibandingkan tanggal lainnya. Oleh karena itu, variabel Day dinilai tidak memberikan informasi prediktif yang cukup, sehingga keikutsertaannya justru menambah kompleksitas model tanpa diimbangi peningkatan akurasi yang berarti.

5.1 Persiapan Data PCA

Sebelum menjalankan PCA, data perlu distandarisasi menggunakan Z-Score karena PCA sangat sensitif terhadap skala. Variabel dengan satuan besar akan mendominasi komponen utama jika tidak distandarisasi terlebih dahulu.

library(dplyr)

# Pakai df_num_2 hasil feature selection dari bab 4
df_pca <- df_num_2 %>%
  select(-`Rented Bike Count`)

# Standardisasi Z-Score
df_scaled <- scale(df_pca)

# PCA
pca_result <- prcomp(df_scaled, center = TRUE, scale. = TRUE)
summary(pca_result)

## Importance of components:
##                           PC1    PC2    PC3     PC4     PC5     PC6     PC7
## Standard deviation     1.5996 1.4689 1.3439 1.23200 1.17465 1.05996 1.03472
## Proportion of Variance 0.1421 0.1199 0.1003 0.08432 0.07666 0.06242 0.05948
## Cumulative Proportion  0.1421 0.2620 0.3623 0.44667 0.52332 0.58574 0.64522
##                            PC8     PC9    PC10    PC11    PC12    PC13    PC14
## Standard deviation     0.99119 0.97130 0.92640 0.91136 0.85938 0.77120 0.68081
## Proportion of Variance 0.05458 0.05241 0.04768 0.04614 0.04103 0.03304 0.02575
## Cumulative Proportion  0.69980 0.75222 0.79989 0.84604 0.88707 0.92011 0.94586
##                           PC15    PC16    PC17    PC18
## Standard deviation     0.60643 0.46684 0.46002 0.42098
## Proportion of Variance 0.02043 0.01211 0.01176 0.00985
## Cumulative Proportion  0.96629 0.97840 0.99015 1.00000

5.2 Penentuan Jumlah Komponen Utama

Penentuan jumlah komponen utama yang optimal dilakukan menggunakan dua kriteria, yaitu Kaiser Criterion (eigenvalue > 1) dan Cumulative Variance (minimal 70–90% variasi terjelaskan).

# Scree Plot
variance <- pca_result$sdev^2 / sum(pca_result$sdev^2) * 100
barplot(variance,
        names.arg = paste0("PC", 1:length(variance)),
        col = "#5B84B1",
        main = "Scree Plot - Proporsi Variasi Tiap Komponen PCA",
        ylab = "Proporsi Variasi (%)",
        xlab = "Komponen Utama")

# Cumulative Variance
cumvar <- cumsum(pca_result$sdev^2 / sum(pca_result$sdev^2)) * 100
plot(cumvar, type = "b",
     col = "#5B84B1",
     pch = 16,
     main = "Cumulative Variance PCA",
     xlab = "Jumlah Komponen Utama",
     ylab = "Cumulative Variance (%)")
abline(h = 80, col = "red", lty = 2)
abline(h = 90, col = "darkgreen", lty = 2)
legend("bottomright", legend = c("80%", "90%"),
       col = c("red", "darkgreen"), lty = 2)

5.3 Loading Factor

Loading factor menunjukkan seberapa besar kontribusi setiap variabel asli terhadap pembentukan masing-masing komponen utama. Nilai loading mendekati 1 atau -1 menunjukkan hubungan yang kuat, sedangkan mendekati 0 menunjukkan hubungan yang lemah.

# Tabel Loading Factor 5 PC pertama
loading_df <- as.data.frame(round(pca_result$rotation[, 1:5], 3))
knitr::kable(loading_df, caption = "Loading Factor 5 Komponen Utama Pertama")

# Kontribusi variabel ke PC1
contrib_pc1 <- abs(pca_result$rotation[,1]) / sum(abs(pca_result$rotation[,1])) * 100
contrib_df1 <- data.frame(
  Variabel = names(contrib_pc1),
  Kontribusi = contrib_pc1) %>%
  arrange(desc(Kontribusi)) %>%
  head(10)

ggplot(contrib_df1, aes(x = reorder(Variabel, Kontribusi), y = Kontribusi)) +
  geom_col(fill = "#5B84B1") +
  coord_flip() +
  labs(title = "Kontribusi Variabel ke PC1",
       x = "", y = "Kontribusi (%)") +
  theme_minimal()

Berdasarkan grafik kontribusi variabel terhadap PC1, terdapat tiga variabel dengan kontribusi tertinggi yaitu Humidity (%), Solar Radiation (MJ/m2), dan Precipitation, masing-masing berkontribusi sekitar 11–12%. Ketiganya merupakan indikator kondisi atmosfer yang saling berkaitan — kelembaban tinggi, radiasi matahari rendah, dan adanya curah hujan/salju cenderung terjadi bersamaan pada kondisi cuaca buruk. Variabel Visibility (10m) dan Snowfall (cm) berada di posisi berikutnya dengan kontribusi sekitar 9%, yang juga mencerminkan kondisi cuaca jarak pandang rendah umumnya terjadi saat cuaca buruk atau bersalju. Variabel Year, Temperature (°C), Wind speed (m/s), Seasons, dan Hour masing-masing berkontribusi sekitar 6–8%, menunjukkan peran yang lebih merata sebagai faktor pendukung. Secara keseluruhan, PC1 didominasi oleh variabel-variabel yang merepresentasikan kondisi cuaca dan atmosfer, sehingga PC1 dapat diinterpretasikan sebagai komponen yang menangkap “intensitas kondisi cuaca” di mana nilai PC1 yang tinggi mencerminkan kondisi cuaca buruk (lembab, hujan/salju, visibilitas rendah, radiasi rendah), dan sebaliknya nilai PC1 rendah mencerminkan cuaca cerah dan kondusif untuk bersepeda.

# Kontribusi variabel ke PC2
contrib_pc2 <- abs(pca_result$rotation[,2]) / sum(abs(pca_result$rotation[,2])) * 100
contrib_df2 <- data.frame(
  Variabel = names(contrib_pc2),
  Kontribusi = contrib_pc2) %>%
  arrange(desc(Kontribusi)) %>%
  head(10)

ggplot(contrib_df2, aes(x = reorder(Variabel, Kontribusi), y = Kontribusi)) +
  geom_col(fill = "#e67e22") +
  coord_flip() +
  labs(title = "Kontribusi Variabel ke PC2",
       x = "", y = "Kontribusi (%)") +
  theme_minimal()

Berdasarkan grafik kontribusi variabel terhadap PC2 menunjukkan bahwa Variabel dengan kontribusi tertinggi pada PC2 adalah Seasons (~16%), jauh melampaui variabel lainnya. Diikuti oleh Temperature (°C), Humidity (%), dan Year yang masing-masing berkontribusi sekitar 13–14%. Keempat variabel teratas ini sangat erat kaitannya dengan siklus musiman musim secara langsung menentukan suhu dan kelembaban, sementara Year mencerminkan tren jangka panjang antar tahun pengamatan. Variabel Visibility (10m) berkontribusi sekitar 9%, sementara Wind speed (m/s), Rainfall (mm), Snowfall (cm), Holiday, dan Month masing-masing berkontribusi lebih kecil di kisaran 4–7%. Secara keseluruhan, PC2 dapat diinterpretasikan sebagai komponen yang menangkap “siklus musiman dan termal” — di mana PC2 membedakan kondisi antar musim secara sistematis, khususnya perbedaan antara musim panas yang hangat dan lembab dengan musim dingin yang dingin dan kering.

5.4 Visualisasi Biplot

Biplot menampilkan dua informasi sekaligus, yaitu arah dan kekuatan hubungan antar variabel melalui panah merah. Biplot ini hanya menampilkan arah variabel tanpa titik observasi agar lebih mudah dibaca mengingat dataset memiliki 8.760 observasi.

# Biplot lebih rapi dengan ggplot
loading_plot <- as.data.frame(pca_result$rotation[, 1:2])
loading_plot$Variabel <- rownames(loading_plot)

ggplot(loading_plot, aes(x = PC1, y = PC2, label = Variabel)) +
  geom_segment(aes(x = 0, y = 0, xend = PC1, yend = PC2),
               arrow = arrow(length = unit(0.3, "cm")),
               color = "#e74c3c") +
  geom_text(hjust = 0.5, vjust = -0.5, size = 3.5, color = "#e74c3c") +
  geom_hline(yintercept = 0, linetype = "dashed", color = "gray") +
  geom_vline(xintercept = 0, linetype = "dashed", color = "gray") +
  labs(title = "Biplot PCA - Seoul Bike Sharing",
       x = "PC1", y = "PC2") +
  theme_minimal()

Berdasarkan biplot di atas didapatkan kesimpulan:

Variabel yang Saling Berdekatan (Berkorelasi) - Temperature (°C) dan Year: panah berhimpitan ke arah bawah kanan, menunjukkan keduanya berkorelasi positif. - Humidity (%) dan Precipitation: panah berdekatan ke kiri bawah, menunjukkan keduanya berkorelasi positif sebagai indikator cuaca basah. - Visibility dan Wind speed: panah berdekatan ke kanan atas, berkorelasi positif pada kondisi cuaca cerah.

Variabel yang Berlawanan Arah (Berkorelasi Negatif) - Solar Radiation vs Humidity/Precipitation: panah berlawanan arah (kanan vs kiri), menunjukkan bahwa radiasi matahari tinggi terjadi saat kelembaban dan curah hujan rendah, dan sebaliknya. - Seasons vs Temperature: panah berlawanan arah pada PC2, mengonfirmasi bahwa musim dingin (Seasons = Winter) identik dengan suhu rendah.

5.5 Ringkasan

# Tabel Ringkasan Proporsi Variasi
prop_var <- data.frame(
  PC = paste0("PC", 1:length(pca_result$sdev)),
  Eigenvalue = round(pca_result$sdev^2, 3),
  PropVariance = round(pca_result$sdev^2 / sum(pca_result$sdev^2) * 100, 2),
  CumVariance = round(cumsum(pca_result$sdev^2 / sum(pca_result$sdev^2)) * 100, 2)
)

knitr::kable(prop_var,
             col.names = c("Komponen", "Eigenvalue", "Proporsi Variasi (%)", "Kumulatif (%)"),
             caption = "Ringkasan Proporsi Variasi Komponen PCA")

Kaiser Criterion (Eigenvalue > 1): Terdapat 7 komponen utama yang memiliki eigenvalue > 1, yaitu PC1 hingga PC7. Artinya, ketujuh komponen ini dianggap signifikan karena masing-masing mampu menjelaskan variasi yang lebih besar dibandingkan satu variabel tunggal yang distandarisasi.

Proporsi Variasi per Komponen: PC1 menjadi komponen paling informatif dengan menjelaskan 14,21% variasi data, diikuti PC2 (11,99%) dan PC3 (10,03%). Penurunan proporsi variasi berjalan cukup merata dari PC1 hingga PC7, yang menunjukkan bahwa informasi dalam dataset ini tersebar secara relatif merata ke banyak komponen. Hal ini terjadi karena 17 variabel input memiliki karakteristik yang sangat beragam mulai dari kondisi cuaca, waktu harian, hingga kalender musiman.

Dengan hanya 7 komponen utama (berdasarkan Kaiser Criterion), sudah dapat merepresentasikan 64,52% informasi dari 17 variabel asli artinya reduksi dimensi dari 17 variabel menjadi 7 komponen berhasil mempertahankan hampir dua pertiga informasi data.