PKL-Model Regresi Data Panel
A. Library dan Setup
knitr::opts_chunk$set(width = Inf)
options(tibble.width = Inf)
library(readxl)
library(dplyr)
library(psych)
library(ggplot2)
library(tidyr)
library(viridis)
library(corrplot)
library(patchwork)
library(scales)
library(reshape2)
library(plm)
library(tseries)
library(lmtest)
library(nortest)
library(sandwich)
library(car)B. Data
2.1. Load Dataset
# Load dataset
data_panel <- read_excel("D:/UNY/MySta/Bismillah-Tuntas/Magang_Lana/Laporan PKL/Dataset_Kependudukan.xlsx")
data_panel## # A tibble: 70 × 7
## Kemantren Tahun Kelahiran Kematian Migrasi_Masuk Migrasi_Keluar
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Tegalrejo 2020 375 250 16.1 18.2
## 2 Jetis 2020 271 235 15.7 17.9
## 3 Gondokusuman 2020 411 342 19.5 20.0
## 4 Danurejan 2020 176 192 15.3 18.3
## 5 Gedong Tengen 2020 168 194 10.9 17.6
## 6 Ngampilan 2020 192 157 14.4 17.0
## 7 Wirobrajan 2020 295 214 14.4 16.9
## 8 Mantrijeron 2020 361 264 14.2 16.3
## 9 Kraton 2020 200 183 13.4 15.6
## 10 Gondomanan 2020 166 141 14.8 16.6
## Pertumbuhan
## <dbl>
## 1 0.57
## 2 0.91
## 3 0.75
## 4 0.26
## 5 -0.03
## 6 -0.36
## 7 0.77
## 8 0.56
## 9 0.34
## 10 0.05
## # ℹ 60 more rows2.2. Struktur Data
## tibble [70 × 7] (S3: tbl_df/tbl/data.frame)
## $ Kemantren : chr [1:70] "Tegalrejo" "Jetis" "Gondokusuman" "Danurejan" ...
## $ Tahun : num [1:70] 2020 2020 2020 2020 2020 2020 2020 2020 2020 2020 ...
## $ Kelahiran : num [1:70] 375 271 411 176 168 192 295 361 200 166 ...
## $ Kematian : num [1:70] 250 235 342 192 194 157 214 264 183 141 ...
## $ Migrasi_Masuk : num [1:70] 16.1 15.7 19.5 15.3 10.9 ...
## $ Migrasi_Keluar: num [1:70] 18.2 17.9 20 18.3 17.6 ...
## $ Pertumbuhan : num [1:70] 0.57 0.91 0.75 0.26 -0.03 -0.36 0.77 0.56 0.34 0.05 ...2.3. Jumlah Observasi
# Cek Jumlah Observasi Tiap Kemantren dan Tahun
data_panel %>%
count(Kemantren) %>%
arrange(desc(n))## # A tibble: 14 × 2
## Kemantren n
## <chr> <int>
## 1 Danurejan 5
## 2 Gedong Tengen 5
## 3 Gondokusuman 5
## 4 Gondomanan 5
## 5 Jetis 5
## 6 Kotagede 5
## 7 Kraton 5
## 8 Mantrijeron 5
## 9 Mergangsan 5
## 10 Ngampilan 5
## 11 Pakualaman 5
## 12 Tegalrejo 5
## 13 Umbulharjo 5
## 14 Wirobrajan 5## # A tibble: 5 × 2
## Tahun n
## <dbl> <int>
## 1 2020 14
## 2 2021 14
## 3 2022 14
## 4 2023 14
## 5 2024 14C. Data Preprocessing
3.1. Missing Value dan Duplicate
## Kemantren Tahun Kelahiran Kematian Migrasi_Masuk
## 0 0 0 0 0
## Migrasi_Keluar Pertumbuhan
## 0 0## [1] 03.2. Outlier
# Cek Outlier
outlier_table <- data_panel %>%
group_by(Kemantren) %>%
summarise(across(
where(is.numeric),
~ sum(.x %in% boxplot.stats(.x)$out),
.names = "{.col}_out"
))
outlier_table## # A tibble: 14 × 7
## Kemantren Tahun_out Kelahiran_out Kematian_out Migrasi_Masuk_out
## <chr> <int> <int> <int> <int>
## 1 Danurejan 0 1 1 0
## 2 Gedong Tengen 0 0 1 1
## 3 Gondokusuman 0 0 1 2
## 4 Gondomanan 0 0 1 2
## 5 Jetis 0 2 1 1
## 6 Kotagede 0 0 2 2
## 7 Kraton 0 2 1 1
## 8 Mantrijeron 0 0 1 0
## 9 Mergangsan 0 2 2 1
## 10 Ngampilan 0 0 1 1
## 11 Pakualaman 0 0 1 0
## 12 Tegalrejo 0 0 1 1
## 13 Umbulharjo 0 1 0 0
## 14 Wirobrajan 0 2 1 0
## Migrasi_Keluar_out Pertumbuhan_out
## <int> <int>
## 1 1 1
## 2 0 0
## 3 1 1
## 4 1 0
## 5 0 0
## 6 1 2
## 7 1 0
## 8 0 0
## 9 1 0
## 10 0 0
## 11 0 1
## 12 0 0
## 13 1 1
## 14 0 0# Menampilkan baris yang mengandung outlier
outlier_rows <- data_panel %>%
group_by(Kemantren) %>%
filter(if_any(
where(is.numeric),
~ .x %in% boxplot.stats(.x)$out
))
outlier_rows## # A tibble: 35 × 7
## # Groups: Kemantren [14]
## Kemantren Tahun Kelahiran Kematian Migrasi_Masuk Migrasi_Keluar
## <chr> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 Jetis 2020 271 235 15.7 17.9
## 2 Gondokusuman 2020 411 342 19.5 20.0
## 3 Wirobrajan 2020 295 214 14.4 16.9
## 4 Kraton 2020 200 183 13.4 15.6
## 5 Mergangsan 2020 324 276 15.0 18.2
## 6 Kotagede 2020 382 237 17.2 16
## 7 Tegalrejo 2021 304 454 18.0 22.6
## 8 Jetis 2021 187 387 17.2 18.4
## 9 Gondokusuman 2021 289 526 23.0 24.1
## 10 Danurejan 2021 106 281 16.9 20.7
## Pertumbuhan
## <dbl>
## 1 0.91
## 2 0.75
## 3 0.77
## 4 0.34
## 5 -0.09
## 6 0.74
## 7 0.04
## 8 -0.35
## 9 0.09
## 10 -0.03
## # ℹ 25 more rowsPada tahun 2020 dan 2021 ditemukan sejumlah nilai ekstrem (outlier) pada variabel kelahiran dan kematian. Outlier tersebut bukan merupakan kesalahan input, melainkan merefleksikan kondisi pandemi COVID-19 yang menyebabkan lonjakan angka kematian serta perubahan dinamika kependudukan di beberapa kemantren.
D. Exploratory Data Analysis
4.1. Statistik Deskriptif
## Kemantren Tahun Kelahiran Kematian
## Length:70 Min. :2020 Min. : 64.0 Min. : 84.0
## Class :character 1st Qu.:2021 1st Qu.:155.8 1st Qu.:184.0
## Mode :character Median :2022 Median :225.0 Median :236.0
## Mean :2022 Mean :246.3 Mean :266.7
## 3rd Qu.:2023 3rd Qu.:310.0 3rd Qu.:308.0
## Max. :2024 Max. :767.0 Max. :758.0
## Migrasi_Masuk Migrasi_Keluar Pertumbuhan
## Min. :10.31 Min. :15.57 Min. :-2.170
## 1st Qu.:15.24 1st Qu.:17.39 1st Qu.:-0.675
## Median :17.17 Median :18.66 Median :-0.060
## Mean :17.98 Mean :20.14 Mean :-0.096
## 3rd Qu.:20.59 3rd Qu.:22.33 3rd Qu.: 0.545
## Max. :29.42 Max. :32.96 Max. : 1.040## Tahun Kelahiran Kematian Migrasi_Masuk Migrasi_Keluar
## 1.4244246 143.8486207 125.8161670 3.8200441 3.8644309
## Pertumbuhan
## 0.7010911## Tahun Kelahiran Kematian Migrasi_Masuk Migrasi_Keluar
## 2.028986e+00 2.069243e+04 1.582971e+04 1.459274e+01 1.493383e+01
## Pertumbuhan
## 4.915287e-01num_vars <- names(data_panel)[sapply(data_panel, is.numeric)]
# tahun saat nilai minimum dan maksimum terjadi
get_min_max_year <- function(df, var){
min_val <- min(df[[var]], na.rm = TRUE)
max_val <- max(df[[var]], na.rm = TRUE)
year_min <- df$Tahun[df[[var]] == min_val]
year_max <- df$Tahun[df[[var]] == max_val]
data.frame(
Variabel = var,
Min = min_val,
Tahun_Min = paste(unique(year_min), collapse = ", "),
Max = max_val,
Tahun_Max = paste(unique(year_max), collapse = ", "),
row.names = NULL
)
}
result <- do.call(rbind, lapply(num_vars, function(v) get_min_max_year(data_panel, v)))
print(result)## Variabel Min Tahun_Min Max Tahun_Max
## 1 Tahun 2020.00 2020 2024.00 2024
## 2 Kelahiran 64.00 2021 767.00 2020
## 3 Kematian 84.00 2020 758.00 2021
## 4 Migrasi_Masuk 10.31 2020 29.42 2022
## 5 Migrasi_Keluar 15.57 2024 32.96 2022
## 6 Pertumbuhan -2.17 2022 1.04 20224.2. Plot Distribusi Variabel
numeric_vars <- names(select(data_panel, where(is.numeric), -Kemantren))
plot_list <- lapply(numeric_vars, \(v)
ggplot(data_panel, aes(.data[[v]])) +
geom_histogram(bins = 15, fill = "#2b8cbe", color = "white") +
labs(title = paste("Distribusi", v), x = v, y = "Frekuensi") +
theme_minimal(base_size = 10) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, size = 10, color = "#045a8d"),
axis.text = element_text(size = 7)
)
)
wrap_plots(plotlist = plot_list, ncol = 2) +
plot_annotation(
title = "Distribusi Variabel Numerik",
theme = theme(plot.title = element_text(hjust = 0.5, size = 12, face = "bold"))
)
4.3. Plot Tren Per Tahun
plot_trend <- lapply(numeric_vars[-1], function(v) {
ggplot(data_panel, aes(x = Tahun, y = .data[[v]],
color = Kemantren, group = Kemantren)) +
geom_line(linewidth = 1) +
geom_point(size = 2) +
scale_color_viridis(discrete = TRUE, option = "C") +
labs(title = paste("Tren", v),
y = v, x = "Tahun") +
theme_minimal(base_size = 11) +
theme(
plot.title = element_text(face = "bold", hjust = 0.5, color = "#045a8d"),
legend.position = "right"
)
})
wrap_plots(plotlist = plot_trend, ncol = 2, guides = "collect") &
theme(legend.position = "bottom")
4.4. Korelasi Global
numeric_vars <- data_panel |>
select(where(is.numeric), -Tahun)
panel.cor <- function(x, y, digits = 4, cex = 1.3) {
r <- cor(x, y, use = "complete.obs")
txt <- formatC(r, format = "f", digits = digits)
par(usr = c(0,1,0,1))
text(0.5, 0.5, txt, cex = cex)
}
pairs(
numeric_vars,
upper.panel = panel.cor,
lower.panel = panel.smooth
)
E. Analisis Regresi Data Panel
5.1. Pembentukan Data Panel
## Panel structure:## Balanced Panel: n = 14, T = 5, N = 70##
## Danurejan Gedong Tengen Gondokusuman Gondomanan Jetis
## 5 5 5 5 5
## Kotagede Kraton Mantrijeron Mergangsan Ngampilan
## 5 5 5 5 5
## Pakualaman Tegalrejo Umbulharjo Wirobrajan
## 5 5 5 5##
## 2020 2021 2022 2023 2024
## 14 14 14 14 145.2. Estimasi Model
a. Common Effect Model (CEM)/ Pooled OLS
model_cem <- plm(Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar, data = pdata, model = "pooling")
summary(model_cem)## Pooling Model
##
## Call:
## plm(formula = Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk +
## Migrasi_Keluar, data = pdata, model = "pooling")
##
## Balanced Panel: n = 14, T = 5, N = 70
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.947130 -0.224430 -0.019293 0.228881 0.773227
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.32871767 0.23348664 5.6908 3.269e-07 ***
## Kelahiran 0.00183118 0.00044772 4.0900 0.0001212 ***
## Kematian -0.00038921 0.00046815 -0.8314 0.4087947
## Migrasi_Masuk 0.08945354 0.01685538 5.3071 1.441e-06 ***
## Migrasi_Keluar -0.16786186 0.01559556 -10.7634 4.402e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.915
## Residual Sum of Squares: 7.5559
## R-Squared: 0.77721
## Adj. R-Squared: 0.7635
## F-statistic: 56.6898 on 4 and 65 DF, p-value: < 2.22e-16b. Fixed Effect Model (FEM)
model_fem <- plm(Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar, data = pdata, model = "within")
summary(model_fem)## Oneway (individual) effect Within Model
##
## Call:
## plm(formula = Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk +
## Migrasi_Keluar, data = pdata, model = "within")
##
## Balanced Panel: n = 14, T = 5, N = 70
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.710737 -0.153707 0.010609 0.159802 0.651482
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## Kelahiran 0.00201161 0.00122195 1.6462 0.10575
## Kematian -0.00021728 0.00077379 -0.2808 0.77998
## Migrasi_Masuk 0.04306617 0.02107689 2.0433 0.04611 *
## Migrasi_Keluar -0.14672417 0.01819332 -8.0647 1.006e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 19.58
## Residual Sum of Squares: 5.1689
## R-Squared: 0.736
## Adj. R-Squared: 0.6497
## F-statistic: 36.2433 on 4 and 52 DF, p-value: 1.8428e-14## Danurejan Gedong Tengen Gondokusuman Gondomanan Jetis
## 1.7990 1.3872 2.0444 1.4072 1.7129
## Kotagede Kraton Mantrijeron Mergangsan Ngampilan
## 1.9430 1.3404 1.4819 1.4835 1.5364
## Pakualaman Tegalrejo Umbulharjo Wirobrajan
## 1.7357 1.7809 1.4835 1.9199c. Random Effect Model (REM)
model_rem <- plm(Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar, data = pdata, model = "random", random.method = "swar")
summary(model_rem)## Oneway (individual) effect Random Effect Model
## (Swamy-Arora's transformation)
##
## Call:
## plm(formula = Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk +
## Migrasi_Keluar, data = pdata, model = "random", random.method = "swar")
##
## Balanced Panel: n = 14, T = 5, N = 70
##
## Effects:
## var std.dev share
## idiosyncratic 0.0994 0.3153 1
## individual 0.0000 0.0000 0
## theta: 0
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.947130 -0.224430 -0.019293 0.228881 0.773227
##
## Coefficients:
## Estimate Std. Error z-value Pr(>|z|)
## (Intercept) 1.32871767 0.23348664 5.6908 1.265e-08 ***
## Kelahiran 0.00183118 0.00044772 4.0900 4.314e-05 ***
## Kematian -0.00038921 0.00046815 -0.8314 0.4058
## Migrasi_Masuk 0.08945354 0.01685538 5.3071 1.114e-07 ***
## Migrasi_Keluar -0.16786186 0.01559556 -10.7634 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.915
## Residual Sum of Squares: 7.5559
## R-Squared: 0.77721
## Adj. R-Squared: 0.7635
## Chisq: 226.759 on 4 DF, p-value: < 2.22e-165.3. Pemilihan Model Terbaik
a. Uji Chow (cem vs fem)
##
## F statistic
##
## data: Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar
## F = 1.8472, df1 = 13, df2 = 52, p-value = 0.05987
## alternative hypothesis: unstabilityPoolability Test (Chow Test)
H0: tidak ada efek individu yang berpengaruh terhadap model dengan model mengikuti CEM
H1 : Terdapat minimal satu efek individu berpengaruh terhadap model dengan FEM
Karena p-value = 0.059, sedikit di atas 0.05, maka H0 tidak ditolak. Oleh karena itu, dapat disimpulkan tidak ada efek individu berpengaruh terhadap model sehingga model mengikuti common effect model.
B. Uji Hausman (fem vs rem)
##
## Hausman Test
##
## data: Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar
## chisq = 5.0137, df = 4, p-value = 0.2859
## alternative hypothesis: one model is inconsistentHausman p = 0.2859 > 0.05, artinya FEM dan REM konsisten. Dengan Demikian, dipilih Random Effect Model (REM)
Kesimpulan Pemilihan Model Terbaik:
Meskipun Uji Hausman menyarankan random effect model (REM) apabila dibandingkan dengan FEM, hasil estimasi model REM pada penelitian ini ternyata identik dengan hasil CEM. Kesamaan tersebut terjadi karena ragam efek individual pada model REM bernilai nol sehingga model REM secara matematis menyatu dengan model common effect. Dengan demikian, common effect model merupakan model terbaik dalam penelitian ini.
F. Uji Asumsi Klasik
6.1. Uji Normalitas
##
## Jarque Bera Test
##
## data: resid(model_cem)
## X-squared = 0.8184, df = 2, p-value = 0.6642Interpretasi:
Karena p-value > 0.05, maka H₀ tidak ditolak, sehingga residual pada model REM berdistribusi normal. Dengan demikian, asumsi normalitas terpenuhi.
6.2. Uji Heteroskedastisitas
##
## studentized Breusch-Pagan test
##
## data: model_cem
## BP = 1.1471, df = 4, p-value = 0.8867Interpretasi:
Karena p-value > 0.05, maka H₀ tidak ditolak, sehingga terjadi homoskedastisitas pada model. Artinya, varians error homogen.
6.3. Uji Autokorelasi
##
## Durbin-Watson test for serial correlation in panel models
##
## data: Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk + Migrasi_Keluar
## DW = 2.0507, p-value = 0.5355
## alternative hypothesis: serial correlation in idiosyncratic errorsInterpretasi:
Nilai p-value > 0.05 menunjukkan bahwa tidak terdapat autokorelasi pada residual model. Dengan demikian, asumsi tidak ada autokorelasi terpenuhi.
6.4. Uji Multikolinearitas
## Kelahiran Kematian Migrasi_Masuk Migrasi_Keluar
## 2.462063 2.059259 2.460874 2.155999Interpretasi:
Berdasarkan kriteria Gujarati (VIF > 10 → terjadi multikolinearitas). Keempat variabel memiliki VIF < 10, maka asumsi tidak adanya multikolinearitas terpenuhi. Uji multikolinearitas dilakukan pada model Pooling karena VIF tidak dapat dihitung pada FEM/REM
G. Uji Signifikansi Parameter
7.1. Uji Serentak (F-Test)
## Pooling Model
##
## Call:
## plm(formula = Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk +
## Migrasi_Keluar, data = pdata, model = "pooling")
##
## Balanced Panel: n = 14, T = 5, N = 70
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.947130 -0.224430 -0.019293 0.228881 0.773227
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.32871767 0.23348664 5.6908 3.269e-07 ***
## Kelahiran 0.00183118 0.00044772 4.0900 0.0001212 ***
## Kematian -0.00038921 0.00046815 -0.8314 0.4087947
## Migrasi_Masuk 0.08945354 0.01685538 5.3071 1.441e-06 ***
## Migrasi_Keluar -0.16786186 0.01559556 -10.7634 4.402e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.915
## Residual Sum of Squares: 7.5559
## R-Squared: 0.77721
## Adj. R-Squared: 0.7635
## F-statistic: 56.6898 on 4 and 65 DF, p-value: < 2.22e-16pvalue < 0.05, maka secara simultan variabel X mempengaruhi variabel Y : Pertumbuhan Penduduk.
7.2. Uji Parsial (T-Test)
##
## t test of coefficients:
##
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.32871767 0.23348664 5.6908 3.269e-07 ***
## Kelahiran 0.00183118 0.00044772 4.0900 0.0001212 ***
## Kematian -0.00038921 0.00046815 -0.8314 0.4087947
## Migrasi_Masuk 0.08945354 0.01685538 5.3071 1.441e-06 ***
## Migrasi_Keluar -0.16786186 0.01559556 -10.7634 4.402e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 17.3. Koefisien Determinasi
## Pooling Model
##
## Call:
## plm(formula = Pertumbuhan ~ Kelahiran + Kematian + Migrasi_Masuk +
## Migrasi_Keluar, data = pdata, model = "pooling")
##
## Balanced Panel: n = 14, T = 5, N = 70
##
## Residuals:
## Min. 1st Qu. Median 3rd Qu. Max.
## -0.947130 -0.224430 -0.019293 0.228881 0.773227
##
## Coefficients:
## Estimate Std. Error t-value Pr(>|t|)
## (Intercept) 1.32871767 0.23348664 5.6908 3.269e-07 ***
## Kelahiran 0.00183118 0.00044772 4.0900 0.0001212 ***
## Kematian -0.00038921 0.00046815 -0.8314 0.4087947
## Migrasi_Masuk 0.08945354 0.01685538 5.3071 1.441e-06 ***
## Migrasi_Keluar -0.16786186 0.01559556 -10.7634 4.402e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Total Sum of Squares: 33.915
## Residual Sum of Squares: 7.5559
## R-Squared: 0.77721
## Adj. R-Squared: 0.7635
## F-statistic: 56.6898 on 4 and 65 DF, p-value: < 2.22e-16