1. Persiapan Data dan Pemuatan Package

Sebelum memulai analisis multivariat, langkah pertama adalah memuat semua package R yang diperlukan untuk pengujian MANOVA, MANCOVA, uji asumsi klasik, serta visualisasi tingkat lanjut.

# Memuat library utama
library(MVN)
library(heplots)
library(car)
library(dplyr)

1.1 Memuat Dataset

Dataset yang digunakan adalah Bike Sharing Dataset (Capital Bikeshare) yang diperoleh dari UCI Machine Learning Repository. Pastikan file day.csv berada di dalam working directory Anda.

# Memuat dataset lokal
bike_data <- read.csv("day.csv")

1.2 Pra-Pemrosesan dan Transformasi Data

Di sini kita akan mengubah variabel prediktor kategorik menjadi factor agar dibaca dengan benar oleh fungsi pemodelan R. Selain itu, transformasi logaritma natural (\(log(x + 1)\)) dilakukan pada variabel dependen untuk mengurangi efek skewness (kemencengan data) dan menjaga pemenuhan asumsi normalitas.

# Konversi variabel menjadi faktor
bike_data$season <- as.factor(bike_data$season)
bike_data$weathersit <- as.factor(bike_data$weathersit)

# Transformasi logaritmik untuk variabel dependen (Y)
bike_data$log_casual <- log(bike_data$casual + 1)
bike_data$log_registered <- log(bike_data$registered + 1)

# PENTING: Membuat matriks Y dengan nama kolom yang spesifik dan bersih 
# agar tidak memicu error mfcol=c(nv,nv) pada visualisasi heplot
Y_matrix <- cbind(
  log_casual = bike_data$log_casual, 
  log_registered = bike_data$log_registered
)

2. Statistika Deskriptif

Analisis deskriptif digunakan untuk melihat karakteristik dasar dari data, termasuk nilai pemusatan data (mean, median) dan ukuran penyebaran data (standar deviasi).

# Menampilkan ringkasan statistik (Min, Q1, Median, Mean, Q3, Max)
summary(bike_data[, c("casual", "registered", "temp")])
##      casual         registered        temp        
##  Min.   :   2.0   Min.   :  20   Min.   :0.05913  
##  1st Qu.: 315.5   1st Qu.:2497   1st Qu.:0.33708  
##  Median : 713.0   Median :3662   Median :0.49833  
##  Mean   : 848.2   Mean   :3656   Mean   :0.49538  
##  3rd Qu.:1096.0   3rd Qu.:4776   3rd Qu.:0.65542  
##  Max.   :3410.0   Max.   :6946   Max.   :0.86167
# Menghitung Standar Deviasi (SD) untuk masing-masing variabel utama
cat("Standar Deviasi Casual (Y1)     : ", sd(bike_data$casual), "\n")
## Standar Deviasi Casual (Y1)     :  686.6225
cat("Standar Deviasi Registered (Y2) : ", sd(bike_data$registered), "\n")
## Standar Deviasi Registered (Y2) :  1560.256
cat("Standar Deviasi Temp (X3)       : ", sd(bike_data$temp), "\n")
## Standar Deviasi Temp (X3)       :  0.183051

3. Pengujian Asumsi Klasik Multivariat

Validitas kesimpulan dari pengujian MANOVA dan MANCOVA sangat bergantung pada dipenuhinya asumsi-asumsi klasik multivariat berikut.

3.1 Uji Normalitas Multivariat (Mardia’s Test)

Uji Mardia digunakan untuk menguji apakah variasi kombinasi linear dari variabel dependen berdistribusi normal multivariat melalui parameter skewness dan kurtosis.

# Memanggil fungsi mvn(). Secara default fungsi ini otomatis menjalankan Uji Mardia, 
# sehingga kita tidak perlu menulis parameter tambahan yang rawan memicu error antar-versi
mardia_test <- mvn(data = bike_data[, c("log_casual", "log_registered")])
print(mardia_test$multivariateNormality)
## NULL

3.2 Uji Homogenitas Matriks Kovarians (Box’s M Test)

Uji ini bertujuan untuk memastikan bahwa matriks varians-kovarians dari variabel dependen bernilai sama (homogen) di setiap kelompok/kombinasi perlakuan.

box_m_result <- boxM(Y_matrix ~ season * weathersit, data = bike_data)
print(box_m_result)
## 
##  Box's M-test for Homogeneity of Covariance Matrices 
## 
## data:  bike_data 
## Chi-Sq (approx.) = 277.8829, df = 33, p-value = < 2.2e-16

3.3 Uji Homogenitas Kemiringan Regresi (Asumsi MANCOVA)

Sebelum menjalankan MANCOVA, kemiringan regresi (regression slopes) antara variabel kovariat (temp) dan variabel dependen harus sama di seluruh kelompok. Hal ini ditunjukkan oleh tidak signifikannya interaksi antara kovariat dan faktor utama.

model_asumsi <- manova(Y_matrix ~ temp * season * weathersit, data = bike_data)
summary(model_asumsi, test = "Pillai")
##                         Df  Pillai approx F num Df den Df    Pr(>F)    
## temp                     1 0.61244   557.82      2    706 < 2.2e-16 ***
## season                   3 0.21265    28.04      6   1414 < 2.2e-16 ***
## weathersit               2 0.29127    60.26      4   1414 < 2.2e-16 ***
## temp:season              3 0.14466    18.38      6   1414 < 2.2e-16 ***
## temp:weathersit          2 0.00590     1.05      4   1414 0.3822299    
## season:weathersit        6 0.02964     1.77     12   1414 0.0477434 *  
## temp:season:weathersit   6 0.04884     2.95     12   1414 0.0004551 ***
## Residuals              707                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4. Pemodelan Utama Analisis Multivariat

Setelah memeriksa asumsi, kita dapat menjalankan pengujian signifikansi kelompok secara multivariat.

4.1 Hasil Uji MANOVA (Tanpa Kovariat)

Pengujian MANOVA digunakan untuk menguji apakah terdapat perbedaan rata-rata kombinasi variabel dependen (casual dan registered) yang signifikan berdasarkan faktor musim (season), cuaca (weathersit), maupun interaksi keduanya.

model_manova <- manova(Y_matrix ~ season * weathersit, data = bike_data)
summary(model_manova, test = "Pillai")
##                    Df  Pillai approx F num Df den Df  Pr(>F)    
## season              3 0.51586   83.304      6   1438 < 2e-16 ***
## weathersit          2 0.23935   48.871      4   1438 < 2e-16 ***
## season:weathersit   6 0.02654    1.612     12   1438 0.08208 .  
## Residuals         719                                           
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.2 Hasil Uji MANCOVA (Dengan Kovariat Suhu)

Pengujian MANCOVA memperluas MANOVA dengan memasukkan variabel kontinu temp (suhu) sebagai kovariat untuk mengontrol pengaruh lingkungan luar secara statistik.

model_mancova <- manova(Y_matrix ~ temp + season * weathersit, data = bike_data)
summary(model_mancova, test = "Pillai")
##                    Df  Pillai approx F num Df den Df    Pr(>F)    
## temp                1 0.57435   483.75      2    717 < 2.2e-16 ***
## season              3 0.18920    25.01      6   1436 < 2.2e-16 ***
## weathersit          2 0.25888    53.38      4   1436 < 2.2e-16 ***
## season:weathersit   6 0.04384     2.68     12   1436  0.001413 ** 
## Residuals         718                                             
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

5. Visualisasi Analisis Multivariat

Visualisasi sangat penting untuk membantu menginterpretasikan hasil uji statistik dan melihat pola data secara langsung.

5.1 Visualisasi Uji Asumsi

A. Chi-Square Q-Q Plot (Normalitas Multivariat)

Plot ini menggambarkan jarak kuadrat Mahalanobis terhadap distribusi Chi-Square untuk mendeteksi deviasi dari normalitas multivariat. Kita menggunakan fungsi cqplot dari package heplots yang sangat stabil untuk membaca objek MANOVA.

# Visualisasi Normalitas Multivariat
cqplot(model_manova, main = "Chi-Square Q-Q Plot (Normalitas Multivariat)")

B. Plot Box’s M (Homogenitas Matriks Kovarians)

Visualisasi berbasis dot plot untuk melihat perbandingan log determinan matriks kovarians antar kelompok kombinasi perlakuan.

plot(box_m_result, main = "Plot Box's M: Homogenitas Matriks Kovarians")

C. Elips Sebaran (Covariance Ellipses)

Penggambaran tingkat kekuatan hubungan varians antar-variabel dependen yang dikelompokkan berdasarkan musim (season).

covEllipses(Y_matrix ~ season, data = bike_data, 
            fill = TRUE, pooled = FALSE, 
            main = "Elips Kovarians Berdasarkan Musim")

5.2 Visualisasi Hasil MANOVA & MANCOVA (HE Plots)

Hypothesis-Error (HE) Plots memproyeksikan kekuatan bukti pengujian menggunakan bentuk geometri elips untuk mewakili variabilitas akibat Efek Hipotesis (H) dibandingkan dengan Variabilitas Error (E).

A. HE Plot pada Model MANOVA

heplot(model_manova, 
       fill = TRUE, fill.alpha = 0.1, 
       main = "HE Plot MANOVA: Efek Musim & Cuaca",
       xlab = "Log Pengguna Kasual (Y1)", 
       ylab = "Log Pengguna Terdaftar (Y2)")

B. HE Plot pada Model MANCOVA (Intervensi Kovariat Suhu)

heplot(model_mancova, 
       fill = TRUE, fill.alpha = 0.1, 
       main = "HE Plot MANCOVA: Intervensi Suhu Lingkungan",
       xlab = "Log Pengguna Kasual (Y1)", 
       ylab = "Log Pengguna Terdaftar (Y2)")