1 Latar Belakang

Analisis ini bertujuan memilih model terbaik untuk menjelaskan Indeks Ketahanan Pangan (IKP) kabupaten/kota di Kalimantan tahun 2021. Pemilihan variabel dilakukan berdasarkan kerangka ketahanan pangan FAO (2006) yang mencakup empat pilar utama: Availability, Access, Utilization, dan Stability.

Variabel Pilar FAO
IPM Access + Utilization
ProdPangan Availability
PDRB Access
Internet Access + Stability
PertPend Stability

AHH dan RLS tidak digunakan karena merupakan komponen pembentuk IPM dan menyebabkan multikolinearitas tinggi pada analisis awal.

2 Library

library(readxl); library(tidyverse); library(leaps); library(MASS)
library(glmnet); library(car); library(lmtest); library(caret)
library(knitr); library(kableExtra); library(ggrepel)
library(corrplot); library(lm.beta)

select  <- dplyr::select; filter  <- dplyr::filter
mutate  <- dplyr::mutate; arrange <- dplyr::arrange
lag     <- dplyr::lag

set.seed(123)

3 Import Data

raw <- read_excel("Data Semifinal NSC 2023 (1).xlsx", sheet = "Sheet2")
colnames(raw) <- c(
  "No","Kabupaten","AHH_2020","AHH_2021","AHH_2022",
  "IKP_2019","IKP_2020","IKP_2021","IKP_2022",
  "IPM_2020","IPM_2021","IPM_2022","Internet_2020","Internet_2021",
  "PDRB_2020","PDRB_2021","PDRB_2022","PertPend_2020","PertPend_2021",
  "PertPend_2022","ProdPangan_2020","ProdPangan_2021","ProdPangan_2022",
  "RLS_2020","RLS_2021","RLS_2022")
raw <- raw[-1,] %>% mutate(across(-Kabupaten,~suppressWarnings(as.numeric(.))))
kalimantan <- raw %>% filter(No >= 315 & No <= 370)
cat("Jumlah awal:", nrow(kalimantan))
## Jumlah awal: 56
## Jumlah awal: 56

4 Data Analisis

data_analisis <- kalimantan %>%
  select(Kabupaten, IKP=IKP_2021, IPM=IPM_2021, Internet=Internet_2021,
         PDRB=PDRB_2021, PertPend=PertPend_2021, ProdPangan=ProdPangan_2021) %>%
  mutate(
    Internet   = log1p(Internet),
    PDRB       = log1p(PDRB),
    ProdPangan = log1p(ProdPangan)
  ) %>% na.omit()
cat("Jumlah observasi final:", nrow(data_analisis))
## Jumlah observasi final: 54
## Jumlah observasi final: 54

5 Statistik Deskriptif

data_analisis %>%
  select(IKP, IPM, Internet,
         PDRB, PertPend, ProdPangan) %>%
  summary()
##       IKP             IPM           Internet          PDRB      
##  Min.   :53.15   Min.   :62.90   Min.   :2.639   Min.   :15.33  
##  1st Qu.:72.28   1st Qu.:67.83   1st Qu.:3.534   1st Qu.:16.01  
##  Median :77.72   Median :70.10   Median :4.119   Median :16.46  
##  Mean   :75.67   Mean   :70.90   Mean   :4.083   Mean   :16.62  
##  3rd Qu.:80.97   3rd Qu.:72.53   3rd Qu.:4.644   3rd Qu.:17.21  
##  Max.   :88.68   Max.   :80.82   Max.   :5.389   Max.   :18.99  
##     PertPend        ProdPangan    
##  Min.   :0.2700   Min.   : 6.993  
##  1st Qu.:0.6725   1st Qu.: 9.265  
##  Median :1.0200   Median :10.424  
##  Mean   :1.2030   Mean   :10.221  
##  3rd Qu.:1.6000   3rd Qu.:11.143  
##  Max.   :4.3600   Max.   :12.693

Interpretasi: IKP memiliki rata-rata 75,67 dengan rentang yang cukup lebar (53,15 – 88,68), mengindikasikan adanya ketimpangan ketahanan pangan yang signifikan antar kabupaten/kota di Kalimantan. Kabupaten dengan IKP terendah (53,15) menggambarkan kondisi ketahanan pangan yang jauh tertinggal dibanding yang tertinggi (88,68). IPM rata-rata 70,90 dengan minimum 62,90, menandakan masih ada daerah dengan kualitas pembangunan manusia yang rendah. PertPend memiliki rata-rata 1,20% namun nilai maksimumnya mencapai 4,36%, menunjukkan beberapa kabupaten mengalami tekanan demografis yang sangat tinggi. Variabel Internet dan PDRB telah ditransformasi log untuk menstabilkan distribusi dan mengurangi pengaruh nilai ekstrem.

6 Korelasi

var_model <- data_analisis %>%
  select(IKP, IPM, Internet,
         PDRB, PertPend, ProdPangan)
cor_matrix <- cor(var_model)
round(cor_matrix, 3)
##               IKP    IPM Internet   PDRB PertPend ProdPangan
## IKP         1.000  0.264    0.225  0.249   -0.373      0.265
## IPM         0.264  1.000   -0.315  0.600   -0.169      0.250
## Internet    0.225 -0.315    1.000  0.214   -0.313      0.267
## PDRB        0.249  0.600    0.214  1.000   -0.108      0.457
## PertPend   -0.373 -0.169   -0.313 -0.108    1.000      0.017
## ProdPangan  0.265  0.250    0.267  0.457    0.017      1.000
corrplot(cor_matrix, method = "color", type = "upper",
         tl.col = "black", addCoef.col = "black", number.cex = 0.7)

Interpretasi: PertPend memiliki korelasi negatif terkuat dengan IKP (r = −0,373), artinya daerah dengan pertumbuhan penduduk lebih tinggi cenderung memiliki ketahanan pangan lebih rendah — mencerminkan tekanan sisi permintaan pangan. IPM, Internet, dan ProdPangan berkorelasi positif dengan IKP (r = 0,264; 0,225; 0,265), namun belum kuat. Korelasi IPM dengan PDRB sebesar 0,600 perlu diwaspadai sebagai indikasi multikolinearitas moderat, meskipun masih perlu dikonfirmasi dengan VIF.

7 Model OLS Full

model_full <- lm(
  IKP ~ IPM + Internet + PDRB +
    PertPend + ProdPangan,
  data = data_analisis)
summary(model_full)
## 
## Call:
## lm(formula = IKP ~ IPM + Internet + PDRB + PertPend + ProdPangan, 
##     data = data_analisis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.256  -2.110   1.655   3.869  11.168 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  33.9873    23.7646   1.430   0.1591  
## IPM           0.5099     0.3909   1.304   0.1984  
## Internet      2.2590     2.1727   1.040   0.3037  
## PDRB         -0.6575     1.8570  -0.354   0.7248  
## PertPend     -3.1034     1.6230  -1.912   0.0618 .
## ProdPangan    1.0733     0.8604   1.247   0.2183  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.572 on 48 degrees of freedom
## Multiple R-squared:  0.2492, Adjusted R-squared:  0.171 
## F-statistic: 3.187 on 5 and 48 DF,  p-value: 0.01453

Interpretasi: Model penuh dengan 5 prediktor secara bersama-sama signifikan (F = 3,187; p = 0,015), namun tidak satu pun variabel signifikan secara individual pada α = 5%. Hanya PertPend yang mendekati signifikan (p = 0,062). Adjusted R² hanya 17,1%, menandakan model ini tidak efisien dan mengindikasikan overfitting akibat terlalu banyak prediktor relatif terhadap ukuran sampel (n = 54). Seleksi variabel diperlukan.

8 Uji Normalitas (Model Full)

res_full <- residuals(model_full)
par(mfrow = c(1, 2))
qqnorm(res_full, pch = 16, col = "steelblue")
qqline(res_full, col = "red", lwd = 2)
hist(res_full, breaks = 10, col = "steelblue", border = "white")

par(mfrow = c(1, 1))
sw_test <- shapiro.test(res_full)
sw_test
## 
##  Shapiro-Wilk normality test
## 
## data:  res_full
## W = 0.90392, p-value = 0.0003879

Interpretasi: Uji Shapiro-Wilk menghasilkan W = 0,904 dan p-value = 0,0004 < 0,05, sehingga H₀ normalitas ditolak. Residual model penuh tidak berdistribusi normal. Hal ini kemungkinan disebabkan oleh adanya beberapa kabupaten dengan nilai ekstrem. Meski dengan n = 54 teorema limit pusat (CLT) mulai berlaku, pelanggaran normalitas ini menandakan perlunya kehati-hatian dalam inferensi.

9 Uji Homoskedastisitas

par(mfrow = c(1, 2))
plot(model_full$fitted.values, res_full, pch = 16, col = "steelblue",
     xlab = "Fitted", ylab = "Residual")
abline(h = 0, col = "red", lwd = 2)
plot(model_full$fitted.values, sqrt(abs(res_full)), pch = 16, col = "steelblue")

par(mfrow = c(1, 1))
bp_test <- bptest(model_full)
bp_test
## 
##  studentized Breusch-Pagan test
## 
## data:  model_full
## BP = 1.6031, df = 5, p-value = 0.9009

Interpretasi: Uji Breusch-Pagan menghasilkan BP = 1,603 dengan p-value = 0,9009 >> 0,05, sehingga H₀ homoskedastisitas gagal ditolak. Varians residual konstan (homoskedastis), artinya asumsi ini terpenuhi. Estimator OLS masih bersifat BLUE untuk komponen ini.

10 VIF

vif_values <- vif(model_full)
data.frame(Variabel = names(vif_values), VIF = round(vif_values, 3)) %>% kable()
Variabel VIF
IPM IPM 2.777
Internet Internet 2.117
PDRB PDRB 2.428
PertPend PertPend 1.352
ProdPangan ProdPangan 1.378

Interpretasi: Semua nilai VIF berada di bawah 3 (jauh di bawah ambang kritis 10), dengan nilai tertinggi pada IPM (VIF = 2,777) dan PDRB (VIF = 2,428). Multikolinearitas tidak menjadi masalah serius dalam model ini. Meski keduanya saling berkorelasi cukup tinggi (r = 0,60), dampaknya terhadap estimasi koefisien masih dalam batas yang dapat diterima.

11 Best Subset Selection

subset_result <- regsubsets(
  IKP ~ IPM + Internet + PDRB +
    PertPend + ProdPangan,
  data = data_analisis, nvmax = 5)
subset_summary <- summary(subset_result)
subset_table <- data.frame(
  n_var    = 1:5,
  Adj_R2   = round(subset_summary$adjr2, 4),
  Cp       = round(subset_summary$cp, 4),
  BIC      = round(subset_summary$bic, 4),
  Variabel = apply(subset_summary$which[,-1], 1,
                   function(x) paste(names(x)[x], collapse = ", ")))
subset_table %>% kable()
n_var Adj_R2 Cp BIC Variabel
1 0.1223 5.0570 -0.0958 PertPend
2 0.1819 2.3347 -0.9493 PertPend, ProdPangan
3 0.1858 3.1104 1.7100 IPM, PertPend, ProdPangan
4 0.1858 4.1254 4.6050 IPM, Internet, PertPend, ProdPangan
5 0.1710 6.0000 8.4531 IPM, Internet, PDRB, PertPend, ProdPangan

Interpretasi: Model dengan 2 variabel (PertPend + ProdPangan) memiliki BIC terkecil (−0,9493), menjadikannya kandidat terbaik menurut kriteria BIC yang menghukum kompleksitas lebih keras. Peningkatan Adj R² dari model 2 variabel (0,1819) ke 3 variabel (0,1858) hanya sebesar 0,004 — manfaat yang terlalu kecil dibanding penalti BIC yang langsung naik ke 1,71. Model 5 variabel (Full) bahkan memiliki Adj R² terendah (0,171), mengkonfirmasi adanya overfitting. Prinsip parsimoni mendukung kuat pemilihan model 2 variabel.

12 Stepwise

model_null <- lm(IKP ~ 1, data = data_analisis)
step_both  <- stepAIC(model_null,
  scope = list(lower = model_null, upper = model_full),
  direction = "both", trace = FALSE)
summary(step_both)
## 
## Call:
## lm(formula = IKP ~ PertPend + ProdPangan, data = data_analisis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.924  -2.898   1.028   4.650  11.520 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.4541     7.6688   8.405  3.4e-11 ***
## PertPend     -4.2117     1.3867  -3.037  0.00376 ** 
## ProdPangan    1.5932     0.7283   2.187  0.03332 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.523 on 51 degrees of freedom
## Multiple R-squared:  0.2127, Adjusted R-squared:  0.1819 
## F-statistic: 6.891 on 2 and 51 DF,  p-value: 0.002244

Interpretasi: Stepwise berbasis AIC arah dua arah (both) secara konsisten memilih model yang identik dengan Best Subset: PertPend dan ProdPangan. Pada model ini:

  • PertPend (β = −4,21; p = 0,004): Signifikan pada α = 1%. Setiap kenaikan 1% pertumbuhan penduduk dikaitkan dengan penurunan IKP sebesar 4,21 poin, ceteris paribus. Ini mencerminkan tekanan demografis nyata terhadap ketahanan pangan.
  • ProdPangan (β = 1,59; p = 0,033): Signifikan pada α = 5%. Setiap kenaikan 1 satuan log produksi pangan dikaitkan dengan kenaikan IKP sebesar 1,59 poin, menunjukkan bahwa kapasitas produksi pangan lokal berperan positif.
  • Adjusted R² = 0,1819: Model menjelaskan ~18,2% variasi IKP. Angka ini tergolong rendah namun wajar mengingat IKP adalah indeks komposit multi-dimensi.
  • F-statistic = 6,891 (p = 0,002): Model secara keseluruhan signifikan.

13 LASSO

x_all <- model.matrix(
  IKP ~ IPM + Internet + PDRB +
    PertPend + ProdPangan, data_analisis)[,-1]
y_all <- data_analisis$IKP

cv_lasso    <- cv.glmnet(x_all, y_all, alpha = 1, nfolds = nrow(data_analisis))
plot(cv_lasso)

lasso_model <- glmnet(x_all, y_all, alpha = 1, lambda = cv_lasso$lambda.1se)
coef_lasso  <- as.matrix(coef(lasso_model))
coef_df <- data.frame(
  Variabel  = rownames(coef_lasso)[-1],
  Koefisien = round(coef_lasso[-1,1], 4),
  Status    = ifelse(coef_lasso[-1,1] == 0, "Dieliminasi", "Aktif"))
coef_df %>% kable(caption = "LASSO - lambda.1se")
LASSO - lambda.1se
Variabel Koefisien Status
IPM IPM 0 Dieliminasi
Internet Internet 0 Dieliminasi
PDRB PDRB 0 Dieliminasi
PertPend PertPend 0 Dieliminasi
ProdPangan ProdPangan 0 Dieliminasi

Interpretasi (lambda.1se): Dengan lambda.1se = 3,0706, LASSO mengeliminasi seluruh 5 variabel, menghasilkan model null. Ini terjadi karena penalti terlalu besar untuk n = 54, bukan karena variabelnya tidak penting. Lambda.1se dirancang konservatif dan pada sampel kecil sering terlalu agresif. Lihat lambda.min sebagai pembanding.

LASSO lambda.min

lasso_min  <- glmnet(x_all, y_all, alpha = 1, lambda = cv_lasso$lambda.min)
coef_min   <- as.matrix(coef(lasso_min))
coef_df_min <- data.frame(
  Variabel  = rownames(coef_min)[-1],
  Koefisien = round(coef_min[-1,1], 4),
  Status    = ifelse(coef_min[-1,1] == 0, "Dieliminasi", "Aktif"))
cat("Lambda.1se :", round(cv_lasso$lambda.1se, 4),
    "\nLambda.min  :", round(cv_lasso$lambda.min, 4), "\n")
## Lambda.1se : 3.0706 
## Lambda.min  : 0.1884
coef_df_min %>%
  kable(caption = "LASSO - lambda.min (lebih longgar)") %>%
  kable_styling(bootstrap_options = c("striped", "hover"))
LASSO - lambda.min (lebih longgar)
Variabel Koefisien Status
IPM IPM 0.3653 Aktif
Internet Internet 1.5838 Aktif
PDRB PDRB 0.0000 Dieliminasi
PertPend PertPend -3.1040 Aktif
ProdPangan ProdPangan 0.9560 Aktif

Interpretasi (lambda.min): Dengan lambda.min = 0,1884 (jauh lebih kecil dari lambda.1se = 3,0706), LASSO mempertahankan 4 dari 5 variabel:

Variabel Koefisien Status
IPM 0,3653 Aktif
Internet 1,5838 Aktif
PDRB 0,0000 Dieliminasi
PertPend −3,1040 Aktif
ProdPangan 0,9560 Aktif

PDRB adalah satu-satunya variabel yang dieliminasi, konsisten dengan koefisiennya yang tidak signifikan di model OLS Full (p = 0,725). Hasil ini mengkonfirmasi bahwa PertPend tetap menjadi prediktor paling penting (koefisien terbesar dalam nilai absolut = −3,104) bahkan setelah regularisasi, selaras dengan temuan Best Subset dan Stepwise.

14 Fungsi LOOCV

loocv_rmse_mae <- function(model_formula, data) {
  n <- nrow(data); errs <- numeric(n)
  for (i in seq_len(n)) {
    fit     <- lm(model_formula, data = data[-i,], na.action = na.exclude)
    pred    <- predict(fit, newdata = data[i,])
    errs[i] <- data$IKP[i] - pred
  }
  c(RMSE = sqrt(mean(errs^2)), MAE = mean(abs(errs)))
}

15 Perbandingan Model

best_n       <- which.min(subset_summary$bic)
best_vars    <- names(which(subset_summary$which[best_n,-1]))
formula_best <- as.formula(paste("IKP ~", paste(best_vars, collapse=" + ")))
vars_step    <- names(coef(step_both))[-1]
formula_step <- as.formula(paste("IKP ~", paste(vars_step, collapse=" + ")))

rmse_full  <- loocv_rmse_mae(formula(model_full), data_analisis)
rmse_best  <- loocv_rmse_mae(formula_best, data_analisis)
rmse_step  <- loocv_rmse_mae(formula_step, data_analisis)
rmse_lasso <- sqrt(min(cv_lasso$cvm))

metrik_df <- data.frame(
  Metode = c("OLS Full","Best Subset","Stepwise","LASSO"),
  RMSE   = c(rmse_full["RMSE"], rmse_best["RMSE"], rmse_step["RMSE"], rmse_lasso),
  MAE    = c(rmse_full["MAE"],  rmse_best["MAE"],  rmse_step["MAE"],  NA))
metrik_df %>% kable()
Metode RMSE MAE
OLS Full 7.951717 5.928998
Best Subset 7.722365 5.754353
Stepwise 7.722365 5.754353
LASSO 7.882631 NA

Interpretasi: Berdasarkan RMSE LOOCV, Best Subset dan Stepwise unggul dengan RMSE = 7,7224 dan MAE = 5,7544 — keduanya lebih kecil dari OLS Full (RMSE 7,9517) dan LASSO (RMSE 7,8826). Ini menunjukkan model 2 variabel bukan hanya lebih sederhana, tetapi juga lebih akurat secara prediktif secara out-of-sample.

Tabel Komprehensif (AIC + BIC + Adj R² + RMSE)

model_best_lm <- lm(formula_best, data = data_analisis)
tabel_komplit <- data.frame(
  Model   = c("OLS Full","Best Subset (2 var)","Stepwise (2 var)"),
  Adj_R2  = round(c(summary(model_full)$adj.r.squared,
                    summary(model_best_lm)$adj.r.squared,
                    summary(step_both)$adj.r.squared), 4),
  AIC     = round(c(AIC(model_full), AIC(model_best_lm), AIC(step_both)), 2),
  BIC     = round(c(BIC(model_full), BIC(model_best_lm), BIC(step_both)), 2),
  RMSE_CV = round(c(rmse_full["RMSE"], rmse_best["RMSE"], rmse_step["RMSE"]), 4),
  MAE_CV  = round(c(rmse_full["MAE"],  rmse_best["MAE"],  rmse_step["MAE"]),  4))

tabel_komplit %>%
  kable(caption = "Perbandingan Komprehensif Semua Model") %>%
  kable_styling(bootstrap_options = c("striped","hover")) %>%
  row_spec(which.min(tabel_komplit$BIC),
           bold = TRUE, background = "#d4edda", color = "black")
Perbandingan Komprehensif Semua Model
Model Adj_R2 AIC BIC RMSE_CV MAE_CV
OLS Full 0.1710 379.53 393.45 7.9517 5.9290
Best Subset (2 var) 0.1819 376.10 384.05 7.7224 5.7544
Stepwise (2 var) 0.1819 376.10 384.05 7.7224 5.7544

Interpretasi: Tabel komprehensif menunjukkan keunggulan konsisten Best Subset / Stepwise (2 var) di seluruh kriteria:

Kriteria OLS Full Best Subset / Stepwise Keunggulan
Adj R² 0,1710 0,1819 +0,011
AIC 379,53 376,10 −3,43
BIC 393,45 384,05 −9,40
RMSE CV 7,9517 7,7224 −0,229
MAE CV 5,9290 5,7544 −0,175

Selisih BIC sebesar 9,4 poin antara OLS Full dan model 2 variabel merupakan bukti kuat (very strong evidence) keunggulan model parsimonious. Model terbaik adalah Best Subset / Stepwise dengan PertPend + ProdPangan.

16 Model Final

rmse_compare <- c(Full = rmse_full["RMSE"], Best = rmse_best["RMSE"],
                  Stepwise = rmse_step["RMSE"])
best_method  <- names(which.min(rmse_compare))
best_method
## [1] "Best.RMSE"
## [1] "Best.RMSE"

if (best_method == "Full") {
  model_final <- model_full
} else if (best_method == "Best") {
  model_final <- lm(formula_best, data = data_analisis)
} else {
  model_final <- step_both
}
summary(model_final)
## 
## Call:
## lm(formula = IKP ~ PertPend + ProdPangan, data = data_analisis)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -21.924  -2.898   1.028   4.650  11.520 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  64.4541     7.6688   8.405  3.4e-11 ***
## PertPend     -4.2117     1.3867  -3.037  0.00376 ** 
## ProdPangan    1.5932     0.7283   2.187  0.03332 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 7.523 on 51 degrees of freedom
## Multiple R-squared:  0.2127, Adjusted R-squared:  0.1819 
## F-statistic: 6.891 on 2 and 51 DF,  p-value: 0.002244

Interpretasi: Kode memilih model dengan RMSE LOOCV terkecil secara otomatis. Output best_method menunjukkan “Best.RMSE”, artinya model yang dipilih adalah hasil Best Subset (PertPend + ProdPangan). Model final identik dengan model Stepwise karena keduanya menghasilkan formula yang sama.

Uji Normalitas Residual Model Final

res_final <- residuals(model_final)
par(mfrow = c(1, 2))
qqnorm(res_final, pch = 16, col = "steelblue", main = "Q-Q Plot - Model Final")
qqline(res_final, col = "red", lwd = 2)
hist(res_final, breaks = 10, col = "steelblue", border = "white",
     main = "Histogram Residual - Model Final", xlab = "Residual")

par(mfrow = c(1, 1))
sw_final <- shapiro.test(res_final)
sw_final
## 
##  Shapiro-Wilk normality test
## 
## data:  res_final
## W = 0.92543, p-value = 0.002414
cat("\nKesimpulan:",
    ifelse(sw_final$p.value > 0.05,
           "Residual NORMAL (H0 gagal ditolak)",
           "Residual TIDAK NORMAL (H0 ditolak)"),
    "\np-value =", round(sw_final$p.value, 4))
## 
## Kesimpulan: Residual TIDAK NORMAL (H0 ditolak) 
## p-value = 0.0024

Interpretasi: Uji Shapiro-Wilk pada residual model final menghasilkan W = 0,925 dan p-value = 0,0024 < 0,05. H₀ normalitas ditolak — residual model final juga tidak normal. Meski demikian, dibanding model full (W = 0,904; p = 0,0004), normalitas model final sedikit membaik (W lebih mendekati 1). Dengan n = 54, pelanggaran ini belum kritis, namun perlu dicatat sebagai keterbatasan. Pemeriksaan Cook’s Distance penting untuk mengidentifikasi apakah outlier adalah penyebab utama.

17 Standardized Coefficient

lm.beta(model_final)
## 
## Call:
## lm(formula = IKP ~ PertPend + ProdPangan, data = data_analisis)
## 
## Standardized Coefficients::
## (Intercept)    PertPend  ProdPangan 
##          NA  -0.3773987   0.2718149

Interpretasi: Koefisien terstandar (β*) memungkinkan perbandingan pengaruh relatif antar prediktor dengan satuan berbeda:

  • PertPend (β* = −0,377): Prediktor dengan pengaruh relatif terbesar. Setiap kenaikan 1 standar deviasi pertumbuhan penduduk dikaitkan dengan penurunan IKP sebesar 0,377 standar deviasi.
  • ProdPangan (β* = 0,272): Pengaruh positif yang lebih lemah. Setiap kenaikan 1 standar deviasi log produksi pangan dikaitkan dengan kenaikan IKP sebesar 0,272 standar deviasi.

Kesimpulan: tekanan demografis (pertumbuhan penduduk) lebih dominan menekan IKP dibanding kontribusi positif produksi pangan lokal.

18 Koefisien Final

coef_final <- summary(model_final)$coefficients
data.frame(
  Variabel  = rownames(coef_final),
  Beta      = round(coef_final[,1], 4),
  Std_Error = round(coef_final[,2], 4),
  t_value   = round(coef_final[,3], 4),
  p_value   = round(coef_final[,4], 4)) %>% kable()
Variabel Beta Std_Error t_value p_value
(Intercept) (Intercept) 64.4541 7.6688 8.4047 0.0000
PertPend PertPend -4.2117 1.3867 -3.0371 0.0038
ProdPangan ProdPangan 1.5932 0.7283 2.1874 0.0333 
confint(model_final, level = 0.95)
##                 2.5 %    97.5 %
## (Intercept) 49.058315 79.849908
## PertPend    -6.995691 -1.427673
## ProdPangan   0.130978  3.055358

Interpretasi: - Intercept (64,45; p < 0,001): Nilai dasar IKP secara matematis ketika kedua prediktor bernilai nol. - PertPend (β = −4,21; CI 95%: [−6,996; −1,428]): Interval kepercayaan seluruhnya negatif, mengkonfirmasi pengaruh negatif yang signifikan dan konsisten. Setiap kenaikan 1% pertumbuhan penduduk dikaitkan dengan penurunan IKP antara 1,43 hingga 7,00 poin. - ProdPangan (β = 1,59; CI 95%: [0,131; 3,055]): Interval kepercayaan seluruhnya positif, mengkonfirmasi pengaruh positif yang signifikan. Peningkatan produksi pangan lokal secara konsisten meningkatkan IKP.

19 Visualisasi Aktual vs Prediksi

pred_final  <- predict(model_final)
resid_final <- residuals(model_final)
df_viz <- data.frame(
  Kabupaten = data_analisis$Kabupaten,
  Aktual    = data_analisis$IKP,
  Prediksi  = pred_final,
  Residual  = resid_final)
Q1 <- quantile(abs(resid_final), 0.25)
Q3 <- quantile(abs(resid_final), 0.75)
outlier_batas <- Q3 + 1.5 * IQR(abs(resid_final))
df_viz$Label  <- ifelse(abs(df_viz$Residual) > outlier_batas, df_viz$Kabupaten, "")
ggplot(df_viz, aes(x = Aktual, y = Prediksi)) +
  geom_point(aes(color = abs(Residual)), size = 3) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed") +
  geom_smooth(method = "lm", se = TRUE) +
  geom_text_repel(aes(label = Label), color = "red") +
  theme_bw() +
  labs(title = "Aktual vs Prediksi IKP 2021", color = "|Residual|")

20 Plot Diagnostik

par(mfrow = c(2, 2))
plot(model_final, which = 1:4, pch = 16, col = "steelblue")

par(mfrow = c(1, 1))

21 Cook’s Distance

cook_cutoff <- 4 / nrow(data_analisis)
cat("Cutoff Cook's Distance:", round(cook_cutoff, 4))
## Cutoff Cook's Distance: 0.0741

Tabel Kabupaten Berpengaruh

cooksd          <- cooks.distance(model_final)
influential_idx <- which(cooksd > cook_cutoff)
influential_df  <- data.frame(
  Kabupaten      = data_analisis$Kabupaten[influential_idx],
  Cooks_D        = round(cooksd[influential_idx], 4),
  IKP_Aktual     = data_analisis$IKP[influential_idx],
  PertPend  = data_analisis$PertPend[influential_idx],
  ProdPangan = round(data_analisis$ProdPangan[influential_idx], 3)
) %>% arrange(desc(Cooks_D))
cat("Jumlah observasi berpengaruh:", nrow(influential_df), "\n")
## Jumlah observasi berpengaruh: 3
influential_df %>%
  kable(caption = "Kabupaten dengan Cook's Distance > Cutoff") %>%
  kable_styling(bootstrap_options = c("striped","hover"))
Kabupaten dengan Cook’s Distance > Cutoff
Kabupaten Cooks_D IKP_Aktual PertPend ProdPangan
14 Kotawaringin Barat 0.2863 60.10 0.59 12.588
24 Gunung Mas 0.1906 53.15 1.68 7.397
52 Tana Tidung 0.0758 56.00 4.36 8.408 
plot(cooksd, type = "h", col = "steelblue", pch = 16,
     main = "Cook's Distance per Observasi",
     ylab = "Cook's Distance", xlab = "Indeks Observasi")
abline(h = cook_cutoff, col = "red", lty = 2, lwd = 2)
text(influential_idx, cooksd[influential_idx],
     labels = data_analisis$Kabupaten[influential_idx],
     pos = 3, cex = 0.7, col = "red")

Interpretasi: Terdapat 3 kabupaten dengan Cook’s Distance melebihi cutoff (4/54 = 0,0741):

Kabupaten Cook’s D IKP PertPend ProdPangan
Kotawaringin Barat 0,2863 60,10 0,59 12,588
Gunung Mas 0,1906 53,15 1,68 7,397
Tana Tidung 0,0758 56,00 4,36 8,408
  • Kotawaringin Barat (Cook’s D = 0,286) adalah observasi paling berpengaruh. Meski pertumbuhannya rendah (0,59%), IKP-nya hanya 60,10 dengan ProdPangan sangat tinggi (12,588) — kombinasi yang anomali dan berdampak besar pada koefisien.
  • Gunung Mas memiliki IKP terendah dalam dataset (53,15), menjadikannya leverage point yang signifikan.
  • Tana Tidung memiliki pertumbuhan penduduk tertinggi (4,36%) dengan IKP rendah (56,00) — konsisten dengan arah model, namun nilainya ekstrem.

Ketiga kabupaten ini perlu diverifikasi datanya. Jika valid, perlu dipertimbangkan analisis sensitivitas dengan dan tanpa observasi berpengaruh tersebut.

22 Ringkasan Hasil

data.frame(
  Aspek  = c("Jumlah observasi","Variabel kandidat","Model final","RMSE","MAE"),
  Hasil  = c(nrow(data_analisis), 5, best_method,
             round(min(rmse_compare), 4),
             round(min(c(rmse_full["MAE"], rmse_best["MAE"], rmse_step["MAE"])), 4))
) %>% kable()
Aspek Hasil
Jumlah observasi 54
Variabel kandidat 5
Model final Best.RMSE
RMSE 7.7224
MAE 5.7544

23 Keterbatasan Penelitian

  • Data bersifat cross-sectional sehingga hubungan yang diperoleh tidak dapat diinterpretasikan sebagai hubungan kausal.
  • Jumlah observasi relatif kecil (n = 54) sehingga pemilihan model dapat sensitif terhadap observasi tertentu — terbukti dari 3 kabupaten dengan Cook’s D tinggi.
  • IKP merupakan indeks komposit sehingga kemungkinan masih terdapat variabel penting lain yang belum dimasukkan (kemiskinan, stunting, infrastruktur, logistik).
  • Residual tidak normal pada model final (p = 0,002) — inferensi koefisien perlu ditafsirkan dengan kehati-hatian.
  • LASSO dengan lambda.1se terlalu konservatif untuk n = 54; hasil lambda.min lebih representatif untuk pembanding seleksi variabel.

24 Kesimpulan & Penentuan Model Terbaik

Model Terbaik: Best Subset / Stepwise (2 Variabel)

Berdasarkan seluruh kriteria evaluasi dari hasil syntax di atas, model terbaik adalah model dengan 2 prediktor: PertPend dan ProdPangan, yang dihasilkan secara konsisten oleh metode Best Subset Selection dan Stepwise Regression.

Dasar Pemilihan

Kriteria OLS Full Best Subset / Stepwise LASSO (1se)
Adj R² 0,1710 0,1819 ✅ — (null)
AIC 379,53 376,10 ✅ —
BIC 393,45 384,05 ✅ —
RMSE LOOCV 7,9517 7,7224 ✅ 7,8826
MAE LOOCV 5,9290 5,7544 ✅ —
Var. signifikan Tidak ada Ya (keduanya) ✅ —
Parsimoni 5 var 2 var ✅ 0 var

Selisih BIC antara OLS Full dan model terpilih sebesar 9,4 poin merupakan bukti sangat kuat (very strong evidence) keunggulan model parsimonious (Raftery, 1995).

Persamaan Model Final

\[\widehat{IKP}_{2021} = 64{,}45 - 4{,}21 \cdot \text{PertPend}_{2021} + 1{,}59 \cdot \log(\text{ProdPangan}_{2021})\]

Interpretasi Substantif

  • Pertumbuhan penduduk (β* = −0,377) adalah prediktor dominan yang menekan ketahanan pangan. Setiap kenaikan 1% pertumbuhan penduduk dikaitkan dengan penurunan IKP rata-rata 4,21 poin (CI 95%: −6,996 s.d. −1,428).
  • Produksi pangan lokal (β* = 0,272) berkontribusi positif yang signifikan. Peningkatan kapasitas produksi lokal terbukti meningkatkan ketahanan pangan (CI 95%: 0,131 s.d. 3,055).
  • Hasil LASSO dengan lambda.min mengkonfirmasi pentingnya PertPend (koefisien −3,104) sebagai prediktor paling robust di antara semua variabel kandidat, selaras dengan temuan Best Subset dan Stepwise.

Implikasi Kebijakan

Peningkatan ketahanan pangan di Kalimantan sebaiknya difokuskan pada dua strategi utama yang didukung oleh model: (1) pengendalian laju pertumbuhan penduduk di kabupaten/kota dengan pertumbuhan tinggi seperti Tana Tidung, dan (2) peningkatan kapasitas produksi pangan lokal terutama di daerah dengan ProdPangan rendah seperti Gunung Mas.