Occupational Safety Incidents on Construction Projects: A Negative Binomial Model for Overdispersed Counts

# Contoh kode R: Negative Binomial Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(MASS)
})
set.seed(42)
# ---- Gaya grafik ----
pal <- c("#1f6feb", "#e8590c", "#2f9e44", "#ae3ec9", "#f59f00")
warna <- list(bg = "#f4f6fa", panel = "#ffffff", teks = "#1f2933", redup = "#5f6b78", grid = "#e3e8ee", bar = "#1d3557")
tata <- 1
theme_set(theme_minimal(base_size = 11) + theme(legend.position = "bottom"))

# ---- Data sintetis ----
n <- 200
d <- data.frame(total_work_hours = rnorm(n, 50, 10),
                safety_training_hours = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$total_work_hours))
z2 <- as.numeric(scale(d$safety_training_hours))
mu0 <- runif(1, 3, 8)                                       # rata-rata hitungan dasar
b1 <- runif(1, 0.3, 0.55)
b2 <- runif(1, 0.1, 0.3) * sample(c(-1, 1), 1)
theta0 <- runif(1, 0.9, 2.5)                                # parameter dispersi: makin kecil makin overdispersi
d$safety_incidents <- rnbinom(n, size = theta0, mu = mu0 * exp(b1 * z1 + b2 * z2))
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: Poisson sebagai pembanding, lalu binomial negatif (MASS::glm.nb) ----
fit_p <- glm(safety_incidents ~ total_work_hours + safety_training_hours, family = poisson, data = d)
disp_p <- sum(residuals(fit_p, type = "pearson")^2) / fit_p$df.residual   # Pearson chi2/df model Poisson
cat(sprintf("Dispersi Poisson (Pearson chi2/df) = %.3f\n", disp_p))
fit <- MASS::glm.nb(safety_incidents ~ total_work_hours + safety_training_hours, data = d)
s <- summary(fit)
print(s)
irr_tab <- exp(cbind(IRR = coef(fit), confint.default(fit)))   # IRR per satuan + CI Wald 95%
print(round(irr_tab, 4))

# Uji LR overdispersi: H0 theta tak hingga (Poisson); p dibagi dua karena di batas ruang parameter
lr_stat <- as.numeric(2 * (logLik(fit) - logLik(fit_p)))
p_od <- 0.5 * pchisq(max(lr_stat, 0), df = 1, lower.tail = FALSE)
aic_tab <- c(Poisson = AIC(fit_p), NegBin = AIC(fit))
print(aic_tab)
cat(sprintf("LR overdispersi = %.2f, p = %.4g; theta = %.3f (SE %.3f)\n", lr_stat, p_od, fit$theta, fit$SE.theta))

# IRR per 1 simpangan baku dari kedua model (untuk forest plot)
sdx <- c(sd(d$total_work_hours), sd(d$safety_training_hours))
irr_sd <- function(m, lab) {
  sm <- summary(m)$coefficients
  data.frame(var = c("total_work_hours", "safety_training_hours"), model = lab,
             irr = exp(sm[2:3, 1] * sdx), lo = exp((sm[2:3, 1] - 1.96 * sm[2:3, 2]) * sdx),
             hi = exp((sm[2:3, 1] + 1.96 * sm[2:3, 2]) * sdx), p = sm[2:3, 4])
}
ir <- rbind(irr_sd(fit_p, "Poisson"), irr_sd(fit, "Negative binomial"))
print(ir)
ir_nb <- ir[ir$model == "Negative binomial", ]


# ---- Grafik ----
# Distribusi hitungan teramati vs prediksi Poisson dan binomial negatif
mu_p <- fitted(fit_p); mu_nb <- fitted(fit)
kmax <- max(8, as.numeric(quantile(d$safety_incidents, 0.97)))
k <- 0:kmax
dk <- data.frame(k = k, teramati = as.numeric(table(factor(pmin(d$safety_incidents, kmax + 1), levels = k))))
dl <- rbind(data.frame(k = k, model = "Poisson", f = sapply(k, function(i) sum(dpois(i, mu_p)))),
            data.frame(k = k, model = "Negative binomial", f = sapply(k, function(i) sum(dnbinom(i, size = fit$theta, mu = mu_nb)))))
p1 <- ggplot() +
  geom_col(data = dk, aes(k, teramati), fill = pal[1], alpha = 0.55, width = 0.85) +
  geom_line(data = dl, aes(k, f, colour = model), linewidth = 1.1) +
  geom_point(data = dl, aes(k, f, colour = model), size = 1.6) +
  scale_colour_manual(values = c(Poisson = pal[2], "Negative binomial" = pal[3]), name = NULL) +
  labs(title = "Observed counts (bars) vs model predictions", x = "safety_incidents", y = "Frequency")
p2 <- ggplot(ir, aes(irr, var, colour = model)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = warna$redup) +
  geom_errorbar(aes(xmin = lo, xmax = hi), width = 0.2, linewidth = 0.9,
                position = position_dodge(width = 0.5), orientation = "y") +
  geom_point(size = 2.8, position = position_dodge(width = 0.5)) +
  scale_colour_manual(values = c(Poisson = pal[2], "Negative binomial" = pal[3]), name = NULL) +
  scale_x_log10() +
  labs(title = "IRR per 1 SD: NB gives wider, honest CIs", x = "IRR (log scale)", y = NULL) +
  theme(legend.position = "none")
dres <- rbind(data.frame(model = "Poisson", fitted = mu_p, r = residuals(fit_p, type = "pearson")),
              data.frame(model = "Negative binomial", fitted = mu_nb, r = residuals(fit, type = "pearson")))
dres$model <- factor(dres$model, levels = c("Poisson", "Negative binomial"))
p3 <- ggplot(dres, aes(fitted, r, colour = model)) +
  geom_hline(yintercept = 0, colour = warna$redup) +
  geom_hline(yintercept = c(-2, 2), linetype = "dashed", colour = warna$redup) +
  geom_point(alpha = 0.45, size = 1.3) +
  facet_wrap(~model) +
  scale_colour_manual(values = c(Poisson = pal[2], "Negative binomial" = pal[3]), guide = "none") +
  scale_x_log10() +
  labs(title = "Pearson residuals vs fitted", x = "Fitted mean (log scale)", y = "Pearson residual")
p_hasil <- if (tata == 1) p1 / (p2 | p3) else (p1 / p2) | p3

print(p_hasil)