poisson gamma bayesian regresi

###############################################
# 4. MODEL BAYESIAN POISSON–GAMMA
###############################################
model_string <- "
model {

  for(i in 1:N){

    # Likelihood: Y_i | lambda_i
    Y[i] ~ dpois(lambda[i])

    # Gamma random effect: lambda_i | theta_i, mu_i
    lambda[i] ~ dgamma(theta[i], theta[i] / mu[i])

    # theta_i ~ Gamma(alpha, beta)
    theta[i] ~ dgamma(alpha, beta)

    # mu_i dengan offset Premium
    mu[i] <- exp(logPremium[i] + a * X1[i] + b * X2[i])
  }

  # Priors untuk hyperparameter
  alpha ~ dgamma(0.001, 0.001)
  beta  ~ dgamma(0.001, 0.001)

  # Priors regresi
  a ~ dnorm(0, 1.0E-6)
  b ~ dnorm(0, 1.0E-6)

}
"

###############################################
# 5. INISIALISASI
###############################################
inits <- function() {
  list(
    alpha  = 1,
    beta   = 1,
    a      = 0,
    b      = 0,
    theta  = rep(1, nrow(data)),
    lambda = pmax(data$Claims, 1)
  )
}

params <- c("a", "b", "alpha", "beta")

###############################################
# 6. JALANKAN MODEL MCMC
###############################################
set.seed(123)

jm <- jags.model(
  textConnection(model_string),
  data = jags_data,
  inits = inits,
  n.chains = 3,
  n.adapt = 2000
)

## Compiling model graph
##    Resolving undeclared variables
##    Allocating nodes
## Graph information:
##    Observed stochastic nodes: 20
##    Unobserved stochastic nodes: 44
##    Total graph size: 197
## 
## Initializing model

update(jm, 5000)   # burn-in

###############################################
# 7. SAMPLING POSTERIOR
###############################################
samples <- coda.samples(
  jm,
  variable.names = params,
  n.iter = 10000
)

print(summary(samples))

## 
## Iterations = 7001:17000
## Thinning interval = 1 
## Number of chains = 3 
## Sample size per chain = 10000 
## 
## 1. Empirical mean and standard deviation for each variable,
##    plus standard error of the mean:
## 
##           Mean       SD  Naive SE Time-series SE
## a      0.01179  0.06198 0.0003579      8.316e-04
## alpha 73.27009 77.05409 0.4448720      1.191e+01
## b      0.27720  0.03653 0.0002109      5.056e-04
## beta   9.39253  9.96552 0.0575359      1.166e+00
## 
## 2. Quantiles for each variable:
## 
##          2.5%      25%      50%       75%    97.5%
## a     -0.1090 -0.02895  0.01136   0.05204   0.1361
## alpha  1.8627 14.95286 41.60658 109.59403 270.9615
## b      0.2070  0.25305  0.27640   0.30073   0.3508
## beta   0.1596  1.77102  5.06890  15.59300  34.6057

###############################################
# 8. SIGNIFIKANSI BAYESIAN
###############################################
# 95% HPD credible interval
print(HPDinterval(samples[, c("a", "b")]))

## [[1]]
##        lower     upper
## a -0.1144820 0.1252550
## b  0.2039914 0.3445428
## attr(,"Probability")
## [1] 0.95
## 
## [[2]]
##        lower     upper
## a -0.1099392 0.1403184
## b  0.2025838 0.3505642
## attr(,"Probability")
## [1] 0.95
## 
## [[3]]
##        lower     upper
## a -0.1110491 0.1304076
## b  0.2069691 0.3492461
## attr(,"Probability")
## [1] 0.95

# Probabilitas parameter positif
s_mat <- as.matrix(samples)
prob_a_pos <- mean(s_mat[, "a"] > 0)
prob_b_pos <- mean(s_mat[, "b"] > 0)

cat("\nProbabilitas a > 0:", prob_a_pos, "\n")

## 
## Probabilitas a > 0: 0.5769

cat("Probabilitas b > 0:", prob_b_pos, "\n")

## Probabilitas b > 0: 1

###############################################
# 9. TRACE PLOT
##############################################

# Tutup semua device grafik (jika ada)
graphics.off()

# Buka jendela grafik baru (Windows)
if (.Platform$OS.type == "windows") {
  windows(width = 12, height = 8)   # ukuran besar agar aman
}

# Atur margin lebih kecil (bawah, kiri, atas, kanan)
###############################################
# TRACE + DENSITY (semua parameter)
###############################################
par(mar = c(4,4,2,1), mfrow = c(1,1))
plot(samples)     # aman untuk knit
par(mar = c(4,4,2,1), mfrow = c(1,1))

traceplot(samples[, "a"],
          main = "Trace plot a (Merit)")
par(mar = c(4,4,2,1), mfrow = c(1,1))

traceplot(samples[, "b"],
          main = "Trace plot b (Class)")
par(mfrow = c(2,1), mar = c(4,4,2,1))

traceplot(samples[, "a"], main = "Trace plot a (Merit)")
traceplot(samples[, "b"], main = "Trace plot b (Class)")

###############################################
# 10. GELMAN–RUBIN (UJI KONVERGENSI)
###############################################
gelman_diag <- gelman.diag(samples, autoburnin = FALSE)
print(gelman_diag)

## Potential scale reduction factors:
## 
##       Point est. Upper C.I.
## a           1.00       1.01
## alpha       2.90       9.08
## b           1.00       1.00
## beta        2.66       6.92
## 
## Multivariate psrf
## 
## 2.3

gelman.plot(samples, autoburnin = FALSE)

###############################################
# 11. RINGKASAN LENGKAP PARAMETER
###############################################
sum_samp <- summary(samples)
sum_samp$statistics

##              Mean          SD     Naive SE Time-series SE
## a      0.01179071  0.06198409 0.0003578653   8.315901e-04
## alpha 73.27008611 77.05408849 0.4448719873   1.191383e+01
## b      0.27719897  0.03653486 0.0002109341   5.055631e-04
## beta   9.39253296  9.96551599 0.0575359334   1.166070e+00

sum_samp$quantiles

##             2.5%         25%         50%          75%       97.5%
## a     -0.1089664 -0.02895203  0.01135655   0.05204045   0.1360682
## alpha  1.8627303 14.95285604 41.60657852 109.59402887 270.9615269
## b      0.2069683  0.25304577  0.27640306   0.30072672   0.3508037
## beta   0.1596330  1.77102350  5.06889689  15.59299784  34.6056669