pkgs <- c("tidyverse","knitr","kableExtra","nnet","MASS","broom","lmtest","AER","readxl")
new_pkg <- pkgs[!pkgs %in% installed.packages()[,"Package"]]
if (length(new_pkg)) install.packages(new_pkg, repos="https://cloud.r-project.org", quiet=TRUE)

library(tidyverse)
library(knitr)
library(kableExtra)
library(nnet)
library(MASS)
library(broom)
library(lmtest)
library(AER)
library(readxl)

Regresi Logistik Biner, Multinomial, Ordinal & Poisson

Dataset: Absenteeism at Work — UCI Machine Learning Repository  |  01 June 2026

1 Deskripsi Data

1.1 Sumber Data

Dataset Absenteeism at Work berasal dari UCI Machine Learning Repository (ID: 445). Data mencatat ketidakhadiran kerja periode Juli 2007–Juli 2010 di perusahaan kurir Brazil. Terdiri dari 740 observasi dan 21 variabel.

Sitasi: Martiniano, A., Ferreira, R. P., Sassi, R. J., & Affonso, C. (2012). Application of a neuro fuzzy network in prediction of absenteeism at work. CISTI 2012. IEEE.

1.2 Memuat Data

# ── Ganti path sesuai lokasi file di komputer kamu ──
file_path <- file.choose()   # <-- GANTI INI

ext <- tolower(tools::file_ext(file_path))
if (ext %in% c("xls","xlsx")) {
  df_raw <- as.data.frame(read_excel(file_path, sheet = 1))
} else {
  df_raw <- read.csv2(file_path, sep = ";", stringsAsFactors = FALSE,
                      fileEncoding = "latin1")
}

cat("Dimensi:", nrow(df_raw), "x", ncol(df_raw), "\n")
#> Dimensi: 740 x 21
names(df_raw) <- trimws(names(df_raw))
cat("Kolom asli:\n"); print(names(df_raw))
#> Kolom asli:
#>  [1] "ID"                              "Reason for absence"             
#>  [3] "Month of absence"                "Day of the week"                
#>  [5] "Seasons"                         "Transportation expense"         
#>  [7] "Distance from Residence to Work" "Service time"                   
#>  [9] "Age"                             "Work load Average/day"          
#> [11] "Hit target"                      "Disciplinary failure"           
#> [13] "Education"                       "Son"                            
#> [15] "Social drinker"                  "Social smoker"                  
#> [17] "Pet"                             "Weight"                         
#> [19] "Height"                          "Body mass index"                
#> [21] "Absenteeism time in hours"

1.3 Pembersihan dan Persiapan Data

df <- df_raw

# Rename ke nama standar (sesuai nama kolom UCI)
names(df)[names(df) == "ID"]                               <- "ID"
names(df)[names(df) == "Reason for absence"]               <- "reason_absence"
names(df)[names(df) == "Month of absence"]                 <- "month_absence"
names(df)[names(df) == "Day of the week"]                  <- "day_week"
names(df)[names(df) == "Seasons"]                          <- "seasons"
names(df)[names(df) == "Transportation expense"]           <- "transport_expense"
names(df)[names(df) == "Distance from Residence to Work"]  <- "distance_residence"
names(df)[names(df) == "Service time"]                     <- "service_time"
names(df)[names(df) == "Age"]                              <- "age"
names(df)[grepl("Work load", names(df), ignore.case=TRUE)] <- "workload_avg"
names(df)[names(df) == "Hit target"]                       <- "hit_target"
names(df)[names(df) == "Disciplinary failure"]             <- "disciplinary"
names(df)[names(df) == "Education"]                        <- "education"
names(df)[names(df) == "Son"]                              <- "children"
names(df)[names(df) == "Social drinker"]                   <- "drinker"
names(df)[names(df) == "Social smoker"]                    <- "smoker"
names(df)[names(df) == "Pet"]                              <- "pet"
names(df)[names(df) == "Weight"]                           <- "weight"
names(df)[names(df) == "Height"]                           <- "height"
names(df)[names(df) == "Body mass index"]                  <- "bmi"
names(df)[names(df) == "Absenteeism time in hours"]        <- "absent_hours"

# Konversi numerik
for (col in names(df)) {
  df[[col]] <- suppressWarnings(as.numeric(as.character(df[[col]])))
}

# Filter valid
df <- df[!is.na(df$absent_hours) & df$absent_hours >= 0, ]

cat("Dimensi bersih:", nrow(df), "x", ncol(df), "\n")
#> Dimensi bersih: 740 x 21
cat("Kolom final:\n"); print(names(df))
#> Kolom final:
#>  [1] "ID"                 "reason_absence"     "month_absence"     
#>  [4] "day_week"           "seasons"            "transport_expense" 
#>  [7] "distance_residence" "service_time"       "age"               
#> [10] "workload_avg"       "hit_target"         "disciplinary"      
#> [13] "education"          "children"           "drinker"           
#> [16] "smoker"             "pet"                "weight"            
#> [19] "height"             "bmi"                "absent_hours"

1.4 Definisi Variabel

Tabel 1.1. Definisi Variabel
No Variabel Tipe
1 ID ID
2 Reason for Absence Nominal
3 Month of Absence Nominal
4 Day of Week Nominal
5 Seasons Nominal
6 Transportation Expense Kontinu
7 Distance to Work Kontinu
8 Service Time Kontinu
9 Age Kontinu
10 Work Load Average/Day Kontinu
11 Hit Target Kontinu
12 Disciplinary Failure Biner
13 Education Ordinal
14 Son (Children) Diskrit
15 Social Drinker Biner
16 Social Smoker Biner
17 Pet Diskrit
18 Weight Kontinu
19 Height Kontinu
20 Body Mass Index Kontinu
21 Absenteeism Time in Hours Count

1.5 Statistik Deskriptif

cont_vars <- c("transport_expense","distance_residence","service_time",
               "age","workload_avg","hit_target","weight","height","bmi","absent_hours")

desc_rows <- lapply(cont_vars, function(v) {
  x <- df[[v]]
  data.frame(
    Variabel = v,
    N    = sum(!is.na(x)),
    Mean = round(mean(x, na.rm=TRUE), 2),
    SD   = round(sd(x, na.rm=TRUE), 2),
    Min  = round(min(x, na.rm=TRUE), 2),
    Med  = round(median(x, na.rm=TRUE), 2),
    Max  = round(max(x, na.rm=TRUE), 2),
    stringsAsFactors = FALSE
  )
})
desc_tbl <- do.call(rbind, desc_rows)

kable(desc_tbl, caption="**Tabel 1.2. Statistik Deskriptif**",
      align=c("l",rep("c",6))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, background="#f0f4ff") %>%
  row_spec(10, background="#fef9c3", bold=TRUE)
Tabel 1.2. Statistik Deskriptif
Variabel N Mean SD Min Med Max
transport_expense 740 221.33 66.95 118 225 388
distance_residence 740 29.63 14.84 5 26 52
service_time 740 12.55 4.38 1 13 29
age 740 36.45 6.48 27 37 58
workload_avg 740 271490.24 39058.12 205917 264249 378884
hit_target 740 94.59 3.78 81 95 100
weight 740 79.04 12.88 56 83 108
height 740 172.11 6.03 163 170 196
bmi 740 26.68 4.29 19 25 38
absent_hours 740 6.92 13.33 0 3 120 
cat_df <- data.frame(
  Variabel = c("Disciplinary Failure = 1","Social Drinker = 1",
               "Social Smoker = 1","Education (median)","Children (mean)"),
  Nilai = c(
    sum(df$disciplinary == 1, na.rm=TRUE),
    sum(df$drinker == 1, na.rm=TRUE),
    sum(df$smoker == 1, na.rm=TRUE),
    median(df$education, na.rm=TRUE),
    round(mean(df$children, na.rm=TRUE), 2)
  ),
  stringsAsFactors = FALSE
)
kable(cat_df, caption="**Tabel 1.3. Ringkasan Variabel Kategori**",
      align=c("l","c")) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=12)
Tabel 1.3. Ringkasan Variabel Kategori
Variabel Nilai
Disciplinary Failure = 1 40.00
Social Drinker = 1 420.00
Social Smoker = 1 54.00
Education (median) 1.00
Children (mean) 1.02

1.6 Visualisasi Eksploratif

par(mfrow=c(2,3), mar=c(4,4,3,1.5), bg="#f8fafc", oma=c(0,0,3,0))

hist(df$absent_hours, breaks=30, col="#2563eb", border="white",
     main="Distribusi Absent Hours", xlab="Jam", ylab="Frekuensi",
     cex.axis=.82, las=1)
box(col="#cbd5e1")

boxplot(absent_hours ~ drinker, data=df, col=c("#2563eb","#dc2626"),
        names=c("Non-Drinker","Drinker"), main="Absent vs Social Drinker",
        ylab="Jam", border="#1b2a4a", cex.axis=.82)
box(col="#cbd5e1")

boxplot(absent_hours ~ education, data=df,
        col=c("#7c3aed","#2563eb","#0d9488","#16a34a"),
        names=c("SD","SMA","S1","S2/S3"), main="Absent vs Pendidikan",
        ylab="Jam", border="#1b2a4a", cex.axis=.82)
box(col="#cbd5e1")

plot(df$age, df$absent_hours, pch=19, cex=.5,
     col=adjustcolor("#2563eb",.4), main="Usia vs Absent Hours",
     xlab="Usia", ylab="Jam", cex.axis=.82, las=1)
abline(lm(absent_hours ~ age, data=df), col="#dc2626", lwd=2, lty=2)
box(col="#cbd5e1")

plot(df$bmi, df$absent_hours, pch=19, cex=.5,
     col=adjustcolor("#0d9488",.4), main="BMI vs Absent Hours",
     xlab="BMI", ylab="Jam", cex.axis=.82, las=1)
abline(lm(absent_hours ~ bmi, data=df), col="#dc2626", lwd=2, lty=2)
box(col="#cbd5e1")

barplot(tapply(df$absent_hours, df$disciplinary, mean, na.rm=TRUE),
        col=c("#2563eb","#dc2626"), names.arg=c("Tidak","Ya"),
        main="Rata-rata Jam vs Disiplin", ylab="Rata-rata Jam",
        border="white", cex.axis=.82)
box(col="#cbd5e1")

mtext("Visualisasi Eksploratif: Absenteeism at Work",
      outer=TRUE, cex=1.05, font=2, col="#1b2a4a")

2 Regresi Logistik Biner

2.1 Pembentukan Variabel Respons Biner

Variabel respons biner dibentuk dari jam ketidakhadiran:

  • Absen Tinggi (1): jam ketidakhadiran \(\geq\) median
  • Absen Rendah (0): jam ketidakhadiran \(<\) median
med_abs <- median(df$absent_hours, na.rm=TRUE)
cat("Median absent hours:", med_abs, "\n")
#> Median absent hours: 3
# Buat semua variabel turunan
df$absent_biner    <- as.integer(df$absent_hours >= med_abs)
df$education_f     <- factor(df$education, levels=1:4,
                              labels=c("SD","SMA","S1","S2/S3"))
df$drinker_f       <- factor(df$drinker, levels=c(0,1),
                              labels=c("Tidak","Ya"))
df$smoker_f        <- factor(df$smoker, levels=c(0,1),
                              labels=c("Tidak","Ya"))
df$disciplinary_f  <- factor(df$disciplinary, levels=c(0,1),
                              labels=c("Tidak","Ya"))

cat("Distribusi biner:\n")
#> Distribusi biner:
print(table(df$absent_biner))
#> 
#>   0   1 
#> 289 451
# Variabel 3 kategori untuk multinomial & ordinal
df$absent_3cat <- ifelse(df$absent_hours == 0, "Tidak Absen",
                  ifelse(df$absent_hours <= med_abs, "Absen Rendah", "Absen Tinggi"))
df$absent_3cat <- factor(df$absent_3cat,
                          levels=c("Tidak Absen","Absen Rendah","Absen Tinggi"))
df$absent_ord  <- ordered(df$absent_3cat,
                           levels=c("Tidak Absen","Absen Rendah","Absen Tinggi"))

cat("\nDistribusi 3 kategori:\n")
#> 
#> Distribusi 3 kategori:
print(table(df$absent_3cat))
#> 
#>  Tidak Absen Absen Rendah Absen Tinggi 
#>           44          357          339

2.2 Dasar Teori dan Model Matematis

Model regresi logistik biner memodelkan log-odds:

\[\log\left(\frac{P(Y_i=1 \mid \mathbf{x}_i)}{P(Y_i=0 \mid \mathbf{x}_i)}\right) = \beta_0 + \beta_1 x_{i1} + \cdots + \beta_p x_{ip}\]

Probabilitas kejadian:

\[P(Y_i=1 \mid \mathbf{x}_i) = \frac{1}{1+\exp(-\eta_i)}, \quad \eta_i = \beta_0 + \mathbf{x}_i^\top\boldsymbol{\beta}\]

2.3 Fungsi Likelihood dan MLE

\[L(\boldsymbol{\beta}) = \prod_{i=1}^{n} \pi_i^{y_i}(1-\pi_i)^{1-y_i}\]

\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^{n} \left[y_i\log\pi_i + (1-y_i)\log(1-\pi_i)\right]\]

Diestimasi menggunakan IRLS (Iteratively Reweighted Least Squares) via glm().

2.4 Estimasi Model

mod_biner <- glm(
  absent_biner ~ age + bmi + service_time + distance_residence +
    transport_expense + workload_avg + hit_target +
    drinker_f + smoker_f + disciplinary_f + education_f + children,
  data=df, family=binomial(link="logit")
)
summary(mod_biner)
#> 
#> Call:
#> glm(formula = absent_biner ~ age + bmi + service_time + distance_residence + 
#>     transport_expense + workload_avg + hit_target + drinker_f + 
#>     smoker_f + disciplinary_f + education_f + children, family = binomial(link = "logit"), 
#>     data = df)
#> 
#> Coefficients:
#>                      Estimate Std. Error z value Pr(>|z|)   
#> (Intercept)        -2.170e-02  2.536e+00  -0.009  0.99317   
#> age                -2.249e-02  1.988e-02  -1.132  0.25777   
#> bmi                 4.653e-02  2.705e-02   1.720  0.08539 . 
#> service_time       -1.890e-02  3.248e-02  -0.582  0.56069   
#> distance_residence -5.835e-03  7.704e-03  -0.757  0.44884   
#> transport_expense   4.317e-03  1.651e-03   2.615  0.00893 **
#> workload_avg       -2.069e-07  2.178e-06  -0.095  0.92431   
#> hit_target         -9.848e-03  2.318e-02  -0.425  0.67099   
#> drinker_fYa         5.938e-01  2.212e-01   2.685  0.00726 **
#> smoker_fYa          2.375e-01  4.073e-01   0.583  0.55980   
#> disciplinary_fYa   -1.870e+01  5.980e+02  -0.031  0.97506   
#> education_fSMA      2.985e-01  4.081e-01   0.731  0.46449   
#> education_fS1       1.992e-01  3.053e-01   0.653  0.51398   
#> education_fS2/S3   -8.263e-02  1.179e+00  -0.070  0.94412   
#> children            3.244e-01  9.872e-02   3.286  0.00102 **
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for binomial family taken to be 1)
#> 
#>     Null deviance: 990.10  on 739  degrees of freedom
#> Residual deviance: 854.59  on 725  degrees of freedom
#> AIC: 884.59
#> 
#> Number of Fisher Scoring iterations: 16

2.5 Uji Signifikansi Model: Likelihood Ratio Test

\(H_0: \beta_1 = \cdots = \beta_p = 0\)

\[G^2 = -2[\ell_{\text{null}} - \ell_{\text{penuh}}] \sim \chi^2_p\]

mod_null <- glm(absent_biner ~ 1, data=df, family=binomial)
G2_b  <- 2*(as.numeric(logLik(mod_biner)) - as.numeric(logLik(mod_null)))
df_b  <- mod_biner$rank - 1
p_G2_b <- pchisq(G2_b, df=df_b, lower.tail=FALSE)

cat("G² =", round(G2_b,4), "| df =", df_b,
    "| p-value =", format(p_G2_b, scientific=TRUE), "\n")
#> G² = 135.5104 | df = 14 | p-value = 5.523651e-22
cat("Keputusan:", ifelse(p_G2_b<0.05,"Tolak H0","Gagal Tolak H0"), "\n")
#> Keputusan: Tolak H0

2.6 Uji Parameter: Wald Test

\[Z_j = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)} \sim N(0,1)\]

sm_b     <- summary(mod_biner)$coefficients
ci_b     <- confint.default(mod_biner)
tbl_b    <- data.frame(
  Parameter = rownames(sm_b),
  Beta      = round(sm_b[,1], 4),
  SE        = round(sm_b[,2], 4),
  z         = round(sm_b[,3], 4),
  p_value   = round(sm_b[,4], 5),
  OR        = round(exp(sm_b[,1]), 4),
  CI_025    = round(exp(ci_b[,1]), 4),
  CI_975    = round(exp(ci_b[,2]), 4),
  Sig       = ifelse(sm_b[,4]<.001,"***",
              ifelse(sm_b[,4]<.01,"**",
              ifelse(sm_b[,4]<.05,"*",
              ifelse(sm_b[,4]<.1,".","")))),
  stringsAsFactors = FALSE, row.names = NULL
)

kable(tbl_b,
      caption="**Tabel 2.1. Uji Wald dan Odds Ratio — Logistik Biner**",
      col.names=c("Parameter","β","SE","z","p-value","OR","CI 2.5%","CI 97.5%","Sig"),
      align=c("l",rep("c",8))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, background="#f0f4ff") %>%
  column_spec(6, bold=TRUE, color="#1b2a4a") %>%
  row_spec(which(sm_b[,4]<0.05), background="#dcfce7")
Tabel 2.1. Uji Wald dan Odds Ratio — Logistik Biner
Parameter β SE z p-value OR CI 2.5% CI 97.5% Sig
(Intercept) -0.0217 2.5360 -0.0086 0.99317 0.9785 0.0068 141.0105
age -0.0225 0.0199 -1.1317 0.25777 0.9778 0.9404 1.0166
bmi 0.0465 0.0270 1.7202 0.08539 1.0476 0.9935 1.1047 .
service_time -0.0189 0.0325 -0.5818 0.56069 0.9813 0.9207 1.0458
distance_residence -0.0058 0.0077 -0.7573 0.44884 0.9942 0.9793 1.0093
transport_expense 0.0043 0.0017 2.6146 0.00893 1.0043 1.0011 1.0076 **
workload_avg 0.0000 0.0000 -0.0950 0.92431 1.0000 1.0000 1.0000
hit_target -0.0098 0.0232 -0.4248 0.67099 0.9902 0.9462 1.0362
drinker_fYa 0.5938 0.2212 2.6848 0.00726 1.8109 1.1739 2.7935 **
smoker_fYa 0.2375 0.4073 0.5831 0.55980 1.2681 0.5708 2.8174
disciplinary_fYa -18.6999 598.0358 -0.0313 0.97506 0.0000 0.0000 Inf
education_fSMA 0.2985 0.4081 0.7315 0.46449 1.3478 0.6057 2.9990
education_fS1 0.1992 0.3053 0.6526 0.51398 1.2205 0.6709 2.2202
education_fS2/S3 -0.0826 1.1789 -0.0701 0.94412 0.9207 0.0913 9.2805
children 0.3244 0.0987 3.2863 0.00102 1.3832 1.1399 1.6785 **

2.7 Forest Plot Odds Ratio

or_plot <- tbl_b[tbl_b$Parameter != "(Intercept)", ]
or_plot <- or_plot[order(or_plot$OR, decreasing=TRUE), ]

# Hilangkan nilai Inf dan NA
finite_vals <- c(or_plot$OR, or_plot$CI_025, or_plot$CI_975)
finite_vals <- finite_vals[is.finite(finite_vals)]

par(mar=c(4,12,3,2), bg="#f8fafc")
y_pos <- seq_len(nrow(or_plot))
cols  <- ifelse(or_plot$p_value < 0.05, "#2563eb", "#94a3b8")

xlim_r <- max(finite_vals, na.rm=TRUE) + 0.5

plot(or_plot$OR, y_pos, xlim=c(0, xlim_r), ylim=c(.5, nrow(or_plot)+.5),
     pch=18, cex=1.6, col=cols, xlab="Odds Ratio (95% CI)", ylab="",
     main="Forest Plot Odds Ratio — Regresi Logistik Biner",
     yaxt="n", las=1, cex.axis=.82)
for (i in y_pos)
  lines(c(or_plot$CI_025[i], or_plot$CI_975[i]), c(i,i), lwd=2, col=cols[i])
abline(v=1, lty=2, col="#dc2626", lwd=1.5)
axis(2, at=y_pos, labels=or_plot$Parameter, las=1, cex.axis=.77)
legend("topright", legend=c("p<0.05","p≥0.05"),
       col=c("#2563eb","#94a3b8"), pch=18, bty="n", cex=.82)
box(col="#cbd5e1")

2.8 Goodness of Fit

ll_null <- as.numeric(logLik(mod_null))
ll_full <- as.numeric(logLik(mod_biner))
n_obs   <- nrow(df)

mcfadden <- 1 - ll_full/ll_null
r2_cox   <- 1 - exp((2/n_obs)*(ll_null - ll_full))
r2_nag   <- r2_cox / (1 - exp((2/n_obs)*ll_null))

pred_p   <- predict(mod_biner, type="response")
pred_cls <- as.integer(pred_p >= 0.5)
akurasi  <- mean(pred_cls == df$absent_biner) * 100

cat("AIC              :", round(AIC(mod_biner),3),"\n")
#> AIC              : 884.594
cat("BIC              :", round(BIC(mod_biner),3),"\n")
#> BIC              : 953.693
cat("McFadden R²      :", round(mcfadden,4),"\n")
#> McFadden R²      : 0.1369
cat("Nagelkerke R²    :", round(r2_nag,4),"\n")
#> Nagelkerke R²    : 0.2269
cat("Akurasi Klasifikasi:", round(akurasi,2),"%\n")
#> Akurasi Klasifikasi: 66.49 %
cat("\nConfusion Matrix:\n")
#> 
#> Confusion Matrix:
print(table(Prediksi=pred_cls, Aktual=df$absent_biner))
#>         Aktual
#> Prediksi   0   1
#>        0 100  59
#>        1 189 392

2.9 Interpretasi

  • OR > 1 menunjukkan prediktor meningkatkan odds absen tinggi; OR < 1 bersifat protektif.
  • Prediktor dengan p < 0.05 secara statistik signifikan mempengaruhi odds ketidakhadiran tinggi setelah mengendalikan variabel lain.
  • Variabel yang tidak signifikan (p ≥ 0.05) seperti hit_target, age, dan service_time tidak memberikan bukti statistik yang cukup untuk menyimpulkan pengaruhnya terhadap odds absensi tinggi setelah mengendalikan variabel lain. Koefisien variabel tidak signifikan tidak boleh diinterpretasikan seolah-olah memiliki efek yang bermakna.

2.10 Kesimpulan Regresi Logistik Biner

Model logistik biner secara keseluruhan signifikan (LRT, p < 0.001). Beberapa prediktor demografis dan pekerjaan terbukti bermakna secara statistik dalam menjelaskan variasi status absensi tinggi/rendah, dengan OR yang dapat digunakan sebagai ukuran risiko relatif.

3 Regresi Logistik Multinomial

3.1 Pembentukan Variabel Respons Nominal

Tiga kategori nominal (tanpa asumsi urutan): Tidak Absen, Absen Rendah, Absen Tinggi. Kategori referensi (baseline): Tidak Absen.

cat("Distribusi 3 kategori:\n")
#> Distribusi 3 kategori:
print(table(df$absent_3cat))
#> 
#>  Tidak Absen Absen Rendah Absen Tinggi 
#>           44          357          339

3.2 Model Baseline-Category Logit

Untuk \(K=3\) kategori dengan referensi \(c\) (Tidak Absen):

\[\log\left(\frac{P(Y_i=k)}{P(Y_i=\text{Tidak Absen})}\right) = \alpha_k + \boldsymbol{\beta}_k^\top\mathbf{x}_i, \quad k \in \{\text{Absen Rendah, Absen Tinggi}\}\]

Probabilitas setiap kategori: \[P(Y_i=k \mid \mathbf{x}_i) = \frac{\exp(\alpha_k+\boldsymbol{\beta}_k^\top\mathbf{x}_i)}{1+\sum_{h}\exp(\alpha_h+\boldsymbol{\beta}_h^\top\mathbf{x}_i)}\]

3.3 Fungsi Likelihood dan MLE

\[L(\boldsymbol{\beta}) = \prod_{i=1}^n\prod_{k=1}^K \pi_{ik}^{y_{ik}}, \quad \ell(\boldsymbol{\beta}) = \sum_i\sum_k y_{ik}\log\pi_{ik}\]

Diestimasi dengan quasi-Newton (BFGS) via nnet::multinom().

3.4 Estimasi Model

# PENTING: Skalasi prediktor kontinu sebelum multinom 
# Tanpa skalasi, prediktor dengan rentang besar
# menghasilkan koefisien sangat besar → model tidak konvergen / probabilitas 0/1 sempurna

df_sc <- df
cont_pred <- c("age","bmi","service_time","distance_residence",
               "transport_expense","workload_avg","hit_target",
               "weight","height","children")
for (v in cont_pred) {
  if (v %in% names(df_sc))
    df_sc[[paste0(v,"_s")]] <- as.numeric(scale(df_sc[[v]]))
}

df_sc$absent_3cat <- relevel(df_sc$absent_3cat, ref="Tidak Absen")

cat("Cek skalasi (harus mean≈0, sd≈1):\n")
#> Cek skalasi (harus mean≈0, sd≈1):
print(round(sapply(paste0(cont_pred,"_s"), function(v)
  c(mean=mean(df_sc[[v]],na.rm=T), sd=sd(df_sc[[v]],na.rm=T))), 3))
#>      age_s bmi_s service_time_s distance_residence_s transport_expense_s
#> mean     0     0              0                    0                   0
#> sd       1     1              1                    1                   1
#>      workload_avg_s hit_target_s weight_s height_s children_s
#> mean              0            0        0        0          0
#> sd                1            1        1        1          1

Catatan: Quasi-Complete Separation pada disciplinary_f

Sebelum estimasi model, dilakukan pemeriksaan distribusi silang antara disciplinary_f dan variabel respons. Hasilnya menunjukkan bahwa seluruh observasi dengan disciplinary_f = "Ya" (nilai = 1) hanya terdapat pada satu kategori respons (Absen Tinggi), sehingga tidak ada variasi di kategori lain.

Kondisi ini disebut quasi-complete separation: prediktor memisahkan kategori respons secara hampir sempurna, menyebabkan:

  • Estimasi koefisien \(\hat{\beta}_{\text{disciplinary\_fYa}}\) menyimpang menuju \(+\infty\) atau \(-\infty\)
  • Standard error menjadi sangat besar (tidak terbatas)
  • Nilai z dan p-value tidak bermakna
  • Algoritma optimasi tidak konvergen dengan benar

Solusi: Variabel disciplinary_f dikeluarkan dari model multinomial dan digantikan dengan catatan ini. Variabel ini tetap digunakan pada model Logistik Biner dan Poisson karena pada model-model tersebut tidak terjadi separation.

#  Verifikasi quasi-complete separation pada disciplinary_f 
cat("Tabel silang: disciplinary_f × absent_3cat\n")
#> Tabel silang: disciplinary_f × absent_3cat
tbl_sep <- table(df_sc$disciplinary_f, df_sc$absent_3cat)
print(tbl_sep)
#>        
#>         Tidak Absen Absen Rendah Absen Tinggi
#>   Tidak           4          357          339
#>   Ya             40            0            0
cat("\nProporsi per baris:\n")
#> 
#> Proporsi per baris:
print(round(prop.table(tbl_sep, margin = 1), 4))
#>        
#>         Tidak Absen Absen Rendah Absen Tinggi
#>   Tidak      0.0057       0.5100       0.4843
#>   Ya         1.0000       0.0000       0.0000
# Identifikasi sel dengan 0 atau sangat sedikit observasi
cat("\nSel dengan frekuensi = 0 (indikasi separation):\n")
#> 
#> Sel dengan frekuensi = 0 (indikasi separation):
print(which(tbl_sep == 0, arr.ind = TRUE))
#>    row col
#> Ya   2   2
#> Ya   2   3
cat("\nKesimpulan: Seluruh observasi disciplinary_f = 'Ya' terkonsentrasi pada\n")
#> 
#> Kesimpulan: Seluruh observasi disciplinary_f = 'Ya' terkonsentrasi pada
cat("satu kategori respons → quasi-complete separation terdeteksi.\n")
#> satu kategori respons → quasi-complete separation terdeteksi.
cat("Variabel disciplinary_f DIKELUARKAN dari model multinomial.\n")
#> Variabel disciplinary_f DIKELUARKAN dari model multinomial.
# Model dengan prediktor terskalasi tanpa disciplinary_f
# disciplinary_f dikeluarkan karena quasi-complete separation
# (semua obs. disciplinary_f = "Ya" berada pada satu kategori respons)
set.seed(123)
mod_multi <- nnet::multinom(
  absent_3cat ~ age_s + bmi_s + service_time_s + distance_residence_s +
    transport_expense_s + workload_avg_s + hit_target_s +
    drinker_f + education_f + children_s,
  data    = df_sc,
  trace   = TRUE,    # tampilkan iterasi untuk cek konvergensi
  maxit   = 500,
  MaxNWts = 5000
)
#> # weights:  42 (26 variable)
#> initial  value 812.973094 
#> iter  10 value 604.802512
#> iter  20 value 578.904411
#> iter  30 value 576.012860
#> iter  40 value 575.907015
#> final  value 575.906484 
#> converged
cat("\n CEK KONVERGENSI \n")
#> 
#>  CEK KONVERGENSI
cat("Residual Deviance :", round(mod_multi$deviance, 3), "\n")
#> Residual Deviance : 1151.813
cat("AIC               :", round(AIC(mod_multi), 3), "\n")
#> AIC               : 1203.813
cat("Jumlah iterasi    :", mod_multi$convergence, "(0 = konvergen)\n")
#> Jumlah iterasi    : 0 (0 = konvergen)
# Cek koefisien masuk akal (tidak ekstrem)
cf_check <- summary(mod_multi)$coefficients
cat("\nRange koefisien (harus wajar, tidak >|10|):\n")
#> 
#> Range koefisien (harus wajar, tidak >|10|):
cat("Min:", round(min(cf_check),3), " Max:", round(max(cf_check),3), "\n")
#> Min: -0.841  Max: 11.541
cat("Ada koefisien ekstrem (|β|>10)?", any(abs(cf_check)>10), "\n")
#> Ada koefisien ekstrem (|β|>10)? TRUE
# Identifikasi koefisien ekstrem jika ada
if (any(abs(cf_check)>10)) {
  cat("\n  PERINGATAN: Koefisien ekstrem terdeteksi!\n")
  idx_ext <- which(abs(cf_check) > 10, arr.ind = TRUE)
  for (i in seq_len(nrow(idx_ext))) {
    cat("  Kategori:", rownames(cf_check)[idx_ext[i,1]],
        "| Prediktor:", colnames(cf_check)[idx_ext[i,2]],
        "| β =", round(cf_check[idx_ext[i,1], idx_ext[i,2]], 3), "\n")
  }
  cat("  → Koefisien ekstrem ini disebabkan oleh sparse data (sedikit observasi)\n")
  cat("    pada kategori prediktor tersebut dan TIDAK dapat diinterpretasikan.\n")
}
#> 
#>   PERINGATAN: Koefisien ekstrem terdeteksi!
#>   Kategori: Absen Rendah | Prediktor: education_fS2/S3 | β = 11.541 
#>   Kategori: Absen Tinggi | Prediktor: education_fS2/S3 | β = 10.774 
#>   → Koefisien ekstrem ini disebabkan oleh sparse data (sedikit observasi)
#>     pada kategori prediktor tersebut dan TIDAK dapat diinterpretasikan.

3.5 Uji Signifikansi Model

mod_multi0 <- nnet::multinom(absent_3cat ~ 1, data=df_sc, trace=FALSE)
G2_m  <- 2*(as.numeric(logLik(mod_multi)) - as.numeric(logLik(mod_multi0)))
df_m  <- mod_multi$edf - mod_multi0$edf
p_m   <- pchisq(G2_m, df=df_m, lower.tail=FALSE)

cat("G² =", round(G2_m,4), "| df =", df_m,
    "| p =", format(p_m, scientific=TRUE), "\n")
#> G² = 146.2892 | df = 24 | p = 1.614882e-19
cat("Keputusan:", ifelse(p_m<0.05,"Tolak H0 — model signifikan",
                         "Gagal Tolak H0"), "\n")
#> Keputusan: Tolak H0 — model signifikan

3.6 Uji Parameter dan RRR

\[z_j = \frac{\hat{\beta}_j}{\text{SE}(\hat{\beta}_j)}, \quad \text{RRR} = \exp(\hat{\beta}_j)\]

Catatan: Koefisien \(\hat{\beta}_j\) di sini merupakan efek per satu standar deviasi perubahan prediktor (karena prediktor telah diskalasi). Interpretasi RRR harus mempertimbangkan hal ini.

sm_m  <- summary(mod_multi)
cf_m  <- sm_m$coefficients
se_m  <- sm_m$standard.errors
z_m   <- cf_m / se_m
p_mat <- 2*(1 - pnorm(abs(z_m)))

# Susun per kategori
make_tbl <- function(cat_idx) {
  data.frame(
    Parameter = colnames(cf_m),
    Beta      = round(cf_m[cat_idx,], 4),
    SE        = round(se_m[cat_idx,], 4),
    z         = round(z_m[cat_idx,], 4),
    p_value   = round(p_mat[cat_idx,], 5),
    RRR       = round(exp(cf_m[cat_idx,]), 4),
    Sig       = ifelse(p_mat[cat_idx,]<.001,"***",
                ifelse(p_mat[cat_idx,]<.01,"**",
                ifelse(p_mat[cat_idx,]<.05,"*",
                ifelse(p_mat[cat_idx,]<.1,".","")))),
    stringsAsFactors=FALSE, row.names=NULL
  )
}

cat(" Kategori: Absen Rendah vs Tidak Absen \n")
#>  Kategori: Absen Rendah vs Tidak Absen
tbl_m1 <- make_tbl(1)
kable(tbl_m1, caption="**Tabel 3.1a. Absen Rendah vs Tidak Absen (Baseline)**",
      align=c("l",rep("c",6))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE, background="#1b2a4a", color="white") %>%
  column_spec(6, bold=TRUE, color="#1b2a4a") %>%
  row_spec(which(tbl_m1$p_value<0.05), background="#dcfce7")
Tabel 3.1a. Absen Rendah vs Tidak Absen (Baseline)
Parameter Beta SE z p_value RRR Sig
(Intercept) 2.8855 0.3804 7.5852 0.00000 17.9116 ***
age_s -0.4695 0.2253 -2.0842 0.03714 0.6253
bmi_s -0.4185 0.1836 -2.2791 0.02266 0.6581
service_time_s 0.4426 0.2887 1.5330 0.12528 1.5568
distance_residence_s 0.6026 0.2261 2.6650 0.00770 1.8269 **
transport_expense_s -0.6641 0.2037 -3.2602 0.00111 0.5147 **
workload_avg_s -0.1797 0.1669 -1.0771 0.28144 0.8355
hit_target_s 0.4279 0.1482 2.8877 0.00388 1.5341 **
drinker_fYa -0.7946 0.4268 -1.8619 0.06261 0.4518 .
education_fSMA -0.8407 0.7609 -1.1049 0.26918 0.4314
education_fS1 0.7404 1.1103 0.6668 0.50488 2.0968
education_fS2/S3 11.5413 0.5122 22.5312 0.00000 102874.9515 ***
children_s -0.4300 0.2015 -2.1341 0.03284 0.6505 
cat(" Kategori: Absen Tinggi vs Tidak Absen \n")
#>  Kategori: Absen Tinggi vs Tidak Absen
tbl_m2 <- make_tbl(2)
kable(tbl_m2, caption="**Tabel 3.1b. Absen Tinggi vs Tidak Absen (Baseline)**",
      align=c("l",rep("c",6))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=11) %>%
  column_spec(1, bold=TRUE, background="#1b2a4a", color="white") %>%
  column_spec(6, bold=TRUE, color="#1b2a4a") %>%
  row_spec(which(tbl_m2$p_value<0.05), background="#dcfce7")
Tabel 3.1b. Absen Tinggi vs Tidak Absen (Baseline)
Parameter Beta SE z p_value RRR Sig
(Intercept) 2.3185 0.3847 6.0260 0.00000 10.1602 ***
age_s -0.5367 0.2188 -2.4528 0.01417 0.5846
bmi_s -0.2493 0.1742 -1.4307 0.15252 0.7794
service_time_s 0.3778 0.2825 1.3375 0.18107 1.4591
distance_residence_s 0.2619 0.2163 1.2108 0.22595 1.2994
transport_expense_s -0.1500 0.1964 -0.7637 0.44506 0.8607
workload_avg_s -0.0755 0.1620 -0.4658 0.64133 0.9273
hit_target_s 0.3451 0.1408 2.4516 0.01422 1.4121
drinker_fYa 0.0538 0.4233 0.1270 0.89897 1.0552
education_fSMA -0.1684 0.7430 -0.2266 0.82071 0.8450
education_fS1 1.0856 1.1100 0.9780 0.32807 2.9612
education_fS2/S3 10.7741 0.5122 21.0336 0.00000 47769.8493 ***
children_s -0.0817 0.1948 -0.4194 0.67496 0.9215

3.7 Predicted Probabilities

# Grid prediksi: variasikan workload_avg_s, prediktor lain di nilai 0 (= rata-rata asli)
wl_s_seq <- seq(min(df_sc$workload_avg_s, na.rm=TRUE),
                max(df_sc$workload_avg_s, na.rm=TRUE), length.out=100)

newdat_m <- data.frame(
  age_s               = 0,
  bmi_s               = 0,
  service_time_s      = 0,
  distance_residence_s= 0,
  transport_expense_s = 0,
  workload_avg_s      = wl_s_seq,
  hit_target_s        = 0,
  drinker_f           = factor("Tidak", levels=levels(df_sc$drinker_f)),
  education_f         = factor("SMA",   levels=levels(df_sc$education_f)),
  children_s          = 0
)
# Catatan: disciplinary_f tidak disertakan karena telah dikeluarkan dari model
# multinomial (quasi-complete separation)

pp_m <- predict(mod_multi, newdata=newdat_m, type="probs")

# Konversi kembali ke skala asli untuk sumbu x
wl_orig <- wl_s_seq * sd(df$workload_avg,na.rm=TRUE) + mean(df$workload_avg,na.rm=TRUE)

par(mar=c(4.5,4.5,3.5,9.5), bg="#f8fafc")
matplot(wl_orig, pp_m, type="l", lwd=2.5, lty=1,
        col=c("#1b2a4a","#2563eb","#dc2626"),
        xlab="Workload Average per Day (skala asli)",
        ylab="Probabilitas Prediksi",
        main="Predicted Probabilities — Multinomial\n(prediktor lain di rata-rata)",
        ylim=c(0,1), las=1, cex.axis=.82)
legend("right", xpd=TRUE, inset=c(-0.32,0),
       legend=levels(df_sc$absent_3cat),
       col=c("#1b2a4a","#2563eb","#dc2626"),
       lwd=2.5, bty="n", cex=.82)
box(col="#cbd5e1")

3.8 Interpretasi

Catatan penting: Prediktor telah diskalasi (z-score) sebelum estimasi untuk mencegah masalah konvergensi numerik. Koefisien \(\hat{\beta}\) merepresentasikan efek per satu standar deviasi perubahan prediktor kontinu.

  • RRR > 1: Kenaikan satu SD prediktor meningkatkan odds relatif masuk kategori tersebut vs Tidak Absen.
  • RRR < 1: Bersifat protektif terhadap kategori tersebut vs baseline.
  • Dua persamaan logit diestimasi secara simultan: (1) Absen Rendah vs Tidak Absen, (2) Absen Tinggi vs Tidak Absen.
  • Prediktor yang signifikan pada keduanya mengindikasikan faktor yang konsisten membedakan kelompok absensi dari kelompok tidak absen.

3.9 Kesimpulan Regresi Logistik Multinomial

Model baseline-category logit diestimasi pada prediktor terskalasi untuk menjamin konvergensi numerik. Variabel disciplinary_f dikeluarkan dari model karena terdeteksi quasi-complete separation (seluruh observasi dengan disciplinary_f = "Ya" terkonsentrasi pada satu kategori respons, menyebabkan estimasi koefisien tidak stabil). Model akhir tanpa disciplinary_f berhasil konvergen (convergence = 0). Tabel 3.1a dan 3.1b menyajikan koefisien dan RRR yang valid. Predicted probabilities menunjukkan bagaimana distribusi kategori absensi berubah seiring variasi beban kerja.

Catatan Koefisien Ekstrem education_fS2/S3: Apabila terdeteksi koefisien ekstrem (|β| > 10) pada kategori education_fS2/S3, hal ini disebabkan oleh masalah sparse data kategori S2/S3 memiliki sangat sedikit observasi sehingga estimasi tidak stabil. RRR yang sangat besar (misalnya > 100) pada kategori ini tidak dapat diinterpretasikan secara substantif dan perlu penanganan lanjutan: penggabungan kategori S1+S2/S3 menjadi satu kategori “Perguruan Tinggi”, atau pengeluaran kategori tersebut dari model.

4 Regresi Logistik Ordinal

4.1 Pembentukan Variabel Respons Ordinal

Variabel respons ordinal menggunakan urutan yang bermakna: Tidak Absen < Absen Rendah < Absen Tinggi. Model memanfaatkan informasi urutan ini — berbeda dengan multinomial yang mengabaikannya.

cat("Distribusi ordinal:\n")
#> Distribusi ordinal:
print(table(df$absent_ord))
#> 
#>  Tidak Absen Absen Rendah Absen Tinggi 
#>           44          357          339
cat("\nProporsi per kategori:\n")
#> 
#> Proporsi per kategori:
print(round(prop.table(table(df$absent_ord)), 4))
#> 
#>  Tidak Absen Absen Rendah Absen Tinggi 
#>       0.0595       0.4824       0.4581
cat("\nSemua kategori ≥ 5 observasi:", all(table(df$absent_ord) >= 5), "\n")
#> 
#> Semua kategori ≥ 5 observasi: TRUE

4.2 Model Cumulative Logit (Proportional Odds)

\[\log\left(\frac{P(Y_i \leq j \mid \mathbf{x}_i)}{P(Y_i > j \mid \mathbf{x}_i)}\right) = \alpha_j - \boldsymbol{\beta}^\top\mathbf{x}_i, \quad j = 1, 2\]

  • \(\alpha_j\) adalah cut-point (intercept) ke-\(j\), dengan \(\alpha_1 < \alpha_2\)
  • \(\boldsymbol{\beta}\) diasumsikan sama untuk semua cut-point (proportional odds)

Cut-point merepresentasikan log-odds kumulatif ketika semua prediktor = 0. Harus selalu dilaporkan bersama koefisien.

4.3 Fungsi Likelihood dan MLE

\[L(\boldsymbol{\alpha},\boldsymbol{\beta}) = \prod_{i=1}^n \prod_{j=1}^J \left[F(\alpha_j - \eta_i) - F(\alpha_{j-1} - \eta_i)\right]^{\mathbf{1}(Y_i=j)}\]

di mana \(F\) adalah fungsi logistik dan \(\eta_i = \boldsymbol{\beta}^\top\mathbf{x}_i\).

Diestimasi dengan Fisher Scoring via MASS::polr().

4.4 Estimasi Model

Catatan: Konsistensi Penanganan disciplinary_f

Seperti pada model multinomial, disciplinary_f juga menunjukkan quasi-complete separation pada model ordinal (β sangat besar dengan SE ≈ 0, menghasilkan z-value ekstrem seperti z = −7.4×10⁹). Kondisi ini identik dengan yang sudah terdeteksi di bagian multinomial. Demi konsistensi dan validitas estimasi, disciplinary_f dikeluarkan juga dari model ordinal. Jika dimasukkan, koefisiennya tidak dapat diinterpretasikan dan akan merusak validitas seluruh model.

# Gunakan prediktor terskalasi (sama dengan multinomial) untuk konsistensi
# disciplinary_f DIKELUARKAN dari model ordinal karena quasi-complete separation
# (sama dengan perlakuan di model multinomial — konsistensi wajib dijaga)
mod_ord <- tryCatch(
  MASS::polr(
    absent_ord ~ age_s + bmi_s + service_time_s + distance_residence_s +
      transport_expense_s + workload_avg_s + hit_target_s +
      drinker_f + education_f + children_s,
    data=df_sc, Hess=TRUE, method="logistic"
  ),
  error = function(e) {
    message("Model penuh gagal. Mencoba model parsimoni...")
    MASS::polr(
      absent_ord ~ age_s + bmi_s + workload_avg_s + drinker_f,
      data=df_sc, Hess=TRUE, method="logistic"
    )
  }
)

cat("Model ordinal berhasil diestimasi (tanpa disciplinary_f).\n\n")
#> Model ordinal berhasil diestimasi (tanpa disciplinary_f).
# Tampilkan summary lengkap termasuk cut-points
print(summary(mod_ord))
#> Call:
#> MASS::polr(formula = absent_ord ~ age_s + bmi_s + service_time_s + 
#>     distance_residence_s + transport_expense_s + workload_avg_s + 
#>     hit_target_s + drinker_f + education_f + children_s, data = df_sc, 
#>     Hess = TRUE, method = "logistic")
#> 
#> Coefficients:
#>                         Value Std. Error t value
#> age_s                -0.20896    0.11562 -1.8072
#> bmi_s                 0.03949    0.10016  0.3943
#> service_time_s        0.06808    0.12142  0.5607
#> distance_residence_s -0.20480    0.09880 -2.0729
#> transport_expense_s   0.33641    0.09785  3.4381
#> workload_avg_s        0.05153    0.07630  0.6754
#> hit_target_s          0.06207    0.07771  0.7988
#> drinker_fYa           0.61158    0.19758  3.0954
#> education_fSMA        0.46526    0.33215  1.4008
#> education_fS1         0.28405    0.28027  1.0135
#> education_fS2/S3     -0.35866    0.98480 -0.3642
#> children_s            0.25811    0.08984  2.8730
#> 
#> Intercepts:
#>                           Value    Std. Error t value 
#> Tidak Absen|Absen Rendah   -2.5088   0.1971   -12.7270
#> Absen Rendah|Absen Tinggi   0.5529   0.1491     3.7091
#> 
#> Residual Deviance: 1240.062 
#> AIC: 1268.062
cat("\n CUT-POINTS (Intercepts/Thresholds) \n")
#> 
#>  CUT-POINTS (Intercepts/Thresholds)
cat("Nilai cut-point merepresentasikan log-odds kumulatif saat semua prediktor = 0\n")
#> Nilai cut-point merepresentasikan log-odds kumulatif saat semua prediktor = 0
print(round(mod_ord$zeta, 4))
#>  Tidak Absen|Absen Rendah Absen Rendah|Absen Tinggi 
#>                   -2.5088                    0.5529
cat("Exp(cut-point):", round(exp(mod_ord$zeta), 4),
    "(odds kumulatif pada nilai rata-rata prediktor)\n")
#> Exp(cut-point): 0.0814 1.7383 (odds kumulatif pada nilai rata-rata prediktor)

4.5 Uji Signifikansi Model (LRT)

mod_ord0 <- MASS::polr(absent_ord ~ 1, data=df_sc, Hess=TRUE, method="logistic")
G2_o  <- 2*(as.numeric(logLik(mod_ord)) - as.numeric(logLik(mod_ord0)))
df_o  <- length(coef(mod_ord))
p_o   <- pchisq(G2_o, df=df_o, lower.tail=FALSE)

cat("Log-likelihood model penuh :", round(as.numeric(logLik(mod_ord)),3),"\n")
#> Log-likelihood model penuh : -620.031
cat("Log-likelihood model null  :", round(as.numeric(logLik(mod_ord0)),3),"\n")
#> Log-likelihood model null  : -649.051
cat("G² =", round(G2_o,4), "| df =", df_o,
    "| p =", format(p_o,scientific=TRUE), "\n")
#> G² = 58.0406 | df = 12 | p = 5.124738e-08
cat("Keputusan:", ifelse(p_o<0.05,"Tolak H0","Gagal Tolak H0"),"\n")
#> Keputusan: Tolak H0
cat("AIC:", round(AIC(mod_ord),3), "| BIC:", round(BIC(mod_ord),3),"\n")
#> AIC: 1268.062 | BIC: 1332.555

4.6 Uji Parameter dan Odds Ratio

cf_o  <- coef(mod_ord)
vc_o  <- vcov(mod_ord)
se_o  <- sqrt(diag(vc_o)[names(cf_o)])
z_o   <- cf_o / se_o
p_o2  <- 2*(1 - pnorm(abs(z_o)))
ci_o  <- confint.default(mod_ord)[names(cf_o), , drop=FALSE]

tbl_o <- data.frame(
  Parameter = names(cf_o),
  Beta      = round(cf_o, 4),
  SE        = round(se_o, 4),
  z         = round(z_o, 4),
  p_value   = round(p_o2, 5),
  OR        = round(exp(cf_o), 4),
  CI_025    = round(exp(ci_o[,1]), 4),
  CI_975    = round(exp(ci_o[,2]), 4),
  Sig       = ifelse(p_o2<.001,"***",
              ifelse(p_o2<.01,"**",
              ifelse(p_o2<.05,"*",
              ifelse(p_o2<.1,".","")))),
  stringsAsFactors=FALSE, row.names=NULL
)

kable(tbl_o,
      caption="**Tabel 4.1. Koefisien, OR, Uji Wald — Logistik Ordinal** (prediktor terskalasi)",
      align=c("l",rep("c",8))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, background="#f0f4ff") %>%
  column_spec(6, bold=TRUE, color="#1b2a4a") %>%
  row_spec(which(p_o2<0.05), background="#dcfce7") %>%
  footnote(general="β dan OR merupakan efek per 1 SD perubahan prediktor kontinu (karena terskalasi).")
Tabel 4.1. Koefisien, OR, Uji Wald — Logistik Ordinal (prediktor terskalasi)
Parameter Beta SE z p_value OR CI_025 CI_975 Sig
age_s -0.2090 0.1156 -1.8072 0.07072 0.8114 0.6469 1.0178 .
bmi_s 0.0395 0.1002 0.3943 0.69335 1.0403 0.8549 1.2659
service_time_s 0.0681 0.1214 0.5607 0.57501 1.0705 0.8437 1.3581
distance_residence_s -0.2048 0.0988 -2.0729 0.03818 0.8148 0.6714 0.9889
transport_expense_s 0.3364 0.0978 3.4381 0.00059 1.3999 1.1556 1.6959 ***
workload_avg_s 0.0515 0.0763 0.6754 0.49944 1.0529 0.9066 1.2227
hit_target_s 0.0621 0.0777 0.7988 0.42441 1.0640 0.9137 1.2391
drinker_fYa 0.6116 0.1976 3.0954 0.00197 1.8433 1.2515 2.7151 **
education_fSMA 0.4653 0.3322 1.4008 0.16129 1.5924 0.8305 3.0535
education_fS1 0.2840 0.2803 1.0135 0.31084 1.3285 0.7670 2.3011
education_fS2/S3 -0.3587 0.9848 -0.3642 0.71571 0.6986 0.1014 4.8139
children_s 0.2581 0.0898 2.8730 0.00407 1.2945 1.0855 1.5437 **
Note:
β dan OR merupakan efek per 1 SD perubahan prediktor kontinu (karena terskalasi).
Tabel 4.2. Cut-Points (Thresholds) Model Ordinal
Cut-Point α (logit) exp(α) Interpretasi
Tidak Absen&#124;Absen Rendah -2.5088 0.0814 Log-odds kumulatif P(Y≤Tidak Absen) saat prediktor=0
Absen Rendah&#124;Absen Tinggi 0.5529 1.7383 Log-odds kumulatif P(Y≤Absen Rendah) saat prediktor=0

4.7 Pemeriksaan Asumsi Proportional Odds

Asumsi Proportional Odds (PO) menyatakan bahwa efek prediktor konstan di seluruh cut-point. Jika dilanggar, OR yang dilaporkan tidak konsisten di berbagai level respons dan inferensi menjadi tidak valid.

Pengujian menggunakan Likelihood Ratio Test: bandingkan model PO (polr) vs model tanpa asumsi PO (multinomial). Jika p < 0.05, asumsi PO dilanggar.

Catatan: Perbandingan LRT antara model ordinal dan multinomial hanya valid jika kedua model menggunakan set prediktor yang sama. Dalam laporan ini, setelah revisi, kedua model menggunakan prediktor yang identik (tanpa disciplinary_f), sehingga perbandingan ini sah secara formal.

# LRT: ordinal (asumsi PO) vs multinomial (tanpa asumsi PO)
ll_ord   <- as.numeric(logLik(mod_ord))
ll_multi <- as.numeric(logLik(mod_multi))
G2_po    <- 2 * (ll_multi - ll_ord)   # multinomial selalu >= ordinal LL

# df = selisih jumlah parameter
p_ord_free  <- length(coef(mod_ord))   + length(mod_ord$zeta)
p_multi_free <- mod_multi$edf
df_po <- max(p_multi_free - p_ord_free, 1)

p_po <- pchisq(G2_po, df=df_po, lower.tail=FALSE)

cat(" UJI ASUMSI PROPORTIONAL ODDS \n")
#>  UJI ASUMSI PROPORTIONAL ODDS
cat("H0: Asumsi proportional odds TERPENUHI\n")
#> H0: Asumsi proportional odds TERPENUHI
cat("H1: Asumsi proportional odds DILANGGAR\n\n")
#> H1: Asumsi proportional odds DILANGGAR
cat("LL model ordinal (PO)       :", round(ll_ord,3), "\n")
#> LL model ordinal (PO)       : -620.031
cat("LL model multinomial (free) :", round(ll_multi,3), "\n")
#> LL model multinomial (free) : -575.906
cat("G² = 2*(LL_multi - LL_ord)  :", round(G2_po,4), "\n")
#> G² = 2*(LL_multi - LL_ord)  : 88.2486
cat("df                          :", df_po, "\n")
#> df                          : 12
cat("P-value                     :", round(p_po,4), "\n\n")
#> P-value                     : 0
if (p_po > 0.05) {
  cat("KEPUTUSAN: Gagal Tolak H0\n")
  cat("Asumsi proportional odds TERPENUHI (p =", round(p_po,4), "> 0.05)\n")
  cat("→ Model polr valid untuk inferensi.\n")
  asumsi_po_ok <- TRUE
} else {
  cat("KEPUTUSAN: Tolak H0\n")
  cat("Asumsi proportional odds TIDAK TERPENUHI (p =", round(p_po,4), "< 0.05)\n")
  cat("→ Model polr TIDAK VALID untuk inferensi langsung.\n")
  cat("→ Alternatif yang direkomendasikan:\n")
  cat("   1. Partial Proportional Odds Model (VGAM::vglm dengan constraint parsial)\n")
  cat("   2. Kembali ke model Multinomial (Bagian 3)\n")
  cat("   3. Model Adjacent-Category Logit\n")
  asumsi_po_ok <- FALSE
}
#> KEPUTUSAN: Tolak H0
#> Asumsi proportional odds TIDAK TERPENUHI (p = 0 < 0.05)
#> → Model polr TIDAK VALID untuk inferensi langsung.
#> → Alternatif yang direkomendasikan:
#>    1. Partial Proportional Odds Model (VGAM::vglm dengan constraint parsial)
#>    2. Kembali ke model Multinomial (Bagian 3)
#>    3. Model Adjacent-Category Logit
# Visualisasi empiris: proporsi kumulatif per cut-point vs prediktor signifikan
# (pendekatan grafis untuk cek PO secara visual)

# Ambil prediktor paling signifikan
sig_preds <- names(p_o2)[p_o2 < 0.1]
if (length(sig_preds) == 0) sig_preds <- names(sort(p_o2))[1:2]

cat("\nPrediktor signifikan untuk visualisasi PO:", sig_preds[1], "\n")
#> 
#> Prediktor signifikan untuk visualisasi PO: age_s
# Gunakan workload_avg_s sebagai contoh
df_sc$wl_qtile <- cut(df_sc$workload_avg_s, breaks=4,
                      labels=c("Q1","Q2","Q3","Q4"))

emp_po <- tapply(as.numeric(df_sc$absent_ord), df_sc$wl_qtile, function(y) {
  c(cp1 = mean(y <= 1), cp2 = mean(y <= 2))
})

cp_mat <- do.call(rbind, emp_po)

par(mar=c(4,4,3,1.5), bg="#f8fafc")
matplot(1:4, cp_mat, type="b", pch=c(19,17), cex=1.2,
        col=c("#2563eb","#dc2626"), lwd=2,
        xlab="Kuartil Workload Average", ylab="Proporsi Kumulatif Empiris",
        main="Cek Visual Proportional Odds\n(garis paralel = asumsi terpenuhi)",
        xaxt="n", ylim=c(0,1), las=1, cex.axis=.82)
axis(1, at=1:4, labels=c("Q1","Q2","Q3","Q4"))
legend("bottomright",
       legend=c("P(Y≤Tidak Absen)","P(Y≤Absen Rendah)"),
       col=c("#2563eb","#dc2626"), pch=c(19,17), lwd=2, bty="n", cex=.85)
abline(h=c(0,1), lty=3, col="#94a3b8")
box(col="#cbd5e1")

Interpretasi plot: Jika asumsi proportional odds terpenuhi, kedua kurva (garis biru dan merah) harus kira-kira paralel. Jika salah satu kurva memotong atau memiliki kemiringan yang sangat berbeda, ini mengindikasikan pelanggaran asumsi untuk cut-point tersebut.

4.8 Predicted Probabilities

# Grid prediksi ordinal — disciplinary_f tidak disertakan (dikeluarkan dari model)
newdat_o <- data.frame(
  age_s               = 0,
  bmi_s               = 0,
  service_time_s      = 0,
  distance_residence_s= 0,
  transport_expense_s = 0,
  workload_avg_s      = wl_s_seq,
  hit_target_s        = 0,
  drinker_f           = factor("Tidak", levels=levels(df_sc$drinker_f)),
  education_f         = factor("SMA",   levels=levels(df_sc$education_f)),
  children_s          = 0
)
pp_o_grid <- predict(mod_ord, newdata=newdat_o, type="probs")

par(mar=c(4.5,4.5,3.5,9.5), bg="#f8fafc")
matplot(wl_orig, pp_o_grid, type="l", lwd=2.5, lty=1,
        col=c("#1b2a4a","#2563eb","#dc2626"),
        xlab="Workload Average per Day (skala asli)",
        ylab="Probabilitas Prediksi",
        main="Predicted Probabilities — Ordinal (Proportional Odds)\n(prediktor lain di rata-rata)",
        ylim=c(0,1), las=1, cex.axis=.82)
legend("right", xpd=TRUE, inset=c(-0.32,0),
       legend=levels(df_sc$absent_ord),
       col=c("#1b2a4a","#2563eb","#dc2626"),
       lwd=2.5, bty="n", cex=.82)
box(col="#cbd5e1")

4.9 Interpretasi

** Asumsi PO TIDAK terpenuhi — interpretasi OR harus dilakukan dengan sangat hati-hati.**

Karena asumsi proportional odds dilanggar:

  1. OR dari polr TIDAK BOLEH diinterpretasikan sebagai efek yang konsisten di seluruh level respons.
  2. Koefisien yang dilaporkan merupakan rata-rata efek tertimbang yang tidak memiliki interpretasi tunggal yang valid.
  3. Rekomendasi: Gunakan hasil dari model Multinomial (Bagian 3) sebagai referensi utama, karena tidak memerlukan asumsi PO.
  4. Jika ingin tetap menggunakan model ordinal, gunakan Partial Proportional Odds (VGAM::vglm) yang hanya membatasi prediktor yang memenuhi asumsi PO.

4.10 Kesimpulan Regresi Logistik Ordinal

** Kesimpulan Kritis:** Model proportional odds (polr) berhasil diestimasi secara numerik, namun asumsi proportional odds dilanggar (p < 0.05 pada uji LRT). Ini berarti efek prediktor tidak konstan di seluruh cut-point, sehingga satu set koefisien tidak dapat merepresentasikan hubungan secara akurat. Hasil di bagian ini tidak boleh dijadikan dasar inferensi. Gunakan model Multinomial (Bagian 3) atau Partial Proportional Odds sebagai alternatif.

5 Regresi Poisson

5.1 Variabel Respons: Jam Ketidakhadiran (Count)

absent_hours merupakan variabel count — cocok dimodelkan dengan distribusi Poisson sebagai titik awal.

par(mfrow=c(1,2), mar=c(4,4,3,1.5), bg="#f8fafc")

hist(df$absent_hours, breaks=30, col="#2563eb", border="white",
     main="Distribusi Absent Hours", xlab="Jam", cex.axis=.82, las=1)
abline(v=mean(df$absent_hours), col="#dc2626", lwd=2, lty=2)
legend("topright", paste("Mean=",round(mean(df$absent_hours),1)),
       col="#dc2626", lwd=2, bty="n", cex=.82)
box(col="#cbd5e1")

plot(log1p(table(df$absent_hours)), type="h", col="#2563eb", lwd=2,
     main="Log Frekuensi per Nilai", xlab="Jam", ylab="log(Frekuensi+1)",
     cex.axis=.82, las=1)
box(col="#cbd5e1")

5.3 Fungsi Likelihood dan MLE

\[\ell(\boldsymbol{\beta}) = \sum_{i=1}^n \left[y_i\log\mu_i - \mu_i - \log(y_i!)\right]\]

5.4 Estimasi Model

mod_poi <- glm(
  absent_hours ~ age + bmi + service_time + distance_residence +
    transport_expense + workload_avg + hit_target +
    drinker_f + smoker_f + disciplinary_f + education_f + children,
  data=df, family=poisson(link="log")
)
summary(mod_poi)
#> 
#> Call:
#> glm(formula = absent_hours ~ age + bmi + service_time + distance_residence + 
#>     transport_expense + workload_avg + hit_target + drinker_f + 
#>     smoker_f + disciplinary_f + education_f + children, family = poisson(link = "log"), 
#>     data = df)
#> 
#> Coefficients:
#>                      Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)         8.640e-01  4.313e-01   2.003  0.04517 *  
#> age                 2.099e-02  2.855e-03   7.354 1.93e-13 ***
#> bmi                -3.836e-02  4.276e-03  -8.971  < 2e-16 ***
#> service_time       -1.586e-03  5.611e-03  -0.283  0.77747    
#> distance_residence -2.016e-02  1.137e-03 -17.736  < 2e-16 ***
#> transport_expense   1.208e-03  2.474e-04   4.883 1.04e-06 ***
#> workload_avg        5.046e-07  3.582e-07   1.409  0.15888    
#> hit_target          1.266e-02  3.926e-03   3.224  0.00127 ** 
#> drinker_fYa         4.301e-01  3.785e-02  11.362  < 2e-16 ***
#> smoker_fYa         -2.841e-01  6.436e-02  -4.415 1.01e-05 ***
#> disciplinary_fYa   -1.729e+01  1.853e+02  -0.093  0.92568    
#> education_fSMA      7.392e-04  6.996e-02   0.011  0.99157    
#> education_fS1      -1.948e-01  6.030e-02  -3.230  0.00124 ** 
#> education_fS2/S3   -9.948e-01  2.207e-01  -4.507 6.58e-06 ***
#> children            1.283e-01  1.380e-02   9.291  < 2e-16 ***
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#> 
#> (Dispersion parameter for poisson family taken to be 1)
#> 
#>     Null deviance: 8092.7  on 739  degrees of freedom
#> Residual deviance: 6664.5  on 725  degrees of freedom
#> AIC: 8991.5
#> 
#> Number of Fisher Scoring iterations: 13

5.5 Uji Signifikansi Model

mod_poi0 <- glm(absent_hours ~ 1, data=df, family=poisson)
G2_p  <- 2*(as.numeric(logLik(mod_poi)) - as.numeric(logLik(mod_poi0)))
df_p  <- mod_poi$rank - 1
p_G2p <- pchisq(G2_p, df=df_p, lower.tail=FALSE)
cat("G² =",round(G2_p,4),"| df =",df_p,"| p =",format(p_G2p,scientific=TRUE),"\n")
#> G² = 1428.159 | df = 14 | p = 1.406109e-296
cat("Keputusan:", ifelse(p_G2p<0.05,"Tolak H0","Gagal Tolak H0"),"\n")
#> Keputusan: Tolak H0

5.6 Uji Parameter dan Incidence Rate Ratio (IRR)

\[\text{IRR}_j = \exp(\hat{\beta}_j)\]

IRR = 1.5 berarti kenaikan satu satuan prediktor meningkatkan rata-rata jam absen 50%, ceteris paribus.

sm_p   <- summary(mod_poi)$coefficients
ci_p   <- confint.default(mod_poi)
tbl_p  <- data.frame(
  Parameter = rownames(sm_p),
  Beta      = round(sm_p[,1], 4),
  SE        = round(sm_p[,2], 4),
  z         = round(sm_p[,3], 4),
  p_value   = round(sm_p[,4], 5),
  IRR       = round(exp(sm_p[,1]), 4),
  CI_025    = round(exp(ci_p[,1]), 4),
  CI_975    = round(exp(ci_p[,2]), 4),
  Sig       = ifelse(sm_p[,4]<.001,"***",
              ifelse(sm_p[,4]<.01,"**",
              ifelse(sm_p[,4]<.05,"*",
              ifelse(sm_p[,4]<.1,".","")))),
  stringsAsFactors=FALSE, row.names=NULL
)

kable(tbl_p,
      caption="**Tabel 5.1. Koefisien, IRR, Uji Wald — Regresi Poisson**",
      col.names=c("Parameter","β","SE","z","p-value","IRR","CI 2.5%","CI 97.5%","Sig"),
      align=c("l",rep("c",8))) %>%
  kable_styling(bootstrap_options=c("striped","hover","bordered","condensed"),
                full_width=TRUE, font_size=12) %>%
  column_spec(1, bold=TRUE, background="#f0f4ff") %>%
  column_spec(6, bold=TRUE, color="#1b2a4a") %>%
  row_spec(which(sm_p[,4]<0.05), background="#dcfce7")
Tabel 5.1. Koefisien, IRR, Uji Wald — Regresi Poisson
Parameter β SE z p-value IRR CI 2.5% CI 97.5% Sig
(Intercept) 0.8640 0.4313 2.0031 0.04517 2.3725 1.0188 5.525200e+00
age 0.0210 0.0029 7.3538 0.00000 1.0212 1.0155 1.026900e+00 ***
bmi -0.0384 0.0043 -8.9708 0.00000 0.9624 0.9543 9.705000e-01 ***
service_time -0.0016 0.0056 -0.2826 0.77747 0.9984 0.9875 1.009500e+00
distance_residence -0.0202 0.0011 -17.7364 0.00000 0.9800 0.9779 9.822000e-01 ***
transport_expense 0.0012 0.0002 4.8835 0.00000 1.0012 1.0007 1.001700e+00 ***
workload_avg 0.0000 0.0000 1.4088 0.15888 1.0000 1.0000 1.000000e+00
hit_target 0.0127 0.0039 3.2237 0.00127 1.0127 1.0050 1.020600e+00 **
drinker_fYa 0.4301 0.0379 11.3622 0.00000 1.5374 1.4275 1.655800e+00 ***
smoker_fYa -0.2841 0.0644 -4.4150 0.00001 0.7527 0.6635 8.538000e-01 ***
disciplinary_fYa -17.2872 185.3160 -0.0933 0.92568 0.0000 0.0000 1.712215e+150
education_fSMA 0.0007 0.0700 0.0106 0.99157 1.0007 0.8725 1.147800e+00
education_fS1 -0.1948 0.0603 -3.2299 0.00124 0.8230 0.7313 9.263000e-01 **
education_fS2/S3 -0.9948 0.2207 -4.5069 0.00001 0.3698 0.2399 5.700000e-01 ***
children 0.1283 0.0138 9.2911 0.00000 1.1368 1.1065 1.168000e+00 *** 
# Deteksi CI ekstrem (indikasi separation/estimasi tidak stabil)
ci_upper_p <- exp(ci_p[,2])
idx_extreme <- which(!is.finite(ci_upper_p) | ci_upper_p > 1e10)
if (length(idx_extreme) > 0) {
  cat("\n  PERINGATAN: CI ekstrem terdeteksi pada parameter berikut:\n")
  for (idx in idx_extreme) {
    cat("  Parameter:", rownames(ci_p)[idx],
        "| CI upper =", formatC(ci_upper_p[idx], format="e", digits=3), "\n")
  }
  cat("  → CI yang sangat besar (mis. > 1×10¹⁰ atau Inf) mengindikasikan\n")
  cat("    quasi-complete separation atau estimasi tidak stabil.\n")
  cat("    IRR parameter ini TIDAK DAPAT diinterpretasikan secara substantif.\n")
  cat("    Pertimbangkan pengeluaran variabel tersebut atau penggabungan kategori.\n")
}
#> 
#>   PERINGATAN: CI ekstrem terdeteksi pada parameter berikut:
#>   Parameter: disciplinary_fYa | CI upper = 1.712e+150 
#>   → CI yang sangat besar (mis. > 1×10¹⁰ atau Inf) mengindikasikan
#>     quasi-complete separation atau estimasi tidak stabil.
#>     IRR parameter ini TIDAK DAPAT diinterpretasikan secara substantif.
#>     Pertimbangkan pengeluaran variabel tersebut atau penggabungan kategori.

5.7 Pemeriksaan Equidispersion dan Overdispersion

Equidispersion: \(E(Y)=\text{Var}(Y)=\mu\). Jika \(\text{Var}(Y)>\mu\) maka terjadi overdispersion.

\[\hat{\phi} = \frac{\text{Devians}}{df_{\text{residual}}}\]

phi_hat <- mod_poi$deviance / mod_poi$df.residual
cat("Mean absent_hours  :", round(mean(df$absent_hours),3),"\n")
#> Mean absent_hours  : 6.924
cat("Var  absent_hours  :", round(var(df$absent_hours),3),"\n")
#> Var  absent_hours  : 177.716
cat("Rasio Var/Mean     :", round(var(df$absent_hours)/mean(df$absent_hours),3),"\n")
#> Rasio Var/Mean     : 25.665
cat("Devians/df (phi)   :", round(phi_hat,4),"\n")
#> Devians/df (phi)   : 9.1924
cat("Status             :", ifelse(phi_hat>1.5,"OVERDISPERSION","Equidispersion OK"),"\n")
#> Status             : OVERDISPERSION
disp_test <- AER::dispersiontest(mod_poi)
cat("\nUji Overdispersion (AER):\n")
#> 
#> Uji Overdispersion (AER):
print(disp_test)
#> 
#>  Overdispersion test
#> 
#> data:  mod_poi
#> z = 4.0746, p-value = 2.305e-05
#> alternative hypothesis: true dispersion is greater than 1
#> sample estimates:
#> dispersion 
#>   18.02666

5.8 Alternatif Model: Quasi-Poisson dan Negative Binomial

# Quasi-Poisson
mod_qpoi <- glm(
  absent_hours ~ age + bmi + service_time + distance_residence +
    transport_expense + workload_avg + hit_target +
    drinker_f + smoker_f + disciplinary_f + education_f + children,
  data=df, family=quasipoisson(link="log")
)
cat("Quasi-Poisson dispersion:", round(summary(mod_qpoi)$dispersion,4),"\n")
#> Quasi-Poisson dispersion: 18.3519
# Negative Binomial
mod_nb <- MASS::glm.nb(
  absent_hours ~ age + bmi + service_time + distance_residence +
    transport_expense + workload_avg + hit_target +
    drinker_f + smoker_f + disciplinary_f + education_f + children,
  data=df
)
cat("NB theta (1/overdispersion):", round(mod_nb$theta,4),"\n")
#> NB theta (1/overdispersion): 1.1805
cat("AIC Poisson :", round(AIC(mod_poi),3),"\n")
#> AIC Poisson : 8991.533
cat("AIC Neg.Bin :", round(AIC(mod_nb),3),"\n")
#> AIC Neg.Bin : 4206.91
cat("NB lebih baik:", AIC(mod_nb) < AIC(mod_poi),"\n")
#> NB lebih baik: TRUE
# Forest plot perbandingan IRR: Poisson vs NB
sm_nb  <- summary(mod_nb)$coefficients
irr_p  <- exp(sm_p[-1, 1])
irr_nb <- exp(sm_nb[-1, 1])
nms    <- rownames(sm_p)[-1]

par(mar=c(4,11,3,2), bg="#f8fafc")
y_p <- seq_len(length(nms))
plot(NULL, xlim=c(0, max(c(irr_p,irr_nb))*1.3), ylim=c(.5,length(nms)+.5),
     xlab="IRR", ylab="", main="IRR: Poisson vs Negative Binomial",
     yaxt="n", las=1, cex.axis=.82)
points(irr_p,  y_p-.15, pch=19, col="#2563eb", cex=1.1)
points(irr_nb, y_p+.15, pch=19, col="#dc2626", cex=1.1)
abline(v=1, lty=2, col="#94a3b8", lwd=1.3)
axis(2, at=y_p, labels=nms, las=1, cex.axis=.75)
legend("topright", legend=c("Poisson","Neg. Binomial"),
       col=c("#2563eb","#dc2626"), pch=19, bty="n", cex=.82)
box(col="#cbd5e1")

5.9 Interpretasi

  • IRR > 1: Prediktor meningkatkan rata-rata jam absen; IRR < 1 bersifat protektif.
  • Jika overdispersion terdeteksi (phi >> 1 atau uji signifikan), Negative Binomial direkomendasikan karena memberikan SE lebih akurat dan inference lebih valid.
  • Quasi-Poisson mengoreksi SE tanpa mengubah estimasi koefisien — pilihan praktis jika distribusi NB tidak konvergen.
  • Variabel dengan CI ekstrem (mis. disciplinary_fYa dengan CI upper hingga 10¹⁵⁰): Ini adalah gejala quasi-complete separation atau sparse data. IRR variabel tersebut tidak dapat diinterpretasikan — hanya dapat dicatat sebagai keterbatasan model. Pada model Negative Binomial, masalah ini umumnya lebih ringan namun tetap perlu diwaspadai.

5.10 Kesimpulan Regresi Poisson

Model Poisson signifikan secara keseluruhan. Pemeriksaan overdispersion menunjukkan bahwa variansi data melebihi mean, sehingga Negative Binomial lebih direkomendasikan. AIC model NB lebih rendah mengkonfirmasi kecocokan yang lebih baik untuk data ketidakhadiran ini.

6 Kesimpulan Umum

6.1 Ringkasan Keempat Model

Tabel 6.1. Ringkasan Perbandingan Keempat Model
Aspek Biner Multinomial Ordinal Poisson
Jenis Respons Biner (0/1) Nominal (3 kat.) Ordinal (3 kat.) Count (jam)
Model Logistik Biner Baseline-Category Logit Proportional Odds Poisson / NB
Link Function Logit Logit/kategori Cumulative Logit Log
Estimasi IRLS Quasi-Newton Fisher Scoring IRLS
Ukuran Efek Odds Ratio RRR = exp(β) OR kumulatif IRR = exp(β)
Asumsi Kunci Linearitas log-odds IIA Proportional Odds Equidispersion
Uji Signifikansi LRT + Wald LRT + Wald LRT + Wald LRT + Wald

6.2 Perbandingan Kegunaan

Model Gunakan Ketika Keunggulan
Logistik Biner Respons 2 kategori Sederhana, OR intuitif
Multinomial Respons nominal ≥ 3 kategori Fleksibel, tanpa asumsi urutan
Ordinal Respons ordinal ≥ 3 kategori Parsimonious, manfaatkan urutan
Poisson/NB Respons count Modelkan count langsung, IRR intuitif

Untuk dataset Absenteeism at Work:

  • Biner → apakah karyawan berisiko absen berlebihan?
  • Multinomial → termasuk kelompok absensi mana? (tanpa asumsi urutan)
  • Ordinal → seberapa tinggi tingkat absensi? (dengan asumsi urutan)
  • Poisson/NB → berapa jam prediksi absen? (skala count)

Referensi: Agresti, A. (2013). Categorical Data Analysis (3rd ed.). Wiley. | Martiniano et al. (2012). CISTI 2012. IEEE. | Cameron & Trivedi (1998). Regression Analysis of Count Data. Cambridge.