2.2.1 Variabel Numerik

num_vars <- df %>% select(-ocean_proximity)
summary(num_vars)

   longitude         latitude     housing_median_age  total_rooms   
 Min.   :-124.3   Min.   :32.54   Min.   : 1.00      Min.   :    2  
 1st Qu.:-121.8   1st Qu.:33.93   1st Qu.:18.00      1st Qu.: 1448  
 Median :-118.5   Median :34.26   Median :29.00      Median : 2127  
 Mean   :-119.6   Mean   :35.63   Mean   :28.64      Mean   : 2636  
 3rd Qu.:-118.0   3rd Qu.:37.71   3rd Qu.:37.00      3rd Qu.: 3148  
 Max.   :-114.3   Max.   :41.95   Max.   :52.00      Max.   :39320  
                                                                    
 total_bedrooms     population      households     median_income    
 Min.   :   1.0   Min.   :    3   Min.   :   1.0   Min.   : 0.4999  
 1st Qu.: 296.0   1st Qu.:  787   1st Qu.: 280.0   1st Qu.: 2.5634  
 Median : 435.0   Median : 1166   Median : 409.0   Median : 3.5348  
 Mean   : 537.9   Mean   : 1425   Mean   : 499.5   Mean   : 3.8707  
 3rd Qu.: 647.0   3rd Qu.: 1725   3rd Qu.: 605.0   3rd Qu.: 4.7432  
 Max.   :6445.0   Max.   :35682   Max.   :6082.0   Max.   :15.0001  
 NA's   :207                                                        
 median_house_value
 Min.   : 14999    
 1st Qu.:119600    
 Median :179700    
 Mean   :206856    
 3rd Qu.:264725    
 Max.   :500001    
                   

desc_table <- num_vars %>%
  summarise(across(everything(), list(
    N        = ~sum(!is.na(.)),
    Mean     = ~round(mean(., na.rm = TRUE), 2),
    Median   = ~round(median(., na.rm = TRUE), 2),
    SD       = ~round(sd(., na.rm = TRUE), 2),
    Min      = ~round(min(., na.rm = TRUE), 2),
    Max      = ~round(max(., na.rm = TRUE), 2),
    Skewness = ~round(mean(((. - mean(.,na.rm=T))/sd(.,na.rm=T))^3,na.rm=T),3)
  ))) %>%
  pivot_longer(everything(),
               names_to  = c("Variabel",".value"),
               names_sep = "_(?=[^_]+$)") %>%
  arrange(desc(abs(Skewness)))

kable(desc_table,
      caption = "Statistik Deskriptif Variabel Numerik",
      format.args = list(big.mark = ",")) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(8, bold = TRUE,
              color = ifelse(abs(desc_table$Skewness) > 1, "red", "darkgreen")) %>%
  row_spec(0, bold = TRUE)

Berdasarkan statistik deskriptif yang diperoleh, variabel population memiliki nilai skewness tertinggi (4,935), menunjukkan distribusi yang sangat menceng ke kanan (strong positive skewness). Hal ini terlihat dari rata-rata (1.425,48) yang lebih besar daripada median (1.166), serta adanya nilai maksimum yang sangat tinggi (35.682). Kondisi tersebut mengindikasikan bahwa sebagian besar blok sensus memiliki jumlah penduduk relatif rendah hingga sedang, tetapi terdapat beberapa blok dengan jumlah penduduk yang sangat besar sehingga menarik distribusi ke arah kanan.

Variabel total_rooms, total_bedrooms, dan households juga menunjukkan skewness positif yang tinggi, masing-masing sebesar 4,147; 3,459; dan 3,410. Nilai rata-rata ketiga variabel tersebut lebih besar daripada mediannya, yang mengindikasikan bahwa sebagian besar blok sensus memiliki jumlah ruangan, kamar tidur, dan rumah tangga yang relatif moderat, namun terdapat sejumlah kecil blok dengan ukuran yang sangat besar sehingga menyebabkan distribusi tidak simetris. Selain itu, variabel total_bedrooms memiliki 207 nilai hilang (missing values) yang perlu ditangani sebelum dilakukan analisis lanjutan.

Variabel median_income memiliki skewness positif sedang (1,646), yang menunjukkan distribusi masih cenderung menceng ke kanan. Nilai median sebesar 3,53 berarti pendapatan median rumah tangga sekitar US$35.300 per tahun karena variabel ini dinyatakan dalam satuan puluhan ribu dolar Amerika. Adanya perbedaan antara rata-rata (3,87) dan median (3,53) menunjukkan bahwa terdapat sejumlah wilayah dengan tingkat pendapatan yang jauh lebih tinggi dibandingkan mayoritas wilayah lainnya.

Sementara itu, variabel geografis latitude dan longitude memiliki nilai skewness yang relatif kecil, masing-masing sebesar 0,466 dan -0,298. Hal ini menunjukkan bahwa distribusi kedua variabel cenderung mendekati simetris meskipun terdapat sedikit kecenderungan ke kanan pada latitude dan ke kiri pada longitude. Rentang nilai yang cukup luas pada kedua variabel tersebut mencerminkan cakupan wilayah pengamatan yang tersebar di berbagai lokasi di negara bagian California.

Variabel housing_median_age memiliki skewness yang sangat dekat dengan nol (0,060), sehingga distribusinya dapat dianggap hampir simetris. Hal ini juga didukung oleh nilai mean (28,64) dan median (29,00) yang hampir sama. Dengan demikian, usia median rumah pada dataset ini relatif tersebar secara seimbang di sekitar nilai tengahnya.

Sebagai variabel target, median_house_value memiliki skewness positif sebesar 0,978. Nilai rata-rata (US$206.856) yang lebih tinggi daripada median (US$179.700) menunjukkan adanya kecenderungan distribusi ke arah kanan akibat keberadaan rumah-rumah bernilai tinggi. Namun, perlu diperhatikan bahwa nilai maksimum sebesar US$500.001 mengindikasikan adanya censoring atau batas atas pada data asli California Housing, sehingga distribusi harga rumah kemungkinan terpotong pada nilai maksimum tersebut.

Secara keseluruhan, sebagian besar variabel yang berkaitan dengan ukuran dan karakteristik blok sensus menunjukkan distribusi yang tidak simetris dengan kecenderungan menceng ke kanan. Kondisi ini mengindikasikan adanya beberapa observasi bernilai sangat besar yang berpotensi memengaruhi analisis. Oleh karena itu, tahap prapengolahan data seperti penanganan nilai hilang, standarisasi variabel, dan eksplorasi lebih lanjut terhadap pencilan perlu dilakukan sebelum menerapkan teknik reduksi dimensi.

2.2.2 Variabel Kategorik

ocean_tbl <- df %>%
  count(ocean_proximity) %>%
  mutate(Persentase = paste0(round(n/sum(n)*100, 1), "%"),
         Rata_Harga = dollar(round(tapply(df$median_house_value,
                                          df$ocean_proximity, mean)[ocean_proximity], 0))) %>%
  rename(Kategori = ocean_proximity, Frekuensi = n)

kable(ocean_tbl, caption = "Distribusi Frekuensi dan Rata-rata Harga per Kategori") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE)

Berdasarkan statistik deskriptif yang diperoleh, Kategori <1H OCEAN (kurang dari 1 jam dari laut) mendominasi dataset dengan 9.136 blok sensus atau sekitar 44,3% dari total observasi. Kategori INLAND menempati urutan kedua dengan 6.551 blok (31,7%), diikuti oleh NEAR OCEAN (12,9%) dan NEAR BAY (11,1%). Sementara itu, kategori ISLAND hanya memiliki 5 observasi sehingga representasinya sangat terbatas dan hasil yang diperoleh pada kategori ini perlu ditafsirkan dengan hati-hati.

Berdasarkan rata-rata harga rumah, terlihat bahwa wilayah yang berdekatan dengan perairan cenderung memiliki nilai rumah yang lebih tinggi dibandingkan wilayah INLAND. Kategori INLAND memiliki rata-rata harga rumah terendah, yaitu sebesar US$124.805. Sebaliknya, kategori NEAR BAY memiliki rata-rata harga tertinggi di antara kategori dengan jumlah observasi yang memadai, yaitu sebesar US$259.212, diikuti oleh NEAR OCEAN (US$249.434) dan <1H OCEAN (US$240.084). Meskipun kategori ISLAND memiliki rata-rata harga rumah tertinggi (US$380.440), jumlah observasinya yang sangat sedikit menyebabkan nilai tersebut kurang representatif untuk menggambarkan kondisi populasi secara umum.

Temuan ini menunjukkan bahwa lokasi geografis, khususnya kedekatan dengan wilayah perairan, berpotensi menjadi faktor penting yang memengaruhi nilai rumah di California. Oleh karena itu, variabel ocean_proximity layak dipertimbangkan sebagai salah satu variabel yang berkontribusi dalam menjelaskan variasi harga rumah pada dataset ini.

2.3.1 Histogram Semua Variabel

num_long <- num_vars %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai")

ggplot(num_long, aes(x = Nilai, fill = Variabel)) +
  geom_histogram(bins = 40, color = "white", alpha = 0.85) +
  facet_wrap(~Variabel, scales = "free", ncol = 3) +
  scale_fill_viridis_d(option = "D") +
  scale_x_continuous(labels = comma) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Distribusi Setiap Variabel Numerik",
    subtitle = "California Housing Dataset",
    x = NULL, y = "Frekuensi"
  ) +
  theme_bw(base_size = 12) +
  theme(
    plot.title       = element_text(face = "bold", size = 15),
    strip.text       = element_text(face = "bold"),
    legend.position  = "none",
    panel.grid.minor = element_blank()
  )

Berdasarkan histrogram di atas, didapatkan interpretasi sebagai berikut.

• housing_median_age: Distribusi paling mendekati seragam (uniform) di antara semua variabel, dengan sebaran yang relatif merata antara 1 hingga 52 tahun. Nilai skewness yang sangat kecil (0,06) menunjukkan bahwa distribusi variabel ini hampir simetris. Selain itu, terlihat beberapa puncak frekuensi pada rentang umur tertentu, terutama di sekitar 35???37 tahun dan mendekati 52 tahun, yang menunjukkan adanya konsentrasi observasi pada kelompok umur tersebut.

• total_rooms: Distribusi sangat menceng ke kanan (strong positive skewness) dengan nilai skewness sebesar 4,15. Sebagian besar blok sensus memiliki total kamar antara 1.000 sampai 5.000, namun terdapat ekor distribusi yang panjang hingga mencapai 39.320 kamar. Perbedaan antara median (2.127) dan mean (2.636) juga mengonfirmasi adanya pengaruh observasi bernilai sangat besar.

• total_bedrooms: Menunjukkan pola yang serupa dengan total_rooms (skewness = 3,46), yaitu distribusi yang menceng ke kanan dengan median 435 kamar tidur dan nilai maksimum mencapai 6.445 kamar tidur. Pola ini konsisten karena total_bedrooms merupakan bagian dari total_rooms.

• population: Merupakan variabel dengan distribusi paling menceng ke kanan (skewness = 4,94). Sebagian besar blok sensus memiliki populasi sekitar 500 sampai 2.500 jiwa, namun terdapat beberapa blok dengan populasi yang sangat besar hingga mencapai 35.682 jiwa. Kehadiran observasi-observasi ekstrem tersebut menyebabkan terbentuknya ekor distribusi yang panjang di sisi kanan.

• households: Memiliki pola distribusi yang mirip dengan population (skewness = 3,41), yang menunjukkan bahwa sebagian besar blok sensus memiliki jumlah rumah tangga relatif rendah hingga sedang, namun terdapat sejumlah kecil blok dengan jumlah rumah tangga yang sangat besar. Hal ini mengindikasikan adanya variasi ukuran blok sensus yang cukup tinggi.

• median_income: Menunjukkan skewness positif yang cukup kuat (skewness = 1,65), sehingga distribusinya masih cenderung menceng ke kanan dan belum dapat dianggap normal. Sebagian besar observasi terkonsentrasi pada tingkat pendapatan menengah, dengan sejumlah kecil wilayah yang memiliki pendapatan jauh lebih tinggi dibandingkan mayoritas wilayah lainnya. Selain itu, terlihat adanya penumpukan kecil pada nilai maksimum 15,0 yang mengindikasikan adanya batas atas (capping) pada data asli.

• median_house_value: Memiliki distribusi menceng ke kanan dengan skewness moderat (0,98). Sebagian besar nilai rumah berada pada kisaran US$100.000???US$300.000, namun masih terdapat sejumlah rumah dengan nilai yang jauh lebih tinggi. Terlihat pula lonjakan frekuensi yang cukup jelas pada nilai maksimum US$500.001 yang menunjukkan adanya capping atau batas atas pada data asli. Kondisi ini perlu diperhatikan dalam analisis lanjutan karena dapat memengaruhi hasil pemodelan.

• longitude dan latitude: Kedua variabel geografis ini menunjukkan pola distribusi yang berbeda dari variabel numerik lainnya. Distribusinya bersifat multimodal, ditandai dengan adanya beberapa puncak frekuensi yang mencerminkan konsentrasi observasi pada wilayah geografis tertentu di California. Pada variabel latitude, puncak utama terlihat di sekitar 34??° dan 37???38??°, yang menunjukkan wilayah dengan kepadatan observasi yang tinggi. Oleh karena itu, pola distribusi kedua variabel ini lebih merefleksikan karakteristik spasial daripada distribusi statistik konvensional.

2.3.2 Distribusi Variabel Target

ggplot(df, aes(x = median_house_value)) +
  geom_histogram(bins = 50, fill = "#2c7fb8", color = "white", alpha = 0.85) +
  geom_vline(xintercept = mean(df$median_house_value, na.rm=TRUE),
             color = "red", linetype = "dashed", size = 1.2) +
  geom_vline(xintercept = median(df$median_house_value, na.rm=TRUE),
             color = "darkorange", linetype = "dashed", size = 1.2) +
  annotate("text", x = mean(df$median_house_value)+15000, y = 1100,
           label = paste0("Mean\n$", comma(round(mean(df$median_house_value)))),
           color = "red", fontface = "bold", size = 3.5) +
  annotate("text", x = median(df$median_house_value)-55000, y = 1100,
           label = paste0("Median\n$", comma(round(median(df$median_house_value)))),
           color = "darkorange", fontface = "bold", size = 3.5) +
  scale_x_continuous(labels = dollar_format(prefix = "$", big.mark = ",")) +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Distribusi Median House Value (Variabel Target)",
    subtitle = "Garis merah = Mean | Garis oranye = Median | Puncak di $500.001 = data ter-cap",
    x = "Median House Value (USD)", y = "Frekuensi"
  ) +
  theme_bw(base_size = 12) +
  theme(plot.title = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 10))

Berdasarkan histogram di atas, dapat diinterpretasikan bahwa distribusi median_house_value menunjukkan pola right-skewed dengan nilai mean (US$206.856) yang lebih tinggi daripada median (US$179.700), dengan selisih sekitar US$27.000. Kondisi ini menunjukkan bahwa sebagian besar blok sensus memiliki harga rumah di bawah rata-rata, sementara sejumlah kecil rumah dengan nilai sangat tinggi menarik rata-rata ke arah kanan.

Selain itu, terlihat lonjakan frekuensi yang sangat mencolok pada nilai sekitar US$500.001 di ujung kanan histogram. Fenomena ini mengindikasikan adanya capping atau batas atas pada data asli California Housing, yaitu observasi dengan nilai rumah yang melebihi batas tertentu dicatat sebagai US$500.001. Berdasarkan data, terdapat sekitar 965 observasi (4,7% dari total data) yang berada pada nilai tersebut. Akibatnya, distribusi harga rumah pada kelompok bernilai tinggi tidak sepenuhnya tercerminkan dalam data dan berpotensi memengaruhi hasil analisis maupun pemodelan prediktif, khususnya pada segmen rumah dengan harga tinggi.

2.3.3 Distribusi Variabel Kategorik

p1 <- ggplot(df, aes(x = reorder(ocean_proximity, -table(ocean_proximity)[ocean_proximity]),
                      fill = ocean_proximity)) +
  geom_bar(color = "white", alpha = 0.9, show.legend = FALSE) +
  geom_text(stat = "count", aes(label = after_stat(count)), vjust = -0.4,
            fontface = "bold", size = 3.8) +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = comma, expand = expansion(mult = c(0, 0.1))) +
  labs(title = "Jumlah Blok per Kategori", x = NULL, y = "Frekuensi") +
  theme_bw(base_size = 11) +
  theme(plot.title = element_text(face = "bold"),
        axis.text.x = element_text(angle = 20, hjust = 1),
        panel.grid.major.x = element_blank())

p2 <- ggplot(df, aes(x = reorder(ocean_proximity, median_house_value,
                                   FUN = median),
                      y = median_house_value, fill = ocean_proximity)) +
  geom_boxplot(alpha = 0.8, color = "gray30", show.legend = FALSE,
               outlier.alpha = 0.2, outlier.size = 0.7) +
  coord_flip() +
  scale_fill_brewer(palette = "Set2") +
  scale_y_continuous(labels = dollar_format(prefix = "$", big.mark = ",")) +
  labs(title    = "Distribusi Harga Rumah per Kategori",
       subtitle = "Diurutkan berdasarkan median harga",
       x = NULL, y = "Median House Value") +
  theme_bw(base_size = 11) +
  theme(plot.title = element_text(face = "bold"),
        plot.subtitle = element_text(size = 9))

grid.arrange(p1, p2, ncol = 2)

Berdasarkan visualisasi tersebut, dapat diinterpertasikan sebagai berikut.

Panel kiri (Frekuensi): Distribusi kategori ocean_proximity tidak seimbang. Kategori <1H OCEAN mendominasi dataset dengan 9.136 blok sensus (44,3%), diikuti oleh INLAND sebanyak 6.551 blok (31,7%). Sementara itu, kategori NEAR OCEAN dan NEAR BAY masing-masing mencakup 12,9% dan 11,1% dari total observasi. Kategori ISLAND hanya terdiri atas 5 observasi sehingga statistik yang dihasilkan pada kategori ini perlu diinterpretasikan dengan sangat hati-hati karena kurang representatif.

Panel kanan (Distribusi harga): Boxplot menunjukkan adanya perbedaan harga rumah yang cukup jelas antar kategori lokasi. Kategori INLAND memiliki median harga rumah terendah (sekitar US$108.500), sedangkan kategori NEAR BAY dan NEAR OCEAN memiliki median harga yang jauh lebih tinggi, masing-masing sekitar US$233.800 dan US$229.450. Kategori <1H OCEAN berada di antara keduanya dengan median sekitar US$214.850. Temuan ini mengindikasikan bahwa kedekatan terhadap wilayah perairan berkaitan dengan peningkatan nilai rumah.

Kategori ISLAND memiliki median harga tertinggi (sekitar US$414.700), namun karena hanya terdiri atas 5 observasi, nilai tersebut tidak dapat digeneralisasikan untuk keseluruhan populasi. Selain itu, sebagian besar kategori menunjukkan keberadaan pencilan (outlier) pada sisi atas distribusi, yang mengindikasikan adanya sejumlah rumah dengan nilai jauh lebih tinggi dibandingkan mayoritas rumah dalam kategori yang sama. Hal ini konsisten dengan hasil analisis sebelumnya yang menunjukkan bahwa variabel median_house_value memiliki distribusi yang menceng ke kanan (right-skewed).

2.4.1 Tabel Missing Values

missing_df <- df %>%
  summarise(across(everything(), ~sum(is.na(.)))) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "N_Missing") %>%
  mutate(
    Pct_Missing = round(N_Missing / nrow(df) * 100, 2),
    Status      = ifelse(N_Missing == 0, "Lengkap", "Ada Missing")
  ) %>%
  arrange(desc(N_Missing))

kable(missing_df,
      caption = "Rekap Missing Values Seluruh Variabel",
      col.names = c("Variabel","Jumlah Missing","% Missing","Status")) %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(missing_df$N_Missing > 0), background = "#fcf8e3", bold = TRUE) %>%
  row_spec(0, bold = TRUE)

2.4.2 Visualisasi Missing Values

ggplot(missing_df, aes(x = reorder(Variabel, N_Missing), y = Pct_Missing,
                        fill = Status)) +
  geom_col(width = 0.6, alpha = 0.85) +
  geom_text(aes(label = ifelse(N_Missing > 0,
                               paste0(N_Missing, " data (", Pct_Missing, "%)"),
                               "Tidak ada")),
            hjust = -0.05, size = 3.5, fontface = "bold") +
  coord_flip(ylim = c(0, 3)) +
  scale_fill_manual(values = c("Lengkap" = "#2ca25f", "Ada Missing" = "#e67e22")) +
  scale_y_continuous(labels = function(x) paste0(x, "%")) +
  labs(
    title    = "Persentase Missing Values per Variabel",
    subtitle = "Hanya variabel total_bedrooms yang memiliki nilai hilang",
    x = NULL, y = "Persentase Missing (%)", fill = "Status"
  ) +
  theme_bw(base_size = 12) +
  theme(plot.title    = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 10),
        panel.grid.major.y = element_blank())

Berdasarkan analisis missing values, secara keseluruhan, kualitas data tergolong sangat baik. Dari 10 variabel, hanya satu variabel yang memiliki missing values yaitu total_bedrooms dengan 207 nilai hilang (1,0% dari total 20.640 observasi). Jumlah ini tergolong kecil dan tidak akan berdampak signifikan pada analisis. Penyebab paling plausibel adalah ketidaklengkapan pencatatan pada beberapa blok di dataset. Untuk penanganan, imputasi menggunakan median atau menggunakan rasio total_bedrooms/total_rooms dari blok-blok yang datanya lengkap adalah pendekatan yang paling tepat untuk konteks dataset ini, mengingat distribusinya yang right-skewed (imputasi mean akan menghasilkan nilai yang bias).

2.5.1 Boxplot Outlier

box_vars <- c("housing_median_age","total_rooms","total_bedrooms",
              "population","households","median_income","median_house_value")

box_long <- df %>%
  select(all_of(box_vars)) %>%
  pivot_longer(everything(), names_to = "Variabel", values_to = "Nilai")

ggplot(box_long, aes(x = Variabel, y = Nilai, fill = Variabel)) +
  geom_boxplot(alpha = 0.75, color = "gray30",
               outlier.color = "red", outlier.alpha = 0.2,
               outlier.size = 0.7, width = 0.55) +
  facet_wrap(~Variabel, scales = "free", ncol = 4) +
  scale_fill_viridis_d(option = "C") +
  scale_y_continuous(labels = comma) +
  labs(
    title    = "Boxplot Deteksi Outlier",
    subtitle = "Titik merah di luar whisker = outlier (lebih dari 1.5 ?? IQR dari Q1/Q3)",
    x = NULL, y = "Nilai"
  ) +
  theme_bw(base_size = 11) +
  theme(
    plot.title       = element_text(face = "bold", size = 14),
    plot.subtitle    = element_text(size = 9),
    strip.text       = element_text(face = "bold"),
    legend.position  = "none",
    axis.text.x      = element_blank(),
    axis.ticks.x     = element_blank()
  )

Berdasarkan boxplot di atas, dapat diinterpretasikan sebagai berikut.

• housing_median_age: Tidak terdapat outlier berdasarkan metode IQR (0% data berada di luar batas whisker). Dengan Q1 = 18 tahun dan Q3 = 37 tahun, seluruh observasi masih berada dalam rentang yang dianggap wajar menurut kriteria IQR. Variabel ini merupakan satu-satunya variabel numerik yang tidak memiliki outlier.

• total_rooms: Terdapat 1.287 outlier (6,2%) di atas batas atas IQR sebesar 5.698,38 kamar. Sebagian besar blok sensus memiliki jumlah kamar yang relatif rendah hingga sedang, namun terdapat sejumlah blok dengan jumlah kamar yang sangat besar hingga mencapai 39.320 kamar. Keberadaan observasi ekstrem ini menyebabkan distribusi menjadi sangat menceng ke kanan.

• total_bedrooms: Terdapat 1.271 outlier (6,2%) dengan batas atas IQR sebesar 1.173,50 kamar tidur. Nilai maksimum mencapai 6.445 kamar tidur, menunjukkan adanya sejumlah blok sensus dengan ukuran yang jauh lebih besar dibandingkan mayoritas observasi. Pola ini sejalan dengan distribusi total_rooms karena kedua variabel menggambarkan karakteristik ukuran hunian.

• population: Memiliki 1.196 outlier (5,8%) dengan batas atas IQR sebesar 3.132 jiwa. Nilai maksimum yang mencapai 35.682 jiwa menunjukkan adanya beberapa blok sensus dengan kepadatan penduduk yang sangat tinggi dibandingkan mayoritas blok lainnya. Variabel ini juga merupakan salah satu variabel dengan tingkat skewness tertinggi dalam dataset.

• households: Terdapat 1.220 outlier (5,9%) di atas batas atas IQR sebesar 1.092,50 rumah tangga. Nilai maksimum mencapai 6.082 rumah tangga, yang menunjukkan adanya sejumlah blok sensus dengan jumlah rumah tangga yang jauh lebih besar dibandingkan sebagian besar observasi lainnya.

• median_income: Memiliki 681 outlier (3,3%), jumlah yang relatif lebih sedikit dibandingkan variabel agregat lainnya. Batas atas IQR sebesar 8,01 menunjukkan bahwa sebagian besar wilayah memiliki tingkat pendapatan yang tidak terlalu ekstrem, meskipun masih terdapat sejumlah wilayah dengan pendapatan yang jauh lebih tinggi dibandingkan mayoritas observasi.

• median_house_value: Terdapat 1.071 outlier (5,2%) berdasarkan metode IQR dengan batas atas sebesar US$482.412,50. Sebagian besar outlier berada di sekitar nilai maksimum US$500.001 yang merupakan batas atas (capping) pada data asli California Housing. Oleh karena itu, sebagian outlier pada variabel ini kemungkinan dipengaruhi oleh proses pembatasan nilai dalam dataset, sehingga interpretasinya perlu dilakukan dengan hati-hati.

Variabel longitude dan latitude tidak disertakan dalam analisis outlier karena keduanya merupakan variabel koordinat geografis. Nilai ekstrem pada kedua variabel tersebut mencerminkan lokasi spasial yang valid, bukan pengamatan yang menyimpang dari populasi, sehingga analisis outlier lebih difokuskan pada variabel karakteristik perumahan, kependudukan, dan ekonomi.

2.5.2 Tabel Rekap Outlier

outlier_stats <- df %>%
  select(all_of(box_vars)) %>%
  summarise(across(everything(), list(
    Q1         = ~round(quantile(.,0.25,na.rm=T), 2),
    Q3         = ~round(quantile(.,0.75,na.rm=T), 2),
    IQRval     = ~round(IQR(.,na.rm=T), 2),
    BatasAtas  = ~round(quantile(.,0.75,na.rm=T) + 1.5*IQR(.,na.rm=T), 2),
    NOutlier   = ~sum(. > quantile(.,0.75,na.rm=T)+1.5*IQR(.,na.rm=T) |
                      . < quantile(.,0.25,na.rm=T)-1.5*IQR(.,na.rm=T), na.rm=T),
    PctOutlier = ~round(sum(. > quantile(.,0.75,na.rm=T)+1.5*IQR(.,na.rm=T) |
                            . < quantile(.,0.25,na.rm=T)-1.5*IQR(.,na.rm=T),
                            na.rm=T)/n()*100, 1)
  ))) %>%
  pivot_longer(everything(),
               names_to  = c("Variabel",".value"),
               names_sep = "_(?=[^_]+$)") %>%
  arrange(desc(PctOutlier))

kable(outlier_stats,
      caption = "Rekap Outlier Berdasarkan Metode IQR",
      col.names = c("Variabel","Q1","Q3","IQR","Batas Atas",
                    "N Outlier","% Outlier"),
      format.args = list(big.mark = ",")) %>%
  kable_styling(bootstrap_options = c("striped","hover","condensed"),
                full_width = TRUE) %>%
  row_spec(which(outlier_stats$PctOutlier > 5), background = "#fdf2f2") %>%
  row_spec(which(outlier_stats$PctOutlier == 0), background = "#eafaf1") %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE)

2.6.1 Heatmap Korelasi

corr_vars <- df %>% select(-ocean_proximity, -longitude, -latitude)
cor_matrix <- cor(corr_vars, use = "pairwise.complete.obs")

corrplot(cor_matrix,
         method      = "color",
         type        = "upper",
         order       = "hclust",
         col         = colorRampPalette(c("#d73027","white","#1a9850"))(200),
         addCoef.col = "black",
         number.cex  = 0.82,
         tl.col      = "black",
         tl.cex      = 0.92,
         tl.srt      = 45,
         diag        = FALSE,
         cl.ratio    = 0.15,
         title       = "Heatmap Korelasi Pearson Antar Variabel Numerik",
         mar         = c(0, 0, 2, 0))

Berdasarkan heatmap korelasi tersebut, dapat diinterpretasikan sebagai berikut.

Korelasi antar variabel prediktor (potensi multikolinearitas):

• total_bedrooms ??? households memiliki korelasi tertinggi (r = 0,980), menunjukkan hubungan linear yang hampir sempurna. Hal ini mengindikasikan bahwa kedua variabel membawa informasi yang sangat mirip mengenai karakteristik ukuran blok sensus.

• total_rooms ??? total_bedrooms (r = 0,930) dan total_rooms ??? households (r = 0,919) juga menunjukkan korelasi yang sangat kuat. Bersama dengan population, variabel-variabel tersebut membentuk kelompok dengan tingkat keterkaitan yang tinggi sehingga berpotensi menimbulkan multikolinearitas.

• population ??? households (r = 0,907), population ??? total_bedrooms (r = 0,878), dan population ??? total_rooms (r = 0,857) menunjukkan bahwa blok sensus dengan jumlah rumah tangga yang besar umumnya juga memiliki populasi dan jumlah ruangan yang lebih besar. Korelasi yang tinggi ini mengindikasikan adanya tumpang tindih informasi antarvariabel.

• housing_median_age memiliki korelasi negatif sedang dengan variabel-variabel agregat, yaitu population (-0,30), households (-0,30), total_bedrooms (-0,32), dan total_rooms (-0,36). Hal ini menunjukkan bahwa blok dengan usia median rumah yang lebih tinggi cenderung memiliki ukuran dan jumlah penduduk yang lebih kecil.

Korelasi variabel prediktor dengan target (median_house_value):

• median_income memiliki korelasi positif paling kuat terhadap harga rumah (r = 0,688). Hasil ini menunjukkan bahwa wilayah dengan pendapatan median yang lebih tinggi cenderung memiliki nilai rumah yang lebih tinggi, sehingga variabel ini merupakan prediktor yang paling informatif terhadap median_house_value.

• housing_median_age memiliki korelasi positif yang lemah terhadap harga rumah (r = 0,106), menunjukkan bahwa hubungan linear antara usia rumah dan harga rumah relatif kecil.

• total_rooms (r = 0,134), households (r = 0,066), dan total_bedrooms (r = 0,050) hanya memiliki korelasi positif yang sangat lemah terhadap harga rumah.

• population memiliki korelasi negatif yang sangat kecil terhadap harga rumah (r = ???0,025). Nilai ini menunjukkan bahwa jumlah penduduk pada suatu blok sensus hampir tidak memiliki hubungan linear dengan harga rumah.

Secara keseluruhan, heatmap menunjukkan adanya kelompok variabel dengan korelasi sangat tinggi (population, households, total_rooms, dan total_bedrooms) yang mengindikasikan potensi multikolinearitas. Temuan ini mendukung perlunya penerapan teknik reduksi dimensi seperti PCA untuk merangkum informasi dari variabel-variabel yang saling berkorelasi kuat tersebut ke dalam sejumlah komponen yang lebih ringkas.

2.6.2 Scatter Plot vs Target

predictors <- c("median_income","housing_median_age",
                "total_rooms","population","households","total_bedrooms")

plots <- lapply(predictors, function(var) {
  r_val <- round(cor(df[[var]], df$median_house_value, use = "complete.obs"), 3)
  ggplot(df, aes_string(x = var, y = "median_house_value")) +
    geom_point(alpha = 0.10, color = "#2c7fb8", size = 0.6) +
    geom_smooth(method = "lm", color = "red", se = FALSE, size = 1) +
    scale_y_continuous(labels = dollar_format(prefix = "$", big.mark = ",")) +
    scale_x_continuous(labels = comma) +
    labs(title = paste0(var, "\n(r = ", r_val, ")"),
         x = var, y = "Harga Rumah") +
    theme_bw(base_size = 9) +
    theme(plot.title = element_text(face = "bold", size = 9, hjust = 0.5))
})

grid.arrange(grobs = plots, ncol = 3)

Berdasarkan scatter plot tersebut, didapat interpretasi hubungan sebagai berikut.

• median_income vs harga rumah: Menunjukkan hubungan linear positif paling kuat di antara seluruh variabel prediktor (r = 0,688). Semakin tinggi pendapatan median suatu blok sensus, semakin tinggi pula nilai rumah pada wilayah tersebut. Pola titik-titik data membentuk tren naik yang cukup jelas, meskipun masih terdapat variasi yang cukup besar pada setiap tingkat pendapatan. Selain itu, terlihat bahwa sebaran harga rumah cenderung semakin lebar pada tingkat pendapatan yang tinggi, yang mengindikasikan peningkatan variasi harga rumah pada kelompok berpendapatan tinggi.

• housing_median_age vs harga rumah: Menunjukkan hubungan positif yang sangat lemah (r = 0,106). Sebaran titik yang sangat menyebar dan garis regresi yang relatif datar menunjukkan bahwa usia median rumah hanya memiliki hubungan linear yang kecil terhadap harga rumah. Meskipun demikian, terdapat kecenderungan bahwa wilayah dengan usia rumah yang lebih tua memiliki nilai rumah yang sedikit lebih tinggi.

• total_rooms vs harga rumah: Memiliki korelasi positif yang lemah (r = 0,134). Meskipun garis regresi menunjukkan kecenderungan meningkat, sebaran titik yang sangat luas menunjukkan bahwa jumlah kamar bukan merupakan prediktor utama harga rumah. Banyaknya kamar pada suatu blok sensus tidak selalu diikuti oleh peningkatan nilai rumah yang signifikan.

• total_bedrooms vs harga rumah: Menunjukkan korelasi positif yang sangat lemah (r = 0,050). Pola sebaran data tidak memperlihatkan hubungan yang kuat, sehingga jumlah kamar tidur secara langsung kurang mampu menjelaskan variasi harga rumah.

• households vs harga rumah: Memiliki korelasi positif yang sangat lemah (r = 0,066). Walaupun terdapat sedikit kecenderungan peningkatan harga seiring bertambahnya jumlah rumah tangga, hubungan tersebut relatif kecil dan tidak menunjukkan pola yang jelas.

• population vs harga rumah: Menunjukkan korelasi negatif yang sangat lemah (r = -0,025). Nilai korelasi yang mendekati nol menunjukkan bahwa jumlah penduduk pada suatu blok sensus hampir tidak memiliki hubungan linear dengan harga rumah. Hal ini juga terlihat dari pola titik yang menyebar tanpa arah yang jelas.

Secara keseluruhan, scatter plot mengonfirmasi hasil analisis korelasi sebelumnya bahwa median_income merupakan variabel yang memiliki hubungan paling kuat dengan median_house_value, sedangkan variabel ukuran blok sensus seperti population, households, total_rooms, dan total_bedrooms memiliki hubungan linear yang relatif lemah terhadap harga rumah. Temuan ini menunjukkan bahwa kondisi ekonomi wilayah lebih berpengaruh terhadap harga rumah dibandingkan ukuran atau kepadatan blok sensus itu sendiri.

2.6.3 Tabel Korelasi dengan Target

cor_target <- cor_matrix["median_house_value",] %>%
  as.data.frame() %>%
  rownames_to_column("Variabel") %>%
  rename(Korelasi = ".") %>%
  filter(Variabel != "median_house_value") %>%
  mutate(
    Korelasi      = round(Korelasi, 4),
    Kekuatan      = case_when(
      abs(Korelasi) >= 0.6 ~ "Kuat",
      abs(Korelasi) >= 0.3 ~ "Sedang",
      abs(Korelasi) >= 0.1 ~ "Lemah",
      TRUE                 ~ "Sangat Lemah"
    ),
    Arah = ifelse(Korelasi > 0, "Positif", "Negatif")
  ) %>%
  arrange(desc(abs(Korelasi)))

kable(cor_target,
      caption = "Korelasi Setiap Variabel terhadap median_house_value") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(1, bold = TRUE, background = "#fcf8e3") %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE)

2.7.1 Heatmap Korelasi Prediktor

pred_vars <- corr_vars %>% select(-median_house_value)
cor_pred  <- cor(pred_vars, use = "pairwise.complete.obs")

ggcorrplot(cor_pred,
           hc.order    = TRUE,
           type        = "lower",
           lab         = TRUE,
           lab_size    = 3.8,
           colors      = c("#d73027","white","#1a9850"),
           outline.col = "white",
           ggtheme     = theme_bw()) +
  labs(title    = "Korelasi Antar Variabel Prediktor",
       subtitle = "Fokus pada korelasi yang mendekati ??±1 (indikasi multikolinearitas)") +
  theme(
    plot.title    = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(size = 10),
    axis.text.x   = element_text(angle = 45, hjust = 1)
  )

Berdasarkan heatmap tersebut, dapat diinterpretasikan bahwa heatmap menunjukkan adanya kelompok variabel yang memiliki korelasi sangat tinggi satu sama lain, yaitu population, total_rooms, total_bedrooms, dan households, dengan koefisien korelasi berkisar antara 0,86 hingga 0,98. Nilai korelasi yang sangat tinggi tersebut mengindikasikan adanya potensi multikolinearitas yang kuat serta menunjukkan bahwa variabel-variabel tersebut kemungkinan mengandung informasi yang saling tumpang tindih mengenai ukuran atau kapasitas suatu blok sensus.

Sebaliknya, variabel median_income dan housing_median_age memiliki korelasi yang relatif rendah terhadap kelompok variabel tersebut, dengan nilai korelasi berkisar antara ???0,36 hingga 0,20. Hal ini menunjukkan bahwa kedua variabel tersebut menangkap aspek yang berbeda dari karakteristik wilayah dibandingkan variabel-variabel yang berkaitan dengan ukuran blok sensus. Temuan ini mengindikasikan bahwa proses reduksi dimensi, seperti PCA, berpotensi efektif dalam merangkum informasi dari kelompok variabel yang saling berkorelasi tinggi tanpa kehilangan terlalu banyak informasi.

2.7.2 VIF (Variance Inflation Factor)

df_vif   <- df %>% select(-ocean_proximity, -longitude, -latitude) %>% drop_na()
model_vif <- lm(median_house_value ~ ., data = df_vif)

vif_df <- data.frame(
  Variabel = names(vif(model_vif)),
  VIF      = round(vif(model_vif), 2)
) %>%
  mutate(
    Kategori = case_when(
      VIF > 10 ~ "PARAH (>10)",
      VIF > 5  ~ "MODERAT (5-10)",
      VIF > 2  ~ "RINGAN (2-5)",
      TRUE     ~ "AMAN (<2)"
    )
  ) %>%
  arrange(desc(VIF))

kable(vif_df, caption = "Variance Inflation Factor (VIF) per Variabel Prediktor") %>%
  kable_styling(bootstrap_options = c("striped","hover"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(vif_df$VIF > 10), background = "#fdf2f2", bold = TRUE) %>%
  row_spec(which(vif_df$VIF > 5 & vif_df$VIF <= 10), background = "#fef9e7") %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE) %>%
  column_spec(3, bold = TRUE,
              color = ifelse(vif_df$VIF > 10, "red",
                      ifelse(vif_df$VIF > 5, "darkorange", "darkgreen")))

2.7.3 Visualisasi VIF

ggplot(vif_df, aes(x = reorder(Variabel, VIF), y = VIF,
                    fill = Kategori)) +
  geom_col(width = 0.6, alpha = 0.88) +
  geom_hline(yintercept = 10, linetype = "dashed", color = "red",     size = 1.0) +
  geom_hline(yintercept = 5,  linetype = "dashed", color = "darkorange", size = 1.0) +
  geom_text(aes(label = round(VIF, 1)), hjust = -0.15, fontface = "bold", size = 4) +
  annotate("text", x = 0.7, y = 11.2, label = "VIF = 10 (batas parah)",
           color = "red", fontface = "italic", size = 3.2, hjust = 0) +
  annotate("text", x = 0.7, y = 6.2,  label = "VIF = 5 (batas moderat)",
           color = "darkorange", fontface = "italic", size = 3.2, hjust = 0) +
  coord_flip(ylim = c(0, max(vif_df$VIF) * 1.15)) +
  scale_fill_manual(values = c("PARAH (>10)"     = "#e74c3c",
                               "MODERAT (5-10)"  = "#e67e22",
                               "RINGAN (2-5)"    = "#f1c40f",
                               "AMAN (<2)"       = "#2ca25f")) +
  labs(
    title    = "Nilai VIF per Variabel Prediktor",
    subtitle = "VIF > 10 = multikolinearitas parah | VIF > 5 = perlu perhatian",
    x = NULL, y = "VIF", fill = "Kategori"
  ) +
  theme_bw(base_size = 12) +
  theme(plot.title    = element_text(face = "bold", size = 14),
        plot.subtitle = element_text(size = 10),
        panel.grid.major.y = element_blank(),
        legend.position = "bottom")

VIF (Variance Inflation Factor) mengukur seberapa besar varians estimasi koefisien suatu variabel meningkat akibat adanya multikolinearitas dengan variabel prediktor lainnya. Secara umum, nilai VIF > 10 menunjukkan adanya multikolinearitas yang tinggi dan perlu mendapat perhatian khusus. Berdasarkan nilai VIF dan visualisasi heatmap di atas, didapat interpretasi sebagai berikut.

• total_bedrooms (VIF = 35,50) memiliki nilai VIF tertinggi dan menunjukkan adanya multikolinearitas yang sangat kuat. Nilai ini mengindikasikan bahwa sebagian besar informasi yang terkandung dalam variabel tersebut dapat dijelaskan oleh kombinasi variabel prediktor lainnya.

• households (VIF = 33,81) juga menunjukkan tingkat multikolinearitas yang sangat tinggi. Hasil ini konsisten dengan analisis korelasi sebelumnya yang menunjukkan hubungan yang sangat kuat antara households, total_bedrooms, total_rooms, dan population.

• total_rooms (VIF = 11,90) masih berada pada kategori multikolinearitas tinggi. Variabel ini memiliki korelasi yang kuat dengan beberapa variabel ukuran blok sensus lainnya sehingga mengandung informasi yang relatif tumpang tindih.

• population (VIF = 6,22) termasuk dalam kategori moderat. Meskipun masih menunjukkan adanya hubungan yang cukup kuat dengan variabel lain, variabel ini tampaknya masih menyimpan informasi unik yang tidak sepenuhnya dijelaskan oleh prediktor lainnya.

• median_income (VIF = 1,51) berada pada kategori aman dan menunjukkan bahwa variabel ini relatif independen dari prediktor lainnya.

• housing_median_age (VIF = 1,16) memiliki nilai VIF terendah, sehingga dapat dianggap hampir tidak mengalami masalah multikolinearitas.

Hasil VIF mengonfirmasi temuan pada analisis korelasi sebelumnya bahwa kelompok variabel population, households, total_rooms, dan total_bedrooms memiliki tingkat keterkaitan yang sangat tinggi. Kondisi ini menunjukkan adanya redundansi informasi yang cukup besar antarvariabel. Oleh karena itu, penerapan teknik reduksi dimensi seperti Principal Component Analysis (PCA) menjadi relevan untuk merangkum informasi dari variabel-variabel tersebut ke dalam sejumlah komponen utama yang lebih ringkas tanpa kehilangan informasi yang signifikan.

5.4.1 Tabel Variansi Komponen

eigenvalues <- pca_result$sdev^2
prop_var    <- eigenvalues / sum(eigenvalues)
cumul_var   <- cumsum(prop_var)

var_table <- data.frame(
  Komponen           = paste0("PC", seq_along(eigenvalues)),
  Eigenvalue         = round(eigenvalues, 4),
  Proporsi_Variansi  = paste0(round(prop_var  * 100, 2), "%"),
  Kumulatif_Variansi = paste0(round(cumul_var * 100, 2), "%"),
  Status_Kaiser      = ifelse(eigenvalues >= 1, "??? Dipertahankan", "??? Diabaikan")
)

kable(var_table,
      caption = "Proporsi Variansi dan Eigenvalue Setiap Komponen Utama",
      col.names = c("Komponen", "Eigenvalue", "Proporsi Variansi",
                    "Kumulatif Variansi", "Status Kaiser (??? 1)")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = FALSE, position = "center") %>%
  row_spec(which(eigenvalues >= 1), background = "#eafaf1", bold = TRUE) %>%
  row_spec(which(eigenvalues <  1), background = "#fdf2f2") %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE)

Berdasarkan tabel di atas, dapat diinterpretasikan sebagai berikut.

• Kaiser Criterion (eigenvalue ??? 1) digunakan sebagai acuan pemilihan komponen. Komponen dengan eigenvalue ??? 1 menjelaskan variasi yang lebih besar dari satu variabel tunggal yang telah distandarisasi, sehingga dianggap layak dipertahankan.
• PC1 menjelaskan porsi terbesar variasi karena merangkum kelompok variabel yang saling berkorelasi tinggi (total_rooms, total_bedrooms, households, population).
• Berdasarkan kriteria kumulatif variansi ??? 80%, sebagian besar informasi data sudah terangkum dalam beberapa komponen pertama saja.

5.4.2 Scree Plot

fviz_eig(pca_result,
         addlabels  = TRUE,
         ylim       = c(0, 60),
         barfill    = "#2c7fb8",
         barcolor   = "white",
         linecolor  = "#e74c3c",
         ggtheme    = theme_bw(base_size = 12)) +
  geom_hline(yintercept = (1 / length(eigenvalues)) * 100,
             linetype = "dashed", color = "darkorange", size = 1) +
  annotate("text",
           x     = length(eigenvalues) - 0.5,
           y     = (1 / length(eigenvalues)) * 100 + 2,
           label = "Batas Kaiser (eigenvalue = 1)",
           color = "darkorange", fontface = "italic", size = 3.5, hjust = 1) +
  labs(
    title    = "Scree Plot ??? Proporsi Variansi Setiap Komponen Utama",
    subtitle = "Batang biru = % variansi | Garis merah = tren | Garis oranye = batas Kaiser",
    x = "Komponen Utama (PC)",
    y = "Persentase Variansi (%)"
  )

Scree plot di atas menunjukkan penurunan yang tajam (elbow) pada komponen-komponen awal, yang menandakan bahwa sebagian besar variasi dalam data terkonsentrasi pada beberapa komponen pertama. Titik elbow yaitu titik di mana kurva mulai mendatar menjadi panduan visual untuk menentukan jumlah komponen yang optimal karena penambahan komponen setelah titik tersebut tidak lagi memberikan peningkatan variansi yang berarti.

5.4.3 Kumulatif Variansi

cumvar_df <- data.frame(
  PC        = factor(paste0("PC", seq_along(cumul_var)),
                     levels = paste0("PC", seq_along(cumul_var))),
  Kumulatif = cumul_var * 100
)

ggplot(cumvar_df, aes(x = PC, y = Kumulatif, group = 1)) +
  geom_line(color = "#2c7fb8", size = 1.2) +
  geom_point(aes(color = Kumulatif >= 80), size = 4) +
  geom_hline(yintercept = 80, linetype = "dashed", color = "red",        size = 1) +
  geom_hline(yintercept = 90, linetype = "dashed", color = "darkorange", size = 1) +
  annotate("text", x = length(cumul_var) - 0.5, y = 81.5,
           label = "80% Kumulatif", color = "red",
           fontface = "italic", size = 3.5, hjust = 1) +
  annotate("text", x = length(cumul_var) - 0.5, y = 91.5,
           label = "90% Kumulatif", color = "darkorange",
           fontface = "italic", size = 3.5, hjust = 1) +
  scale_color_manual(values = c("FALSE" = "#7f8c8d", "TRUE" = "#27ae60"),
                     labels = c("FALSE" = "< 80%", "TRUE" = "??? 80%")) +
  scale_y_continuous(labels = function(x) paste0(x, "%"), limits = c(0, 103)) +
  labs(
    title    = "Kumulatif Proporsi Variansi PCA",
    subtitle = "Titik hijau = komponen yang sudah mencapai kumulatif ??? 80%",
    x = "Komponen Utama", y = "Kumulatif Variansi (%)", color = "Status"
  ) +
  theme_bw(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold", size = 14),
    plot.subtitle   = element_text(size = 10),
    axis.text.x     = element_text(angle = 30, hjust = 1),
    legend.position = "bottom"
  )

5.5.1 Tabel Loading Factor

loading_mat <- pca_result$rotation[, 1:5]
loading_df  <- as.data.frame(round(loading_mat, 4))

kable(loading_df,
      caption = "Loading Factor: Kontribusi Variabel terhadap 5 PC Pertama") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"),
                full_width = TRUE) %>%
  row_spec(0, bold = TRUE) %>%
  column_spec(1, bold = TRUE)

5.5.2 Heatmap Loading Factor

load_long <- as.data.frame(loading_mat) %>%
  rownames_to_column("Variabel") %>%
  pivot_longer(-Variabel, names_to = "PC", values_to = "Loading")

ggplot(load_long, aes(x = PC, y = Variabel, fill = Loading)) +
  geom_tile(color = "white", size = 0.5) +
  geom_text(aes(label = round(Loading, 2),
                color = abs(Loading) > 0.30),
            size = 3.8, fontface = "bold") +
  scale_fill_gradient2(low  = "#d73027", mid = "white", high = "#1a9850",
                       midpoint = 0, limits = c(-1, 1), name = "Loading") +
  scale_color_manual(values = c("FALSE" = "gray65", "TRUE" = "black"),
                     guide = "none") +
  scale_x_discrete(position = "top") +
  labs(
    title    = "Heatmap Loading Factor ??? 5 Komponen Utama Pertama",
    subtitle = "Hijau = kontribusi positif kuat | Merah = kontribusi negatif kuat | Putih = lemah",
    x = NULL, y = NULL
  ) +
  theme_bw(base_size = 12) +
  theme(
    plot.title      = element_text(face = "bold", size = 14),
    plot.subtitle   = element_text(size = 10),
    axis.text.y     = element_text(face = "bold"),
    legend.position = "right"
  )

5.5.3 Interpretasi Loading Factor

Berdasarkan tabel dan heatmap loading factor, setiap komponen utama dapat diinterpretasikan maknanya secara substantif sebagai berikut.

PC1 ??? “Ukuran dan Kapasitas Blok Sensus”

PC1 menunjukkan loading positif yang kuat pada variabel-variabel yang mencerminkan ukuran fisik dan kapasitas blok sensus: total_rooms, total_bedrooms, households, dan population. Blok dengan nilai PC1 tinggi adalah blok besar yang dihuni banyak rumah tangga dan penduduk. Komponen ini pada dasarnya merangkum seberapa besar ukuran suatu blok sensus sebuah dimensi tunggal yang sebelumnya tersebar di empat variabel yang saling berkorelasi sangat tinggi (VIF 6???35).

PC2 ??? “Kemakmuran dan Kualitas Hunian”

PC2 menangkap kontras antara variabel ekonomi (median_income dan rooms_per_household dengan loading positif) versus variabel kepadatan seperti bedrooms_ratio (loading negatif). Komponen ini merepresentasikan tingkat kemakmuran dan kualitas relatif hunian kawasan dengan PC2 tinggi adalah kawasan berpendapatan tinggi dengan hunian yang luas (sedikit proporsi kamar tidur dari total ruangan), mencirikan hunian berkelas atau mewah.

PC3 ??? “Posisi Geografis Utara-Selatan”

PC3 didominasi oleh variabel latitude dengan loading yang besar. Komponen ini merepresentasikan posisi geografis utara-selatan di California kawasan Bay Area dan Silicon Valley di utara versus kawasan Los Angeles dan San Diego di selatan.

PC4 ??? “Posisi Geografis Timur-Barat dan Usia Bangunan”

PC4 menangkap variasi pada longitude dan sedikit dari housing_median_age. Komponen ini merepresentasikan posisi geografis timur-barat (pesisir vs. pedalaman) yang juga berkaitan dengan pola perkembangan historis bangunan di California.

PC5 ??? “Kepadatan Penghuni per Rumah Tangga”

PC5 menunjukkan loading yang menonjol pada population_per_household dan merepresentasikan kepadatan penghuni per unit rumah tangga apakah rumah tangga di suatu blok cenderung kecil atau besar jumlah penghuninya.

5.6.1 Biplot PC1 vs PC2

fviz_pca_biplot(pca_result,
                axes          = c(1, 2),
                geom.ind      = "point",
                col.ind       = housing_clean$median_house_value,
                gradient.cols = c("#d73027", "#fee090", "#1a9850"),
                col.var       = "black",
                repel         = TRUE,
                alpha.ind     = 0.4,
                pointsize     = 1.5,
                arrowsize     = 0.8,
                labelsize     = 3.8,
                ggtheme       = theme_bw(base_size = 12)) +
  scale_color_gradient(low    = "#d73027",
                       high   = "#1a9850",
                       name   = "Harga Rumah (USD)",
                       labels = dollar_format(prefix = "$", big.mark = ",")) +
  labs(
    title    = "Biplot PCA ??? PC1 vs PC2",
    subtitle = "Titik = blok sensus (warna = harga rumah) | Panah = arah & kekuatan loading variabel",
    x = paste0("PC1 (", round(prop_var[1] * 100, 1), "% Variansi)"),
    y = paste0("PC2 (", round(prop_var[2] * 100, 1), "% Variansi)")
  ) +
  theme(
    plot.title      = element_text(face = "bold", size = 14),
    plot.subtitle   = element_text(size = 9),
    legend.position = "right"
  )

Biplot di atas memvisualisasikan dua informasi sekaligus. Titik merepresentasikan setiap blok sensus yang diwarnai berdasarkan median_house_valuetitik hijau menandakan harga rumah tinggi, titik merah menandakan harga rendah. Panah menunjukkan arah loading setiap variabel asli: variabel dengan panah searah saling berkorelasi positif, sedangkan panah berlawanan arah berkorelasi negatif.

Terlihat bahwa median_income dan rooms_per_household memiliki panah yang cenderung menuju area titik-titik hijau (harga tinggi), mengonfirmasi bahwa kedua variabel ini adalah prediktor paling kuat harga rumah. Sebaliknya, kelompok total_rooms, total_bedrooms, households, population memiliki panah yang searah satu sama lain membuktikan multikolinearitas yang telah terdeteksi sebelumnya namun tidak mengarah secara dominan ke area harga tinggi, selaras dengan korelasi lemah keempat variabel tersebut terhadap median_house_value.

5.6.2 Kontribusi Variabel ke PC1 & PC2

p_c1 <- fviz_contrib(pca_result, choice = "var", axes = 1,
                      fill = "#2c7fb8", color = "white",
                      ggtheme = theme_bw(base_size = 11)) +
  labs(title = "Kontribusi Variabel ke PC1",
       x = NULL, y = "Kontribusi (%)") +
  theme(plot.title  = element_text(face = "bold"),
        axis.text.x = element_text(angle = 40, hjust = 1))

p_c2 <- fviz_contrib(pca_result, choice = "var", axes = 2,
                      fill = "#e74c3c", color = "white",
                      ggtheme = theme_bw(base_size = 11)) +
  labs(title = "Kontribusi Variabel ke PC2",
       x = NULL, y = "Kontribusi (%)") +
  theme(plot.title  = element_text(face = "bold"),
        axis.text.x = element_text(angle = 40, hjust = 1))

grid.arrange(p_c1, p_c2, ncol = 2)

Garis putus-putus pada grafik kontribusi menunjukkan nilai rata-rata kontribusi jika setiap variabel berkontribusi secara merata (~9,1% per variabel untuk 11 variabel). Variabel yang berada di atas garis putus-putus merupakan penyusun utama komponen tersebut. Pada PC1, kontributor utama adalah kelompok variabel ukuran blok, sedangkan pada PC2 kontributor utama adalah median_income, rooms_per_household, dan bedrooms_ratio sesuai dengan interpretasi substantif yang telah dijelaskan sebelumnya.

5.6.3 Cos² ??? Kualitas Representasi Variabel

fviz_cos2(pca_result, choice = "var", axes = 1:2,
          fill    = "#27ae60",
          color   = "white",
          ggtheme = theme_bw(base_size = 11)) +
  labs(title    = "Cos² Variabel pada PC1 + PC2",
       subtitle = "Seberapa baik setiap variabel terwakili oleh dua komponen pertama",
       x = NULL, y = "Cos² (Kualitas Representasi)") +
  theme(
    plot.title    = element_text(face = "bold", size = 13),
    plot.subtitle = element_text(size = 9),
    axis.text.x   = element_text(angle = 40, hjust = 1)
  )

Nilai cos² mengukur kualitas representasi suatu variabel pada ruang yang dibentuk oleh komponen-komponen terpilih. Nilai mendekati 1 berarti variabel sangat baik terwakili oleh dua komponen pertama, sedangkan nilai mendekati 0 berarti variabel lebih banyak terwakili oleh komponen-komponen selanjutnya.