library(readxl)    
library(tidyverse) 
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## āœ” dplyr     1.1.4     āœ” readr     2.1.5
## āœ” forcats   1.0.0     āœ” stringr   1.5.1
## āœ” ggplot2   4.0.0     āœ” tibble    3.3.0
## āœ” lubridate 1.9.4     āœ” tidyr     1.3.1
## āœ” purrr     1.1.0     
## ── 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
library(plm)       
## 
## Attaching package: 'plm'
## 
## The following objects are masked from 'package:dplyr':
## 
##     between, lag, lead
library(gplots)    
## 
## Attaching package: 'gplots'
## 
## The following object is masked from 'package:stats':
## 
##     lowess
library(stargazer) 
## 
## Please cite as: 
## 
##  Hlavac, Marek (2022). stargazer: Well-Formatted Regression and Summary Statistics Tables.
##  R package version 5.2.3. https://CRAN.R-project.org/package=stargazer
library(dplyr)
library(ggplot2)
data = read_excel("D:\\Semester 6\\ADP\\DATA PANEL_ (1).xlsx")
## New names:
## • `` -> `...10`
## • `` -> `...11`
## • `` -> `...12`
str(data)
## tibble [190 Ɨ 12] (S3: tbl_df/tbl/data.frame)
##  $ Provinsi: chr [1:190] "Aceh" NA NA NA ...
##  $ Tahun   : num [1:190] 2020 2021 2022 2023 2024 ...
##  $ Y       : chr [1:190] "31.8" "30.96" "26.6" "23.68" ...
##  $ X1      : num [1:190] 41.7 42.1 36.6 42.6 43.7 ...
##  $ X2      : chr [1:190] "0.51400000000000001" "0.503" "0.504" "0.48899999999999999" ...
##  $ X3      : chr [1:190] "9.1300000000000008" "9.2200000000000006" "9.24" "9.4" ...
##  $ X4      : chr [1:190] "131580.97" "135274.04" "140971.72" "146932.42000000001" ...
##  $ X5      : chr [1:190] "0.34" "0.35" "0.33" "0.33" ...
##  $ X6      : chr [1:190] "0.74" "0.56999999999999995" "0.5" "0.38" ...
##  $ ...10   : logi [1:190] NA NA NA NA NA NA ...
##  $ ...11   : chr [1:190] NA NA NA NA ...
##  $ ...12   : chr [1:190] NA NA NA NA ...

Persiapan Data

# isi kolom provinsi yang NA degan nilai di atasnya 
data = data %>% fill(Provinsi, .direction = "down")
head(data)
## # A tibble: 6 Ɨ 12
##   Provinsi     Tahun Y        X1 X2    X3    X4    X5    X6    ...10 ...11 ...12
##   <chr>        <dbl> <chr> <dbl> <chr> <chr> <chr> <chr> <chr> <lgl> <chr> <chr>
## 1 Aceh          2020 31.8   41.7 0.51… 9.13… 1315… 0.34  0.74  NA    <NA>  <NA> 
## 2 Aceh          2021 30.96  42.1 0.503 9.22… 1352… 0.35  0.56… NA    <NA>  <NA> 
## 3 Aceh          2022 26.6   36.6 0.504 9.24  1409… 0.33  0.5   NA    <NA>  <NA> 
## 4 Aceh          2023 23.68  42.6 0.48… 9.4   1469… 0.33  0.38  NA    <NA>  <NA> 
## 5 Aceh          2024 21.59  43.7 0.45… 9.49  1537… 0.33  0.42  NA    <NA>  <NA> 
## 6 Sumatera Ut…  2020 33.03  51.7 0.46… 9.27… 5337… 0.56… 1.09… NA    <NA>  <NA>
data = data %>% 
  select(-starts_with("...")) %>%
  mutate(
    `Y` = as.numeric(`Y`),
    `X1` = as.numeric(`X1`),
    `X2` = as.numeric(`X2`),
    `X3` = as.numeric(`X3`), 
    `X4` = as.numeric(`X4`),
    `X5` = as.numeric(`X5`),
    `X6` = as.numeric(`X6`)
  )
## Warning: There were 6 warnings in `mutate()`.
## The first warning was:
## ℹ In argument: `Y = as.numeric(Y)`.
## Caused by warning:
## ! NAs introduced by coercion
## ℹ Run `dplyr::last_dplyr_warnings()` to see the 5 remaining warnings.
str(data)
## tibble [190 Ɨ 9] (S3: tbl_df/tbl/data.frame)
##  $ Provinsi: chr [1:190] "Aceh" "Aceh" "Aceh" "Aceh" ...
##  $ Tahun   : num [1:190] 2020 2021 2022 2023 2024 ...
##  $ Y       : num [1:190] 31.8 31 26.6 23.7 21.6 ...
##  $ X1      : num [1:190] 41.7 42.1 36.6 42.6 43.7 ...
##  $ X2      : num [1:190] 0.514 0.503 0.504 0.489 0.459 0.468 0.445 0.442 0.425 0.399 ...
##  $ X3      : num [1:190] 9.13 9.22 9.24 9.4 9.49 9.28 9.33 9.5 9.63 9.77 ...
##  $ X4      : num [1:190] 131581 135274 140972 146932 153780 ...
##  $ X5      : num [1:190] 0.34 0.35 0.33 0.33 0.33 0.57 0.58 0.11 0 0 ...
##  $ X6      : num [1:190] 0.74 0.57 0.5 0.38 0.42 1.09 0.93 0.82 0.71 1.06 ...

Eksplorasi Data

head(data)
## # A tibble: 6 Ɨ 9
##   Provinsi       Tahun     Y    X1    X2    X3      X4    X5    X6
##   <chr>          <dbl> <dbl> <dbl> <dbl> <dbl>   <dbl> <dbl> <dbl>
## 1 Aceh            2020  31.8  41.7 0.514  9.13 131581.  0.34  0.74
## 2 Aceh            2021  31.0  42.1 0.503  9.22 135274.  0.35  0.57
## 3 Aceh            2022  26.6  36.6 0.504  9.24 140972.  0.33  0.5 
## 4 Aceh            2023  23.7  42.6 0.489  9.4  146932.  0.33  0.38
## 5 Aceh            2024  21.6  43.7 0.459  9.49 153780.  0.33  0.42
## 6 Sumatera Utara  2020  33.0  51.7 0.468  9.28 533746.  0.57  1.09
# Statistika Deskriptif
summary(data)
##    Provinsi             Tahun            Y               X1       
##  Length:190         Min.   :2020   Min.   :16.04   Min.   :34.54  
##  Class :character   1st Qu.:2021   1st Qu.:29.72   1st Qu.:44.15  
##  Mode  :character   Median :2022   Median :33.62   Median :47.37  
##                     Mean   :2022   Mean   :33.61   Mean   :48.47  
##                     3rd Qu.:2023   3rd Qu.:38.29   3rd Qu.:51.73  
##                     Max.   :2024   Max.   :46.28   Max.   :81.34  
##                                    NA's   :16      NA's   :16     
##        X2               X3               X4                X5        
##  Min.   :0.1420   Min.   : 3.400   Min.   :  13172   Min.   :0.0000  
##  1st Qu.:0.4235   1st Qu.: 7.872   1st Qu.:  64008   1st Qu.:0.3425  
##  Median :0.4745   Median : 8.550   Median : 141776   Median :0.4650  
##  Mean   :0.4568   Mean   : 8.523   Mean   : 333834   Mean   :0.5018  
##  3rd Qu.:0.5218   3rd Qu.: 9.195   3rd Qu.: 356561   3rd Qu.:0.5900  
##  Max.   :0.6650   Max.   :11.190   Max.   :2151041   Max.   :1.3300  
##  NA's   :16       NA's   :16       NA's   :12        NA's   :20      
##        X6       
##  Min.   :0.380  
##  1st Qu.:0.930  
##  Median :1.160  
##  Mean   :1.183  
##  3rd Qu.:1.410  
##  Max.   :2.360  
##  NA's   :17
pdata = pdata.frame(data, index = c ("Provinsi", "Tahun"))
class(pdata)
## [1] "pdata.frame" "data.frame"
pdim(pdata)
## Balanced Panel: n = 38, T = 5, N = 190
# Cek Heterogentitas
plotmeans(Y ~ Provinsi, data = data,
          xlab = "Provinsi", ylab = "Y",
          main = "Heterogenitas Antar Provinsi",
          bars = TRUE, p = 0.95, minbar = TRUE, mean.labels = FALSE, 
          xaxt = "n") 
## Warning in qt((1 + p)/2, ns - 1): NaNs produced
# Label sumbu X dibuat miring agar terbaca
axis(1, at = 1:length(unique(data$Provinsi)), 
     labels = unique(data$Provinsi), las = 2, cex.axis = 0.7)

# Cek jumlah data valid (bukan NA) per Provinsi
cek_data = data %>%
  filter(!is.na(Y)) %>%        # Hanya ambil data Y yang ada isinya
  group_by(Provinsi) %>%
  summarise(
    Jumlah_Data = n(),         # Hitung berapa tahun datanya
    Variasi = var(Y)           # Hitung varians (keragaman)
  ) %>%
  arrange(Jumlah_Data)         # Urutkan dari yang paling sedikit

# Tampilkan hasil
print(cek_data, n = 34)
## # A tibble: 38 Ɨ 3
##    Provinsi             Jumlah_Data Variasi
##    <chr>                      <int>   <dbl>
##  1 Papua Barat Daya               1   NA   
##  2 Papua Pegunungan               1   NA   
##  3 Papua Selatan                  1   NA   
##  4 Papua Tengah                   1   NA   
##  5 Aceh                           5   19.8 
##  6 Bali                           5   11.4 
##  7 Banten                         5   45.8 
##  8 Bengkulu                       5   14.6 
##  9 DI Yogyakarta                  5    5.17
## 10 DKI Jakarta                    5   42.4 
## 11 Gorontalo                      5   21.6 
## 12 Jambi                          5   12.2 
## 13 Jawa Barat                     5   34.3 
## 14 Jawa Tengah                    5    7.11
## 15 Jawa Timur                     5    5.73
## 16 Kalimantan Barat               5   13.4 
## 17 Kalimantan Selatan             5   26.2 
## 18 Kalimantan Tengah              5   24.3 
## 19 Kalimantan Utara               5   27.4 
## 20 Kalimantan timur               5   20.7 
## 21 Kep. Bangka Belitung           5   28.6 
## 22 Kep. Riau                      5   10.9 
## 23 Lampung                        5   17.5 
## 24 Maluku                         5   12.6 
## 25 Maluku Utara                   5   19.1 
## 26 NTB                            5   23.9 
## 27 NTT                            5   11.8 
## 28 Papua                          5   91.8 
## 29 Papua Barat                    5   14.9 
## 30 Riau                           5   35.0 
## 31 Sulawesi Barat                 5   19.7 
## 32 Sulawesi Selatan               5   16.3 
## 33 Sulawesi Tengah                5   11.5 
## 34 Sulawesi Tenggara              5   15.6 
## # ℹ 4 more rows
# Visualisasi Tren Waktu
ggplot(data, aes(x = Tahun, y = Y, group = Provinsi)) +
  geom_line(color = "grey", alpha = 0.6) + 
  geom_point(color = "grey", size = 1, alpha = 0.6) + 
  stat_summary(fun =  mean, geom = "line", aes(group = 1), color = "red", size = 1.2) +
  theme_minimal() +
  labs(
    title = "Tren Persentase Pemuda Kawin (20220-2024)",
    subtitle ="Garis Abu: per provinsi | Garis Merah: Rataan Nasional",
    y = "Persentase Pemuda Kawin",
    x = "Tahun"
  )
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_summary()`).
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_line()`).
## Warning: Removed 16 rows containing missing values or values outside the scale range
## (`geom_point()`).

  • garis merah menandakan bahwa rata-rata presentase pemuda kawin di tiap provinsi di Indonesia cenderung menurun dari waktu ke waktu.

  • ada beberapa garis abu yang menaik tajam/menurun tajam, bisa jadi perhatian lebih untuk didalami.

    # 1. Ubah struktur data X menjadi memanjang ke bawah (Long Format)
    # Ini trik agar kita bisa mem-plot semuanya sekaligus
    data_long = data %>%
      select(Y, X1, X2, X3, X4, X5, X6) %>%
      pivot_longer(
        cols = starts_with("X"),      # Pilih semua kolom yang diawali huruf X
        names_to = "Variabel_X",      # Nama kolom baru untuk label (X1, X2, dst)
        values_to = "Nilai_X"         # Nama kolom baru untuk angkanya
      )
    
    # 2. Buat Plot Gabungan (Facet Wrap)
    ggplot(data_long, aes(x = Nilai_X, y = Y)) +
      geom_point(alpha = 0.4, color = "darkblue", size = 1) + # Titik data
      geom_smooth(method = "lm", color = "red", se = FALSE) + # Garis regresi linier
      facet_wrap(~Variabel_X, scales = "free_x", ncol = 3) +  # Memecah gambar per variabel X
      theme_minimal() +
      labs(
        title = "Hubungan Y (Pemuda Kawin) dengan Semua Variabel X",
        subtitle = "Garis Merah: Tren Regresi Linier",
        y = "Persen Pemuda Kawin (Y)",
        x = "Nilai Variabel X"
      ) +
      theme(strip.text = element_text(face = "bold", size = 10))  
    ## `geom_smooth()` using formula = 'y ~ x'
    ## Warning: Removed 101 rows containing non-finite outside the scale range
    ## (`stat_smooth()`).
    ## Warning: Removed 101 rows containing missing values or values outside the scale range
    ## (`geom_point()`).

Garis merah merupakan slope yang menunjukkan hubungan Y dengan variabel peubah

A. Hubungan Negatif

  • Y vs X1 = Y (Persen Pemuda Kawin) meningkat, maka X1 (TPAK Perempuan) menurun.

  • Y vs X3 = Y (Persen Pemuda Kawin) meningkat, maka X3 (RLS Perempuan) menurun.

  • Y vs X4 = Y (Persen Pemuda Kawin) meningkat, maka X4 (PDRB Konstan) menurun.

  • Y vs X5 = Y (Persen Pemuda Kawin) meningkat, maka X5 (Urbanisasi) menurun.

B. Hubungan Positif

  • Y vs X2 = Y (Persen Pemuda Kawin) menurun, maka X2 (IKG) meningkat.

  • Y vs X6 = Y (Persen Pemuda Kawin) menurun, maka X6 (Cerai) meningkat.

# Membuang Provinsi yang datanya < 2 Tahun 
data_clean = data %>% 
  group_by(Provinsi) %>%
  filter(n() > 1) %>% # hanya ambil provinsi dengan data > 1
  ungroup()
table(data_clean$Provinsi)
## 
##                 Aceh                 Bali               Banten 
##                    5                    5                    5 
##             Bengkulu        DI Yogyakarta          DKI Jakarta 
##                    5                    5                    5 
##            Gorontalo                Jambi           Jawa Barat 
##                    5                    5                    5 
##          Jawa Tengah           Jawa Timur     Kalimantan Barat 
##                    5                    5                    5 
##   Kalimantan Selatan    Kalimantan Tengah     Kalimantan timur 
##                    5                    5                    5 
##     Kalimantan Utara Kep. Bangka Belitung            Kep. Riau 
##                    5                    5                    5 
##              Lampung               Maluku         Maluku Utara 
##                    5                    5                    5 
##                  NTB                  NTT                Papua 
##                    5                    5                    5 
##          Papua Barat     Papua Barat Daya     Papua Pegunungan 
##                    5                    5                    5 
##        Papua Selatan         Papua Tengah                 Riau 
##                    5                    5                    5 
##       Sulawesi Barat     Sulawesi Selatan      Sulawesi Tengah 
##                    5                    5                    5 
##    Sulawesi Tenggara       Sulawesi Utara       Sumatera Barat 
##                    5                    5                    5 
##     Sumatera Selatan       Sumatera Utara 
##                    5                    5
# Plot Korelasi: Melihat hubungan antar variabel X
library(corrplot)
## corrplot 0.95 loaded
# Ambil hanya kolom numerik
corr_data = data_clean %>%
  select(Y, X1, X2, X3, X4, X5, X6) %>%
  na.omit() # Buang NA sebentar untuk hitung korelasi

# Hitung matriks
M = cor(corr_data)

# Visualisasi yang cantik
corrplot(M, method = "number", type = "upper", 
         tl.col = "black", tl.srt = 45, 
         title = "Matriks Korelasi", mar=c(0,0,2,0))

# Buat ranking Provinsi
ranking_provinsi = data_clean %>% 
  group_by(Provinsi) %>%
  summarise(Rata_Rata_Y = mean(Y, na.rm = TRUE)) %>%
  arrange(desc(Rata_Rata_Y))

# Barplot 5 Provinsi dengan Persentase Pemuda Kawin tertinggi dan Terendah
top_bottom = rbind(
  head(ranking_provinsi, 5) %>% mutate(Status = "Tertinggi"),
  tail(ranking_provinsi, 5) %>% mutate(Status = "Terendah")
)

# 2. Buat Plot
ggplot(top_bottom, aes(x = reorder(Provinsi, Rata_Rata_Y), y = Rata_Rata_Y, fill = Status)) +
  geom_col() +
  coord_flip() +
  scale_fill_manual(values = c("Tertinggi" = "red", "Terendah" = "darkgreen")) +
  geom_text(aes(label = round(Rata_Rata_Y, 1)), hjust = -0.2) +
  theme_minimal() +
  labs(
    title = "Kontras Ekstrem: Provinsi Tertinggi vs Terendah",
    y = "Persentase Pemuda Kawin",
    x = ""
  )

5 Provinsi dengan Persentase Pemuda Kawin Tertinggi

  1. NTB
  2. Kalimantan Tengah
  3. Jawa Timur
  4. Bengkulu
  5. Kep. Bangka Belitung

5 Provinsi dengan Persentase Pemuda Kawin Terendah

  1. Aceh
  2. Papua Barat Daya
  3. DKI Jakarta
  4. Papua Tengah
  5. Papua Pegunungan
# Plot Korelasi: Melihat hubungan antar variabel X
library(corrplot)

# Ambil hanya kolom numerik
corr_data = data_clean %>%
  select(X1, X2, X3, X4, X5, X6) %>%
  na.omit() # Buang NA sebentar untuk hitung korelasi

# Hitung matriks
M = cor(corr_data)

# Visualisasi yang cantik
corrplot(M, method = "number", type = "upper", 
         tl.col = "black", tl.srt = 45, 
         title = "Matriks Korelasi", mar=c(0,0,2,0))

  • TPAK Perempuan (\(X1\)) vs IKG (\(X2\))= Korelasi Negatif kuat

  • Pendidikan (\(X3\) ) vs IKG (\(X2\)) = Korelasi Negatif cukup kuat.

  • IKG (\(X2\)) vs Urbanisasi (\(X5\)) = Korelasi Negatif kuat.

  • Cerai (\(X6\)) vs Pendidikan (\(X3\)) = Korelasi Negatif cukup kuat.

  • Ekonomi (\(X4\)) vs Urbanisasi (\(X5\)) = Korelasi Positif cukup kuat.

  • Pendidikan (\(X3\)) vs Urbanisasi (\(X5\)) = Korelasi Positif cukup kuat.

  • Cerai (\(X6\)) vs IKG (\(X2\)) = Korelasi Positif cukup kuat.

# 1. Tentukan variabel mana saja yang mau dicek
# (X5 kita skip karena isinya 0, tidak bisa di-log)
vars_x =  c("X1", "X2", "X3", "X4", "X5", "X6")
par(mfrow = c(1, 2)) # 1 Baris, 2 Kolom
for (var in vars_x) {
  nilai = data[[var]]

  # Histogram
  hist(nilai, 
       main = paste("Distribusi ASLI:", var),
       xlab = var, 
       col = "lightblue", 
       border = "white",
       breaks = 15) # Breaks biar batangnya lebih detail
  
  # Periksa pakah semua angkanya > 0? (Log butuh angka positif)
  if (all(nilai > 0, na.rm = TRUE)) {
    hist(log(nilai), 
         main = paste("Distribusi LOG:", var),
         xlab = paste("Log(", var, ")"), 
         col = "lightgreen", 
         border = "white",
         breaks = 15)
  } else {
    # Jika ada angka 0 atau negatif, tampilkan pesan error di gambar
    plot.new()
    text(0.5, 0.5, "Data mengandung 0/ Negatif\nTidak bisa di-Log", col = "red", cex = 1.2)
  }
  
  # Beri jeda sedikit biar kamu sempat lihat (opsional, bisa dihapus)
  Sys.sleep(0.5) 
}

# 4. Kembalikan tata letak gambar ke normal (1 gambar per layar)
par(mfrow = c(1, 1))
# Buat kolom baru "Wilayah"
data2 = data_clean %>%
  mutate(Wilayah = ifelse(Provinsi %in% c("DKI Jakarta", "Jawa Barat", "Jawa Tengah", 
                                          "DI Yogyakarta", "Jawa Timur", "Banten"), 
                          "Jawa", "Luar Jawa"))

# Plot Boxplot membandingkan Jawa vs Luar Jawa
ggplot(data2, aes(x = Wilayah, y = Y, fill = Wilayah)) +
  geom_boxplot() +
  labs(title = "Perbandingan Persentase Pemuda Kawin: Jawa vs Luar Jawa",
       y = "Persentase Pemuda Kawin") +
  theme_minimal()
## Warning: Removed 16 rows containing non-finite outside the scale range
## (`stat_boxplot()`).