1 FUNCTION 1 – METODE OPTIMASI

1.1 Paket yang Dibutuhkan

ip <- rownames(installed.packages())
pkg_need <- c("MASS","ggplot2","dplyr","mvtnorm","coda","stableGR")
for(p in pkg_need){
  if(!(p %in% ip)) install.packages(p, dependencies = TRUE)
}
library(MASS)
library(ggplot2)
library(dplyr)
library(mvtnorm)
library(coda)
library(stableGR)

Interpretasi: Bagian ini mengecek apakah paket MASS, ggplot2, dplyr, mvtnorm, coda, dan stableGR sudah terinstal. Jika belum, paket akan di-install otomatis, lalu dipanggil dengan library(). Paket-paket ini digunakan untuk:

  • MASS: fungsi ginv() dan mvrnorm()
  • ggplot2: visualisasi (histogram Monte Carlo)
  • dplyr: manipulasi data (bind_rows, dll.)
  • mvtnorm: fungsi dmvnorm() untuk prior multivariat normal
  • coda: analisis output MCMC (traceplot, densplot, cumuplot, acf)
  • stableGR: menghitung Gelman–Rubin yang stabil (R-hat)

1.2 IMPORT DATA

df <- read.csv(file.choose(), header = TRUE, sep = ",", dec = ".")
y <- as.numeric(df$Pangan)
X <- model.matrix(~ IPM + LPP + PDRB + TPT, data = df)
n <- nrow(X)
p <- ncol(X)

Interpretasi Tahap 1 – Import Data:

  • read.csv(...): Data diambil dari file .csv yang dipilih pengguna. Data disimpan sebagai df. Di sini diasumsikan ada kolom: Pangan, IPM, LPP, PDRB, TPT.
  • y <- as.numeric(df$Pangan): y adalah variabel respon (Pangan) dalam bentuk vektor numerik.
  • X <- model.matrix(~ IPM + LPP + PDRB + TPT, ...): Membentuk matriks desain yang memuat:
    • Kolom intercept (konstanta)
    • Kolom untuk IPM, LPP, PDRB, dan TPT
  • n dan p:
    • n: jumlah observasi
    • p: jumlah parameter (intercept + 4 kovariat). Nilai ini dipakai di seluruh metode (Newton, Fisher, BFGS, bootstrap, Monte Carlo, MCMC).

1.3 FIT GAMMA REGRESSION (GLM)

fit_glm <- glm(Pangan ~ IPM + LPP + PDRB + TPT,
               data = df, family = Gamma(link = "log"))

beta_glm <- coef(fit_glm)
phi_glm  <- summary(fit_glm)$dispersion

cat("=== ESTIMASI MLE (GLM) ===
")
## === ESTIMASI MLE (GLM) ===
print(beta_glm)
##   (Intercept)           IPM           LPP          PDRB           TPT 
##  2.544493e+00  2.118425e-02  4.999287e-02  2.623636e-09 -1.259747e-02
cat("Dispersion =", phi_glm, "

")
## Dispersion = 0.00603871

Interpretasi Tahap 2 – GLM Gamma:

  • glm(..., family = Gamma(link = "log")) mengestimasi model: \[ \log(\mu_i) = \beta_0 + \beta_1 \mathrm{IPM}_i + \beta_2 \mathrm{LPP}_i + \beta_3 \mathrm{PDRB}_i + \beta_4 \mathrm{TPT}_i \] dengan asumsi \(Y_i \mid X_i \sim \text{Gamma}(\mu_i, \phi)\). Estimasi dilakukan dengan metode MLE (Fisher scoring).
  • beta_glm berisi estimasi \(\hat\beta_0, \hat\beta_1, \hat\beta_2, \hat\beta_3, \hat\beta_4\).
    Tanda positif/negatif menunjukkan arah pengaruh terhadap log(rata-rata Pangan).
  • phi_glm adalah estimasi parameter dispersi \(\phi\). Pada Gamma: \[ \mathrm{Var}(Y_i \mid X_i) = \phi\,\mu_i^2 \] Semakin besar \(\phi\), semakin besar variasi data di sekitar mean.

1.4 LOG-LIKELIHOOD, GRADIENT, HESSIAN, MATRIX INFORMATION

loglik_gamma <- function(beta, X, y, phi){
  eta <- X %*% beta
  eta <- pmin(pmax(eta, -30), 30)  # stabilisasi
  mu  <- exp(eta)

  ll <- -(1/phi) * sum(eta + y/mu)
  return(ll)
}

grad_gamma <- function(beta, X, y, phi){
  eta <- X %*% beta
  eta <- pmin(pmax(eta, -30), 30)
  mu  <- exp(eta)

  g <- t(X) %*% ((y/mu) - 1) / phi
  return(as.vector(g))
}

hess_gamma <- function(beta, X, y, phi){
  eta <- X %*% beta
  eta <- pmin(pmax(eta, -30), 30)
  mu  <- exp(eta)

  w <- as.vector(y/mu)      # bobot
  W <- diag(w)

  H <- (t(X) %*% W %*% X) / phi   # Hessian log-likelihood (tanda negatif di inf_gamma)
  return(H)
}

inf_gamma <- function(beta, X, y, phi){
  H <- hess_gamma(beta, X, y, phi)
  I <- -H                           # matriks informasi Fisher = -H
  return(I)
}

Interpretasi Tahap 3 – Struktur Likelihood & Turunan:

  • loglik_gamma: menghitung log-likelihood (dalam bentuk yang sebanding) untuk model Gamma dengan link log.
    eta = X %*% beta dipotong di [-30,30] untuk mencegah overflow/underflow saat exp(eta).
  • Clipping eta di [-30, 30]:
    • Jika eta terlalu besar → exp(eta) bisa Inf.
    • Jika terlalu kecil → exp(eta) sangat dekat 0.
    • Rentang [-30, 30] cukup aman untuk komputasi floating point.
  • grad_gamma: turunan pertama log-likelihood terhadap vektor parameter \(\beta\); dipakai untuk Newton–Raphson dan BFGS.
  • hess_gamma: matriks turunan kedua log-likelihood (Hessian). W adalah matriks bobot dengan elemen \(y/\mu\).
  • inf_gamma: matriks informasi Fisher \(I = -H\); dipakai di Fisher Scoring.

1.5 METODE NEWTON–RHAPSON

newton_gamma <- function(X, y, phi, beta_init = NULL,
                         maxit = 1000, tol = 1e-8){

  p <- ncol(X)
  if(is.null(beta_init)) beta <- rep(0, p) else beta <- beta_init

  for(iter in 1:maxit){
    g <- grad_gamma(beta, X, y, phi)
    H <- hess_gamma(beta, X, y, phi)

    # REGULARIZER untuk menghindari singular
    H_reg <- H + diag(1e-6, p)

    # Cek apakah matriks numerically valid
    if(rcond(H_reg) < 1e-12){
      step <- MASS::ginv(H_reg) %*% g
    } else {
      step <- solve(H_reg, g)
    }

    beta_new <- beta - step

    if(max(abs(beta_new - beta)) < tol){
      return(list(beta = as.vector(beta_new),
                  iter = iter,
                  converged = TRUE))
    }

    beta <- beta_new
  }

  list(beta = as.vector(beta), iter = maxit, converged = FALSE)
}

start_time <- Sys.time()
nr_res <- newton_gamma(X, y, phi_glm, beta_init = beta_glm)
end_time <- Sys.time()

cat("=== NEWTON–RHAPSON ===
")
## === NEWTON–RHAPSON ===
print(nr_res$beta)
## [1]  2.544493e+00  2.118425e-02  4.999287e-02  2.623630e-09 -1.259747e-02
cat("Iterasi:", nr_res$iter, "
")
## Iterasi: 1
cat("Konvergen:", nr_res$converged, "
")
## Konvergen: TRUE
cat("Waktu (detik):", round(end_time - start_time, 4), "

")
## Waktu (detik): 0.0276

Interpretasi Tahap 4 – Newton–Raphson:

  • Inisialisasi beta: default vektor nol, tetapi di sini dimulai dari beta_glm sehingga lebih dekat ke optimum.
  • Regularisasi H_reg = H + 1e-6 I mencegah error matrix is singular.
  • Jika rcond(H_reg) sangat kecil → gunakan pseudo-inverse MASS::ginv().
  • Update Newton–Raphson: \[ \beta_{baru} = \beta_{lama} - H^{-1} g \]
  • Kriteria konvergensi: perubahan maksimum antar komponen \(< \text{tol}\).
  • Hasil: nr_res$beta hampir identik dengan beta_glm → implementasi Newton–Raphson konsisten dengan GLM.

1.6 METODE FISHER SCORING

fisher_gamma <- function(X, y, phi, beta_init = NULL,
                         maxit = 100, tol = 1e-8){

  p <- ncol(X)
  if(is.null(beta_init)) beta <- rep(0, p) else beta <- beta_init

  for(iter in 1:maxit){

    g     <- grad_gamma(beta, X, y, phi)
    I_mat <- inf_gamma(beta, X, y, phi)

    # regularizer (hindari singular)
    I_reg <- I_mat + diag(1e-6, p)

    if(rcond(I_reg) < 1e-12){
      step <- MASS::ginv(I_reg) %*% g
    } else {
      step <- solve(I_reg, g)
    }

    beta_new <- beta + step     # fisher scoring pakai "+"

    if(max(abs(beta_new - beta)) < tol){
      return(list(beta = as.vector(beta_new),
                  iter = iter,
                  converged = TRUE))
    }

    beta <- beta_new
  }

  list(beta = as.vector(beta), iter = maxit, converged = FALSE)
}

start_time <- Sys.time()
fs_res <- fisher_gamma(X, y, phi_glm, beta_init = beta_glm)
end_time <- Sys.time()

cat("=== FISHER SCORING ===
")
## === FISHER SCORING ===
print(fs_res$beta)
## [1]  2.544493e+00  2.118425e-02  4.999287e-02  2.623630e-09 -1.259747e-02
cat("Iterasi:", fs_res$iter, "
")
## Iterasi: 1
cat("Konvergen:", fs_res$converged, "
")
## Konvergen: TRUE
cat("Waktu (detik):", round(end_time - start_time, 4), "

")
## Waktu (detik): 0.0547

Interpretasi Tahap 5 – Fisher Scoring:

  • Berbeda dengan Newton yang memakai Hessian \(H\), Fisher memakai matriks informasi \(I = -H\).
  • Karena inf_gamma mengembalikan \(-H\), langkah menjadi: \[ \beta_{baru} = \beta_{lama} + I^{-1} g \]
  • Regularisasi I_reg memastikan kestabilan numerik.
  • Karena dimulai dari beta_glm, algoritma langsung konvergen dan hasil \(\hat\beta\) konsisten dengan GLM dan Newton.

1.7 OPTIM BUILT-IN (BFGS)

neglog  <- function(b) -loglik_gamma(b, X, y, phi_glm)
neggrad <- function(b) -grad_gamma(b, X, y, phi_glm)

start_time <- Sys.time()
opt_res <- optim(beta_glm, neglog, neggrad, method = "BFGS")
end_time <- Sys.time()

cat("=== optim() BFGS ===
")
## === optim() BFGS ===
print(opt_res$par)
##   (Intercept)           IPM           LPP          PDRB           TPT 
##  2.544493e+00  2.118425e-02  4.999287e-02  2.623655e-09 -1.259747e-02
cat("Iterasi:", opt_res$counts[1], "
")
## Iterasi: 24
cat("Waktu (detik):", round(end_time - start_time, 4), "

")
## Waktu (detik): 0.0202

Interpretasi Tahap 6 – BFGS dengan optim():

  • optim() meminimalkan fungsi, sehingga digunakan negatif log-likelihood (neglog) dan negatif gradient (neggrad).
  • method = "BFGS" adalah metode quasi-Newton yang mengaproksimasi Hessian.
  • Hasil opt_res$par sangat dekat dengan beta_glm, Newton, dan Fisher → kembali menegaskan konsistensi antar-metode optimasi.

1.8 PERBANDINGAN HASIL

cat("=== PERBANDINGAN ===
")
## === PERBANDINGAN ===
cmp <- cbind(
  GLM      = beta_glm,
  Newton   = nr_res$beta,
  Fisher   = as.vector(fs_res$beta),
  BFGS     = opt_res$par
)
print(cmp)
##                       GLM        Newton        Fisher          BFGS
## (Intercept)  2.544493e+00  2.544493e+00  2.544493e+00  2.544493e+00
## IPM          2.118425e-02  2.118425e-02  2.118425e-02  2.118425e-02
## LPP          4.999287e-02  4.999287e-02  4.999287e-02  4.999287e-02
## PDRB         2.623636e-09  2.623630e-09  2.623630e-09  2.623655e-09
## TPT         -1.259747e-02 -1.259747e-02 -1.259747e-02 -1.259747e-02
cat("
")

Interpretasi Tahap 7 – Perbandingan Estimator:

  • Tabel cmp memperlihatkan bahwa:
    • GLM, Newton, Fisher, dan BFGS memberi estimasi parameter yang hampir identik (perbedaan hanya di digit desimal sangat kecil).
  • Artinya:
    • Implementasi log-likelihood, gradient, dan Hessian sudah benar.
    • Semua algoritma optimasi konsisten dan saling mengkonfirmasi.

1.9 BOOTSTRAPPING (50–200 Sampel)

bootstrap_gamma <- function(X, y, phi, B = 100){
  out <- matrix(NA, B, p)
  for(b in 1:B){
    idx <- sample(1:n, n, replace = TRUE)
    est <- newton_gamma(X[idx, ], y[idx], phi, beta_init = beta_glm)
    out[b, ] <- est$beta
  }
  colnames(out) <- names(beta_glm)
  return(out)
}

set.seed(123)
boot <- bootstrap_gamma(X, y, phi_glm, B = 100)

cat("=== HASIL BOOTSTRAP (RINGKASAN SD PARAMETER) ===
")
## === HASIL BOOTSTRAP (RINGKASAN SD PARAMETER) ===
print(apply(boot, 2, sd))
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 0.0358356985 2.3582237855 0.0511848321 0.0002872976 0.1230500329
cat("
Bootstrapping selesai.
")
## 
## Bootstrapping selesai.

Interpretasi Tahap 8 – Bootstrapping:

  • Setiap bootstrap \(b = 1,\dots,B\) mengambil sampel berukuran n dari data asli dengan pengembalian, lalu mengestimasi \(\beta\) dengan Newton.
  • Matriks boot menyimpan seluruh \(\hat\beta^{(b)}\).
  • apply(..., sd) menghitung standar deviasi bootstrap tiap parameter, yang dapat dipakai sebagai pendekatan standard error.
  • SD kecil → estimator stabil terhadap perubahan sampel;
    SD besar → estimator sensitif terhadap variasi sampel.

2 FUNCTION 2 – SIMULASI MONTE CARLO

2.1 PARAMETER MLE DARI FUNCTION 1

beta_mle <- coef(fit_glm)
phi_mle  <- summary(fit_glm)$dispersion

Interpretasi Tahap 1 – Parameter Acuan:

  • beta_mle dan phi_mle dianggap sebagai “parameter benar” (true parameter) untuk simulasi Monte Carlo.
  • Data sintetis akan disimulasikan dari model Gamma dengan parameter ini, lalu dilihat apakah estimator mampu “kembali” ke nilai tersebut (bias kecil, SE mengecil saat n membesar).

2.2 FUNGSI SIMULASI DATA DARI DISTRIBUSI GAMMA

simulate_gamma <- function(X, beta, phi){
  eta <- X %*% beta
  eta <- pmin(pmax(eta, -30), 30)  # stabilisasi numerik
  mu  <- exp(eta)

  shape <- 1 / phi
  scale <- phi * mu

  rgamma(nrow(X), shape = shape, scale = scale)
}

Interpretasi Tahap 2 – Simulasi Data:

  • eta dan mu dihitung sebagaimana pada model Gamma dengan link log.
  • Parametrisasi Gamma di R:
    • shape = 1/phi
    • scale = phi * mu
  • rgamma(...) menghasilkan data y sintetis yang mengikuti model Gamma dengan parameter beta dan phi yang ditentukan.

2.3 FUNGSI MONTE CARLO (B = 3000)

monte_carlo <- function(N, B = 3000){

  bhat <- matrix(NA, B, p)

  for(b in 1:B){
    idx <- sample(1:n, N, replace = TRUE)
    Xb  <- X[idx, ]
    yb  <- simulate_gamma(Xb, beta_mle, phi_mle)

    est <- newton_gamma(Xb, yb, phi_mle, beta_init = beta_mle)
    bhat[b, ] <- est$beta
  }

  colnames(bhat) <- names(beta_mle)

  mean_hat <- colMeans(bhat)
  bias     <- mean_hat - beta_mle
  se       <- apply(bhat, 2, sd)
  mse      <- bias^2 + se^2

  list(est = bhat, mean = mean_hat, bias = bias, se = se, mse = mse)
}

Interpretasi Tahap 3 – Monte Carlo:

  • Untuk tiap simulasi:
    • Dipilih N baris dari X (dengan pengembalian).
    • Dibangkitkan yb dari model Gamma dengan parameter “benar” beta_mle dan phi_mle.
    • Diestimasi kembali \(\beta\) dengan Newton.
  • Statistik Monte Carlo:
    • mean_hat: rata-rata estimator.
    • bias: mean_hat - beta_mle.
    • se: simpangan baku estimator antar-simulasi.
    • mse = bias^2 + se^2: ukuran kualitas estimator (mengkombinasikan bias dan varians).

2.4 SIMULASI UNTUK n = 30, 50, 100, 500

n_vec     <- c(30, 50, 100, 500)
mc_result <- list()

set.seed(123)
for(N in n_vec){
  cat("\n=== SIMULASI UNTUK n =", N, "===\n")
  mc_result[[as.character(N)]] <- monte_carlo(N, B = 3000)

  cat("Bias:\n")
  print(mc_result[[as.character(N)]]$bias)
  cat("SE:\n")
  print(mc_result[[as.character(N)]]$se)
  cat("MSE:\n")
  print(mc_result[[as.character(N)]]$mse)
}
## 
## === SIMULASI UNTUK n = 30 ===
## Bias:
##   (Intercept)           IPM           LPP          PDRB           TPT 
##   51847629.54 1774089300.04  116842836.63     -31362.56   76405461.97 
## SE:
## (Intercept)         IPM         LPP        PDRB         TPT 
##  2077000918 68700456553  4955388802     1221009  4638730823 
## MSE:
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 4.316621e+18 4.722900e+21 2.456953e+19 1.491846e+12 2.152366e+19 
## 
## === SIMULASI UNTUK n = 50 ===
## Bias:
##   (Intercept)           IPM           LPP          PDRB           TPT 
## -2.154624e-02 -1.459478e+00 -3.101864e-02 -4.419033e-05 -8.420054e-02 
## SE:
## (Intercept)         IPM         LPP        PDRB         TPT 
## 0.056594271 3.723228949 0.082643169 0.001560027 0.218792522 
## MSE:
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 3.667152e-03 1.599251e+01 7.792049e-03 2.435638e-06 5.495990e-02 
## 
## === SIMULASI UNTUK n = 100 ===
## Bias:
##   (Intercept)           IPM           LPP          PDRB           TPT 
## -1.158744e-02 -7.905591e-01 -1.653789e-02 -4.324259e-05 -4.433430e-02 
## SE:
## (Intercept)         IPM         LPP        PDRB         TPT 
##  0.03660044  2.50109803  0.05197184  0.00016144  0.13949271 
## MSE:
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 1.473861e-03 6.880475e+00 2.974574e-03 2.793278e-08 2.142375e-02 
## 
## === SIMULASI UNTUK n = 500 ===
## Bias:
##   (Intercept)           IPM           LPP          PDRB           TPT 
## -2.294178e-04 -1.556486e-02 -3.283242e-04 -1.046110e-06 -8.475388e-04 
## SE:
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 5.160462e-03 3.499583e-01 7.382141e-03 2.316225e-05 1.904729e-02 
## MSE:
##  (Intercept)          IPM          LPP         PDRB          TPT 
## 2.668300e-05 1.227131e-01 5.460381e-05 5.375843e-10 3.635177e-04

Interpretasi Tahap 4 – Perbandingan Ukuran Sampel:

  • n_vec mendefinisikan empat skenario ukuran sampel: 30, 50, 100, 500.
  • Untuk tiap N:
    • Dilakukan 3000 simulasi.
    • Dicetak Bias, SE, dan MSE.
  • Pola umum:
    • Untuk n = 30, bias dan SE cenderung besar → estimator belum stabil.
    • Mulai n = 50, 100, 500, Bias dan SE mengecil.
    • Semakin besar n, MSE cenderung menurun → menunjukkan konsistensi estimator.

2.5 RINGKASAN BIAS, SE, MSE DALAM TABEL

summ_list <- list()

for(N in n_vec){
  res <- mc_result[[as.character(N)]]
  tmp <- data.frame(
    n     = N,
    param = names(beta_mle),
    bias  = as.numeric(res$bias),
    se    = as.numeric(res$se),
    mse   = as.numeric(res$mse)
  )
  summ_list[[as.character(N)]] <- tmp
}

summary_df <- dplyr::bind_rows(summ_list)
print(summary_df)
##      n       param          bias           se          mse
## 1   30 (Intercept)  5.184763e+07 2.077001e+09 4.316621e+18
## 2   30         IPM  1.774089e+09 6.870046e+10 4.722900e+21
## 3   30         LPP  1.168428e+08 4.955389e+09 2.456953e+19
## 4   30        PDRB -3.136256e+04 1.221009e+06 1.491846e+12
## 5   30         TPT  7.640546e+07 4.638731e+09 2.152366e+19
## 6   50 (Intercept) -2.154624e-02 5.659427e-02 3.667152e-03
## 7   50         IPM -1.459478e+00 3.723229e+00 1.599251e+01
## 8   50         LPP -3.101864e-02 8.264317e-02 7.792049e-03
## 9   50        PDRB -4.419033e-05 1.560027e-03 2.435638e-06
## 10  50         TPT -8.420054e-02 2.187925e-01 5.495990e-02
## 11 100 (Intercept) -1.158744e-02 3.660044e-02 1.473861e-03
## 12 100         IPM -7.905591e-01 2.501098e+00 6.880475e+00
## 13 100         LPP -1.653789e-02 5.197184e-02 2.974574e-03
## 14 100        PDRB -4.324259e-05 1.614400e-04 2.793278e-08
## 15 100         TPT -4.433430e-02 1.394927e-01 2.142375e-02
## 16 500 (Intercept) -2.294178e-04 5.160462e-03 2.668300e-05
## 17 500         IPM -1.556486e-02 3.499583e-01 1.227131e-01
## 18 500         LPP -3.283242e-04 7.382141e-03 5.460381e-05
## 19 500        PDRB -1.046110e-06 2.316225e-05 5.375843e-10
## 20 500         TPT -8.475388e-04 1.904729e-02 3.635177e-04

Interpretasi Tahap 5 – Tabel Ringkasan:

  • summary_df merangkum semua hasil Monte Carlo dalam satu tabel dengan kolom:
    • n : ukuran sampel.
    • param : nama parameter (\(\beta_0, \beta_{IPM}, \beta_{LPP}, \beta_{PDRB}, \beta_{TPT}\))
    • bias, se, mse.
  • Dari tabel ini bisa dilihat pola bahwa:
    • Bias mengecil seiring bertambahnya n.
    • SE menurun.
    • MSE menurun → estimator semakin baik untuk sampel besar.

2.6 HISTOGRAM DISTRIBUSI SAMPLING (PARAMETER IPM)

param_name <- "IPM"

plot_list <- list()
for(N in n_vec){
  est_mat <- mc_result[[as.character(N)]]$est
  df_plot <- data.frame(
    est = est_mat[, param_name],
    n   = as.factor(N)
  )
  plot_list[[as.character(N)]] <- df_plot
}

df_plot_all <- dplyr::bind_rows(plot_list)

ggplot(df_plot_all, aes(x = est)) +
  geom_histogram(bins = 40, alpha = 0.6) +
  geom_vline(xintercept = beta_mle[param_name], linetype = "dashed") +
  facet_wrap(~ n, scales = "free") +
  labs(
    title = paste("Distribusi Sampling Estimator untuk Parameter", param_name),
    x     = paste0("beta_hat_", param_name),
    y     = "Frekuensi"
  )

Interpretasi Tahap 6 – Histogram Estimator \(\beta_{IPM}\):

  • Histogram menunjukkan distribusi sampling estimator \(\beta_{IPM}\) untuk berbagai ukuran sampel (n = 30, 50, 100, 500).
  • Garis putus-putus adalah nilai “benar” beta_mle["IPM"].
  • Untuk n kecil, histogram lebih lebar dan bisa bergeser dari nilai benar.
  • Untuk n besar, histogram menyempit dan terpusat di sekitar nilai benar → estimator makin presisi.

3 FUNCTION 3 – MCMC GAMMA REGRESSION

3.1 Paket MCMC

Bagian ini menggunakan paket:

  • mvtnorm: dmvnorm() untuk prior multivariat normal.
  • MASS: mvrnorm() untuk proposal \(\beta\).
  • coda: objek MCMC dan fungsi diagnostik (traceplot, densplot, cumuplot, acf).
  • stableGR: Gelman–Rubin R-hat stabil.

3.2 PRIOR (Multivariate Normal untuk \(\beta\), Inverse-Gamma untuk \(\phi\))

logprior <- function(beta, logphi){
  # prior beta ~ N_p(0, τ^2 I)
  tau2 <- 100
  lp_beta <- mvtnorm::dmvnorm(
    beta,
    mean  = rep(0, p),
    sigma = diag(tau2, p),
    log   = TRUE
  )

  # prior phi ~ IG(a,b) → logphi prior memakai Jacobian
  a <- 2; b <- 2
  phi <- exp(logphi)
  lp_phi <- (a * log(b) - lgamma(a)) + (-a - 1)*logphi - b/phi + logphi

  lp_beta + lp_phi
}

Interpretasi Tahap 2 – Prior:

  • Prior untuk \(\beta\):
    • \(\beta \sim N_p(0, \tau^2 I)\) dengan \(\tau^2 = 100\) → prior lemah (varians besar), tidak terlalu mengikat.
    • Digunakan multivariate normal melalui dmvnorm(), sehingga korelasi awal antar-koefisien diambil nol (kovarian diagonal).
  • Prior untuk \(\phi\):
    • \(\phi \sim \text{Inverse-Gamma}(a,b)\) dengan \(a = b = 2\).
    • Bekerja di skala logphi supaya \(\phi > 0\) terjamin.
    • Term tambahan + logphi adalah Jacobian transformasi dari \(\phi\) ke logphi.

3.3 LOG-POSTERIOR

# --- LOG-LIKELIHOOD ---
loglik <- function(beta, logphi, X, y){
  phi <- exp(logphi)
  eta <- X %*% beta
  eta <- pmin(pmax(eta, -30), 30)
  mu  <- exp(eta)

  shape <- 1/phi
  scale <- phi * mu

  sum(dgamma(y, shape = shape, scale = scale, log = TRUE))
}

# --- LOG-POSTERIOR ---
logpost <- function(beta, logphi, X, y){
  loglik(beta, logphi, X, y) + logprior(beta, logphi)
}

Interpretasi Tahap 3 – Likelihood & Posterior:

  • loglik: log-likelihood Gamma dengan parametrisasi shape dan scale yang konsisten dengan fungsi simulasi.
  • Stabilisasi numerik dengan memotong eta di [-30, 30] menghindari overflow/underflow.
  • logpost: menjumlahkan log-likelihood dan log-prior → menghasilkan log-posterior \[ \log p(\beta, \phi \mid y) = \log L(y \mid \beta, \phi) + \log p(\beta, \phi) \] yang menjadi target algoritma Metropolis–Hastings.

3.4 METROPOLIS–HASTINGS SAMPLER

mcmc_chain <- function(X, y,
                       iter = 20000,
                       beta_init,
                       phi_init,
                       sigma_beta = 0.0001,
                       sigma_phi  = 0.05){

  beta   <- beta_init
  logphi <- log(phi_init)

  BETA   <- matrix(NA, iter, length(beta))
  LOGPHI <- numeric(iter)

  lp <- logpost(beta, logphi, X, y)

  for(i in 1:iter){

    # ---- UPDATE β (Random-Walk MVN) ----
    beta_prop <- as.vector(MASS::mvrnorm(
      1, beta,
      Sigma = diag(rep(sigma_beta, length(beta)))
    ))

    lp_prop <- logpost(beta_prop, logphi, X, y)

    acc <- lp_prop - lp
    if(is.na(acc) || !is.finite(acc)) acc <- -Inf

    if(log(runif(1)) < acc){
      beta <- beta_prop
      lp   <- lp_prop
    }

    # ---- UPDATE log(φ) (Random-Walk Normal) ----
    logphi_prop <- rnorm(1, logphi, sigma_phi)
    lp_prop <- logpost(beta, logphi_prop, X, y)

    acc <- lp_prop - lp
    if(is.na(acc) || !is.finite(acc)) acc <- -Inf

    if(log(runif(1)) < acc){
      logphi <- logphi_prop
      lp     <- lp_prop
    }

    BETA[i, ]  <- beta
    LOGPHI[i]  <- logphi
  }

  list(beta = BETA, logphi = LOGPHI)
}

Interpretasi Tahap 4 – Algoritma MH:

  • Algoritma melakukan Metropolis–Hastings random-walk:
    • Untuk \(\beta\): proposal Normal multivariat \(N(\beta, \sigma_\beta^2 I)\).
    • Untuk logphi: proposal Normal 1-dimensi \(N(\log\phi, \sigma_\phi^2)\).
  • Setiap langkah:
    1. Hitung log-posterior pada proposal (lp_prop).
    2. Hitung selisih acc = lp_prop - lp.
    3. Jika log(runif(1)) < acc → terima proposal, jika tidak → tolak (tetap pada nilai lama).
  • Nilai yang disimpan:
    • BETA[i, ]: sampel posterior \(\beta\) ke-i.
    • LOGPHI[i]: sampel posterior logphi ke-i (nanti dikonversi ke \(\phi\)).

3.5 GENERATE ≥ 2 MCMC CHAINS (20.000 Iterasi)

set.seed(100)
chain1 <- mcmc_chain(X, y, iter = 20000,
                     beta_init = rep(0, p),
                     phi_init  = 1)

set.seed(200)
chain2 <- mcmc_chain(X, y, iter = 20000,
                     beta_init = rnorm(p, 0, 0.1),
                     phi_init  = 0.5)

Interpretasi Tahap 5 – Multi Chain:

  • Dibuat dua rantai dengan titik awal berbeda:
    • Chain 1: \(\beta\) mulai dari nol, \(\phi\) dari 1.
    • Chain 2: \(\beta\) mulai dari nilai acak dekat nol, \(\phi\) dari 0.5.
  • Tujuan:
    • Memeriksa apakah rantai yang mulai dari titik berbeda akhirnya konvergen ke distribusi posterior yang sama.
    • Ini penting untuk evaluasi konvergensi dengan Gelman–Rubin R-hat.

3.6 KONVERSI KE coda::mcmc.list

mcmc1 <- coda::mcmc(cbind(chain1$beta, phi = exp(chain1$logphi)))
mcmc2 <- coda::mcmc(cbind(chain2$beta, phi = exp(chain2$logphi)))

mcmc_list <- coda::mcmc.list(mcmc1, mcmc2)

Interpretasi Tahap 6 – Struktur MCMC:

  • mcmc() mengubah matriks sampel menjadi objek MCMC.
  • mcmc.list menggabungkan beberapa rantai (chain1 dan chain2) ke dalam satu objek.
  • Format ini diperlukan oleh fungsi-fungsi diagnostik coda:
    • traceplot, densplot, cumuplot, acf, dan perhitungan R-hat.

3.7 TRACEPLOT

par(mfrow = c(3, 2), mar = c(3, 3, 2, 1))

traceplot(mcmc_list[, 1], main = "Trace β0")
traceplot(mcmc_list[, 2], main = "Trace β_IPM")
traceplot(mcmc_list[, 3], main = "Trace β_LPP")
traceplot(mcmc_list[, 4], main = "Trace β_PDRB")
traceplot(mcmc_list[, 5], main = "Trace β_TPT")
traceplot(mcmc_list[, 6], main = "Trace φ")

Interpretasi Tahap 7 – Traceplot:

  • Traceplot menampilkan evolusi nilai parameter sepanjang iterasi MCMC untuk kedua rantai.
  • Konvergensi yang baik terlihat saat:
    • Rantai bergerak naik-turun seperti “ulat gemuk” (tidak datar, tidak trending).
    • Tidak ada tren naik atau turun yang jelas.
    • Kedua rantai untuk satu parameter saling overlap dalam rentang yang sama.
  • Jika pola ini muncul, rantai dianggap sudah mencapai distribusi stasioner.

3.8 DENSITY PLOT

par(mfrow = c(3, 2), mar = c(3, 3, 2, 1))

densplot(mcmc_list[, 1], main = "Density β0")
densplot(mcmc_list[, 2], main = "Density β_IPM")
densplot(mcmc_list[, 3], main = "Density β_LPP")
densplot(mcmc_list[, 4], main = "Density β_PDRB")
densplot(mcmc_list[, 5], main = "Density β_TPT")
densplot(mcmc_list[, 6], main = "Density φ")

Interpretasi Tahap 8 – Density Posterior:

  • Density plot menggambarkan bentuk distribusi posterior parameter.
  • Hal yang dicek:
    • Apakah distribusinya tampak wajar (misal unimodal, tidak terlalu tajam/aneh).
    • Apakah kedua rantai memberikan bentuk densitas yang mirip.
  • Dari sini bisa diambil ringkasan posterior (mean, median, credible interval) jika diperlukan.

3.9 CUMULATIVE MEAN PLOT

par(mfrow = c(3, 2), mar = c(3, 3, 2, 1))

cumuplot(mcmc_list[, 1], main = "Mean β0")

cumuplot(mcmc_list[, 2], main = "Mean β_IPM")

cumuplot(mcmc_list[, 3], main = "Mean β_LPP")

cumuplot(mcmc_list[, 4], main = "Mean β_PDRB")

cumuplot(mcmc_list[, 5], main = "Mean β_TPT")

cumuplot(mcmc_list[, 6], main = "Mean φ")

Interpretasi Tahap 9 – Cumulative Mean:

  • Cumulative mean plot menampilkan rata-rata kumulatif parameter terhadap jumlah iterasi.
  • Pada awal iterasi, garis cumulative mean biasanya berfluktuasi tajam.
  • Jika rantai sudah konvergen, garis akan mendatar (stabil) setelah burn-in.
  • Ini menunjukkan bahwa tambahan iterasi tidak banyak mengubah mean → rantai sudah sampling dari distribusi target.

3.10 AUTOCORRELATION (ACF)

par(mfrow = c(3, 2), mar = c(3, 3, 4, 1))

acf(mcmc_list[[1]][, 1], main = "ACF β0")
acf(mcmc_list[[1]][, 2], main = "ACF β_IPM")
acf(mcmc_list[[1]][, 3], main = "ACF β_LPP")
acf(mcmc_list[[1]][, 4], main = "ACF β_PDRB")
acf(mcmc_list[[1]][, 5], main = "ACF β_TPT")
acf(mcmc_list[[1]][, 6], main = "ACF φ")

Interpretasi Tahap 10 – Autocorrelation:

  • Plot ACF menunjukkan korelasi antar-sampel yang dipisah beberapa lag.
  • Idealnya:
    • Autokorelasi turun cepat ke sekitar nol saat lag meningkat.
    • Jika autokorelasi tinggi untuk banyak lag, rantai “lengket” dan efektif sampelnya kecil → mungkin perlu thinning atau tuning proposal (sigma_beta, sigma_phi).

3.11 GELMAN–RUBIN R-hat

burn <- 5000
thin <- 5

chain1_post <- mcmc1[(burn + 1):20000, ][seq(1, 20000 - burn, by = thin), ]
chain2_post <- mcmc2[(burn + 1):20000, ][seq(1, 20000 - burn, by = thin), ]

# tambahkan jitter kecil jika varians terlalu kecil (hindari masalah numerik)
if(any(apply(chain1_post, 2, var) < 1e-10)){
  chain1_post <- chain1_post + matrix(
    rnorm(length(chain1_post), 0, 1e-6),
    nrow = nrow(chain1_post)
  )
}
if(any(apply(chain2_post, 2, var) < 1e-10)){
  chain2_post <- chain2_post + matrix(
    rnorm(length(chain2_post), 0, 1e-6),
    nrow = nrow(chain2_post)
  )
}

mcmc_list_fixed <- coda::mcmc.list(
  coda::mcmc(chain1_post),
  coda::mcmc(chain2_post)
)

cat("=== Gelman–Rubin R-hat (Stabil) ===\n")
## === Gelman–Rubin R-hat (Stabil) ===
print(stableGR::stable.GR(mcmc_list_fixed))
## $psrf
## [1] 1.012282 1.012269 1.012281 1.001400 1.012223 1.009157
## 
## $mpsrf
## [1] 1.001159
## 
## $means
##                                                                       
##  9.433104e-05  6.126312e-02 -1.735143e-02  2.820303e-06 -1.850231e-02 
##           phi 
##  2.818256e+00 
## 
## $n.eff
## [1] 751.4501
## 
## $blather
## [1] FALSE

Interpretasi Tahap 11 – Gelman–Rubin R-hat:

  • Dilakukan burn-in (buang 5000 iterasi awal) dan thinning (ambil tiap 5 iterasi).
  • Jika varians beberapa parameter terlalu kecil, ditambahkan jitter kecil untuk menghindari masalah numerik saat menghitung R-hat.
  • stable.GR() menghitung nilai R-hat stabil:
    • Nilai R-hat mendekati 1 (misal < 1.1) → menunjukkan konvergensi antar-rantai yang baik.
    • Nilai R-hat jauh di atas 1 → rantai belum benar-benar menyatu, perlu iterasi tambahan atau tuning proposal.

4 KESIMPULAN

Function 1 (Optimasi)
Menunjukkan bahwa GLM, Newton–Raphson, Fisher Scoring, dan BFGS memberi estimasi parameter yang konsisten. Bootstrapping memberikan gambaran ketidakpastian estimator (standard error bootstrap).

Function 2 (Simulasi Monte Carlo)
Membuktikan sifat estimator MLE Gamma:

  • Untuk n kecil (misal 30), bias dan SE besar → hasil tidak stabil.
  • Untuk n menengah–besar (n ≥ 50), bias kecil dan SE mengecil.
  • MSE turun seiring pertambahan n → estimator konsisten.

Function 3 (MCMC – Metropolis–Hastings)
Memberikan pendekatan Bayesian untuk model regresi Gamma:

  • Prior: \(\beta \sim N_p(0, 100 I)\), \(\phi \sim IG(2,2)\).
  • Posterior disample dengan MH random-walk untuk \(\beta\) dan logphi.
  • Evaluasi konvergensi:
    • Traceplot (pola “ulat gemuk” dan overlap antar-rantai),
    • Density plot (bentuk posterior yang wajar),
    • Cumulative mean plot (mean kumulatif stabil setelah burn-in),
    • ACF (autokorelasi turun cepat),
    • Gelman–Rubin R-hat (mendekati 1).

Setelah konvergen, posterior mean dan credible interval dapat digunakan sebagai dasar inferensi probabilistik untuk parameter regresi Gamma.