1 Library

# import data
library(readxl)
library(readr)

# eksplorasi dan manipulasi data
library(skimr)
library(rlang)
library(tidyr)
library(dplyr)
library(DescTools)

# analisis regresi data panel
library(plm)

# pengujian asumsi dan inferensi statistik
library(lmtest)
library(sandwich)
library(car)
library(tseries)

# visualisasi data
library(ggplot2)

2 Input Data

# -------------------------
# Tetapkan Variabel 
# -------------------------

file_path <- "AngkasaPura.csv"   
id_var   <- "Bandara"               
time_var <- "Bulan"                 
y_var  <- "Bagasi_kg"                     
x_var <- c("Penumpang_pax")
alpha <- 0.05                       

# -------------------------
# Jadikan Dataframe
# -------------------------
df <- read.csv(file_path, header = TRUE, sep = ",")
df <- as.data.frame(df)


# -------------------------
# Bentuk sebagai Data Panel
# -------------------------
pdata <- pdata.frame(df, index = c(id_var, time_var))

2.1 Preview Data

DataBandara <- read.csv("AngkasaPura.csv", header = TRUE, sep=",")

Rows: 156
Columns: 4
$ Bandara       <chr> "JOG", "JOG", "JOG", "JOG", "JOG", "JOG", "JOG", "JOG", …
$ Bulan         <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5, 6,…
$ Bagasi_kg     <int> 7752894, 6642019, 6868575, 6941794, 8351969, 8069393, 11…
$ Penumpang_pax <int> 1039287, 958092, 1025520, 1038947, 1210433, 1085312, 130…

pdata_DataBandara <- pdata.frame(DataBandara, index = c("Bandara", "Bulan"))

index(pdata_DataBandara)

pdim(pdata_DataBandara)

Balanced Panel: n = 13, T = 12, N = 156

3 Statistik Deskriptif

   Bandara              Bulan         Bagasi_kg        Penumpang_pax    
 Length:156         Min.   : 1.00   Min.   :  344691   Min.   :  38806  
 Class :character   1st Qu.: 3.75   1st Qu.: 2563269   1st Qu.: 312718  
 Mode  :character   Median : 6.50   Median : 3978169   Median : 582034  
                    Mean   : 6.50   Mean   : 8194582   Mean   : 958164  
                    3rd Qu.: 9.25   3rd Qu.:10160061   3rd Qu.:1204518  
                    Max.   :12.00   Max.   :37696730   Max.   :3518057

skim(DataBandara)

Data summary
Name	DataBandara
Number of rows	156
Number of columns	4
_______________________
Column type frequency:
character	1
numeric	3
________________________
Group variables	None

Variable type: character

skim_variable	n_missing	complete_rate	min	max	empty	n_unique	whitespace
Bandara	0	1	3	3	0	13	0

Variable type: numeric

skim_variable	complete_rate	mean	sd	p0	p25	p50	p75	p100	hist
Bulan	1	6.5	3.46	1	3.75	6.5	9.25	12	▇▅▅▅▇
Bagasi_kg	1	8194581.5	9005086.64	344691	2563269.25	3978169.0	10160061.00	37696730	▇▂▁▁▁
Penumpang_pax	1	958163.8	948867.59	38806	312718.25	582034.0	1204518.50	3518057	▇▃▁▁▁

3.0.1 Per Bandara

DataBandara %>%
  group_by(Bandara) %>%
  summarise(
    mean_bagasi = mean(Bagasi_kg),
    sd_bagasi = sd(Bagasi_kg),
    mean_penumpang = mean(Penumpang_pax),
    sd_penumpang = sd(Penumpang_pax),
    min_penumpang = min(Penumpang_pax),
    max_penumpang = max(Penumpang_pax)
  )

3.0.2 Per Bulan

DataBandara %>%
  group_by(Bulan) %>%
  summarise(
    avg_penumpang = mean(Penumpang_pax),
    avg_bagasi = mean(Bagasi_kg)
  )

4 Visualisasi Data Eksplorasi

4.0.1 Tren Penumpang per Bulan

ggplot(DataBandara, aes(x = Bulan, y = Penumpang_pax, color = Bandara)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(
    title = "Tren Jumlah Penumpang per Bulan",
    x = "Bulan",
    y = "Jumlah Penumpang"
  ) +
  theme_minimal()

4.0.2 Tren Bagasi per Bulan

ggplot(DataBandara, aes(x = Bulan, y = Bagasi_kg, color = Bandara)) +
  geom_line(linewidth = 1) +
  geom_point(size = 2) +
  labs(
    title = "Tren Volume Bagasi per Bulan",
    x = "Bulan",
    y = "Bagasi (kg)"
  ) +
  theme_minimal()

4.0.3 Scatter Plot : Bagasi v Penumpang

ggplot(DataBandara, aes(x = Penumpang_pax, y = Bagasi_kg, color = Bandara)) +
  geom_point(size = 3, alpha = 0.8) +
  geom_smooth(method = "lm", se = FALSE, linewidth = 1) +
  labs(
    title = "Hubungan Jumlah Penumpang dan Volume Bagasi",
    x = "Jumlah Penumpang",
    y = "Bagasi (kg)"
  ) +
  theme_minimal()

`geom_smooth()` using formula = 'y ~ x'

4.0.4 Heatmap bulan x bandara

ggplot(DataBandara, aes(x = Bulan, y = Bandara, fill = Penumpang_pax)) +
  geom_tile() +
  labs(
    title = "Heatmap Jumlah Penumpang per Bulan dan Bandara",
    x = "Bulan",
    y = "Bandara",
    fill = "Penumpang"
  ) +
  theme_minimal()

5 Pemilihan Model Regresi Panel

# -------------------------
# BENTUK FORMULA
# -------------------------
formula_panel <- as.formula(
  paste(y_var, "~", paste(x_var, collapse = " + "))
)

# -------------------------
# ESTIMASI MODEL
# -------------------------
model_cem <- plm(formula_panel, data = pdata, model = "pooling")
model_fem <- plm(formula_panel, data = pdata, model = "within")
model_rem <- plm(formula_panel, data = pdata, model = "random")

# -------------------------
# UJI PEMILIHAN MODEL
# -------------------------

# Chow Test: CEM vs FEM
chow_test <- pFtest(model_fem, model_cem)

# LM Test: CEM vs REM
lm_test <- plmtest(model_cem, type = "bp")

# Hausman Test: FEM vs REM
hausman_test <- phtest(model_fem, model_rem)

# -------------------------
# LOGIKA PEMILIHAN MODEL
# -------------------------

# - Chow dulu
#   - jika tidak signifikan -> CEM vs REM pakai LM
#   - jika signifikan -> FEM vs REM pakai Hausman

if (chow_test$p.value > alpha) {
  
  # CEM cukup, lanjut cek apakah REM lebih baik dari CEM
  if (lm_test$p.value > alpha) {
    selected_model <- model_cem
    selected_name  <- "CEM / Pooled OLS"
    decision_note  <- "Chow tidak signifikan -> CEM lebih tepat daripada FEM; LM tidak signifikan -> CEM lebih tepat daripada REM."
  } else {
    selected_model <- model_rem
    selected_name  <- "REM / Random Effect Model"
    decision_note  <- "Chow tidak signifikan -> CEM lebih tepat daripada FEM; LM signifikan -> REM lebih tepat daripada CEM."
  }
  
} else {
  
  # FEM lebih baik dari CEM, lanjut cek FEM vs REM
  if (hausman_test$p.value > alpha) {
    selected_model <- model_rem
    selected_name  <- "REM / Random Effect Model"
    decision_note  <- "Chow signifikan -> FEM lebih tepat daripada CEM; Hausman tidak signifikan -> REM lebih tepat daripada FEM."
  } else {
    selected_model <- model_fem
    selected_name  <- "FEM / Fixed Effect Model"
    decision_note  <- "Chow signifikan -> FEM lebih tepat daripada CEM; Hausman signifikan -> FEM lebih tepat daripada REM."
  }
}


==============================
HASIL UJI PEMILIHAN MODEL
==============================

1. Chow Test (CEM vs FEM)

    F test for individual effects

data:  formula_panel
F = 142.49, df1 = 12, df2 = 142, p-value < 2.2e-16
alternative hypothesis: significant effects

p-value Chow Test = 8.091963e-73 

2. LM Test (CEM vs REM)

    Lagrange Multiplier Test - (Breusch-Pagan)

data:  formula_panel
chisq = 684.86, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

p-value LM Test = 5.85414e-151 

3. Hausman Test (FEM vs REM)

    Hausman Test

data:  formula_panel
chisq = 13.51, df = 1, p-value = 0.0002373
alternative hypothesis: one model is inconsistent

p-value Hausman Test = 0.0002372822 

==============================
MODEL TERPILIH : FEM / Fixed Effect Model 
==============================
Chow signifikan -> FEM lebih tepat daripada CEM; Hausman signifikan -> FEM lebih tepat daripada REM. 

==============================
SUMMARY MODEL TERPILIH
==============================
Oneway (individual) effect Within Model

Call:
plm(formula = formula_panel, data = pdata, model = "within")

Balanced Panel: n = 13, T = 12, N = 156

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1234731  -207994    -7163        0   160032  2055122 

Coefficients:
              Estimate Std. Error t-value  Pr(>|t|)    
Penumpang_pax  11.2771     0.3123   36.11 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    3.3112e+14
Residual Sum of Squares: 3.2518e+13
R-Squared:      0.90179
Adj. R-Squared: 0.8928
F-statistic: 1303.94 on 1 and 142 DF, p-value: < 2.22e-16

6 Uji Asumsi (Homoskedastisitas, Non-Autokorelasi, dan Cross-Sectional Dependence)

cat("\n==============================\n")
cat("UJI ASUMSI MODEL PANEL TERPILIH\n")
cat("==============================\n")

# ------------------------------
# A. Ambil residual model
# ------------------------------
resid_model <- residuals(selected_model)

# ------------------------------
# B. UJI HETEROSKEDASTISITAS
# ------------------------------
hetero_test <- bptest(selected_model)

cat("\n------------------------------\n")
cat("1. UJI HETEROSKEDASTISITAS\n")
cat("------------------------------\n")
print(hetero_test)

if (hetero_test$p.value > alpha) {
  keputusan_hetero <- "Tidak terjadi heteroskedastisitas (p-value > alpha)."
  hetero_flag <- FALSE
} else {
  keputusan_hetero <- "Terjadi heteroskedastisitas (p-value <= alpha)."
  hetero_flag <- TRUE
}
cat(keputusan_hetero, "\n")

# ------------------------------
# C. UJI AUTOKORELASI PANEL
# ------------------------------
auto_test <- pbgtest(selected_model)

cat("\n------------------------------\n")
cat("2. UJI AUTOKORELASI PANEL\n")
cat("------------------------------\n")
print(auto_test)

if (auto_test$p.value > alpha) {
  keputusan_auto <- "Tidak terjadi autokorelasi panel (p-value > alpha)."
  auto_flag <- FALSE
} else {
  keputusan_auto <- "Terjadi autokorelasi panel (p-value <= alpha)."
  auto_flag <- TRUE
}
cat(keputusan_auto, "\n")

# ------------------------------
# D. UJI CROSS-SECTIONAL DEPENDENCE
# ------------------------------
csd_test <- pcdtest(selected_model, test = "cd")

cat("\n-------------------------------------------\n")
cat("3. UJI CROSS-SECTIONAL DEPENDENCE\n")
cat("-------------------------------------------\n")
print(csd_test)

if (csd_test$p.value > alpha) {
  keputusan_csd <- "Tidak terdapat cross-sectional dependence (p-value > alpha)."
  csd_flag <- FALSE
} else {
  keputusan_csd <- "Terdapat cross-sectional dependence (p-value <= alpha)."
  csd_flag <- TRUE
}
cat(keputusan_csd, "\n")

# ==============================
# RINGKASAN KEPUTUSAN UJI ASUMSI
# ==============================
cat("\n==============================\n")
cat("RINGKASAN KEPUTUSAN UJI ASUMSI\n")
cat("==============================\n")
cat("Heteroskedastisitas : ", keputusan_hetero, "\n")
cat("Autokorelasi panel  : ", keputusan_auto, "\n")
cat("Cross-sectional dep : ", keputusan_csd, "\n")


==============================
UJI ASUMSI MODEL PANEL TERPILIH
==============================

------------------------------
1. UJI HETEROSKEDASTISITAS
------------------------------

    studentized Breusch-Pagan test

data:  selected_model
BP = 99.408, df = 1, p-value < 2.2e-16

Terjadi heteroskedastisitas (p-value <= alpha). 

------------------------------
2. UJI AUTOKORELASI PANEL
------------------------------

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  formula_panel
chisq = 46.925, df = 12, p-value = 4.804e-06
alternative hypothesis: serial correlation in idiosyncratic errors

Terjadi autokorelasi panel (p-value <= alpha). 

-------------------------------------------
3. UJI CROSS-SECTIONAL DEPENDENCE
-------------------------------------------

    Pesaran CD test for cross-sectional dependence in panels

data:  Bagasi_kg ~ Penumpang_pax
z = 15.832, p-value < 2.2e-16
alternative hypothesis: cross-sectional dependence

Terdapat cross-sectional dependence (p-value <= alpha). 

==============================
RINGKASAN KEPUTUSAN UJI ASUMSI
==============================
Heteroskedastisitas :  Terjadi heteroskedastisitas (p-value <= alpha). 
Autokorelasi panel  :  Terjadi autokorelasi panel (p-value <= alpha). 
Cross-sectional dep :  Terdapat cross-sectional dependence (p-value <= alpha).

7 Robust Method

cat("\n==============================\n")
cat("HASIL AKHIR ESTIMASI\n")
cat("==============================\n")
cat('\n')
cat("Model terpilih:", selected_name, "\n")
cat("Terdapat heteroskedastisitas, autokorelasi panel, dan cross-sectional dependence.\n")
cat("Oleh karena itu, inferensi akhir menggunakan Driscoll-Kraay robust standard error.\n\n")

final_result <- coeftest(
  selected_model,
  vcov = vcovSCC(selected_model, type = "HC1")
)

print(final_result)


==============================
HASIL AKHIR ESTIMASI
==============================

Model terpilih: FEM / Fixed Effect Model 
Terdapat heteroskedastisitas, autokorelasi panel, dan cross-sectional dependence.
Oleh karena itu, inferensi akhir menggunakan Driscoll-Kraay robust standard error.


t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
Penumpang_pax 11.27706    0.74315  15.175 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

8 Ringkasan

Hasil uji asumsi pada model panel terpilih menunjukkan adanya heteroskedastisitas, autokorelasi panel, dan cross-sectional dependence. Oleh karena itu, pengujian signifikansi parameter dilakukan menggunakan Driscoll-Kraay robust standard error agar inferensi statistik tetap valid.

Karena ketiga pelanggaran asumsi tersebut terdeteksi secara simultan, standard error konvensional tidak lagi digunakan sebagai dasar inferensi. Sebagai gantinya, penelitian ini menggunakan Driscoll-Kraay robust standard error yang robust terhadap heteroskedastisitas, autokorelasi, dan cross-sectional dependence.

9 Pemisah

10 Pemisah

11 Pemisah

12 Pemisah

13 (Shortcut) Otomatis Semua Mua Nya

# =========================
# AUTO SELECT PANEL MODEL
# =========================

# -------------------------
# 1) INPUT USER
# -------------------------
file_path <- "AngkasaPura.csv"   

id_var   <- "Bandara"               
time_var <- "Bulan"                 

y_var  <- "Bagasi_kg"                     
x_var <- c("Penumpang_pax")

alpha <- 0.05                       

# -------------------------
# 2) READ DATA
# -------------------------
df <- read.csv(file_path, header = TRUE, sep = ",")
df <- as.data.frame(df)


# -------------------------
# 3) BENTUK DATA PANEL
# -------------------------
pdata <- pdata.frame(df, index = c(id_var, time_var))

# 4) BENTUK FORMULA
# -------------------------
formula_panel <- as.formula(
  paste(y_var, "~", paste(x_var, collapse = " + "))
)

# -------------------------
# 5) ESTIMASI MODEL
# -------------------------
model_cem <- plm(formula_panel, data = pdata, model = "pooling")
model_fem <- plm(formula_panel, data = pdata, model = "within")
model_rem <- plm(formula_panel, data = pdata, model = "random")

# -------------------------
# 6) UJI PEMILIHAN MODEL
# -------------------------

# Chow Test: CEM vs FEM
chow_test <- pFtest(model_fem, model_cem)

# LM Test: CEM vs REM
lm_test <- plmtest(model_cem, type = "bp")

# Hausman Test: FEM vs REM
hausman_test <- phtest(model_fem, model_rem)

# -------------------------
# 7) LOGIKA PEMILIHAN MODEL
# -------------------------
# Urutan:
# - Chow dulu
#   - jika tidak signifikan -> CEM vs REM pakai LM
#   - jika signifikan -> FEM vs REM pakai Hausman

if (chow_test$p.value > alpha) {
  
  # CEM cukup, lanjut cek apakah REM lebih baik dari CEM
  if (lm_test$p.value > alpha) {
    selected_model <- model_cem
    selected_name  <- "CEM / Pooled OLS"
    decision_note  <- "Chow tidak signifikan -> CEM lebih tepat daripada FEM; LM tidak signifikan -> CEM lebih tepat daripada REM."
  } else {
    selected_model <- model_rem
    selected_name  <- "REM / Random Effect Model"
    decision_note  <- "Chow tidak signifikan -> CEM lebih tepat daripada FEM; LM signifikan -> REM lebih tepat daripada CEM."
  }
  
} else {
  
  # FEM lebih baik dari CEM, lanjut cek FEM vs REM
  if (hausman_test$p.value > alpha) {
    selected_model <- model_rem
    selected_name  <- "REM / Random Effect Model"
    decision_note  <- "Chow signifikan -> FEM lebih tepat daripada CEM; Hausman tidak signifikan -> REM lebih tepat daripada FEM."
  } else {
    selected_model <- model_fem
    selected_name  <- "FEM / Fixed Effect Model"
    decision_note  <- "Chow signifikan -> FEM lebih tepat daripada CEM; Hausman signifikan -> FEM lebih tepat daripada REM."
  }
}

# -------------------------
# 8) OUTPUT HASIL
# -------------------------
cat("\n==============================\n")
cat("HASIL UJI PEMILIHAN MODEL\n")
cat("==============================\n")

cat("\n1. Chow Test (CEM vs FEM)\n")
print(chow_test)
cat("p-value Chow Test =", chow_test$p.value, "\n")

cat("\n2. LM Test (CEM vs REM)\n")
print(lm_test)
cat("p-value LM Test =", lm_test$p.value, "\n")

cat("\n3. Hausman Test (FEM vs REM)\n")
print(hausman_test)
cat("p-value Hausman Test =", hausman_test$p.value, "\n")

cat("\n==============================\n")
cat("MODEL TERPILIH :", selected_name, "\n")
cat("==============================\n")
cat(decision_note, "\n")

cat("\n==============================\n")
cat("SUMMARY MODEL TERPILIH\n")
cat("==============================\n")
print(summary(selected_model))


==============================
HASIL UJI PEMILIHAN MODEL
==============================

1. Chow Test (CEM vs FEM)

    F test for individual effects

data:  formula_panel
F = 142.49, df1 = 12, df2 = 142, p-value < 2.2e-16
alternative hypothesis: significant effects

p-value Chow Test = 8.091963e-73 

2. LM Test (CEM vs REM)

    Lagrange Multiplier Test - (Breusch-Pagan)

data:  formula_panel
chisq = 684.86, df = 1, p-value < 2.2e-16
alternative hypothesis: significant effects

p-value LM Test = 5.85414e-151 

3. Hausman Test (FEM vs REM)

    Hausman Test

data:  formula_panel
chisq = 13.51, df = 1, p-value = 0.0002373
alternative hypothesis: one model is inconsistent

p-value Hausman Test = 0.0002372822 

==============================
MODEL TERPILIH : FEM / Fixed Effect Model 
==============================
Chow signifikan -> FEM lebih tepat daripada CEM; Hausman signifikan -> FEM lebih tepat daripada REM. 

==============================
SUMMARY MODEL TERPILIH
==============================
Oneway (individual) effect Within Model

Call:
plm(formula = formula_panel, data = pdata, model = "within")

Balanced Panel: n = 13, T = 12, N = 156

Residuals:
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-1234731  -207994    -7163        0   160032  2055122 

Coefficients:
              Estimate Std. Error t-value  Pr(>|t|)    
Penumpang_pax  11.2771     0.3123   36.11 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Total Sum of Squares:    3.3112e+14
Residual Sum of Squares: 3.2518e+13
R-Squared:      0.90179
Adj. R-Squared: 0.8928
F-statistic: 1303.94 on 1 and 142 DF, p-value: < 2.22e-16

# ==============================
# 9) UJI ASUMSI MODEL TERPILIH
# ==============================

cat("\n==============================\n")
cat("UJI ASUMSI MODEL TERPILIH\n")
cat("==============================\n")

# ------------------------------
# A. Ambil residual model
# ------------------------------
resid_model <- residuals(selected_model)


# ------------------------------
# B. HETEROSKEDASTISITAS
# ------------------------------
# bptest pada model panel terpilih
hetero_test <- bptest(selected_model)

cat("\n------------------------------\n")
cat("1. UJI HETEROSKEDASTISITAS\n")
cat("------------------------------\n")
print(hetero_test)

if (hetero_test$p.value > alpha) {
  keputusan_hetero <- "Tidak terjadi heteroskedastisitas (p-value > alpha)."
  hetero_flag <- FALSE
} else {
  keputusan_hetero <- "Terjadi heteroskedastisitas (p-value <= alpha)."
  hetero_flag <- TRUE
}
cat(keputusan_hetero, "\n")

# ------------------------------
# C. AUTOKORELASI
# ------------------------------
# Uji Breusch-Godfrey/Wooldridge style untuk panel
auto_test <- pbgtest(selected_model)

cat("\n------------------------------\n")
cat("2. UJI AUTOKORELASI\n")
cat("------------------------------\n")
print(auto_test)

if (auto_test$p.value > alpha) {
  keputusan_auto <- "Tidak terjadi autokorelasi (p-value > alpha)."
  auto_flag <- FALSE
} else {
  keputusan_auto <- "Terjadi autokorelasi (p-value <= alpha)."
  auto_flag <- TRUE
}
cat(keputusan_auto, "\n")

# ------------------------------
# E. CROSS-SECTIONAL DEPENDENCE
# ------------------------------
csd_test <- pcdtest(selected_model, test = "cd")

cat("\n-------------------------------------------\n")
cat("3. UJI CROSS-SECTIONAL DEPENDENCE\n")
cat("-------------------------------------------\n")
print(csd_test)

if (csd_test$p.value > alpha) {
  keputusan_csd <- "Tidak terdapat cross-sectional dependence (p-value > alpha)."
  csd_flag <- FALSE
} else {
  keputusan_csd <- "Terdapat cross-sectional dependence (p-value <= alpha)."
  csd_flag <- TRUE
}
cat(keputusan_csd, "\n")

# ------------------------------
# F. NORMALITAS RESIDUAL
# ------------------------------
normal_test <- jarque.bera.test(resid_model)

cat("\n------------------------------\n")
cat("4. UJI NORMALITAS RESIDUAL\n")
cat("------------------------------\n")
print(normal_test)

if (normal_test$p.value > alpha) {
  keputusan_normal <- "Residual berdistribusi normal (p-value > alpha)."
  normal_flag <- FALSE
} else {
  keputusan_normal <- "Residual tidak berdistribusi normal (p-value <= alpha)."
  normal_flag <- TRUE
}
cat(keputusan_normal, "\n")

# ==============================
# 10) KEPUTUSAN AKHIR UJI ASUMSI
# ==============================
cat("\n==============================\n")
cat("RINGKASAN KEPUTUSAN UJI ASUMSI\n")
cat("==============================\n")
cat("Heteroskedastisitas: ", keputusan_hetero, "\n")
cat("Autokorelasi       : ", keputusan_auto, "\n")
cat("Cross-sectional dep: ", keputusan_csd, "\n")
cat("Normalitas residual: ", keputusan_normal, "\n")

# ==============================
# 11) ROBUST STANDARD ERROR JIKA DIPERLUKAN
# ==============================
# Kalau ada heteroskedastisitas / autokorelasi / cross-sectional dependence,
# maka tampilkan koefisien dengan robust standard error

if (hetero_flag == FALSE && auto_flag == FALSE && csd_flag == FALSE) {
  
  cat("\n==============================\n")
  cat("KESIMPULAN AKHIR\n")
  cat("==============================\n")
  cat("Model terpilih:", selected_name, "\n")
  cat("Model lolos asumsi utama, sehingga summary biasa dapat digunakan.\n")
  
  cat("\nSUMMARY MODEL FINAL\n")
  print(summary(selected_model))
  
} else {
  
  cat("\n==============================\n")
  cat("KESIMPULAN AKHIR\n")
  cat("==============================\n")
  cat("Model terpilih:", selected_name, "\n")
  cat("Terdapat pelanggaran asumsi pada model, sehingga digunakan robust standard error.\n")
  
  # Robust vcov:
  # - HC1 untuk heteroskedastisitas
  # - method = arellano cukup aman untuk panel
  robust_result <- coeftest(
    selected_model,
    vcov = vcovHC(selected_model, method = "arellano", type = "HC1", cluster = "group")
  )
  
  cat("\nHASIL ESTIMASI DENGAN ROBUST STANDARD ERROR\n")
  print(robust_result)
}


==============================
UJI ASUMSI MODEL TERPILIH
==============================

------------------------------
1. UJI HETEROSKEDASTISITAS
------------------------------

    studentized Breusch-Pagan test

data:  selected_model
BP = 99.408, df = 1, p-value < 2.2e-16

Terjadi heteroskedastisitas (p-value <= alpha). 

------------------------------
2. UJI AUTOKORELASI
------------------------------

    Breusch-Godfrey/Wooldridge test for serial correlation in panel models

data:  formula_panel
chisq = 46.925, df = 12, p-value = 4.804e-06
alternative hypothesis: serial correlation in idiosyncratic errors

Terjadi autokorelasi (p-value <= alpha). 

-------------------------------------------
3. UJI CROSS-SECTIONAL DEPENDENCE
-------------------------------------------

    Pesaran CD test for cross-sectional dependence in panels

data:  Bagasi_kg ~ Penumpang_pax
z = 15.832, p-value < 2.2e-16
alternative hypothesis: cross-sectional dependence

Terdapat cross-sectional dependence (p-value <= alpha). 

------------------------------
4. UJI NORMALITAS RESIDUAL
------------------------------

    Jarque Bera Test

data:  resid_model
X-squared = 78.256, df = 2, p-value < 2.2e-16

Residual tidak berdistribusi normal (p-value <= alpha). 

==============================
RINGKASAN KEPUTUSAN UJI ASUMSI
==============================
Heteroskedastisitas:  Terjadi heteroskedastisitas (p-value <= alpha). 
Autokorelasi       :  Terjadi autokorelasi (p-value <= alpha). 
Cross-sectional dep:  Terdapat cross-sectional dependence (p-value <= alpha). 
Normalitas residual:  Residual tidak berdistribusi normal (p-value <= alpha). 

==============================
KESIMPULAN AKHIR
==============================
Model terpilih: FEM / Fixed Effect Model 
Terdapat pelanggaran asumsi pada model, sehingga digunakan robust standard error.

HASIL ESTIMASI DENGAN ROBUST STANDARD ERROR

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
Penumpang_pax 11.27706    0.87636  12.868 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

if (hetero_flag == FALSE && auto_flag == FALSE && csd_flag == FALSE) {
  
  cat("\n==============================\n")
  cat("KESIMPULAN AKHIR\n")
  cat("==============================\n")
  cat("Model terpilih:", selected_name, "\n")
  cat("Model lolos asumsi utama, sehingga summary biasa dapat digunakan.\n")
  
  print(summary(selected_model))
  
} else {
  
  cat("\n==============================\n")
  cat("KESIMPULAN AKHIR\n")
  cat("==============================\n")
  cat("Model terpilih:", selected_name, "\n")
  cat("Terdapat pelanggaran asumsi pada model, sehingga digunakan robust standard error.\n")
  
  if (csd_flag == TRUE) {
    final_result <- coeftest(
      selected_model,
      vcov = vcovSCC(selected_model, type = "HC1")
    )
    
    cat("\nHASIL ESTIMASI DENGAN DRISCOLL-KRAAY STANDARD ERROR\n")
    print(final_result)
    
  } else {
    final_result <- coeftest(
      selected_model,
      vcov = vcovHC(selected_model, method = "arellano", type = "HC1", cluster = "group")
    )
    
    cat("\nHASIL ESTIMASI DENGAN ROBUST STANDARD ERROR\n")
    print(final_result)
  }
}


==============================
KESIMPULAN AKHIR
==============================
Model terpilih: FEM / Fixed Effect Model 
Terdapat pelanggaran asumsi pada model, sehingga digunakan robust standard error.

HASIL ESTIMASI DENGAN DRISCOLL-KRAAY STANDARD ERROR

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
Penumpang_pax 11.27706    0.74315  15.175 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

14 Uji Asumsi Custom

14.1 (Glesjer/White, DurWatson, CSD pakai pcd)

# ==============================
# 9) UJI ASUMSI (VERSI LENGKAP / EKSPLORATIF)
# ==============================

cat("\n==============================\n")
cat("UJI ASUMSI MODEL TERPILIH\n")
cat("==============================\n")

# Pastikan data terurut berdasarkan id dan waktu
df <- df[order(df[[id_var]], df[[time_var]]), ]
pdata <- pdata.frame(df, index = c(id_var, time_var))

# Ambil residual dari model terpilih
resid_model <- as.numeric(residuals(selected_model))

# Simpan residual ke data frame agar mudah dipakai
df$resid_model     <- resid_model
df$abs_resid_model <- abs(resid_model)

# =========================================
# A. UJI HETEROSKEDASTISITAS - GLEJSER
# =========================================
# Regresi |residual| terhadap variabel independen
formula_glejser <- as.formula(
  paste("abs_resid_model ~", paste(x_var, collapse = " + "))
)

glejser_model <- lm(formula_glejser, data = df)
glejser_summary <- summary(glejser_model)
glejser_fstat <- glejser_summary$fstatistic
glejser_pval <- pf(glejser_fstat[1], glejser_fstat[2], glejser_fstat[3], lower.tail = FALSE)

cat("\n------------------------------\n")
cat("1. UJI HETEROSKEDASTISITAS - GLEJSER\n")
cat("------------------------------\n")
print(glejser_summary)

cat("\nP-value uji F Glejser =", glejser_pval, "\n")

if (glejser_pval > alpha) {
  keputusan_glejser <- "Tidak terdapat indikasi heteroskedastisitas menurut uji Glejser (p-value > alpha)."
  glejser_flag <- FALSE
} else {
  keputusan_glejser <- "Terdapat indikasi heteroskedastisitas menurut uji Glejser (p-value <= alpha)."
  glejser_flag <- TRUE
}
cat(keputusan_glejser, "\n")


# =========================================
# B. UJI HETEROSKEDASTISITAS - WHITE
# =========================================
# White test eksploratif dengan bptest pada model OLS bantu
# Dibuat dari formula utama + kuadrat variabel X
# (untuk 1 variabel X seperti kasus kamu ini aman dan sederhana)

white_formula <- as.formula(
  paste(
    y_var, "~",
    paste(c(x_var, paste0("I(", x_var, "^2)")), collapse = " + ")
  )
)

white_lm <- lm(white_formula, data = df)
white_test <- bptest(white_lm)

cat("\n------------------------------\n")
cat("2. UJI HETEROSKEDASTISITAS - WHITE\n")
cat("------------------------------\n")
print(white_test)

if (white_test$p.value > alpha) {
  keputusan_white <- "Tidak terdapat indikasi heteroskedastisitas menurut uji White (p-value > alpha)."
  white_flag <- FALSE
} else {
  keputusan_white <- "Terdapat indikasi heteroskedastisitas menurut uji White (p-value <= alpha)."
  white_flag <- TRUE
}
cat(keputusan_white, "\n")


# =========================================
# C. UJI AUTOKORELASI - DURBIN WATSON
# =========================================
# DW dibuat secara eksploratif pada model OLS bantu
dw_lm <- lm(formula_panel, data = df)
dw_test <- dwtest(dw_lm)

cat("\n------------------------------\n")
cat("3. UJI AUTOKORELASI - DURBIN WATSON\n")
cat("------------------------------\n")
print(dw_test)

if (dw_test$p.value > alpha) {
  keputusan_dw <- "Tidak terdapat indikasi autokorelasi menurut uji Durbin-Watson (p-value > alpha)."
  dw_flag <- FALSE
} else {
  keputusan_dw <- "Terdapat indikasi autokorelasi menurut uji Durbin-Watson (p-value <= alpha)."
  dw_flag <- TRUE
}
cat(keputusan_dw, "\n")


# =========================================
# D. UJI CROSS-SECTIONAL DEPENDENCE - PCD
# =========================================
csd_test <- pcdtest(selected_model, test = "cd")

cat("\n-------------------------------------------\n")
cat("4. UJI CROSS-SECTIONAL DEPENDENCE - PESARAN CD\n")
cat("-------------------------------------------\n")
print(csd_test)

if (csd_test$p.value > alpha) {
  keputusan_csd <- "Tidak terdapat cross-sectional dependence (p-value > alpha)."
  csd_flag <- FALSE
} else {
  keputusan_csd <- "Terdapat cross-sectional dependence (p-value <= alpha)."
  csd_flag <- TRUE
}
cat(keputusan_csd, "\n")


# =========================================
# E. RINGKASAN UJI ASUMSI
# =========================================
cat("\n==============================\n")
cat("RINGKASAN KEPUTUSAN UJI ASUMSI\n")
cat("==============================\n")
cat("Glejser (hetero)    : ", keputusan_glejser, "\n")
cat("White (hetero)      : ", keputusan_white, "\n")
cat("Durbin-Watson (auto): ", keputusan_dw, "\n")
cat("PCD Test (CSD)      : ", keputusan_csd, "\n")


# =========================================
# F. HASIL AKHIR ROBUST SE (OPSIONAL)
# =========================================
# Kalau ada CSD -> vcovSCC
# Kalau tidak ada CSD tapi ada indikasi hetero/auto -> vcovHC Arellano
# Kalau semua aman -> summary biasa

if (csd_flag == TRUE) {
  
  cat("\n==============================\n")
  cat("HASIL AKHIR\n")
  cat("==============================\n")
  cat("Karena terdapat cross-sectional dependence, digunakan Driscoll-Kraay robust standard error.\n")
  
  final_result <- coeftest(
    selected_model,
    vcov = vcovSCC(selected_model, type = "HC1")
  )
  
  print(final_result)
  
} else if (glejser_flag == TRUE || white_flag == TRUE || dw_flag == TRUE) {
  
  cat("\n==============================\n")
  cat("HASIL AKHIR\n")
  cat("==============================\n")
  cat("Karena terdapat indikasi heteroskedastisitas/autokorelasi, digunakan robust standard error Arellano.\n")
  
  final_result <- coeftest(
    selected_model,
    vcov = vcovHC(selected_model, method = "arellano", type = "HC1", cluster = "group")
  )
  
  print(final_result)
  
} else {
  
  cat("\n==============================\n")
  cat("HASIL AKHIR\n")
  cat("==============================\n")
  cat("Tidak ada indikasi pelanggaran serius pada uji eksploratif, sehingga summary model biasa dapat digunakan.\n")
  
  final_result <- summary(selected_model)
  print(final_result)
}


==============================
UJI ASUMSI MODEL TERPILIH
==============================

------------------------------
1. UJI HETEROSKEDASTISITAS - GLEJSER
------------------------------

Call:
lm(formula = formula_glejser, data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-816076 -128391  -59089  115307 1343213 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.986e+04  2.970e+04   2.689  0.00797 ** 
Penumpang_pax 2.339e-01  2.206e-02  10.604  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 260600 on 154 degrees of freedom
Multiple R-squared:  0.422, Adjusted R-squared:  0.4183 
F-statistic: 112.5 on 1 and 154 DF,  p-value: < 2.2e-16


P-value uji F Glejser = 4.543495e-20 
Terdapat indikasi heteroskedastisitas menurut uji Glejser (p-value <= alpha). 

------------------------------
2. UJI HETEROSKEDASTISITAS - WHITE
------------------------------

    studentized Breusch-Pagan test

data:  white_lm
BP = 114.85, df = 2, p-value < 2.2e-16

Terdapat indikasi heteroskedastisitas menurut uji White (p-value <= alpha). 

------------------------------
3. UJI AUTOKORELASI - DURBIN WATSON
------------------------------

    Durbin-Watson test

data:  dw_lm
DW = 0.26696, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0

Terdapat indikasi autokorelasi menurut uji Durbin-Watson (p-value <= alpha). 

-------------------------------------------
4. UJI CROSS-SECTIONAL DEPENDENCE - PESARAN CD
-------------------------------------------

    Pesaran CD test for cross-sectional dependence in panels

data:  Bagasi_kg ~ Penumpang_pax
z = 15.832, p-value < 2.2e-16
alternative hypothesis: cross-sectional dependence

Terdapat cross-sectional dependence (p-value <= alpha). 

==============================
RINGKASAN KEPUTUSAN UJI ASUMSI
==============================
Glejser (hetero)    :  Terdapat indikasi heteroskedastisitas menurut uji Glejser (p-value <= alpha). 
White (hetero)      :  Terdapat indikasi heteroskedastisitas menurut uji White (p-value <= alpha). 
Durbin-Watson (auto):  Terdapat indikasi autokorelasi menurut uji Durbin-Watson (p-value <= alpha). 
PCD Test (CSD)      :  Terdapat cross-sectional dependence (p-value <= alpha). 

==============================
HASIL AKHIR
==============================
Karena terdapat cross-sectional dependence, digunakan Driscoll-Kraay robust standard error.

t test of coefficients:

              Estimate Std. Error t value  Pr(>|t|)    
Penumpang_pax 11.27706    0.74315  15.175 < 2.2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Ekomet Lanjut : Regresi Data Panel

Mangunzonkkk

2026-03-25

1 Library

2 Input Data

2.1 Preview Data

3 Statistik Deskriptif

3.0.1 Per Bandara

3.0.2 Per Bulan

4 Visualisasi Data Eksplorasi

4.0.1 Tren Penumpang per Bulan

4.0.2 Tren Bagasi per Bulan

4.0.3 Scatter Plot : Bagasi v Penumpang

4.0.4 Heatmap bulan x bandara

5 Pemilihan Model Regresi Panel

6 Uji Asumsi (Homoskedastisitas, Non-Autokorelasi, dan Cross-Sectional Dependence)

7 Robust Method

8 Ringkasan

9 Pemisah

10 Pemisah

11 Pemisah

12 Pemisah

13 (Shortcut) Otomatis Semua Mua Nya

14 Uji Asumsi Custom

14.1 (Glesjer/White, DurWatson, CSD pakai pcd)