1 Pendahuluan

Penelitian ini menganalisis hubungan antara kebutuhan pangan protein (KPP), pembangunan manusia (IPM), laju pertumbuhan penduduk (LPP), Produk Domestik Regional Bruto (PDRB), dan tingkat pengangguran (TPT) di provinsi-provinsi Indonesia tahun 2024 menggunakan regresi Gamma.

Regresi Gamma merupakan bagian dari Generalized Linear Model (GLM) yang cocok untuk data respon:

kontinu,
positif (Y > 0),
dan berdistribusi condong ke kanan (right-skewed).

Selain estimasi klasik menggunakan Maximum Likelihood Estimation (MLE), laporan ini juga membahas:

Metode optimasi (Newton-Raphson) untuk regresi Gamma,
Simulasi Monte Carlo untuk mengkaji sifat penaksir MLE,
Markov Chain Monte Carlo (MCMC) berbasis Metropolis–Hastings untuk pendekatan Bayesian.

2 Data

Data yang digunakan berisi 38 provinsi dengan variabel:

KPP (Y) : Kebutuhan Pangan Protein (gram),
IPM (X₁): Indeks Pembangunan Manusia,
LPP (X₂): Laju Pertumbuhan Penduduk,
PDRB (X₃): Produk Domestik Regional Bruto,
TPT (X₄): Tingkat Pengangguran Agustus (%).

# Membaca data dari file CSV (pilih manual, misal 'Data.csv')
data <- read.csv(file.choose(), header = TRUE, sep = ",", dec = ".")

# Pastikan nama kolom sesuai: Provinsi, Pangan, IPM, LPP, PDRB, TPT
str(data)

## 'data.frame':    38 obs. of  6 variables:
##  $ Provinsi: chr  "Aceh" "Sumatera Utara" "Sumatera Barat" "Riau" ...
##  $ Pangan  : num  60.6 62.1 57.8 59.2 58.7 ...
##  $ IPM     : num  74 74 74.5 74.8 73.4 ...
##  $ LPP     : num  1.39 1.4 1.43 1.37 1.3 1.15 1.32 1.2 1.36 1.5 ...
##  $ PDRB    : num  243202 1146920 332936 1112482 322976 ...
##  $ TPT     : num  5.75 5.6 5.75 3.7 4.48 3.86 3.11 4.19 4.63 6.39 ...

# Definisikan respon dan kovariat
data$Y <- as.numeric(data$Pangan)

Interpretasi:
Chunk ini membaca data dari file CSV, menampilkan struktur data, dan mendefinisikan variabel respon Y sebagai KPP dalam bentuk numerik.

3 Teori Regresi Gamma

3.1 Distribusi Gamma

Jika \(Y\) memiliki distribusi Gamma dengan mean \(\mu\) dan parameter dispersi \(\phi\), maka salah satu bentuk densitasnya adalah:

\[ f_Y(y \mid \mu, \phi) = \frac{1}{\Gamma(1/\phi)} \left(\frac{1}{\phi \mu}\right)^{1/\phi} y^{\frac{1}{\phi}-1} \exp\!\left(-\frac{y}{\phi \mu}\right), \quad y > 0, \]

dengan:

\(\mathbb{E}(Y \mid X) = \mu\),
\(\mathrm{Var}(Y \mid X) = \phi \mu^2\).

3.2 Struktur GLM Gamma dengan Link Log

Bentuk umum model regresi Gamma:

\[ Y_i \mid X_i \sim \text{Gamma}(\mu_i, \phi), \quad i = 1,\dots,n. \]

Dengan fungsi link log:

\[ g(\mu_i) = \log(\mu_i) = \eta_i = X_i^\top \beta, \]

maka:

\[ \mu_i = \exp(X_i^\top \beta), \]

dan model regresi Gamma dengan link log dapat ditulis:

\[ \log(\mu_i) = \beta_0 + \beta_1 \text{IPM}_i + \beta_2 \text{LPP}_i + \beta_3 \text{PDRB}_i + \beta_4 \text{TPT}_i. \]

3.3 Asumsi Regresi Gamma

Variabel respon positif: \(Y_i > 0\).
Varians proporsional terhadap kuadrat mean: \(\mathrm{Var}(Y_i \mid X_i) \propto \mu_i^2\).
Distribusi respon condong ke kanan (right-skewed).
Hubungan mean dan prediktor linear melalui \(g(\mu_i) = \log(\mu_i)\).
Observasi \(Y_i\) saling independen.
Kovariat tidak memiliki multikolinearitas serius.

4 Uji Asumsi Regresi Gamma

# Model Gamma awal untuk uji asumsi
fit_uji <- glm(
  Pangan ~ IPM + LPP + PDRB + TPT,
  family = Gamma(link = "log"),
  data   = data
)
summary(fit_uji)

## 
## Call:
## glm(formula = Pangan ~ IPM + LPP + PDRB + TPT, family = Gamma(link = "log"), 
##     data = data)
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  2.544e+00  2.517e-01  10.109 1.23e-11 ***
## IPM          2.118e-02  3.202e-03   6.616 1.59e-07 ***
## LPP          4.999e-02  5.364e-02   0.932    0.358    
## PDRB         2.624e-09  1.891e-08   0.139    0.890    
## TPT         -1.260e-02  1.073e-02  -1.174    0.249    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for Gamma family taken to be 0.00603871)
## 
##     Null deviance: 0.50303  on 37  degrees of freedom
## Residual deviance: 0.19190  on 33  degrees of freedom
## AIC: 229.49
## 
## Number of Fisher Scoring iterations: 4

Interpretasi:
Chunk ini mengestimasi model regresi Gamma dengan KPP sebagai respon dan IPM, LPP, PDRB, TPT sebagai kovariat menggunakan link log. Output summary(fit_uji) memberikan:

estimasi koefisien \(\hat\beta\),
standard error,
nilai t dan p-value,
serta estimasi parameter dispersi \(\hat\phi\).

4.1 Respon Positif (Y > 0)

sum(data$Y <= 0)

## [1] 0

if(sum(data$Y <= 0) == 0){
  cat("Semua Y > 0 (terpenuhi)\n")
} else {
  cat("Ada Y <= 0 (tidak terpenuhi)\n")
}

## Semua Y > 0 (terpenuhi)

Interpretasi:
Untuk regresi Gamma, semua nilai Y harus positif. Jika jumlah Y <= 0 adalah 0, maka asumsi ini terpenuhi.

4.2 Var(Y) ∝ mean(Y)^2

mu_uji  <- fitted(fit_uji)
plot(mu_uji, (data$Y - mu_uji)^2,
     main = "Varian vs Mean^2",
     xlab = "Mean (μ)",
     ylab = "(Y - μ)^2")

Interpretasi:
Plot \(\mu\) versus \((Y-\mu)^2\) menggambarkan apakah varian cenderung meningkat ketika mean meningkat. Pola menyebar naik seiring naiknya \(\mu\) mendukung asumsi Gamma bahwa \(\mathrm{Var}(Y \mid X) \propto \mu^2\).

4.3 Distribusi Gamma (Right-Skewed)

hist(data$Y,
     main = "Histogram Y (KPP)",
     xlab = "Y",
     col  = "lightblue",
     border = "white")

library(e1071)
skewness(data$Y)

## [1] -0.7736714

Interpretasi:
Histogram Y yang condong ke kanan dan nilai skewness > 0 menunjukkan distribusi right-skewed yang konsisten dengan asumsi distribusi Gamma.

4.4 Linearitas log(Y) vs η

eta_uji <- fit_uji$linear.predictors
correlation <- cor(log(data$Y), eta_uji)
correlation

## [1] 0.8019885

plot(eta_uji, log(data$Y),
     pch = 19,
     col = "blue",
     xlab = "Linear Predictor η = Xβ",
     ylab = "log(Y)",
     main = "log(Y) vs Linear Predictor η")
abline(lm(log(data$Y) ~ eta_uji), col = "red", lwd = 2)

Interpretasi:
- Nilai korelasi correlation yang cukup besar (mendekati 1) menunjukkan hubungan hampir linier antara log(Y) dan η.
- Plot log(Y) vs η yang mengikuti garis lurus (garis merah) mendukung asumsi linearitas pada skala link log.

4.5 Independensi Pengamatan

acf(data$Y, main = "ACF Respon Y (KPP)")

Interpretasi:
Plot ACF digunakan untuk mendeteksi autokorelasi. Jika spike-spike ACF pada berbagai lag berada dalam batas kepercayaan, maka tidak ada autokorelasi kuat dan asumsi independensi dapat dianggap wajar.

4.6 Multikolinearitas

library(car)
vif(fit_uji)

##      IPM      LPP     PDRB      TPT 
## 1.666369 1.812292 1.683272 1.411736

Interpretasi:
Nilai VIF yang relatif kecil (misal < 5) menunjukkan tidak ada multikolinearitas serius antara IPM, LPP, PDRB, dan TPT. Jika semua VIF jauh di bawah 10, asumsi ini dianggap terpenuhi.

5 Metode Penaksiran: Newton-Raphson

Untuk data \((Y_i, X_i)\) dan model Gamma dengan link log:

\[ \log(\mu_i) = X_i^\top \beta, \quad \mu_i = \exp(X_i^\top \beta), \]

log-likelihood (tanpa konstanta) dapat ditulis:

\[ \ell(\beta) \propto -\frac{1}{\phi} \sum_{i=1}^n \left( \log \mu_i + \frac{Y_i}{\mu_i} \right). \]

Gradien (turunan pertama) terhadap \(\beta\):

\[ g(\beta) = \frac{\partial \ell}{\partial \beta} = \frac{1}{\phi} X^\top \bigg(\frac{Y}{\mu} - 1\bigg). \]

Hessian (turunan kedua) terhadap \(\beta\):

\[ H(\beta) = \frac{\partial^2 \ell}{\partial \beta \partial \beta^\top} = \frac{1}{\phi} X^\top W X, \]

dengan \(W = \mathrm{diag}(Y_i / \mu_i)\).

Iterasi Newton-Raphson:

\[ \beta^{(t+1)} = \beta^{(t)} - H(\beta^{(t)})^{-1} g(\beta^{(t)}). \]

6 Simulasi Monte Carlo

Untuk mempelajari sifat penaksir \(\hat\beta\) (MLE) digunakan simulasi Monte Carlo.

Untuk setiap simulasi \(j = 1,\dots,M\):

Tentukan ukuran sampel \(n\) dan parameter “benar” \(\beta^{\text{true}}\), \(\phi\).
Hitung \(\mu_i = \exp(X_i^\top \beta^{\text{true}})\).
Bangkitkan respon: \[ Y_i^{(j)} \sim \text{Gamma}(\text{shape} = 1/\phi, \text{scale} = \phi \mu_i). \]
Estimasi \(\hat\beta^{(j)} = \arg\max_\beta \ell(\beta \mid Y^{(j)}, X)\) dengan Newton-Raphson (atau GLM).
Simpan \(\hat\beta^{(j)}\) dalam matriks: \[ B = [\hat\beta^{(j)}]_{j=1}^M \in \mathbb{R}^{M \times p}. \]

Dari matriks \(B\) dihitung:

Bias: \[ \text{bias}(\hat\beta_k) = \bar\beta_k - \beta_k^{\text{true}}, \quad \bar\beta_k = \frac{1}{M} \sum_{j=1}^M \hat\beta_k^{(j)}. \]
Empirical standard error (SE): \[ \widehat{\text{SE}}_{\text{emp}}(\hat\beta_k) = \sqrt{\frac{1}{M-1} \sum_{j=1}^M (\hat\beta_k^{(j)} - \bar\beta_k)^2 }. \]
Mean Squared Error (MSE): \[ \text{MSE}(\hat\beta_k) = \frac{1}{M} \sum_{j=1}^M (\hat\beta_k^{(j)} - \beta_k^{\text{true}})^2. \]

Semakin besar \(n\), biasanya bias dan SE cenderung menurun, sehingga MSE juga menurun dan penaksir menjadi konsisten.

7 Markov Chain Monte Carlo (MCMC)

Kerangka Bayesian Regresi Gamma

Dalam pendekatan Bayesian, tujuan utama adalah melakukan penarikan sampel dari distribusi posterior bersama:

\[ p(\beta, \phi \mid Y, X), \]

dengan \(\beta = (\beta_0, \beta_1, \dots, \beta_p)^\top\) dan \(\phi\) parameter dispersi.

Misalkan prior:

\[ \beta \sim N_p(0, \tau^2 I_p), \quad \phi \sim \text{Inverse-Gamma}(a, b). \]

Posterior (tidak terstandarisasi):

\[ p(\beta, \phi \mid Y, X) \propto L(Y \mid \beta, \phi, X) \, p(\beta) \, p(\phi). \]

Karena bentuk posterior tidak dapat dihitung secara tertutup, digunakan algoritma Metropolis–Hastings (MH) dalam kerangka MCMC untuk menghasilkan deret sampel yang, setelah konvergen, mengikuti distribusi posterior tersebut.

Skema MCMC (Metropolis–Hastings)

Untuk setiap iterasi \(t = 1,\dots,T\):

A. Update \(\beta\) (random-walk multivariate normal)

Usulkan kandidat: \[ \beta^ \sim N_p(\beta^{(t-1)}, \Sigma_\beta), \quad \Sigma_\beta = \sigma_\beta^2 I_p. \]
Hitung rasio penerimaan: \[ \alpha_\beta = \min\Bigg\{ 1,\; \frac{L(Y \mid \beta^, \phi^{(t-1)}) p(\beta^)} {L(Y \mid \beta^{(t-1)}, \phi^{(t-1)}) p(\beta^{(t-1)})} \Bigg\}. \]
Terima dengan probabilitas \(\alpha_\beta\).

B. Update \(\phi\) (random-walk pada skala log)

Usulkan kandidat pada skala log: \[ \log \phi^ \sim N(\log \phi^{(t-1)}, \sigma_\phi^2), \quad \phi^ = \exp(\log \phi^). \]
Hitung rasio penerimaan: \[ \alpha_\phi = \min\Bigg\{ 1,\; \frac{L(Y \mid \beta^{(t)}, \phi^) p(\phi^)} {L(Y \mid \beta^{(t)}, \phi^{(t-1)}) p(\phi^{(t-1)})} \Bigg\}. \]
Terima dengan probabilitas \(\alpha_\phi\).

Set sampel \(\{(\beta^{(t)}, \phi^{(t)})\}_{t=1}^T\) membentuk rantai MCMC yang (setelah konvergen) mengikuti distribusi posterior.

Ringkasan Inferensial Posterior

Setelah burn-in dan thinning, diperoleh \(S\) sampel posterior \(\{(\beta^{(s)}, \phi^{(s)})\}_{s=1}^S\).

Posterior mean: \[ \bar\beta_k = \frac{1}{S} \sum_{s=1}^S \beta_k^{(s)}. \]
Posterior standard deviation: \[ \text{SD}(\beta_k) = \sqrt{\frac{1}{S-1} \sum_{s=1}^S (\beta_k^{(s)} - \bar\beta_k)^2 }. \]
Credible interval 95%: \[ \text{CI}_{95\%}(\beta_k) = [\beta_{k,(0.025)}, \beta_{k,(0.975)}], \] yaitu kuantil ke-2.5% dan 97.5% dari sampel posterior.

Diagnostik Konvergensi MCMC

Untuk menilai apakah rantai sudah konvergen digunakan:

Traceplot: pola “ulat gemuk” tanpa tren.
Density plot: distribusi posterior stabil, unimodal.
Cumulative mean plot: rata-rata kumulatif mendatar setelah burn-in.
Autocorrelation Function (ACF): autokorelasi turun cepat ke nol.
Gelman–Rubin R-hat (untuk ≥ 2 rantai): nilai R-hat < 1.1 indikasi baik.

8 Kesimpulan

Regresi Gamma dengan link log cocok digunakan untuk data KPP yang positif dan cenderung miring ke kanan.
Uji asumsi menunjukkan bahwa:
- semua Y > 0,
- pola varian vs mean konsisten dengan \(\mathrm{Var}(Y) \propto \mu^2\),
- tidak terdapat multikolinearitas serius,
- dan tidak ada autokorelasi kuat yang mengganggu.

Uji Asumsi Regresi Gamma

Ega Saherti

07 December 2025