# Library Bawaan & Manipulasi Data
library(readxl)
library(readr)
library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

library(stringr)

# Library Analisis Missing Value
library(naniar)
library(VIM)

## Loading required package: colorspace

## Loading required package: grid

## VIM is ready to use.

## Suggestions and bug-reports can be submitted at: https://github.com/statistikat/VIM/issues

## 
## Attaching package: 'VIM'

## The following object is masked from 'package:datasets':
## 
##     sleep

# Library Visualisasi & Tabel Interaktif
library(ggplot2)
library(ggcorrplot)
library(DT)

Meteorologi

Input data

# Input Data  
data_met= read_xlsx("D:\\deska_baru2\\Data Meteotologi.xlsx")
str(data_met)

## tibble [12,053 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Tanggal                 : POSIXct[1:12053], format: "1990-01-01" "1990-01-02" ...
##  $ Tn (suhu minimum)       : num [1:12053] 19 19 20 19 19 20 19 19 20 20 ...
##  $ Tx ( suhu maksimum)     : num [1:12053] 28.4 27 26.8 28.4 28.6 22 27.8 28.2 28.2 24.8 ...
##  $ Tavg  (rata rata suhu)  : num [1:12053] 22.2 22.5 23.1 22.5 21.5 20.3 22.3 23.3 23.2 21.8 ...
##  $ RH_avg(kelembaban)      : num [1:12053] 86 83 85 85 86 94 84 72 79 87 ...
##  $ RR (curah hujan )       : num [1:12053] 6 47 0 12 16.5 1 13 1.5 2 0 ...
##  $ ss ( penyinaran)        : num [1:12053] 4.1 0.4 6.2 4.3 NA 0.1 5.2 2.8 5.1 0.4 ...
##  $ ff_avg (angin rata rata): num [1:12053] 0 0 1 2 2 1 1 2 2 2 ...

head(data_met)

## # A tibble: 6 × 8
##   Tanggal             `Tn (suhu minimum)` `Tx ( suhu maksimum)`
##   <dttm>                            <dbl>                 <dbl>
## 1 1990-01-01 00:00:00                  19                  28.4
## 2 1990-01-02 00:00:00                  19                  27  
## 3 1990-01-03 00:00:00                  20                  26.8
## 4 1990-01-04 00:00:00                  19                  28.4
## 5 1990-01-05 00:00:00                  19                  28.6
## 6 1990-01-06 00:00:00                  20                  22  
## # ℹ 5 more variables: `Tavg  (rata rata suhu)` <dbl>,
## #   `RH_avg(kelembaban)` <dbl>, `RR (curah hujan )` <dbl>,
## #   `ss ( penyinaran)` <dbl>, `ff_avg (angin rata rata)` <dbl>

pre-processing Data

# Ganti nama kolom 
data_met = data_met %>%
  rename_with(~ str_remove(., "\\s*\\(.*\\)"))
str(data_met)

## tibble [12,053 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Tanggal: POSIXct[1:12053], format: "1990-01-01" "1990-01-02" ...
##  $ Tn     : num [1:12053] 19 19 20 19 19 20 19 19 20 20 ...
##  $ Tx     : num [1:12053] 28.4 27 26.8 28.4 28.6 22 27.8 28.2 28.2 24.8 ...
##  $ Tavg   : num [1:12053] 22.2 22.5 23.1 22.5 21.5 20.3 22.3 23.3 23.2 21.8 ...
##  $ RH_avg : num [1:12053] 86 83 85 85 86 94 84 72 79 87 ...
##  $ RR     : num [1:12053] 6 47 0 12 16.5 1 13 1.5 2 0 ...
##  $ ss     : num [1:12053] 4.1 0.4 6.2 4.3 NA 0.1 5.2 2.8 5.1 0.4 ...
##  $ ff_avg : num [1:12053] 0 0 1 2 2 1 1 2 2 2 ...

# Mengubah kode 99.9 dalam data menjadi NA
data_met[data_met==8888] = NA

EDA

Nilai Statistik

summary(data_met)

##     Tanggal                 Tn              Tx             Tavg      
##  Min.   :1990-01-01   Min.   :13.00   Min.   :18.80   Min.   :18.40  
##  1st Qu.:1998-04-02   1st Qu.:19.00   1st Qu.:28.10   1st Qu.:22.80  
##  Median :2006-07-02   Median :20.00   Median :29.00   Median :23.40  
##  Mean   :2006-07-02   Mean   :19.38   Mean   :29.01   Mean   :23.39  
##  3rd Qu.:2014-10-01   3rd Qu.:20.00   3rd Qu.:30.00   3rd Qu.:24.00  
##  Max.   :2022-12-31   Max.   :29.00   Max.   :36.00   Max.   :28.90  
##                       NAs    :578     NAs    :326     NAs    :156    
##      RH_avg            RR                ss             ff_avg      
##  Min.   :42.00   Min.   :  0.000   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:73.00   1st Qu.:  0.000   1st Qu.: 3.300   1st Qu.: 1.000  
##  Median :79.00   Median :  0.800   Median : 5.100   Median : 2.000  
##  Mean   :77.99   Mean   :  6.738   Mean   : 4.924   Mean   : 1.736  
##  3rd Qu.:84.00   3rd Qu.:  7.800   3rd Qu.: 6.800   3rd Qu.: 2.000  
##  Max.   :96.00   Max.   :160.000   Max.   :10.000   Max.   :15.000  
##  NAs    :330     NAs    :1342      NAs    :492      NAs    :114

Missing Value

# Missing Value 
data_met = data_met %>%
  mutate(tahun = year(Tanggal))
missing_tahunan = data_met %>%
  group_by(tahun) %>%
  summarise(
    RR_na = sum(is.na(RR)),
    Tn_na = sum(is.na(Tn)),
    Tx_na = sum(is.na(Tx))
  ) 

library(DT)

# Membuat tabel interaktif
datatable(missing_tahunan, 
          options = list(pageLength = 33), 
          caption = 'Tabel Missing Value Tahunan')

# Visualisasi
# A. Menggunakan naniar: Melihat peta sebaran data kosong
# Hitam = Data ada, Abu-abu = Data kosong (NA)
vis_miss(data_met) + 
  labs(title = "Peta Sebaran Missing Value Data Meteorologi")

# B. Menggunakan VIM: Melihat kombinasi variabel yang sering kosong bersamaan
aggr(data_met, 
     col = c('navyblue','red'), 
     numbers = TRUE, 
     sortVars = TRUE, 
     labels = names(data_met), 
     cex.axis = .7, 
     gap = 3, 
     ylab = c("Proporsi Missing Data", "Pola Kombinasi Kosong"))

## Warning in plot.aggr(res, ...): not enough vertical space to display
## frequencies (too many combinations)

## 
##  Variables sorted by number of missings: 
##  Variable       Count
##        RR 0.111341575
##        Tn 0.047954866
##        ss 0.040819713
##    RH_avg 0.027379076
##        Tx 0.027047208
##      Tavg 0.012942836
##    ff_avg 0.009458226
##   Tanggal 0.000000000
##     tahun 0.000000000

Korelasi

# Korelasi 
data_numeric = na.omit(data_met[, 2:8]) # Buang kolom Tanggal dan hapus NA sementara khusus untuk hitung korelasi
data_numeric = na.omit(data_met[, 2:8]) 
matriks_korelasi = cor(data_numeric)

# Visualisasi 
ggcorrplot(matriks_korelasi, 
           method = "square", 
           type = "lower",      # Hanya menampilkan separuh bawah matriks
           lab = TRUE,          # Menampilkan angka korelasi
           colors = c("red", "white", "blue"),
           title = "Matriks Korelasi Variabel Meteorologi Bandung")

## 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 per session.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

matriks_korelasi = cor(data_numeric)

Distribusi Data

# Histogram 
data_long = data_met %>%
  select(Tn, Tx, RR) %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")

# Visualisasi
ggplot(data_long, aes(x = Variabel, y = Nilai, fill = Variabel)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 16) +
  coord_flip() +  
  theme_minimal() +
  labs(title = "Distribusi Suhu (Tn, Tx) & Curah Hujan (RR)",
       x = "Variabel Iklim", 
       y = "Nilai Pengukuran")

## Warning: Removed 2246 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

Tren Curah hujan

# Line chart untuk Curah Hujan (RR)
ggplot(data_met, aes(x = Tanggal, y = RR)) +
  geom_line(color = "dodgerblue", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Tren Runtun Waktu Curah Hujan (RR) Kota Bandung 1990-2022",
       x = "Tahun", 
       y = "Curah Hujan (mm)")

# Line chart untuk Suhu Rata-rata (Tavg)
ggplot(data_met, aes(x = Tanggal, y = Tavg)) +
  geom_line(color = "firebrick", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Tren Runtun Waktu Suhu Rata-rata (Tavg) Kota Bandung 1990-2022",
       x = "Tahun", 
       y = "Suhu Rata-rata (°C)")

Uji Algoritma MF: Skema Penghapusan Data

# Data Ground Truth
gt_met = na.omit(data_met[,2:8])
anyNA(gt_met) # Cek data sudah bersih NA

## [1] FALSE

Skema Rendah (10%)

Duplikasi data dan hilangkan 10% nilai dalam data

set.seed(123)
uji_1 = gt_met
for (col in names(uji_1)) {
  n_hilang = round(0.1 * nrow(uji_1))
  idx = sample(1:nrow(uji_1), n_hilang)
  uji_1[idx, col] = NA
}

# Simpan nilai asli yang nanti akan kita tebak
mask_1 = is.na(uji_1)
nilai_asli_1 = gt_met[mask_1]

Skema Sedang (40%)

Duplikasi data dan hilangkan 40% nilai dalam data

set.seed(123)
uji_2 = gt_met
for (col in names(uji_2)) {
  n_hilang = round(0.4 * nrow(uji_2))
  idx = sample(1:nrow(uji_2), n_hilang)
  uji_2[idx, col] = NA
}

# Simpan nilai asli yang nanti akan kita tebak
mask_2 = is.na(uji_2)
nilai_asli_2 = gt_met[mask_2]

Skema Tinggi (70%)

Duplikasi data dan hilangkan 40% nilai dalam data

set.seed(123)
uji_3 = gt_met
for (col in names(uji_3)) {
  n_hilang = round(0.4 * nrow(uji_3))
  idx = sample(1:nrow(uji_3), n_hilang)
  uji_3[idx, col] = NA
}

# Simpan nilai asli yang nanti akan kita tebak
mask_3 = is.na(uji_3)
nilai_asli_3 = gt_met[mask_3]

Miss Forest

library(missForest)

## 
## Attaching package: 'missForest'

## The following object is masked from 'package:VIM':
## 
##     nrmse

library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ forcats 1.0.1     ✔ tibble  3.3.1
## ✔ purrr   1.2.2

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

MF_1 = missForest(as.matrix(uji_1))
MF_1 = as.data.frame(MF_1$ximp)
MF_2 = missForest(as.matrix(uji_2))
MF_2 = as.data.frame(MF_2$ximp)
MF_3 = missForest(as.matrix(uji_3))
MF_3 = as.data.frame(MF_3$ximp)

prediksi_mf_1 = MF_1[mask_1]
prediksi_mf_2 = MF_2[mask_2]
prediksi_mf_3 = MF_3[mask_3]

Evaluasi

RMSE

# RMSE 
rmse = function(asli, prediksi) sqrt(mean((asli - prediksi)^2, na.rm = TRUE))
rmse.mf_1 = rmse(nilai_asli_1, prediksi_mf_1)
rmse.mf_2 = rmse(nilai_asli_2, prediksi_mf_2)
rmse.mf_3 = rmse(nilai_asli_3, prediksi_mf_3)

Standar Deviasi

# standar deviasi
sd_asli_1 = sd(nilai_asli_1, na.rm = TRUE)
sd_asli_2 = sd(nilai_asli_2, na.rm = TRUE)
sd_asli_3 = sd(nilai_asli_3, na.rm = TRUE)

# Output
results = data.frame(
  Skema = c("Rendah", "Sedang", "Tinggi"),
  missForest = c(rmse.mf_1, rmse.mf_2, rmse.mf_3),
  SD = c(sd_asli_1, sd_asli_2, sd_asli_3)
)

results

##    Skema missForest       SD
## 1 Rendah   4.603724 24.88453
## 2 Sedang   5.362874 24.88799
## 3 Tinggi   5.365027 24.88799

Per Variabel

R-Square

# Dataframe 
eval_met_r2 = data.frame(
  Skema = character(),
  Variabel = character(),
  R_Squared = numeric(),
  stringsAsFactors = FALSE
)

# Daftar skema (Stasiun)
skema_met = list(
  "10%" = list(uji = uji_1, mf = MF_1),
  "40%" = list(uji = uji_2, mf = MF_2),
  "70%" = list(uji = uji_3, mf = MF_3)
)

# Looping perhitungan khusus R-Squared
for (nama_skema in names(skema_met)) {
  d_uji = skema_met[[nama_skema]]$uji
  d_mf = skema_met[[nama_skema]]$mf
  
  for (var in names(d_uji)) {
    idx_na = which(is.na(d_uji[[var]])) 
    
    if(length(idx_na) > 0) {
      aktual = gt_met[[var]][idx_na]    
      prediksi = d_mf[[var]][idx_na]    
      
      eval_met_r2 = rbind(eval_met_r2, data.frame(
        Skema = nama_skema,
        Variabel = var,
        R_Squared = round(cor(aktual, prediksi)^2, 4) # Dibulatkan 4 desimal
      ))
    }
  }
}

tabel_r2_stasiun = eval_met_r2 %>%
  pivot_wider(names_from = Skema, values_from = R_Squared, names_prefix = "R2_Skema_")

# Tampilkan hasil
print(tabel_r2_stasiun)

## # A tibble: 7 × 4
##   Variabel `R2_Skema_10%` `R2_Skema_40%` `R2_Skema_70%`
##   <chr>             <dbl>          <dbl>          <dbl>
## 1 Tn               0.509          0.290          0.280 
## 2 Tx               0.617          0.406          0.400 
## 3 Tavg             0.569          0.380          0.371 
## 4 RH_avg           0.665          0.488          0.480 
## 5 RR               0.144          0.0844         0.0824
## 6 ss               0.457          0.321          0.314 
## 7 ff_avg           0.0504         0.0284         0.0293

RMSE

# EVALUASI RMSE SPESIFIK PER VARIABEL 

hitung_rmse_per_var = function(data_asli, data_uji, data_imputasi) {
  rmse_list = list() # Buat list kosong untuk menyimpan hasil
  
  # Looping untuk setiap nama kolom di dataset
  for (kolom in names(data_asli)) {
    # Cari posisi baris yang nilainya sengaja dihilangkan (NA) HANYA untuk kolom ini
    mask_kolom = is.na(data_uji[[kolom]])
    
    # Cek apakah ada data yang dihilangkan di kolom tersebut
    if (sum(mask_kolom) > 0) {
      # Ambil nilai asli dan tebakan MissForest khusus di posisi yang di-NA-kan
      asli_kolom = data_asli[[kolom]][mask_kolom]
      pred_kolom = data_imputasi[[kolom]][mask_kolom]
      
      # Hitung RMSE khusus kolom ini
      rmse_val = sqrt(mean((asli_kolom - pred_kolom)^2, na.rm = TRUE))
      rmse_list[[kolom]] = rmse_val
    } else {
      # Jika tidak ada yang di-NA-kan, kembalikan NA
      rmse_list[[kolom]] = NA
    }
  }
  
  # Ubah list menjadi dataframe agar rapi
  return(as.data.frame(rmse_list))
}

# EVALUASI RMSE PER VARIABEL (Meteorologi)
rmse_met_10 = hitung_rmse_per_var(gt_met, uji_1, MF_1)
rmse_met_40 = hitung_rmse_per_var(gt_met, uji_2, MF_2)
rmse_met_70 = hitung_rmse_per_var(gt_met, uji_3, MF_3)

# Gabungkan hasil
tabel_rmse_met = rbind(rmse_met_10, rmse_met_40, rmse_met_70)

# Tabel
tabel_rmse_met$Skema = c("Rendah (10%)", "Sedang (40%)", "Tinggi (70%)")
tabel_rmse_met = tabel_rmse_met %>% select(Skema, everything())

# Hasil
print(tabel_rmse_met)

##          Skema       Tn        Tx      Tavg   RH_avg       RR       ss   ff_avg
## 1 Rendah (10%) 1.081451 0.9613954 0.5632190 4.416211 11.05712 1.717957 1.111184
## 2 Sedang (40%) 1.351834 1.1759131 0.7474947 5.643128 12.66877 1.949687 1.187436
## 3 Tinggi (70%) 1.363749 1.1817202 0.7557681 5.683354 12.65274 1.964099 1.186187

Gabungan (Meteorologi dan CHIRPS)

Input Data

# Data Meteorologi (variabel-variabel lain)
data_met = data_met %>% 
  select(-tahun)

str(data_met)

## tibble [12,053 × 8] (S3: tbl_df/tbl/data.frame)
##  $ Tanggal: POSIXct[1:12053], format: "1990-01-01" "1990-01-02" ...
##  $ Tn     : num [1:12053] 19 19 20 19 19 20 19 19 20 20 ...
##  $ Tx     : num [1:12053] 28.4 27 26.8 28.4 28.6 22 27.8 28.2 28.2 24.8 ...
##  $ Tavg   : num [1:12053] 22.2 22.5 23.1 22.5 21.5 20.3 22.3 23.3 23.2 21.8 ...
##  $ RH_avg : num [1:12053] 86 83 85 85 86 94 84 72 79 87 ...
##  $ RR     : num [1:12053] 6 47 0 12 16.5 1 13 1.5 2 0 ...
##  $ ss     : num [1:12053] 4.1 0.4 6.2 4.3 NA 0.1 5.2 2.8 5.1 0.4 ...
##  $ ff_avg : num [1:12053] 0 0 1 2 2 1 1 2 2 2 ...

# Data Chirps (satelit)
data_chrp = read_xlsx("D:\\deska_baru2\\DATA CH Bandung CHRIPS.xlsx")

# Ubah nama kolom ke-2 agar lebih mudah 
data_chrp = data_chrp %>%
  rename(RR_chrp = 2)

str(data_chrp)

## tibble [12,053 × 2] (S3: tbl_df/tbl/data.frame)
##  $ Tanggal: POSIXct[1:12053], format: "1990-01-01" "1990-01-02" ...
##  $ RR_chrp: num [1:12053] 13.22 3.86 4.41 11.57 3.86 ...

head(data_chrp)

## # A tibble: 6 × 2
##   Tanggal             RR_chrp
##   <dttm>                <dbl>
## 1 1990-01-01 00:00:00   13.2 
## 2 1990-01-02 00:00:00    3.86
## 3 1990-01-03 00:00:00    4.41
## 4 1990-01-04 00:00:00   11.6 
## 5 1990-01-05 00:00:00    3.86
## 6 1990-01-06 00:00:00    9.17

Pre-Proecessing Data

# Gabungkan data Meteorologi dengan data Satelit 
data_gab = data_met %>% 
  left_join(data_chrp, by = "Tanggal")

# Ubah kode 8888 jadi NA 
data_gab[data_gab==8888] = NA

str(data_gab)

## tibble [12,053 × 9] (S3: tbl_df/tbl/data.frame)
##  $ Tanggal: POSIXct[1:12053], format: "1990-01-01" "1990-01-02" ...
##  $ Tn     : num [1:12053] 19 19 20 19 19 20 19 19 20 20 ...
##  $ Tx     : num [1:12053] 28.4 27 26.8 28.4 28.6 22 27.8 28.2 28.2 24.8 ...
##  $ Tavg   : num [1:12053] 22.2 22.5 23.1 22.5 21.5 20.3 22.3 23.3 23.2 21.8 ...
##  $ RH_avg : num [1:12053] 86 83 85 85 86 94 84 72 79 87 ...
##  $ RR     : num [1:12053] 6 47 0 12 16.5 1 13 1.5 2 0 ...
##  $ ss     : num [1:12053] 4.1 0.4 6.2 4.3 NA 0.1 5.2 2.8 5.1 0.4 ...
##  $ ff_avg : num [1:12053] 0 0 1 2 2 1 1 2 2 2 ...
##  $ RR_chrp: num [1:12053] 13.22 3.86 4.41 11.57 3.86 ...

ncol(data_gab)

## [1] 9

EDA

Nilai Statistik

summary(data_gab)

##     Tanggal                 Tn              Tx             Tavg      
##  Min.   :1990-01-01   Min.   :13.00   Min.   :18.80   Min.   :18.40  
##  1st Qu.:1998-04-02   1st Qu.:19.00   1st Qu.:28.10   1st Qu.:22.80  
##  Median :2006-07-02   Median :20.00   Median :29.00   Median :23.40  
##  Mean   :2006-07-02   Mean   :19.38   Mean   :29.01   Mean   :23.39  
##  3rd Qu.:2014-10-01   3rd Qu.:20.00   3rd Qu.:30.00   3rd Qu.:24.00  
##  Max.   :2022-12-31   Max.   :29.00   Max.   :36.00   Max.   :28.90  
##                       NAs    :578     NAs    :326     NAs    :156    
##      RH_avg            RR                ss             ff_avg      
##  Min.   :42.00   Min.   :  0.000   Min.   : 0.000   Min.   : 0.000  
##  1st Qu.:73.00   1st Qu.:  0.000   1st Qu.: 3.300   1st Qu.: 1.000  
##  Median :79.00   Median :  0.800   Median : 5.100   Median : 2.000  
##  Mean   :77.99   Mean   :  6.738   Mean   : 4.924   Mean   : 1.736  
##  3rd Qu.:84.00   3rd Qu.:  7.800   3rd Qu.: 6.800   3rd Qu.: 2.000  
##  Max.   :96.00   Max.   :160.000   Max.   :10.000   Max.   :15.000  
##  NAs    :330     NAs    :1342      NAs    :492      NAs    :114     
##     RR_chrp      
##  Min.   : 0.000  
##  1st Qu.: 0.000  
##  Median : 3.079  
##  Mean   : 6.063  
##  3rd Qu.: 9.865  
##  Max.   :90.572  
##

Missing Value

# Missing Value 
data_gab = data_gab %>%
  mutate(tahun = year(Tanggal))
missing_tahunan2 = data_gab %>%
  group_by(tahun) %>%
  summarise(
    RR_na = sum(is.na(RR)),
    Tn_na = sum(is.na(Tn)),
    Tx_na = sum(is.na(Tx))
  ) 

library(DT)

# Membuat tabel interaktif
datatable(missing_tahunan2, 
          options = list(pageLength = 33), 
          caption = 'Tabel Missing Value Tahunan')

# Visualisasi
# A. Menggunakan naniar: Melihat peta sebaran data kosong
# Hitam = Data ada, Abu-abu = Data kosong (NA)
vis_miss(data_gab) + 
  labs(title = "Peta Sebaran Missing Value Data Meteorologi")

# B. Menggunakan VIM: Melihat kombinasi variabel yang sering kosong bersamaan
aggr(data_gab, 
     col = c('navyblue','red'), 
     numbers = TRUE, 
     sortVars = TRUE, 
     labels = names(data_gab), 
     cex.axis = .7, 
     gap = 3, 
     ylab = c("Proporsi Missing Data", "Pola Kombinasi Kosong"))

## Warning in plot.aggr(res, ...): not enough vertical space to display
## frequencies (too many combinations)

## 
##  Variables sorted by number of missings: 
##  Variable       Count
##        RR 0.111341575
##        Tn 0.047954866
##        ss 0.040819713
##    RH_avg 0.027379076
##        Tx 0.027047208
##      Tavg 0.012942836
##    ff_avg 0.009458226
##   Tanggal 0.000000000
##   RR_chrp 0.000000000
##     tahun 0.000000000

Korelasi

# Korelasi 
data_numeric2 = na.omit(data_gab[, 2:8]) # Buang kolom Tanggal dan hapus NA sementara khusus untuk hitung korelasi
data_numeric2 = na.omit(data_gab[, 2:8]) 
matriks_korelasi = cor(data_numeric2)

# Visualisasi 
ggcorrplot(matriks_korelasi, 
           method = "square", 
           type = "lower",      # Hanya menampilkan separuh bawah matriks
           lab = TRUE,          # Menampilkan angka korelasi
           colors = c("red", "white", "blue"),
           title = "Matriks Korelasi Variabel Meteorologi Bandung")

matriks_korelasi = cor(data_numeric2)

Distribusi Data

# Histogram 
data_long2 = data_gab %>%
  select(Tn, Tx, RR) %>%
  pivot_longer(cols = everything(), names_to = "Variabel", values_to = "Nilai")

# Visualisasi
ggplot(data_long2, aes(x = Variabel, y = Nilai, fill = Variabel)) +
  geom_boxplot(alpha = 0.7, outlier.colour = "red", outlier.shape = 16) +
  coord_flip() +  
  theme_minimal() +
  labs(title = "Distribusi Suhu (Tn, Tx) & Curah Hujan (RR)",
       x = "Variabel Iklim", 
       y = "Nilai Pengukuran")

## Warning: Removed 2246 rows containing non-finite outside the scale range
## (`stat_boxplot()`).

Tren Curah Hujan

# Line chart untuk Curah Hujan (RR)
ggplot(data_gab, aes(x = Tanggal, y = RR)) +
  geom_line(color = "dodgerblue", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Tren Runtun Waktu Curah Hujan (RR) Kota Bandung 1990-2022",
       x = "Tahun", 
       y = "Curah Hujan (mm)")

# Line chart untuk Suhu Rata-rata (Tavg)
ggplot(data_gab, aes(x = Tanggal, y = Tavg)) +
  geom_line(color = "firebrick", alpha = 0.8) +
  theme_minimal() +
  labs(title = "Tren Runtun Waktu Suhu Rata-rata (Tavg) Kota Bandung 1990-2022",
       x = "Tahun", 
       y = "Suhu Rata-rata (°C)")

Uji Algoritma

# Data Ground Truth
gt_gab = na.omit(data_gab[,2:8])
anyNA(gt_gab) # Cek data sudah bersih NA

## [1] FALSE

Skema Rendah (10%)

Duplikasi data dan hilangkan 10% nilai dalam data

set.seed(123)
uji1 = gt_gab

for (col in names(uji1)) {
  n_hilang = round(0.1 * nrow(uji1))
  idx = sample(1:nrow(uji1), n_hilang)
  uji1[idx, col] = NA
}

# Simpan nilai asli 10%
mask1 = is.na(uji1)
nilai_asli1 = gt_gab[mask1]

Skema Rendah (40%)

Duplikasi data dan hilangkan 40% nilai dalam data

set.seed(123)
uji2 = gt_gab

for (col in names(uji2)) {
  n_hilang = round(0.4 * nrow(uji2))
  idx = sample(1:nrow(uji2), n_hilang)
  uji2[idx, col] = NA
}

# Simpan nilai asli 40%
mask2 = is.na(uji2)
nilai_asli2 = gt_gab[mask2]

Skema Rendah (70%)

Duplikasi data dan hilangkan 70% nilai dalam data

set.seed(123)
uji3 = gt_gab

for (col in names(uji3)) {
  n_hilang = round(0.7 * nrow(uji3))
  idx = sample(1:nrow(uji3), n_hilang)
  uji3[idx, col] = NA
}

# Simpan nilai asli 70%
mask3 = is.na(uji3)
nilai_asli3 = gt_gab[mask3]

Miss Forest

MF1 = missForest(as.matrix(uji1))
MF1 = as.data.frame(MF1$ximp)
MF2 = missForest(as.matrix(uji2))
MF2 = as.data.frame(MF2$ximp)
MF3 = missForest(as.matrix(uji3))
MF3 = as.data.frame(MF3$ximp)

prediksi_mf1 = MF1[mask1]
prediksi_mf2 = MF2[mask2]
prediksi_mf3 = MF3[mask3]

Evaluasi

RMSE

# RMSE 
rmse = function(asli, prediksi) sqrt(mean((asli - prediksi)^2, na.rm = TRUE))
rmse.mf1 = rmse(nilai_asli1, prediksi_mf1)
rmse.mf2 = rmse(nilai_asli2, prediksi_mf2)
rmse.mf3 = rmse(nilai_asli3, prediksi_mf3)

Standar Deviasi

# standar deviasi
sd_asli1 = sd(nilai_asli1, na.rm = TRUE)
sd_asli2 = sd(nilai_asli2, na.rm = TRUE)
sd_asli3 = sd(nilai_asli3, na.rm = TRUE)

# Output
results = data.frame(
  Skema = c("Rendah", "Sedang", "Tinggi"),
  RMSE = c(rmse.mf1, rmse.mf2, rmse.mf3),
  SD = c(sd_asli1, sd_asli2, sd_asli3)
)

results

##    Skema     RMSE       SD
## 1 Rendah 4.551962 24.88453
## 2 Sedang 5.288929 24.88799
## 3 Tinggi 6.110610 24.91973

Per Variabel

R-square

# Dataframe
eval_gab = data.frame(
  Skema = character(),
  Variabel = character(),
  RMSE = numeric(),
  R_Squared = numeric(),
  stringsAsFactors = FALSE
)

# Daftar skema untuk dilooping (Gabungan)
skema_gab = list(
  "Rendah (10%)" = list(uji = uji1, mf = MF1),
  "Sedang (40%)" = list(uji = uji2, mf = MF2),
  "Tinggi (70%)" = list(uji = uji3, mf = MF3)
)

# Looping perhitungan
for (nama_skema in names(skema_gab)) {
  d_uji = skema_gab[[nama_skema]]$uji
  d_mf = skema_gab[[nama_skema]]$mf
  
  for (var in names(d_uji)) {
    idx_na = which(is.na(d_uji[[var]])) # Cari baris NA
    
    if(length(idx_na) > 0) {
      aktual = gt_gab[[var]][idx_na]    # Kunci jawaban dari gt_gab
      prediksi = d_mf[[var]][idx_na]    # Tebakan dari MissForest
      
      eval_gab = rbind(eval_gab, data.frame(
        Skema = nama_skema,
        Variabel = var,
        RMSE = rmse(aktual, prediksi),
        R_Squared = cor(aktual, prediksi)^2
      ))
    }
  }
}

# Tampilkan hasil
print(eval_gab)

##           Skema Variabel       RMSE   R_Squared
## 1  Rendah (10%)       Tn  1.0788798 0.511606005
## 2  Rendah (10%)       Tx  0.9628569 0.615128035
## 3  Rendah (10%)     Tavg  0.5669648 0.563665329
## 4  Rendah (10%)   RH_avg  4.4037469 0.666674149
## 5  Rendah (10%)       RR 10.9120820 0.160591771
## 6  Rendah (10%)       ss  1.7131708 0.459781748
## 7  Rendah (10%)   ff_avg  1.1083663 0.055447006
## 8  Sedang (40%)       Tn  1.3379208 0.293441960
## 9  Sedang (40%)       Tx  1.1620587 0.412671853
## 10 Sedang (40%)     Tavg  0.7440837 0.380426821
## 11 Sedang (40%)   RH_avg  5.6316101 0.485988087
## 12 Sedang (40%)       RR 12.4620256 0.086327351
## 13 Sedang (40%)       ss  1.9207178 0.332129533
## 14 Sedang (40%)   ff_avg  1.1870419 0.026645267
## 15 Tinggi (70%)       Tn  1.5095452 0.138112994
## 16 Tinggi (70%)       Tx  1.3135139 0.234542563
## 17 Tinggi (70%)     Tavg  0.8585830 0.183102883
## 18 Tinggi (70%)   RH_avg  6.7394621 0.259530964
## 19 Tinggi (70%)       RR 14.3110894 0.027322828
## 20 Tinggi (70%)       ss  2.1887322 0.174929811
## 21 Tinggi (70%)   ff_avg  1.2717971 0.007016194

RMSE

# Data 
rmse_var_10 = hitung_rmse_per_var(gt_gab, uji1, MF1)
rmse_var_40 = hitung_rmse_per_var(gt_gab, uji2, MF2)
rmse_var_70 = hitung_rmse_per_var(gt_gab, uji3, MF3)

# Gabungkan Hasil
tabel_rmse_variabel = rbind(rmse_var_10, rmse_var_40, rmse_var_70)

# Tabel
tabel_rmse_variabel$Skema = c("Rendah (10%)", "Sedang (40%)", "Tinggi (70%)")
tabel_rmse_variabel = tabel_rmse_variabel %>% select(Skema, everything())

# Hasil
print(tabel_rmse_variabel)

##          Skema       Tn        Tx      Tavg   RH_avg       RR       ss   ff_avg
## 1 Rendah (10%) 1.078880 0.9628569 0.5669648 4.403747 10.91208 1.713171 1.108366
## 2 Sedang (40%) 1.337921 1.1620587 0.7440837 5.631610 12.46203 1.920718 1.187042
## 3 Tinggi (70%) 1.509545 1.3135139 0.8585830 6.739462 14.31109 2.188732 1.271797

Imputasi

# Imputasi Meteorologi
#data1 = data_met %>% 
  #select(-Tanggal)
#MF.1 = missForest(as.matrix(data1))
#MF.1 = as.data.frame(MF.1$ximp)


# Imputasi Gabungan (Meteorologi + Chirps)
#data2 = data_gab %>% 
  #select(-Tanggal, -tahun)
#MF.2 = missForest(as.matrix(data2))
#MF.2 = as.data.frame(MF.2$ximp)

Data Final

# Meteorologi 
#final_met = data.frame(
  #Tanggal = data_met$Tanggal, 
  #CH = MF.1$RR,
  #Tn = MF.1$Tn,
  #Tx = MF.1$Tx
#)
#head(final_met)
#anyNA(final_met)

# Gabungan (Meteorologi + Satelit)
#final_gab = data.frame(
  #Tanggal = data_gab$Tanggal, 
  #CH = MF.2$RR,
  #Tn = MF.2$Tn,
  #Tx = MF.2$Tx
#)
#head(final_gab)
#anyNA(final_gab)

# Excel
#library(writexl)
#path_met = "D:\\deska_baru2\\Meteorologi.xlsx"
#path_gab = "D:\\deska_baru2\\Gabungan.xlsx"

# Simpan hasil imputasi stasiun saja ke folder Downloads
#write_xlsx(final_met, path_met)

# Simpan hasil imputasi gabungan (stasiun + satelit) ke folder Downloads
#write_xlsx(final_gab, path_gab)

# CSV
#library(readr)
#path_met_csv = "D:\\deska_baru2\\Meteorologi.csv"
#path_gab_csv = "D:\\deska_baru2\\Gabungan.csv"

# Meteorologi
#write_csv(final_met, path_met_csv)

# Gabungan
#write_csv(final_gab, path_gab_csv)

Imputasi Curah Hujan dengan Miss Forest

Muthi’ah Iffa Karimah

2026-04-10

Meteorologi

Input data

pre-processing Data

EDA

Nilai Statistik

Missing Value

Korelasi

Distribusi Data

Tren Curah hujan

Uji Algoritma MF: Skema Penghapusan Data

Skema Rendah (10%)

Skema Sedang (40%)

Skema Tinggi (70%)

Miss Forest

Evaluasi

RMSE

Standar Deviasi

Per Variabel

R-Square

RMSE

Gabungan (Meteorologi dan CHIRPS)

Input Data

Pre-Proecessing Data

EDA

Nilai Statistik

Missing Value

Korelasi

Distribusi Data

Tren Curah Hujan

Uji Algoritma

Skema Rendah (10%)

Skema Rendah (40%)

Skema Rendah (70%)

Miss Forest

Evaluasi

RMSE

Standar Deviasi

Per Variabel

R-square

RMSE

Imputasi

Data Final