Poisson Modelling of Annual Employee Absence Days

# Contoh kode R: Poisson Regression
suppressPackageStartupMessages({
  library(ggplot2)
  library(patchwork)
  library(lmtest)
})
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(commute_distance = rnorm(n, 50, 10),
                monthly_overtime = rnorm(n, 30, 6))
z1 <- as.numeric(scale(d$commute_distance))
z2 <- as.numeric(scale(d$monthly_overtime))
mu0 <- runif(1, 2.5, 7)                                   # rata-rata hitungan dasar
b1 <- runif(1, 0.25, 0.5)
b2 <- runif(1, 0.1, 0.3) * sample(c(-1, 1), 1)
d$absence_days <- rpois(n, mu0 * exp(b1 * z1 + b2 * z2))         # hitungan Poisson (tanpa overdispersi)
write.csv(d, "data_contoh.csv", row.names = FALSE)

# ---- Analisis: regresi Poisson (glm, family poisson, link log) ----
fit <- glm(absence_days ~ commute_distance + monthly_overtime, family = poisson, data = d)
s <- summary(fit)
print(s)
irr_tab <- exp(cbind(IRR = coef(fit), confint.default(fit)))   # incidence rate ratio per satuan + CI 95%
print(round(irr_tab, 4))

# IRR per 1 simpangan baku prediktor (agar sebanding)
sdx <- c(sd(d$commute_distance), sd(d$monthly_overtime))
se_b <- s$coefficients[2:3, 2]
irr_sd <- data.frame(var = c("commute_distance", "monthly_overtime"),
                     irr = exp(coef(fit)[2:3] * sdx),
                     lo = exp((coef(fit)[2:3] - 1.96 * se_b) * sdx),
                     hi = exp((coef(fit)[2:3] + 1.96 * se_b) * sdx),
                     p = s$coefficients[2:3, 4])
print(irr_sd)

# ---- Uji model dan dispersi ----
fit0 <- glm(absence_days ~ 1, family = poisson, data = d)
lr <- lmtest::lrtest(fit0, fit)                          # uji likelihood ratio terhadap model null
print(lr)
p_lr <- lr[2, "Pr(>Chisq)"]
disp_dev <- fit$deviance / fit$df.residual               # deviance/df: dekat 1 berarti tidak overdispersi
disp_pear <- sum(residuals(fit, type = "pearson")^2) / fit$df.residual
p_gof <- pchisq(fit$deviance, fit$df.residual, lower.tail = FALSE)   # uji kecocokan deviance
r2_dev <- 1 - fit$deviance / fit$null.deviance
cat(sprintf("Deviance/df = %.3f; Pearson chi2/df = %.3f; GOF p = %.3f\n", disp_dev, disp_pear, p_gof))


# ---- Grafik ----
# Frekuensi teramati vs harapan model untuk tiap nilai hitungan
mu_hat <- fitted(fit)
k <- 0:max(d$absence_days)
dk <- data.frame(k = k,
                 teramati = as.numeric(table(factor(d$absence_days, levels = k))),
                 harapan = sapply(k, function(i) sum(dpois(i, mu_hat))))
p1 <- ggplot(dk, aes(k)) +
  geom_col(aes(y = teramati, fill = "Observed"), alpha = 0.75, width = 0.8) +
  geom_line(aes(y = harapan, colour = "Expected (Poisson)"), linewidth = 1) +
  geom_point(aes(y = harapan, colour = "Expected (Poisson)"), size = 1.8) +
  scale_fill_manual(values = c(Observed = pal[1]), name = NULL) +
  scale_colour_manual(values = c("Expected (Poisson)" = pal[2]), name = NULL) +
  labs(title = "Observed vs expected count frequencies", x = "absence_days", y = "Frequency")
# Laju prediksi menurut X1 (X2 pada rata-rata) dengan pita 95%
g <- data.frame(commute_distance = seq(min(d$commute_distance), max(d$commute_distance), length.out = 120), monthly_overtime = mean(d$monthly_overtime))
pr <- predict(fit, newdata = g, type = "link", se.fit = TRUE)
g$m <- exp(pr$fit); g$lo <- exp(pr$fit - 1.96 * pr$se.fit); g$hi <- exp(pr$fit + 1.96 * pr$se.fit)
p2 <- ggplot() +
  geom_jitter(data = d, aes(commute_distance, absence_days), width = 0, height = 0.15, colour = pal[3], alpha = 0.45, size = 1.4) +
  geom_ribbon(data = g, aes(commute_distance, ymin = lo, ymax = hi), fill = pal[2], alpha = 0.22) +
  geom_line(data = g, aes(commute_distance, m), colour = pal[2], linewidth = 1.2) +
  labs(title = "Predicted rate (95% CI)", x = "commute_distance", y = "absence_days")
p3 <- ggplot(irr_sd, aes(irr, var)) +
  geom_vline(xintercept = 1, linetype = "dashed", colour = warna$redup) +
  geom_segment(aes(x = lo, xend = hi, yend = var), colour = pal[1], linewidth = 1.1) +
  geom_point(colour = pal[2], size = 3.2) +
  geom_text(aes(label = sprintf("%.3f", irr)), vjust = -1.1, colour = warna$teks, size = 3.2) +
  scale_x_log10() +
  labs(title = "Incidence rate ratio per 1 SD (95% CI)", x = "IRR (log scale)", y = NULL)
p_hasil <- if (tata == 1) p1 / (p2 | p3) else (p1 / p3) | p2

print(p_hasil)