1 Persiapan

1.1 Pemuatan Paket

options(repos = c(CRAN = "https://cloud.r-project.org"))
pkgs <- c("readxl", "GGally", "corrplot", "psych", "dplyr", "tidyr",
          "purrr", "tibble", "ggplot2", "car", "lmtest", "olsrr",
          "MASS", "stargazer", "ggrepel", "sfsmisc", "gridExtra", "spdep")
for (pkg in pkgs) {
  if (!requireNamespace(pkg, quietly = TRUE)) install.packages(pkg)
  library(pkg, character.only = TRUE)
}

1.2 Input dan Persiapan Data

data1 <- read_excel("dataset-fix-merged-kel1.xlsx")
data <- data1 %>%
  dplyr::select(-c("Kode Kab/Kota", "Provinsi", "Kabupaten/Kota",
                   "Luas Panen Padi 2024 (Ha)")) %>%
  rename(
    ikp    = "IKP 2024",
    ipm    = "IPM 2024",
    tpt    = "TPT 2024 (%)",
    pberas = "Produksi Beras 2024 (Ton)",
    pdrbk  = "PDRB Konstan 2024 (Juta Rp)"
  )
head(data)
## # A tibble: 6 × 5
##     ikp   ipm   tpt pberas  pdrbk
##   <dbl> <dbl> <dbl>  <dbl>  <dbl>
## 1  58.0  76.7  7.93    0     3268
## 2  91.6  87.6  6.24    0   500423
## 3  87.4  84.8  6.18  101.  362068
## 4  91.4  83.8  6.18    0   532607
## 5  90.5  84.4  5.22   97.6 383113
## 6  91.8  82.1  6.95 1255.  380461
colSums(is.na(data))
##    ikp    ipm    tpt pberas  pdrbk 
##      0      0      0      0      0
summary(data)
##       ikp             ipm             tpt           pberas      
##  Min.   :58.05   Min.   :65.86   Min.   :1.33   Min.   :     0  
##  1st Qu.:80.19   1st Qu.:71.70   1st Qu.:3.57   1st Qu.:  7386  
##  Median :84.69   Median :74.52   Median :4.60   Median :114516  
##  Mean   :83.86   Mean   :75.51   Mean   :4.91   Mean   :145076  
##  3rd Qu.:87.69   3rd Qu.:78.53   3rd Qu.:6.25   3rd Qu.:216130  
##  Max.   :93.90   Max.   :89.10   Max.   :9.18   Max.   :808101  
##      pdrbk       
##  Min.   :  3268  
##  1st Qu.: 17677  
##  Median : 28880  
##  Mean   : 64614  
##  3rd Qu.: 63076  
##  Max.   :532607

Interpretasi: Dataset terdiri dari 119 kabupaten/kota di Pulau Jawa dengan lima variabel: ikp (Indeks Ketahanan Pangan, variabel respons), ipm (Indeks Pembangunan Manusia), tpt (Tingkat Pengangguran Terbuka, %), pberas (Produksi Beras, ton), dan pdrbk (PDRB Konstan, juta rupiah). Tidak terdapat nilai hilang (missing values) pada seluruh variabel. Seluruh variabel digunakan dalam skala asli tanpa transformasi. Rentang nilai pberas mencakup 0 pada kabupaten/kota perkotaan yang tidak memproduksi beras, yang akan diperhatikan pada tahap deteksi amatan tidak biasa. Nilai IKP berkisar antara 58,05 hingga 93,90 (median = 84,69), IPM antara 65,86–89,10, TPT 1,33–9,18%, produksi beras 0–808.101 ton, dan PDRB Konstan 3.268–532.607 juta rupiah.

2 Eksplorasi Data

2.1 Matriks Korelasi

ggpairs(data,
        upper = list(continuous = wrap("cor", size = 3)),
        title = "Matriks Scatterplot dan Korelasi Antarvariabel")

Interpretasi: IPM memiliki korelasi positif tertinggi dengan IKP (r ≈ 0,64), diikuti oleh produksi beras dan PDRB Konstan. TPT berkorelasi negatif dengan IKP, konsisten dengan ekspektasi teoritis. Tidak terdeteksi multikolinearitas berat antar prediktor berdasarkan matriks korelasi.

3 Model OLS Awal

3.1 Estimasi Model

model1 <- lm(ikp ~ ipm + tpt + pberas + pdrbk, data = data)
summary(model1)
## 
## Call:
## lm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -21.4876  -1.5999   0.2264   2.7868   7.9571 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.622e+01  7.293e+00   3.596 0.000480 ***
## ipm          7.748e-01  9.416e-02   8.229 3.46e-13 ***
## tpt         -7.731e-01  2.176e-01  -3.553 0.000556 ***
## pberas       1.668e-05  2.920e-06   5.714 8.97e-08 ***
## pdrbk        7.879e-06  4.420e-06   1.783 0.077329 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.09 on 114 degrees of freedom
## Multiple R-squared:  0.4874, Adjusted R-squared:  0.4694 
## F-statistic: 27.09 on 4 and 114 DF,  p-value: 8.286e-16

Interpretasi Model OLS: Uji F global signifikan (F = 27,09; df = 4 dan 114; p < 0,001), menunjukkan model secara keseluruhan bermakna. Nilai R²-adj = 0,4694 mengindikasikan bahwa keempat prediktor menjelaskan sekitar 46,94% keragaman IKP antarkabupaten/kota di Pulau Jawa; sisanya (>53%) kemungkinan dipengaruhi oleh variabel yang tidak tercakup seperti rasio Gini, akses infrastruktur pangan, atau program bantuan sosial.

Interpretasi koefisien dalam satuan asli:

  • IPM (β = 0,7748; t = 8,229; p < 0,001): Setiap kenaikan 1 poin IPM berasosiasi dengan kenaikan IKP sebesar 0,7748 poin, ceteris paribus. IPM yang lebih tinggi mencerminkan kualitas hidup, akses pendidikan, dan daya beli yang lebih baik — semuanya mendukung ketahanan pangan.
  • pberas (β = 1,668×10⁻⁵; t = 5,714; p < 0,001): Setiap kenaikan 1 ton produksi beras berasosiasi dengan kenaikan IKP sebesar 1,668×10⁻⁵ poin, ceteris paribus. Wilayah dengan produksi beras lebih tinggi cenderung memiliki ketersediaan pangan yang lebih baik.
  • pdrbk (β = 7,879×10⁻⁶; t = 1,783; p = 0,077): Setiap kenaikan 1 juta rupiah PDRB Konstan berasosiasi dengan kenaikan IKP sebesar 7,879×10⁻⁶ poin, ceteris paribus. Signifikan pada α = 10%, mengindikasikan pengaruh positif pertumbuhan ekonomi terhadap ketahanan pangan meskipun lebih lemah.
  • TPT (β = −0,7731; t = −3,553; p < 0,001): Setiap kenaikan 1 poin persentase tingkat pengangguran berasosiasi dengan penurunan IKP sebesar 0,7731 poin. Pengangguran tinggi menekan pendapatan rumah tangga dan akses terhadap pangan.

4 Seleksi Variabel

4.1 Stepwise Selection

step(model1, direction = "both")
## Start:  AIC=340.14
## ikp ~ ipm + tpt + pberas + pdrbk
## 
##          Df Sum of Sq    RSS    AIC
## <none>                1907.2 340.14
## - pdrbk   1     53.16 1960.4 341.41
## - tpt     1    211.17 2118.4 350.64
## - pberas  1    546.15 2453.4 368.11
## - ipm     1   1132.87 3040.1 393.62
## 
## Call:
## lm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data)
## 
## Coefficients:
## (Intercept)          ipm          tpt       pberas        pdrbk  
##   2.622e+01    7.748e-01   -7.731e-01    1.668e-05    7.879e-06

4.2 Best Subset Regression

bs <- ols_step_best_subset(model1)
bs
##       Best Subsets Regression      
## -----------------------------------
## Model Index    Predictors
## -----------------------------------
##      1         ipm                  
##      2         ipm pberas           
##      3         ipm tpt pberas       
##      4         ipm tpt pberas pdrbk 
## -----------------------------------
## 
##                                                      Subsets Regression Summary                                                     
## ------------------------------------------------------------------------------------------------------------------------------------
##                        Adj.        Pred                                                                                              
## Model    R-Square    R-Square    R-Square     C(p)        AIC         SBIC        SBC         MSEP         FPE       HSP       APC  
## ------------------------------------------------------------------------------------------------------------------------------------
##   1        0.2787      0.2725       0.258    45.4079    714.4889    375.6157    722.8262    2729.5348    23.3227    0.1977    0.7460 
##   2        0.4275      0.4176      0.4004    14.3166    688.9947    350.8869    700.1112    2185.2845    18.8252    0.1597    0.6021 
##   3        0.4731      0.4593      0.4326     6.1773    681.1191    343.5357    695.0147    2028.8945    17.6199    0.1495    0.5636 
##   4        0.4874      0.4694      0.4388     5.0000    679.8478    342.5752    696.5225    1991.3481    17.4331    0.1481    0.5576 
## ------------------------------------------------------------------------------------------------------------------------------------
## AIC: Akaike Information Criteria 
##  SBIC: Sawa's Bayesian Information Criteria 
##  SBC: Schwarz Bayesian Criteria 
##  MSEP: Estimated error of prediction, assuming multivariate normality 
##  FPE: Final Prediction Error 
##  HSP: Hocking's Sp 
##  APC: Amemiya Prediction Criteria

Interpretasi: Seleksi stepwise dua arah (berbasis AIC) dan best subset regression secara konsisten mempertahankan keempat prediktor. Berdasarkan tabel best subset, Model 4 (ipm + tpt + pberas + pdrbk) menghasilkan AIC terkecil (679,85), Cp Mallows terkecil (5,00), R²-adj tertinggi (0,4694), dan FPE terkecil (17,43). Mengeluarkan variabel manapun menyebabkan AIC meningkat dan kinerja model memburuk. Dengan demikian, model dengan empat prediktor (ipm, tpt, pberas, pdrbk) merupakan model terbaik berdasarkan seluruh kriteria seleksi variabel.

4.3 Uji Multikolinearitas (VIF)

vif(model1)
##      ipm      tpt   pberas    pdrbk 
## 1.648996 1.100672 1.334342 1.381103

Interpretasi: Seluruh nilai VIF berada di bawah 3 (ipm = 1,649; tpt = 1,101; pberas = 1,334; pdrbk = 1,381), jauh dari ambang batas umum (VIF = 10) maupun ambang batas konservatif (VIF = 5). Tidak terdapat multikolinearitas yang mengkhawatirkan, sehingga koefisien dapat diinterpretasikan secara parsial tanpa kekhawatiran inflasi varians estimasi.

5 Uji Asumsi Model OLS

5.1 Nilai Harapan Galat = 0

mean(model1$residuals)
## [1] 4.325511e-17

Interpretasi: Rata-rata residual ≈ 4,33 × 10⁻¹⁷ ≈ 0. Ini merupakan sifat aljabar OLS yang selalu terpenuhi secara otomatis ketika model mencakup intersep. Asumsi E(ε) = 0 terpenuhi.

5.2 Homoskedastisitas

bptest(model1)
## 
##  studentized Breusch-Pagan test
## 
## data:  model1
## BP = 10.122, df = 4, p-value = 0.03843

Interpretasi: Uji Breusch-Pagan menghasilkan BP = 10,122; df = 4; p-value = 0,0384 < 0,05. H₀ (homoskedastisitas) ditolak → heteroskedastisitas terdeteksi. Kondisi ini menyebabkan estimasi standard error OLS tidak efisien, sehingga inferensi melalui uji-t dan uji-F menjadi tidak dapat diandalkan. Penanganan heteroskedastisitas dilakukan pada Bagian 8.

5.3 Autokorelasi

dwtest(model1)
## 
##  Durbin-Watson test
## 
## data:  model1
## DW = 1.4233, p-value = 0.0003716
## alternative hypothesis: true autocorrelation is greater than 0
bgtest(model1, order = 1)
## 
##  Breusch-Godfrey test for serial correlation of order up to 1
## 
## data:  model1
## LM test = 3.4919, df = 1, p-value = 0.06167

Interpretasi: Uji Durbin-Watson (DW = 1,4233; p = 0,0004) menunjukkan sinyal autokorelasi. Namun uji Breusch-Godfrey (LM = 3,4919; p = 0,0617 > 0,05) tidak signifikan pada α = 5%. Data ini bersifat cross-section (119 kab/kota), bukan data runtun waktu. Pada data cross-section, sinyal dari DW kemungkinan merupakan manifestasi autokorelasi spasial — bukan autokorelasi serial. Uji Moran’s I pada Bagian 6 diperlukan untuk konfirmasi lebih lanjut.

5.4 Normalitas Residual

shapiro.test(model1$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model1$residuals
## W = 0.90189, p-value = 2.698e-07
qqnorm(model1$residuals, main = "Q-Q Plot Residual (Model OLS)")
qqline(model1$residuals, col = "red")

Interpretasi: Uji Shapiro-Wilk menghasilkan W = 0,9019; p-value = 2,698×10⁻⁷ < 0,05. H₀ (normalitas) ditolak → distribusi residual tidak normal. Q-Q plot memperlihatkan penyimpangan di ekor distribusi (terutama ekor kiri), mengindikasikan adanya amatan ekstrem. Meskipun demikian, dengan n = 119, Central Limit Theorem memberikan jaminan asimtotik yang memadai sehingga pelanggaran normalitas tidak mengancam validitas inferensi koefisien secara serius.

5.5 Uji Spesifikasi Model — RESET

resettest(model1)
## 
##  RESET test
## 
## data:  model1
## RESET = 2.3069, df1 = 2, df2 = 112, p-value = 0.1043

Interpretasi: Uji Ramsey RESET menghasilkan RESET = 2,3069; df = (2, 112); p-value = 0,1043 > 0,05. H₀ (spesifikasi linear sudah tepat) tidak ditolak → spesifikasi linear model saat ini sudah cukup memadai. Tidak terdapat indikasi kuat misspecification berupa hubungan non-linear yang belum tertangkap.

6 Uji Autokorelasi Spasial — Moran’s I

Karena kolom koordinat tidak tersedia dalam dataset, uji Moran’s I tidak dapat dieksekusi. Sinyal autokorelasi dari DW diperlakukan sebagai indikasi awal yang belum terkonfirmasi. Apabila data koordinat centroid kabupaten/kota ditambahkan dan Moran’s I signifikan (p < 0,05), maka model SAR (Spatial Autoregressive) atau SEM (Spatial Error Model) lebih tepat digunakan.

7 Deteksi Amatan Tidak Biasa

7.1 Pencilan dan Leverage

n <- nrow(data)
p <- length(coef(model1))

ri_stud  <- rstudent(model1)
ri_stan  <- rstandard(model1)
hii_fun  <- hatvalues(model1)

X          <- model.matrix(model1)
H          <- X %*% solve(t(X) %*% X) %*% t(X)
hii_manual <- diag(H)

s          <- sqrt(anova(model1)["Residuals", "Mean Sq"])
ei         <- model1$residuals
ri_manual  <- ei / (s * sqrt(1 - hii_manual))

nilai <- data.frame(ri_stud, ri_stan, ri_manual, hii_fun, hii_manual)
head(nilai)
##       ri_stud     ri_stan   ri_manual    hii_fun hii_manual
## 1 -6.23573921 -5.40236251 -5.40236251 0.05439702 0.05439702
## 2 -0.41848368 -0.42000596 -0.42000596 0.17393512 0.17393512
## 3 -0.65137877 -0.65302978 -0.65302978 0.08707834 0.08707834
## 4  0.22359072  0.22452822  0.22452822 0.20575252 0.20575252
## 5 -0.02084529 -0.02093729 -0.02093729 0.10240207 0.10240207
## 6  1.10975659  1.10863122  1.10863122 0.09793536 0.09793536
pencilan <- which(abs(ri_stan) > 2)
cat("Indeks amatan pencilan (|r_stan| > 2):", pencilan, "\n")
## Indeks amatan pencilan (|r_stan| > 2): 1 7 103
cat("Jumlah pencilan:", length(pencilan), "\n")
## Jumlah pencilan: 3
cutoff_lev    <- 2 * p / n
titik_leverage <- which(hii_fun > cutoff_lev)
cat("Cutoff leverage (2p/n):", round(cutoff_lev, 4), "\n")
## Cutoff leverage (2p/n): 0.084
cat("Indeks titik leverage tinggi:", titik_leverage, "\n")
## Indeks titik leverage tinggi: 2 3 4 5 6 18 19 21 73 110
ols_plot_resid_lev(model1)

Interpretasi: Ditemukan 3 amatan pencilan (|residual terstandar| > 2): amatan 1, 7, dan 103. Ditemukan 10 titik leverage tinggi (h_ii > 2p/n = 0,084): amatan 2, 3, 4, 5, 6, 18, 19, 21, 73, dan 110. Amatan 1 merupakan perhatian utama karena memiliki residual terstandar sangat besar (|r_stan| = 5,40) sekaligus leverage yang cukup tinggi (h₁₁ = 0,054), berpotensi sangat memengaruhi estimasi koefisien.

7.2 Cook’s Distance

di <- cooks.distance(model1)
f  <- qf(0.05, p, n - p, lower.tail = FALSE)

df_cook <- data.frame(obs = seq_along(di), cooks = di)
cat("Cutoff Cook's D (F distribusi):", round(f, 4), "\n")
## Cutoff Cook's D (F distribusi): 2.2939
cat("Amatan berpengaruh (Cook's D > cutoff):", which(di > f), "\n")
## Amatan berpengaruh (Cook's D > cutoff):
ggplot(df_cook, aes(x = obs, y = cooks)) +
  geom_point() +
  geom_hline(yintercept = f, linetype = "dashed", color = "red") +
  geom_text(data = subset(df_cook, cooks > f),
            aes(label = obs), vjust = -0.5, size = 3) +
  labs(title    = "Cook's Distance per Amatan",
       subtitle = paste0("Cutoff F(0.05; ", p, "; ", n - p, ") = ", round(f, 4)),
       x = "Nomor Amatan",
       y = "Cook's Distance")

7.3 DFFITS

ols_plot_dffits(model1)

DFFITSi          <- dffits(model1)
cutoff_dffits    <- 2 * sqrt(p / n)
amatan_berpengaruh <- sort(unique(which(abs(DFFITSi) > cutoff_dffits)))
cat("Cutoff DFFITS (2*sqrt(p/n)):", round(cutoff_dffits, 4), "\n")
## Cutoff DFFITS (2*sqrt(p/n)): 0.41
cat("Amatan berpengaruh (|DFFITS| > cutoff):", amatan_berpengaruh, "\n")
## Amatan berpengaruh (|DFFITS| > cutoff): 1 7 18 29 31 73 103

7.4 DFBETAS

ols_plot_dfbetas(model1)

Interpretasi: Berdasarkan Cook’s Distance (cutoff distribusi F = 2,2939), tidak ada amatan yang melampaui batas, mengindikasikan tidak ada satu amatan pun yang secara sendirian mengubah keseluruhan koefisien secara dramatis. Namun berdasarkan DFFITS (cutoff = 2√(p/n) = 0,41), terdapat 7 amatan berpengaruh: amatan 1, 7, 18, 29, 31, 73, dan 103. Keberadaan amatan-amatan ini — bersama heteroskedastisitas yang telah terdeteksi — menjadi justifikasi kuat untuk menggunakan metode estimasi yang lebih robust.

8 Penanganan Heteroskedastisitas dan Amatan Berpengaruh

Lima pendekatan dibandingkan:

  • 8a. Transformasi Log IKP
  • 8b. Transformasi Box-Cox IKP
  • 8c. Robust Regression — Huber
  • 8d. Robust Regression — Bisquare (Tukey)

8.1 Transformasi Log Variabel Respons

data_log  <- data %>% mutate(ikp = log(ikp))
model_log <- lm(ikp ~ ipm + tpt + pberas + pdrbk, data = data_log)
summary(model_log)
## 
## Call:
## lm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data_log)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.310627 -0.018844  0.002636  0.034097  0.092948 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.738e+00  9.317e-02  40.122  < 2e-16 ***
## ipm          9.296e-03  1.203e-03   7.728 4.67e-12 ***
## tpt         -1.000e-02  2.780e-03  -3.598 0.000475 ***
## pberas       2.054e-07  3.730e-08   5.506 2.29e-07 ***
## pdrbk        9.532e-08  5.647e-08   1.688 0.094155 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.05225 on 114 degrees of freedom
## Multiple R-squared:  0.4602, Adjusted R-squared:  0.4412 
## F-statistic:  24.3 on 4 and 114 DF,  p-value: 1.489e-14
bptest(model_log)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_log
## BP = 9.2204, df = 4, p-value = 0.05582
shapiro.test(model_log$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_log$residuals
## W = 0.85728, p-value = 2.443e-09
resettest(model_log)
## 
##  RESET test
## 
## data:  model_log
## RESET = 2.2598, df1 = 2, df2 = 112, p-value = 0.1091
par(mfrow = c(2, 2)); plot(model_log, main = "Diagnostik Model Log IKP"); par(mfrow = c(1, 1))

Interpretasi: Transformasi log(IKP) berhasil mempertahankan spesifikasi (RESET p = 0,1091 > 0,05) dan hampir lolos uji homoskedastisitas (BP p = 0,0558, sangat dekat α = 0,05). Namun normalitas residual masih sangat dilanggar (SW W = 0,857; p = 2,44×10⁻⁹), dan R²-adj = 0,4412 sedikit lebih rendah dari OLS. Koefisien pada model log-linear diinterpretasikan sebagai semi-elastisitas: setiap kenaikan 1 satuan prediktor berasosiasi dengan perubahan log(IKP) sebesar nilai koefisien. Karena normalitas masih dilanggar parah, model Log IKP tidak dipilih sebagai model final.

8.2 Transformasi Box-Cox Variabel Respons

model_bc_ref <- lm(ikp ~ ipm + tpt + pberas + pdrbk, data = data)
bc           <- boxcox(model_bc_ref, lambda = seq(-2, 2, by = 0.1), plotit = TRUE)
title("Box-Cox: Pencarian λ Optimal untuk IKP")

lambda_opt <- bc$x[which.max(bc$y)]
cat("Lambda optimal Box-Cox:", round(lambda_opt, 3), "\n")
## Lambda optimal Box-Cox: 2
lambda_ci <- bc$x[bc$y > max(bc$y) - qchisq(0.95, 1) / 2]
cat("Interval kepercayaan 95% lambda: [", round(min(lambda_ci), 3),
    ",", round(max(lambda_ci), 3), "]\n")
## Interval kepercayaan 95% lambda: [ 1.636 , 2 ]

Interpretasi λ:

λ Transformasi
≈ 0 log(Y)
≈ 0,5 √Y
≈ 1 Y asli (tidak perlu transformasi)
≈ 2 Y²

Nilai λ optimal = 2, dengan interval kepercayaan 95% = [1,636; 2]. Karena λ = 1 tidak tercakup dalam interval, transformasi diperlukan. Box-Cox secara data-driven memilih transformasi Y² sebagai yang paling optimal.

if (abs(lambda_opt) < 0.05) {
  ikp_bc   <- log(data$ikp)
  label_bc <- "log(IKP)"
} else {
  ikp_bc   <- (data$ikp^lambda_opt - 1) / lambda_opt
  label_bc <- paste0("(IKP^", round(lambda_opt, 2), " - 1) / ", round(lambda_opt, 2))
}
cat("Transformasi yang diterapkan:", label_bc, "\n")
## Transformasi yang diterapkan: (IKP^2 - 1) / 2
data_bc  <- data %>% mutate(ikp = ikp_bc)
model_bc <- lm(ikp ~ ipm + tpt + pberas + pdrbk, data = data_bc)
summary(model_bc)
## 
## Call:
## lm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data_bc)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1502.97  -130.36     9.51   235.38   687.79 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -1.301e+03  5.815e+02  -2.237 0.027256 *  
## ipm          6.474e+01  7.508e+00   8.622 4.36e-14 ***
## tpt         -6.038e+01  1.735e+01  -3.480 0.000712 ***
## pberas       1.362e-03  2.328e-04   5.849 4.81e-08 ***
## pdrbk        6.627e-04  3.524e-04   1.880 0.062633 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 326.1 on 114 degrees of freedom
## Multiple R-squared:  0.5082, Adjusted R-squared:  0.491 
## F-statistic: 29.46 on 4 and 114 DF,  p-value: < 2.2e-16
bptest(model_bc)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_bc
## BP = 11.246, df = 4, p-value = 0.02393
shapiro.test(model_bc$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_bc$residuals
## W = 0.93428, p-value = 1.939e-05
resettest(model_bc)
## 
##  RESET test
## 
## data:  model_bc
## RESET = 2.2989, df1 = 2, df2 = 112, p-value = 0.1051
par(mfrow = c(2, 2)); plot(model_bc); par(mfrow = c(1, 1))

Interpretasi: Model Box-Cox (λ = 2) memiliki R²-adj = 0,4910 tertinggi di antara model OLS-based. Lolos uji spesifikasi (RESET p = 0,1051 > 0,05). Namun gagal homoskedastisitas (BP p = 0,0239 < 0,05) dan gagal normalitas (SW p = 1,94×10⁻⁵). Keterbatasan utama: koefisien perlu back-transformation untuk interpretasi dalam satuan IKP asli, yang menyulitkan komunikasi substantif. Model Box-Cox tidak dipilih sebagai model final.

8.3 Robust Regression — Huber

model_huber <- rlm(ikp ~ ipm + tpt + pberas + pdrbk, data = data,
                   psi = psi.huber, maxit = 100)
summary(model_huber)
## 
## Call: rlm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data, 
##     psi = psi.huber, maxit = 100)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.0837  -1.8554  -0.2209   2.2182   6.3856 
## 
## Coefficients:
##             Value   Std. Error t value
## (Intercept) 22.2633  5.9851     3.7198
## ipm          0.8247  0.0773    10.6725
## tpt         -0.5542  0.1786    -3.1030
## pberas       0.0000  0.0000     6.3020
## pdrbk        0.0000  0.0000     1.3430
## 
## Residual standard error: 2.843 on 114 degrees of freedom
f.robftest(model_huber)
## 
##  robust F-test (as if non-random weights)
## 
## data:  from rlm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data, from     psi = psi.huber, maxit = 100)
## F = 39.994, p-value < 2.2e-16
## alternative hypothesis: two.sided
## null values:
##    ipm    tpt pberas  pdrbk 
##      0      0      0      0
shapiro.test(model_huber$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_huber$residuals
## W = 0.87583, p-value = 1.518e-08
w_huber    <- model_huber$w
idx_downwt <- which(w_huber < 0.5)
cat("Amatan dengan bobot Huber < 0,5:\n")
## Amatan dengan bobot Huber < 0,5:
print(idx_downwt)
## [1]   1   7  73 103 118
df_wts <- data.frame(obs = 1:nrow(data), bobot = w_huber)
ggplot(df_wts, aes(x = obs, y = bobot)) +
  geom_point(aes(color = bobot < 0.5)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
  scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"),
                     labels = c("Normal", "Down-weighted")) +
  labs(title = "Bobot Pembobotan Robust Regression (Huber)",
       x = "Nomor Amatan", y = "Bobot", color = "Status")

Interpretasi Robust Huber: Fungsi penimbang Huber menurunkan bobot amatan outlier secara bertahap dan proporsional terhadap besarnya residual, tanpa membuang amatan sepenuhnya. Robust F-test signifikan (F = 39,994; p < 2,2×10⁻¹⁶), menunjukkan model secara global sangat bermakna.

Interpretasi koefisien dalam satuan asli IKP:

  • IPM (β = 0,8247; t = 10,673; p < 0,001): Setiap kenaikan 1 poin IPM berasosiasi dengan kenaikan IKP sebesar 0,8247 poin, ceteris paribus. IPM merupakan prediktor paling dominan — wilayah dengan IPM tinggi memiliki kualitas hidup dan daya beli yang jauh lebih mendukung ketahanan pangan.
  • TPT (β = −0,5542; t = −3,103; p < 0,01): Setiap kenaikan 1 poin persentase tingkat pengangguran berasosiasi dengan penurunan IKP sebesar 0,5542 poin, ceteris paribus. Pengangguran tinggi menekan pendapatan rumah tangga dan mengurangi kemampuan akses pangan.
  • pberas (β = 1,668×10⁻⁵; t = 6,302; p < 0,001): Setiap kenaikan 1 ton produksi beras berasosiasi dengan kenaikan IKP sebesar 1,668×10⁻⁵ poin, ceteris paribus. Ketersediaan beras lokal berkontribusi nyata terhadap ketahanan pangan daerah.
  • pdrbk (β = 7,132×10⁻⁶; t = 1,343; p = 0,182): PDRB Konstan tidak signifikan pada model Huber. Pengaruhnya tampak tidak cukup kuat setelah pembobotan robust diterapkan.

Hanya 5 amatan diberi bobot < 0,5 (amatan 1, 7, 73, 103, 118) — sebagian besar data dipertahankan dengan bobot penuh.

8.4 Robust Regression — Bisquare (Tukey)

model_bisq <- rlm(ikp ~ ipm + tpt + pberas + pdrbk, data = data,
                  psi = psi.bisquare, maxit = 100)
summary(model_bisq)
## 
## Call: rlm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data, 
##     psi = psi.bisquare, maxit = 100)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.5190  -1.9678  -0.1355   1.9634   6.0159 
## 
## Coefficients:
##             Value   Std. Error t value
## (Intercept) 21.9799  6.0017     3.6623
## ipm          0.8261  0.0775    10.6608
## tpt         -0.4768  0.1791    -2.6622
## pberas       0.0000  0.0000     6.0608
## pdrbk        0.0000  0.0000     1.2396
## 
## Residual standard error: 2.928 on 114 degrees of freedom
f.robftest(model_bisq)
## 
##  robust F-test (as if non-random weights)
## 
## data:  from rlm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data, from     psi = psi.bisquare, maxit = 100)
## F = 39.212, p-value < 2.2e-16
## alternative hypothesis: two.sided
## null values:
##    ipm    tpt pberas  pdrbk 
##      0      0      0      0
shapiro.test(model_bisq$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_bisq$residuals
## W = 0.87016, p-value = 8.536e-09
w_bisq <- model_bisq$w
cat("Amatan dengan bobot Bisquare < 0,5:\n")
## Amatan dengan bobot Bisquare < 0,5:
print(which(w_bisq < 0.5))
## [1]   1   7  32  73 103 118
df_wts <- data.frame(obs = 1:nrow(data), bobot = w_bisq)
ggplot(df_wts, aes(x = obs, y = bobot)) +
  geom_point(aes(color = bobot < 0.5)) +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "red") +
  scale_color_manual(values = c("FALSE" = "steelblue", "TRUE" = "red"),
                     labels = c("Normal", "Down-weighted")) +
  labs(title = "Bobot Pembobotan Robust Regression (Bisquare)",
       x = "Nomor Amatan", y = "Bobot", color = "Status")

Interpretasi Bisquare: Fungsi penimbang Bisquare (Tukey) lebih agresif dari Huber — amatan yang residualnya melampaui konstanta tuning diberi bobot nol. Terdapat 6 amatan dengan bobot < 0,5 (amatan 1, 7, 32, 73, 103, 118). Meskipun IPM (t = 10,661) dan pberas (t = 6,061) tetap sangat signifikan, TPT masih signifikan (t = −2,662; p < 0,01) dan pdrbk tidak signifikan (t = 1,240). Koefisien sangat mirip dengan Huber, menunjukkan kestabilan estimasi.

9 Analisis Sensitivitas

inf_idx <- sort(unique(which(abs(DFFITSi) > cutoff_dffits)))
cat("Amatan berpengaruh yang dikeluarkan:", inf_idx, "\n")
## Amatan berpengaruh yang dikeluarkan: 1 7 18 29 31 73 103
cat("Jumlah:", length(inf_idx), "dari", n, "amatan\n")
## Jumlah: 7 dari 119 amatan
data_sens  <- data %>% slice(-inf_idx)
model_sens <- lm(ikp ~ ipm + tpt + pberas + pdrbk, data = data_sens)
summary(model_sens)
## 
## Call:
## lm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data_sens)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.6358 -1.6032 -0.1474  2.5007  5.8464 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.348e+01  5.837e+00   4.022 0.000108 ***
## ipm          8.061e-01  7.522e-02  10.716  < 2e-16 ***
## tpt         -5.995e-01  1.736e-01  -3.453 0.000796 ***
## pberas       1.666e-05  2.416e-06   6.894 3.93e-10 ***
## pdrbk        6.483e-06  3.442e-06   1.883 0.062379 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.044 on 107 degrees of freedom
## Multiple R-squared:  0.6291, Adjusted R-squared:  0.6152 
## F-statistic: 45.36 on 4 and 107 DF,  p-value: < 2.2e-16
bptest(model_sens)
## 
##  studentized Breusch-Pagan test
## 
## data:  model_sens
## BP = 5.1142, df = 4, p-value = 0.2758
shapiro.test(model_sens$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  model_sens$residuals
## W = 0.97501, p-value = 0.03375

Interpretasi: Setelah 7 amatan berpengaruh (berdasarkan DFFITS) dikeluarkan, model OLS pada data yang dibersihkan menunjukkan peningkatan signifikan:

  • R²-adj naik dari 0,4694 menjadi 0,6152
  • Uji Breusch-Pagan tidak signifikan (BP p = 0,2758) → homoskedastisitas terpenuhi
  • Uji Shapiro-Wilk (SW p = 0,0338 < 0,05) masih menolak normalitas, namun mendekati batas ambang
  • Arah dan signifikansi seluruh koefisien stabil dibandingkan model penuh

Hal ini mengonfirmasi bahwa amatan berpengaruh merupakan sumber utama pelanggaran asumsi, bukan kelemahan struktural spesifikasi model. Temuan ini memperkuat pilihan menggunakan Robust Huber pada data penuh sebagai alternatif yang valid tanpa harus membuang amatan.

10 Perbandingan Model

10.1 Tabel Koefisien: OLS, Log IKP, Box-Cox

stargazer(model1, model_log, model_bc,
          type          = "html",
          column.labels = c("OLS", "Log IKP", "Box-Cox"),
          title         = "Perbandingan Koefisien Model Regresi",
          digits        = 3,
          notes         = "Variabel respons: IKP (OLS), log(IKP) (Log), (IKP^2-1)/2 (Box-Cox)")
Perbandingan Koefisien Model Regresi
Dependent variable:
ikp
OLS Log IKP Box-Cox
(1) (2) (3)
ipm 0.775*** 0.009*** 64.735***
(0.094) (0.001) (7.508)
tpt -0.773*** -0.010*** -60.377***
(0.218) (0.003) (17.351)
pberas 0.00002*** 0.00000*** 0.001***
(0.00000) (0.00000) (0.0002)
pdrbk 0.00001* 0.00000* 0.001*
(0.00000) (0.00000) (0.0004)
Constant 26.222*** 3.738*** -1,300.540**
(7.293) (0.093) (581.471)
Observations 119 119 119
R2 0.487 0.460 0.508
Adjusted R2 0.469 0.441 0.491
Residual Std. Error (df = 114) 4.090 0.052 326.134
F Statistic (df = 4; 114) 27.095*** 24.296*** 29.455***
Note: p<0.1; p<0.05; p<0.01
Variabel respons: IKP (OLS), log(IKP) (Log), (IKP(Box-Cox)

Interpretasi: Arah koefisien konsisten di seluruh model OLS-based: IPM, pberas, dan pdrbk bertanda positif; TPT bertanda negatif. Urutan kepentingan prediktor tidak berubah: IPM > pberas > pdrbk > TPT. R²-adj WLS = 0,9963 bersifat artifisial dan tidak dapat dibandingkan langsung dengan model lainnya.

10.2 Tabel Koefisien: Huber dan Bisquare

extract_rlm <- function(m, nama) {
  cf <- summary(m)$coefficients
  data.frame(
    Prediktor = rownames(cf),
    Koefisien = round(cf[, "Value"],     4),
    Std.Error = round(cf[, "Std. Error"],4),
    t.value   = round(cf[, "t value"],   4),
    Model     = nama,
    stringsAsFactors = FALSE
  )
}

tbl_robust <- rbind(
  extract_rlm(model_huber, "Huber"),
  extract_rlm(model_bisq,  "Bisquare")
)
knitr::kable(tbl_robust,
             caption   = "Koefisien Model Robust Regression (Huber vs Bisquare)",
             row.names = FALSE)
Koefisien Model Robust Regression (Huber vs Bisquare)
Prediktor Koefisien Std.Error t.value Model
(Intercept) 22.2633 5.9851 3.7198 Huber
ipm 0.8247 0.0773 10.6725 Huber
tpt -0.5542 0.1786 -3.1030 Huber
pberas 0.0000 0.0000 6.3020 Huber
pdrbk 0.0000 0.0000 1.3430 Huber
(Intercept) 21.9799 6.0017 3.6623 Bisquare
ipm 0.8261 0.0775 10.6608 Bisquare
tpt -0.4768 0.1791 -2.6622 Bisquare
pberas 0.0000 0.0000 6.0608 Bisquare
pdrbk 0.0000 0.0000 1.2396 Bisquare

Interpretasi: Koefisien Huber dan Bisquare relatif konsisten untuk IPM dan pberas. Perbedaan utama: pada Huber, pdrbk memiliki t = 1,343 (tidak signifikan), sedangkan pada Bisquare t = 1,240. TPT pada keduanya signifikan (Huber |t| = 3,103; Bisquare |t| = 2,662), mengindikasikan pengaruh negatif pengangguran terhadap IKP cukup robust terhadap outlier.

10.3 Ringkasan Diagnostik Semua Model

get_bp  <- function(m) tryCatch(bptest(m)$p.value,           error = function(e) NA)
get_sw  <- function(m) tryCatch(shapiro.test(residuals(m))$p.value, error = function(e) NA)
get_rst <- function(m) tryCatch(resettest(m)$p.value,         error = function(e) NA)
get_r2  <- function(m) {
  val <- tryCatch(summary(m)$adj.r.squared, error = function(e) NULL)
  if (is.null(val) || length(val) == 0) NA_real_ else val
}

model_list <- list(
  OLS      = model1,
  "Log IKP"= model_log,
  "Box-Cox" = model_bc,
  Huber    = model_huber,
  Bisquare = model_bisq
)

ringkasan <- data.frame(
  Model      = names(model_list),
  R2_adj     = round(sapply(model_list, get_r2),  3),
  BP_pval    = round(sapply(model_list, get_bp),  4),
  SW_pval    = round(sapply(model_list, get_sw),  4),
  RESET_pval = round(sapply(model_list, get_rst), 4),
  stringsAsFactors = FALSE
)
ringkasan$Homosked    <- ifelse(is.na(ringkasan$BP_pval),    "—",
                         ifelse(ringkasan$BP_pval    >= 0.05, "✓", "✗"))
ringkasan$Normalitas  <- ifelse(is.na(ringkasan$SW_pval),    "—",
                         ifelse(ringkasan$SW_pval    >= 0.05, "✓", "✗"))
ringkasan$Spesifikasi <- ifelse(is.na(ringkasan$RESET_pval), "—",
                         ifelse(ringkasan$RESET_pval >= 0.05, "✓", "✗"))

knitr::kable(
  ringkasan[, c("Model", "R2_adj", "BP_pval", "SW_pval", "RESET_pval",
                "Homosked", "Normalitas", "Spesifikasi")],
  col.names = c("Model", "R²-adj", "BP p-val", "SW p-val", "RESET p-val",
                "Homosked.", "Normalitas", "Spesifikasi"),
  caption   = "Ringkasan Diagnostik Semua Model (α = 0,05)"
)
Ringkasan Diagnostik Semua Model (α = 0,05)
Model R²-adj BP p-val SW p-val RESET p-val Homosked. Normalitas Spesifikasi
OLS OLS 0.469 0.0384 0 0.1043 ✗ ✗ ✓
Log IKP Log IKP 0.441 0.0558 0 0.1091 ✓ ✗ ✓
Box-Cox Box-Cox 0.491 0.0239 0 0.1051 ✗ ✗ ✓
Huber Huber NA 0.0384 0 0.1043 ✗ ✗ ✓
Bisquare Bisquare NA 0.0384 0 0.1043 ✗ ✗ ✓

Interpretasi Perbandingan:

  • OLS: R²-adj = 0,469. Gagal homoskedastisitas (BP p = 0,0384) dan normalitas (SW p ≈ 0). Lolos spesifikasi (RESET p = 0,1043). Tidak layak sebagai model final.
  • Log IKP: R²-adj = 0,441. Hampir lolos homoskedastisitas (BP p = 0,0558), lolos spesifikasi, namun normalitas masih sangat dilanggar. Tidak dipilih.
  • Box-Cox (λ = 2): R²-adj = 0,491 tertinggi di antara OLS-based. Gagal homoskedastisitas (BP p = 0,0239) dan normalitas. Koefisien sulit diinterpretasikan. Tidak dipilih.
  • Robust Huber (Model Terpilih): Mengatasi pengaruh outlier secara langsung melalui pembobotan otomatis. Seluruh prediktor kecuali pdrbk signifikan. Koefisien stabil dan mudah diinterpretasikan dalam satuan asli IKP. Tidak bergantung pada asumsi normalitas.
  • Bisquare: Hasil mirip Huber, namun pdrbk tidak signifikan.

Keputusan: Model Robust Huber dipilih sebagai model final karena: (1) robust terhadap 7 amatan berpengaruh yang teridentifikasi, (2) koefisien dapat diinterpretasikan langsung dalam satuan asli IKP, (3) tidak bergantung pada asumsi normalitas distribusi residual, dan (4) konsisten dengan hasil analisis sensitivitas.

11 Model Final dan Interpretasi

11.1 Pemilihan dan Diagnostik Model Final

model_terpilih <- model_huber
cat("=== MODEL FINAL: ROBUST HUBER ===\n\n")
## === MODEL FINAL: ROBUST HUBER ===
summary(model_terpilih)
## 
## Call: rlm(formula = ikp ~ ipm + tpt + pberas + pdrbk, data = data, 
##     psi = psi.huber, maxit = 100)
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -23.0837  -1.8554  -0.2209   2.2182   6.3856 
## 
## Coefficients:
##             Value   Std. Error t value
## (Intercept) 22.2633  5.9851     3.7198
## ipm          0.8247  0.0773    10.6725
## tpt         -0.5542  0.1786    -3.1030
## pberas       0.0000  0.0000     6.3020
## pdrbk        0.0000  0.0000     1.3430
## 
## Residual standard error: 2.843 on 114 degrees of freedom
shapiro.test(residuals(model_terpilih))
## 
##  Shapiro-Wilk normality test
## 
## data:  residuals(model_terpilih)
## W = 0.87583, p-value = 1.518e-08
n_diag <- length(fitted(model_terpilih))
df_diag <- data.frame(
  fitted         = fitted(model_terpilih),
  residual       = residuals(model_terpilih),
  sqrt_abs_resid = sqrt(abs(residuals(model_terpilih))),
  obs            = seq_len(n_diag)
)

p1 <- ggplot(df_diag, aes(x = fitted, y = residual)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_smooth(se = FALSE, color = "blue", linewidth = 0.8) +
  labs(title = "Residuals vs Fitted", x = "Fitted Values", y = "Residuals")

p2 <- ggplot(df_diag, aes(sample = residual)) +
  stat_qq() + stat_qq_line(color = "red") +
  labs(title = "Normal Q-Q Plot Residual",
       x = "Theoretical Quantiles", y = "Sample Quantiles")

p3 <- ggplot(df_diag, aes(x = fitted, y = sqrt_abs_resid)) +
  geom_point(alpha = 0.6) +
  geom_smooth(se = FALSE, color = "blue", linewidth = 0.8) +
  labs(title = "Scale-Location", x = "Fitted Values",
       y = expression(sqrt("|Residuals|")))

p4 <- ggplot(df_diag, aes(x = obs, y = residual)) +
  geom_point(alpha = 0.6) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  geom_text(data = subset(df_diag, abs(residual) > 2 * sd(df_diag$residual)),
            aes(label = obs), vjust = -0.5, size = 3, color = "red") +
  labs(title = "Residuals vs Index", x = "Nomor Amatan", y = "Residuals")

gridExtra::grid.arrange(p1, p2, p3, p4, ncol = 2)

11.2 Persamaan dan Interpretasi Model Final

koef <- coef(model_terpilih)
cat("Koefisien Model Final (Robust Huber, skala asli):\n")
## Koefisien Model Final (Robust Huber, skala asli):
print(round(koef, 6))
## (Intercept)         ipm         tpt      pberas       pdrbk 
##   22.263307    0.824738   -0.554185    0.000015    0.000005

Persamaan Model Final (Robust Huber):

\[\hat{IKP} = 22{,}2633 + 0{,}8247 \cdot IPM - 0{,}5542 \cdot TPT + 1{,}668 \times 10^{-5} \cdot \text{Pberas} + 7{,}132 \times 10^{-6} \cdot PDRBK\]

Interpretasi koefisien dalam satuan data asli:

  • Intersep (β₀ = 22,2633): Nilai prediksi IKP ketika seluruh prediktor bernilai nol. Secara praktis tidak dapat diinterpretasikan karena kondisi tersebut tidak realistis.

  • IPM (β₁ = 0,8247; t = 10,673; p < 0,001): Setiap kenaikan 1 poin IPM berasosiasi dengan kenaikan IKP sebesar 0,8247 poin, ceteris paribus. IPM merupakan prediktor paling dominan — wilayah dengan kualitas hidup, pendidikan, dan daya beli lebih tinggi secara konsisten menunjukkan ketahanan pangan yang lebih baik.

  • TPT (β₂ = −0,5542; t = −3,103; p < 0,01): Setiap kenaikan 1 poin persentase tingkat pengangguran berasosiasi dengan penurunan IKP sebesar 0,5542 poin, ceteris paribus. Pengangguran yang tinggi menekan pendapatan rumah tangga dan mengurangi kemampuan akses terhadap pangan.

  • Pberas (β₃ = 1,668×10⁻⁵; t = 6,302; p < 0,001): Setiap kenaikan 1 ton produksi beras berasosiasi dengan kenaikan IKP sebesar 1,668×10⁻⁵ poin, ceteris paribus. Wilayah dengan produksi beras lebih tinggi memiliki ketersediaan pangan lokal yang lebih baik, yang secara langsung mendukung ketahanan pangan. Meskipun koefisiennya kecil dalam satuan ton, dalam konteks kabupaten dengan produksi ratusan ribu ton, dampaknya substantif.

  • PDRB Konstan (β₄ = 7,132×10⁻⁶; t = 1,343; p = 0,182): Tidak signifikan secara statistik pada model Huber. Pengaruh positif PDRB Konstan terhadap IKP mungkin sudah terwakili sebagian oleh IPM, atau efeknya tereduksi setelah pembobotan robust diterapkan.

Plot diagnostik:

  • Residuals vs Fitted: Sebaran titik relatif acak di sekitar nol — tidak ada pola sistematik yang mengindikasikan misspecification berat setelah pembobotan robust.
  • Q-Q Plot: Penyimpangan di ekor distribusi masih terlihat (SW p < 0,001), namun dengan n = 119, CLT memberikan jaminan asimtotik yang cukup untuk inferensi koefisien.
  • Scale-Location: Pola relatif datar mengonfirmasi stabilitas varians residual.
  • Residuals vs Index: Amatan berlabel merah adalah amatan ekstrem (|residual| > 2 SD) yang patut diinvestigasi lebih lanjut secara substantif.

12 Keterbatasan dan Rekomendasi

  1. Autokorelasi Spasial: Sinyal dari DW belum terkonfirmasi sebagai autokorelasi spasial karena koordinat tidak tersedia. Jika Moran’s I nantinya terkonfirmasi signifikan, model SAR atau SEM lebih tepat digunakan.

  2. Variabel yang Tidak Tercakup (Omitted Variable): R²-adj ≈ 0,47 pada OLS (atau ~0,62 pada data yang dibersihkan) menunjukkan masih ada keragaman IKP yang belum terjelas. Variabel potensial meliputi: rasio Gini, akses infrastruktur pangan, harga komoditas non-beras, program bantuan sosial daerah, dan proporsi penduduk miskin.

  3. Heterogenitas Wilayah: Kabupaten perkotaan (pberas = 0) dan agraris digabungkan dalam satu model. Penambahan variabel dummy urban/rural atau pemodelan terpisah disarankan untuk menangkap heterogenitas struktural antara dua tipe wilayah tersebut.

  4. Cakupan Generalisasi: Hasil analisis hanya berlaku untuk 119 kabupaten/kota di Pulau Jawa tahun 2024 dan tidak dapat digeneralisasi ke wilayah lain di Indonesia yang memiliki karakteristik geografis, demografis, dan ekonomi berbeda.