Analisis Regresi Logistik: Biner, Multinomial, Ordinal, dan Poisson

# install.packages(c("MASS","nnet","ordinal","broom","knitr","kableExtra",
#                    "scales","tibble","purrr","car","AER","lmtest",
#                    "ResourceSelection","tidyverse"))

suppressPackageStartupMessages({
  library(MASS)              # Pima.tr, housing, quine, polr, stepAIC
  library(nnet)              # multinom
  library(ordinal)           # clm, nominal_test
  library(broom)
  library(knitr)
  library(kableExtra)
  library(scales)
  library(tibble)
  library(purrr)
  library(car)               # vif
  library(AER)               # dispersiontest
  library(lmtest)            # lrtest
  library(ResourceSelection) # Hosmer-Lemeshow
  library(tidyverse)         # load paling akhir
})

theme_analisis <- function(base_size = 12) {
  theme_minimal(base_size = base_size) +
    theme(
      plot.title = element_text(face = "bold", size = base_size + 4, color = "#0B4F2F"),
      plot.subtitle = element_text(size = base_size - 1, color = "#2E7D32"),
      axis.title = element_text(face = "bold", color = "#1B5E20"),
      axis.text = element_text(color = "#263238"),
      panel.grid.minor = element_blank(),
      panel.grid.major = element_line(color = "#DDEEDD", linewidth = 0.35),
      legend.position = "bottom",
      legend.title = element_text(face = "bold", color = "#1B5E20"),
      strip.text = element_text(face = "bold", color = "#0B4F2F"),
      plot.background = element_rect(fill = "#ffffff", color = NA),
      panel.background = element_rect(fill = "#ffffff", color = NA)
    )
}

pal <- c("#0B4F2F", "#1B5E20", "#2E7D32", "#43A047", "#66BB6A", "#A5D6A7", "#C8E6C9")

kable_rapi <- function(x, caption = NULL, digits = 3) {
  knitr::kable(x, caption = caption, digits = digits) %>%
    kableExtra::kable_styling(
      bootstrap_options = c("striped", "hover", "condensed", "responsive"),
      full_width = TRUE,
      position = "center"
    ) %>%
    kableExtra::row_spec(0, background = "#1B5E20", color = "white", bold = TRUE)
}

stratified_split <- function(y, prop = 0.8, seed = 123) {
  set.seed(seed)
  idx_by_class <- split(seq_along(y), y)
  train_idx <- lapply(
    idx_by_class,
    function(idx) sample(idx, size = floor(length(idx) * prop))
  )
  unlist(train_idx, use.names = FALSE)
}

nagelkerke_r2 <- function(fit, null, n) {
  as.numeric((1 - exp((logLik(null) - logLik(fit)) * (2 / n))) /
               (1 - exp(logLik(null) * (2 / n))))
}

mcfadden_r2 <- function(fit, null) {
  as.numeric(1 - (logLik(fit) / logLik(null)))
}

akurasi <- function(actual, pred) {
  actual <- as.character(actual)
  pred <- as.character(pred)
  mean(actual == pred, na.rm = TRUE)
}

1 Regresi Logistik Biner

1.1 Deskripsi Data

Dataset yang digunakan pada bagian ini adalah Pima.tr dari package MASS. Dataset ini merupakan data nyata mengenai perempuan keturunan Pima Indian yang digunakan untuk mempelajari status diabetes berdasarkan beberapa variabel klinis. Variabel respons yang digunakan adalah type, yaitu status diabetes dengan dua kategori, No dan Yes. Karena variabel respons hanya memiliki dua kategori, model yang digunakan adalah regresi logistik biner.

Variabel dependen: type

• No = tidak diabetes
• Yes = diabetes

Variabel prediktor:

data(Pima.tr, package = "MASS")

tibble(
  Variabel = names(Pima.tr),
  Keterangan = c(
    "Jumlah kehamilan",
    "Konsentrasi glukosa plasma",
    "Tekanan darah diastolik",
    "Ketebalan lipatan kulit trisep",
    "Indeks massa tubuh",
    "Riwayat/pedigree diabetes",
    "Usia",
    "Status diabetes"
  )
) %>%
  kable_rapi(caption = "Keterangan variabel pada dataset Pima.tr")

Keterangan variabel pada dataset Pima.tr

Variabel	Keterangan
npreg	Jumlah kehamilan
glu	Konsentrasi glukosa plasma
bp	Tekanan darah diastolik
skin	Ketebalan lipatan kulit trisep
bmi	Indeks massa tubuh
ped	Riwayat/pedigree diabetes
age	Usia
type	Status diabetes

1.2 Load Library dan Data

raw_pima <- Pima.tr
glimpse(raw_pima)

## Rows: 200
## Columns: 8
## $ npreg <int> 5, 7, 5, 0, 0, 5, 3, 1, 3, 2, 0, 9, 1, 12, 1, 4, 1, 11, 1, 0, 2,…
## $ glu   <int> 86, 195, 77, 165, 107, 97, 83, 193, 142, 128, 137, 154, 189, 92,…
## $ bp    <int> 68, 70, 82, 76, 60, 76, 58, 50, 80, 78, 40, 78, 60, 62, 66, 76, …
## $ skin  <int> 28, 33, 41, 43, 25, 27, 31, 16, 15, 37, 35, 30, 23, 7, 52, 15, 8…
## $ bmi   <dbl> 30.2, 25.1, 35.8, 47.9, 26.4, 35.6, 34.3, 25.9, 32.4, 43.3, 43.1…
## $ ped   <dbl> 0.364, 0.163, 0.156, 0.259, 0.133, 0.378, 0.336, 0.655, 0.200, 1…
## $ age   <int> 24, 55, 35, 26, 23, 52, 25, 24, 63, 31, 33, 45, 59, 44, 29, 21, …
## $ type  <fct> No, Yes, No, No, No, Yes, No, No, No, Yes, Yes, No, Yes, Yes, No…

Interpretasi output: Dataset Pima.tr terdiri dari 200 observasi dan 8 variabel. Variabel type merupakan variabel kategorik dua kelas, sedangkan prediktor lainnya berupa variabel numerik.

1.3 Persiapan Data

1.3.1 Preprocessing dan Pemberian Label

pima <- raw_pima %>%
  mutate(
    type = factor(type, levels = c("No", "Yes")),
    diabetes_bin = ifelse(type == "Yes", 1, 0)
  )

glimpse(pima)

## Rows: 200
## Columns: 9
## $ npreg        <int> 5, 7, 5, 0, 0, 5, 3, 1, 3, 2, 0, 9, 1, 12, 1, 4, 1, 11, 1…
## $ glu          <int> 86, 195, 77, 165, 107, 97, 83, 193, 142, 128, 137, 154, 1…
## $ bp           <int> 68, 70, 82, 76, 60, 76, 58, 50, 80, 78, 40, 78, 60, 62, 6…
## $ skin         <int> 28, 33, 41, 43, 25, 27, 31, 16, 15, 37, 35, 30, 23, 7, 52…
## $ bmi          <dbl> 30.2, 25.1, 35.8, 47.9, 26.4, 35.6, 34.3, 25.9, 32.4, 43.…
## $ ped          <dbl> 0.364, 0.163, 0.156, 0.259, 0.133, 0.378, 0.336, 0.655, 0…
## $ age          <int> 24, 55, 35, 26, 23, 52, 25, 24, 63, 31, 33, 45, 59, 44, 2…
## $ type         <fct> No, Yes, No, No, No, Yes, No, No, No, Yes, Yes, No, Yes, …
## $ diabetes_bin <dbl> 0, 1, 0, 0, 0, 1, 0, 0, 0, 1, 1, 0, 1, 1, 0, 0, 0, 1, 1, …

1.3.2 Pengecekan Data

cat("Missing values per kolom:\n")

## Missing values per kolom:

colSums(is.na(pima))

##        npreg          glu           bp         skin          bmi          ped 
##            0            0            0            0            0            0 
##          age         type diabetes_bin 
##            0            0            0

dist_pima <- pima %>%
  count(type) %>%
  mutate(Proporsi = round(n / sum(n), 4))

dist_pima %>%
  kable_rapi(caption = "Distribusi variabel dependen")

Distribusi variabel dependen

type	n	Proporsi
No	132	0.66
Yes	68	0.34

Interpretasi output: Pengecekan missing value dilakukan untuk memastikan bahwa data sudah siap dianalisis. Distribusi status diabetes digunakan untuk melihat keseimbangan kelas antara kelompok No dan Yes.

1.4 Statistik Deskriptif

1.4.1 Distribusi Variabel Prediktor

pima %>%
  dplyr::select(npreg, glu, bp, skin, bmi, ped, age) %>%
  summary()

##      npreg            glu              bp              skin      
##  Min.   : 0.00   Min.   : 56.0   Min.   : 38.00   Min.   : 7.00  
##  1st Qu.: 1.00   1st Qu.:100.0   1st Qu.: 64.00   1st Qu.:20.75  
##  Median : 2.00   Median :120.5   Median : 70.00   Median :29.00  
##  Mean   : 3.57   Mean   :124.0   Mean   : 71.26   Mean   :29.21  
##  3rd Qu.: 6.00   3rd Qu.:144.0   3rd Qu.: 78.00   3rd Qu.:36.00  
##  Max.   :14.00   Max.   :199.0   Max.   :110.00   Max.   :99.00  
##       bmi             ped              age       
##  Min.   :18.20   Min.   :0.0850   Min.   :21.00  
##  1st Qu.:27.57   1st Qu.:0.2535   1st Qu.:23.00  
##  Median :32.80   Median :0.3725   Median :28.00  
##  Mean   :32.31   Mean   :0.4608   Mean   :32.11  
##  3rd Qu.:36.50   3rd Qu.:0.6160   3rd Qu.:39.25  
##  Max.   :47.90   Max.   :2.2880   Max.   :63.00

Berdasarkan statistik deskriptif, karakteristik responden pada dataset Pima.tr menunjukkan variasi yang cukup besar pada setiap variabel prediktor. Rata-rata jumlah kehamilan (npreg) adalah 3,57 kali, sedangkan rata-rata kadar glukosa plasma (glu) sebesar 124 mg/dL dengan rentang 56–199 mg/dL. Rata-rata tekanan darah diastolik (bp) tercatat 71,26 mmHg, ketebalan lipatan kulit (skin) 29,21 mm, dan indeks massa tubuh (bmi) 32,31 kg/m² yang mengindikasikan sebagian responden berada pada kategori kelebihan berat badan hingga obesitas. Selain itu, nilai rata-rata Diabetes Pedigree Function (ped) sebesar 0,4608 menunjukkan adanya variasi tingkat risiko diabetes berdasarkan riwayat keluarga, sementara rata-rata usia responden adalah 32,11 tahun dengan rentang 21–63 tahun. Secara umum, rentang nilai yang cukup lebar pada variabel glukosa, BMI, dan usia menunjukkan heterogenitas karakteristik responden sehingga variabel-variabel tersebut berpotensi memberikan kontribusi penting dalam memprediksi status diabetes pada model regresi logistik biner.

1.4.2 Tabulasi Silang: Proporsi Diabetes per Kategori

Agar tabulasi silang dapat dibuat untuk prediktor numerik, setiap prediktor dikelompokkan menjadi tiga kategori berdasarkan kuartil.

pima_cat <- pima %>%
  mutate(
    glu_cat  = cut(glu, breaks = quantile(glu, probs = c(0, .33, .67, 1), na.rm = TRUE),
                   include.lowest = TRUE, labels = c("Rendah", "Sedang", "Tinggi")),
    bmi_cat  = cut(bmi, breaks = quantile(bmi, probs = c(0, .33, .67, 1), na.rm = TRUE),
                   include.lowest = TRUE, labels = c("Rendah", "Sedang", "Tinggi")),
    age_cat  = cut(age, breaks = quantile(age, probs = c(0, .33, .67, 1), na.rm = TRUE),
                   include.lowest = TRUE, labels = c("Muda", "Dewasa", "Lebih Tua"))
  )

tab_biner <- function(var, label) {
  pima_cat %>%
    count(!!sym(var), type) %>%
    group_by(!!sym(var)) %>%
    mutate(Proporsi = round(n / sum(n), 3)) %>%
    filter(type == "Yes") %>%
    dplyr::select(Kategori = !!sym(var), n_Diabetes = n, Proporsi_Diabetes = Proporsi) %>%
    kable_rapi(caption = paste("Proporsi diabetes menurut", label))
}

tab_biner("glu_cat", "kategori glukosa")

Proporsi diabetes menurut kategori glukosa

Kategori	n_Diabetes	Proporsi_Diabetes
Rendah	8	0.121
Sedang	20	0.286
Tinggi	40	0.625

tab_biner("bmi_cat", "kategori BMI")

Proporsi diabetes menurut kategori BMI

Kategori	n_Diabetes	Proporsi_Diabetes
Rendah	9	0.136
Sedang	28	0.412
Tinggi	31	0.470

tab_biner("age_cat", "kategori usia")

Proporsi diabetes menurut kategori usia

Kategori	n_Diabetes	Proporsi_Diabetes
Muda	11	0.162
Dewasa	21	0.309
Lebih Tua	36	0.562

Proporsi penderita diabetes cenderung meningkat seiring bertambahnya usia. Kelompok Muda memiliki proporsi diabetes sebesar 16,2% (11 orang), kelompok Dewasa sebesar 30,9% (21 orang), sedangkan kelompok Lebih Tua memiliki proporsi tertinggi yaitu 56,2% (36 orang). Pola ini menunjukkan adanya hubungan positif antara usia dan kejadian diabetes, di mana individu yang lebih tua cenderung memiliki risiko diabetes yang lebih tinggi dibandingkan kelompok usia yang lebih muda. Dengan demikian, variabel usia berpotensi menjadi salah satu prediktor penting dalam model regresi logistik biner untuk memprediksi status diabetes.

1.5 Eksplorasi Visual

ggplot(pima, aes(x = type, fill = type)) +
  geom_bar(alpha = 0.85) +
  scale_fill_manual(values = pal[1:2]) +
  labs(title = "Distribusi Status Diabetes",
       x = "Status Diabetes", y = "Frekuensi", fill = "Status") +
  theme_analisis()

Grafik menunjukkan distribusi status diabetes pada dataset Pima.tr. Jumlah responden yang tidak menderita diabetes (No) lebih banyak dibandingkan responden yang menderita diabetes (Yes). Terlihat bahwa kelompok No memiliki frekuensi sekitar 130 observasi, sedangkan kelompok Yes sekitar 70 observasi. Hal ini menunjukkan bahwa data memiliki ketidakseimbangan kelas (class imbalance) ringan, di mana proporsi penderita diabetes lebih kecil dibandingkan non-diabetes. Meskipun demikian, jumlah kasus diabetes masih cukup besar sehingga model regresi logistik biner tetap dapat dibangun dan digunakan untuk mengidentifikasi faktor-faktor yang berpengaruh terhadap status diabetes.

ggplot(pima, aes(x = type, y = glu, fill = type)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = pal[1:2]) +
  labs(title = "Perbandingan Glukosa berdasarkan Status Diabetes",
       x = "Status Diabetes", y = "Glukosa", fill = "Status") +
  theme_analisis()

Kelompok responden yang menderita diabetes (Yes) memiliki kadar glukosa yang cenderung lebih tinggi dibandingkan kelompok non-diabetes (No). Perbedaan median yang cukup jelas menunjukkan bahwa variabel glukosa (glu) berpotensi menjadi prediktor penting dalam memprediksi status diabetes.

ggplot(pima, aes(x = type, y = bmi, fill = type)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = pal[3:4]) +
  labs(title = "Perbandingan BMI berdasarkan Status Diabetes",
       x = "Status Diabetes", y = "BMI", fill = "Status") +
  theme_analisis()

Kelompok responden yang menderita diabetes (Yes) memiliki nilai BMI yang cenderung lebih tinggi dibandingkan kelompok non-diabetes (No). Hal ini menunjukkan bahwa BMI berpotensi menjadi salah satu faktor yang berhubungan dengan risiko diabetes.

1.6 Pembagian Training dan Testing

Data dibagi menjadi 80% training dan 20% testing secara stratified agar proporsi kelas diabetes tetap seimbang.

set.seed(42)
train_id_pima <- stratified_split(pima$type, prop = 0.8, seed = 42)
train_pima <- pima[train_id_pima, ]
test_pima  <- pima[-train_id_pima, ]

split_pima <- bind_rows(
  train_pima %>% count(type) %>% mutate(Data = "Training"),
  test_pima %>% count(type) %>% mutate(Data = "Testing")
) %>%
  group_by(Data) %>%
  mutate(Proporsi = scales::percent(n / sum(n), accuracy = 0.1)) %>%
  ungroup() %>%
  dplyr::select(Data, Status = type, Jumlah = n, Proporsi)

split_pima %>%
  kable_rapi(caption = "Distribusi kelas pada training dan testing")

Distribusi kelas pada training dan testing

Data	Status	Jumlah	Proporsi
Training	No	105	66.0%
Training	Yes	54	34.0%
Testing	No	27	65.9%
Testing	Yes	14	34.1%

1.7 Model Regresi Logistik Biner

1.7.1 Rumus Model Logistik

Misalkan \(Y_i\) adalah status diabetes pada individu ke-\(i\), dengan:

\[ Y_i = \begin{cases} 1, & \text{jika diabetes} \\ 0, & \text{jika tidak diabetes} \end{cases} \]

Peluang diabetes dinotasikan sebagai:

\[ p_i = P(Y_i = 1 \mid X_i) \]

Model regresi logistik biner dituliskan sebagai:

\[ \log \left(\frac{p_i}{1-p_i}\right) = \beta_0 + \beta_1X_{1i} + \beta_2X_{2i} + \cdots + \beta_kX_{ki} \]

1.7.2 Interpretasi Odds Ratio

Odds Ratio (OR) dihitung dengan:

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

Jika \(OR_j > 1\), maka peningkatan prediktor berkaitan dengan kenaikan odds diabetes. Jika \(OR_j < 1\), maka peningkatan prediktor berkaitan dengan penurunan odds diabetes.

1.7.3 Model Penuh

fit_bin_full <- glm(
  diabetes_bin ~ npreg + glu + bp + skin + bmi + ped + age,
  data = train_pima,
  family = binomial(link = "logit")
)

ringkasan_bin <- data.frame(
  Keterangan = c("Jumlah observasi training", "Null deviance",
                 "Residual deviance", "Derajat bebas residual", "AIC"),
  Nilai = c(
    nobs(fit_bin_full),
    round(fit_bin_full$null.deviance, 3),
    round(fit_bin_full$deviance, 3),
    fit_bin_full$df.residual,
    round(AIC(fit_bin_full), 3)
  )
)

ringkasan_bin %>%
  kable_rapi(caption = "Ringkasan kecocokan model penuh")

Ringkasan kecocokan model penuh

Keterangan	Nilai
Jumlah observasi training	159.000
Null deviance	203.770
Residual deviance	139.801
Derajat bebas residual	151.000
AIC	155.801

Model regresi logistik biner diestimasi menggunakan 159 observasi data training. Nilai residual deviance sebesar 139,801 lebih kecil dibandingkan null deviance sebesar 203,770, yang menunjukkan bahwa variabel-variabel prediktor yang digunakan mampu meningkatkan kecocokan model dibandingkan model tanpa prediktor (model null). Selain itu, nilai AIC sebesar 155,801 dapat digunakan sebagai acuan dalam pemilihan model, di mana nilai AIC yang lebih kecil menunjukkan model yang lebih baik dalam menyeimbangkan antara kecocokan model dan kompleksitas parameter yang digunakan.

or_bin_full <- broom::tidy(fit_bin_full) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_low = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error)
  ) %>%
  arrange(p.value) %>%
  transmute(
    Variabel = term,
    `Odds Ratio` = round(odds_ratio, 3),
    `IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
    `p-value` = signif(p.value, 3)
  )

or_bin_full %>%
  kable_rapi(caption = "Odds ratio dan p-value — Model Penuh")

Odds ratio dan p-value — Model Penuh

Variabel	Odds Ratio	IK 95%	p-value
glu	1.035	1.02 - 1.051	0.000
bmi	1.135	1.031 - 1.25	0.010
age	1.064	1.011 - 1.121	0.018
ped	3.073	0.662 - 14.255	0.152
skin	0.969	0.924 - 1.015	0.183
npreg	1.066	0.925 - 1.229	0.375
bp	0.992	0.954 - 1.031	0.687

Berdasarkan hasil estimasi model regresi logistik biner, variabel glukosa (glu), BMI (bmi), dan usia (age) berpengaruh signifikan terhadap status diabetes karena memiliki p-value < 0,05. Variabel glukosa memiliki pengaruh paling signifikan (OR = 1,035; p < 0,001), yang berarti setiap kenaikan 1 satuan kadar glukosa meningkatkan odds terjadinya diabetes sebesar 3,5%, dengan asumsi variabel lain konstan. Variabel BMI (OR = 1,135; p = 0,010) menunjukkan bahwa setiap kenaikan 1 satuan BMI meningkatkan odds diabetes sebesar 13,5%, sedangkan usia (OR = 1,064; p = 0,018) meningkatkan odds diabetes sebesar 6,4% untuk setiap tambahan satu tahun usia. Sementara itu, variabel ped, skin, npreg, dan bp tidak berpengaruh signifikan secara statistik karena memiliki p-value lebih besar dari 0,05. Dengan demikian, glukosa, BMI, dan usia merupakan prediktor utama yang berkontribusi terhadap kemungkinan seseorang menderita diabetes pada model ini.

1.8 Seleksi Variabel: Stepwise Backward (AIC)

cat("== Stepwise Backward Selection (AIC) ==\n")

## == Stepwise Backward Selection (AIC) ==

fit_bin_red <- stepAIC(fit_bin_full, direction = "backward", trace = TRUE)

## Start:  AIC=155.8
## diabetes_bin ~ npreg + glu + bp + skin + bmi + ped + age
## 
##         Df Deviance    AIC
## - bp     1   139.96 153.96
## - npreg  1   140.59 154.59
## - skin   1   141.45 155.45
## <none>       139.80 155.80
## - ped    1   141.83 155.83
## - age    1   145.74 159.74
## - bmi    1   146.69 160.69
## - glu    1   165.77 179.77
## 
## Step:  AIC=153.96
## diabetes_bin ~ npreg + glu + skin + bmi + ped + age
## 
##         Df Deviance    AIC
## - npreg  1   140.78 152.78
## - skin   1   141.72 153.72
## <none>       139.96 153.96
## - ped    1   141.99 153.99
## - age    1   145.78 157.78
## - bmi    1   146.73 158.73
## - glu    1   165.84 177.84
## 
## Step:  AIC=152.78
## diabetes_bin ~ glu + skin + bmi + ped + age
## 
##        Df Deviance    AIC
## - ped   1   142.55 152.55
## <none>      140.78 152.78
## - skin  1   142.80 152.80
## - bmi   1   147.60 157.60
## - age   1   153.45 163.45
## - glu   1   166.35 176.35
## 
## Step:  AIC=152.55
## diabetes_bin ~ glu + skin + bmi + age
## 
##        Df Deviance    AIC
## <none>      142.55 152.55
## - skin  1   144.61 152.61
## - bmi   1   150.46 158.46
## - age   1   154.56 162.56
## - glu   1   169.15 177.15

cat("\n== Formula Model Reduksi (hasil stepAIC) ==\n")

## 
## == Formula Model Reduksi (hasil stepAIC) ==

print(formula(fit_bin_red))

## diabetes_bin ~ glu + skin + bmi + age

ringkasan_bin_2 <- data.frame(
  Keterangan = c("Null deviance", "Residual deviance",
                 "Derajat bebas residual", "AIC"),
  `Model Penuh` = c(
    round(fit_bin_full$null.deviance, 3),
    round(fit_bin_full$deviance, 3),
    fit_bin_full$df.residual,
    round(AIC(fit_bin_full), 3)
  ),
  `Model Reduksi` = c(
    round(fit_bin_red$null.deviance, 3),
    round(fit_bin_red$deviance, 3),
    fit_bin_red$df.residual,
    round(AIC(fit_bin_red), 3)
  ),
  check.names = FALSE
)

ringkasan_bin_2 %>%
  kable_rapi(caption = "Perbandingan ringkasan Model Penuh vs Reduksi")

Perbandingan ringkasan Model Penuh vs Reduksi

Keterangan	Model Penuh	Model Reduksi
Null deviance	203.770	203.770
Residual deviance	139.801	142.549
Derajat bebas residual	151.000	154.000
AIC	155.801	152.549

Hasil stepwise backward selection menghasilkan model reduksi dengan nilai AIC sebesar 152,549, lebih kecil dibandingkan model penuh yang memiliki AIC sebesar 155,801. Meskipun nilai residual deviance model reduksi sedikit meningkat dari 139,801 menjadi 142,549, peningkatan tersebut relatif kecil dibandingkan pengurangan jumlah parameter yang digunakan. Selain itu, kedua model memiliki null deviance yang sama karena berasal dari data yang sama. Dengan nilai AIC yang lebih rendah, model reduksi dianggap lebih parsimonious, yaitu mampu menjelaskan data dengan baik menggunakan variabel yang lebih sedikit. Oleh karena itu, model reduksi lebih layak dipilih sebagai model final karena memberikan keseimbangan yang lebih baik antara kecocokan model dan kompleksitas model

1.9 Uji Asumsi

1.9.1 Multikolinearitas — VIF

cat("== VIF Model Penuh ==\n")

## == VIF Model Penuh ==

vif(fit_bin_full)

##    npreg      glu       bp     skin      bmi      ped      age 
## 1.504630 1.084280 1.205725 1.846400 1.784990 1.060416 1.725049

cat("\n== VIF Model Reduksi ==\n")

## 
## == VIF Model Reduksi ==

vif(fit_bin_red)

##      glu     skin      bmi      age 
## 1.040253 1.799830 1.737625 1.097124

Hasil pengujian multikolinearitas menunjukkan bahwa seluruh nilai VIF pada model penuh maupun model reduksi berada jauh di bawah batas umum 5. Pada model penuh, nilai VIF berkisar antara 1,060 hingga 1,846, sedangkan pada model reduksi berkisar antara 1,040 hingga 1,800. Nilai VIF yang rendah ini menunjukkan bahwa tidak terdapat hubungan linear yang kuat antar variabel prediktor, sehingga masalah multikolinearitas tidak terdeteksi dalam model. Dengan demikian, setiap variabel prediktor memberikan informasi yang relatif independen dan estimasi koefisien model dapat dianggap stabil serta dapat diinterpretasikan dengan baik.

1.9.2 Likelihood Ratio Test (LRT)

null_bin <- glm(diabetes_bin ~ 1, data = train_pima, family = binomial)

cat("== LRT: Model Penuh vs Null ==\n")

## == LRT: Model Penuh vs Null ==

lrtest(null_bin, fit_bin_full)

cat("\n== LRT: Model Reduksi vs Null ==\n")

## 
## == LRT: Model Reduksi vs Null ==

lrtest(null_bin, fit_bin_red)

cat("\n== LRT: Model Reduksi vs Model Penuh ==\n")

## 
## == LRT: Model Reduksi vs Model Penuh ==

lrtest(fit_bin_red, fit_bin_full)

Hasil LRT menunjukkan bahwa model penuh dan model reduksi sama-sama lebih baik dibandingkan model null (p-value < 0,001). Selain itu, perbandingan model reduksi dengan model penuh menghasilkan p-value = 0,4321 (> 0,05), sehingga variabel yang dihilangkan tidak memberikan peningkatan yang signifikan. Oleh karena itu, model reduksi dapat dipilih karena lebih sederhana namun tetap memiliki performa yang baik.

1.9.3 Hosmer-Lemeshow Goodness-of-Fit

cat("== Hosmer-Lemeshow — Model Penuh ==\n")

## == Hosmer-Lemeshow — Model Penuh ==

hoslem.test(train_pima$diabetes_bin, fitted(fit_bin_full), g = 10)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  train_pima$diabetes_bin, fitted(fit_bin_full)
## X-squared = 11.985, df = 8, p-value = 0.1519

cat("\n== Hosmer-Lemeshow — Model Reduksi ==\n")

## 
## == Hosmer-Lemeshow — Model Reduksi ==

hoslem.test(train_pima$diabetes_bin, fitted(fit_bin_red), g = 10)

## 
##  Hosmer and Lemeshow goodness of fit (GOF) test
## 
## data:  train_pima$diabetes_bin, fitted(fit_bin_red)
## X-squared = 14.513, df = 8, p-value = 0.06934

Hasil uji Hosmer-Lemeshow pada model penuh menghasilkan p-value = 0,1519 dan pada model reduksi p-value = 0,0693. Karena kedua nilai p-value lebih besar dari 0,05, maka tidak terdapat perbedaan yang signifikan antara nilai observasi dan prediksi model. Dengan demikian, baik model penuh maupun model reduksi dapat dikatakan fit dan sesuai untuk digunakan dalam memodelkan status diabetes.

1.9.4 Linearitas Log-Odds (Box-Tidwell)

# Box-Tidwell dilakukan pada prediktor kontinu.
# Untuk menghindari log(0), setiap prediktor ditambah konstanta kecil.
bt_data <- train_pima %>%
  mutate(across(c(npreg, glu, bp, skin, bmi, ped, age), ~ .x + 0.01, .names = "{.col}_pos")) %>%
  mutate(
    npreg_log = npreg_pos * log(npreg_pos),
    glu_log   = glu_pos   * log(glu_pos),
    bp_log    = bp_pos    * log(bp_pos),
    skin_log  = skin_pos  * log(skin_pos),
    bmi_log   = bmi_pos   * log(bmi_pos),
    ped_log   = ped_pos   * log(ped_pos),
    age_log   = age_pos   * log(age_pos)
  )

fit_bt <- glm(
  diabetes_bin ~ npreg + glu + bp + skin + bmi + ped + age +
    npreg_log + glu_log + bp_log + skin_log + bmi_log + ped_log + age_log,
  data = bt_data,
  family = binomial
)

broom::tidy(fit_bt) %>%
  filter(str_detect(term, "_log")) %>%
  transmute(Variabel = term, `p-value` = signif(p.value, 3)) %>%
  kable_rapi(caption = "Uji Box-Tidwell untuk linearitas log-odds")

Uji Box-Tidwell untuk linearitas log-odds

Variabel	p-value
npreg_log	0.201
glu_log	0.623
bp_log	0.366
skin_log	0.680
bmi_log	0.029
ped_log	0.508
age_log	0.508

Hasil uji Box-Tidwell menunjukkan bahwa sebagian besar variabel memiliki p-value > 0,05, yaitu npreg, glu, bp, skin, ped, dan age, sehingga asumsi linearitas antara prediktor kontinu dan logit dapat dianggap terpenuhi. Namun, variabel bmi memiliki p-value = 0,029 (< 0,05), yang mengindikasikan adanya penyimpangan dari asumsi linearitas logit. Meskipun demikian, karena hanya satu variabel yang menunjukkan pelanggaran dan nilainya tidak terlalu jauh dari batas signifikansi, model masih dapat digunakan dengan catatan bahwa hubungan BMI terhadap log-odds diabetes mungkin tidak sepenuhnya linear.

1.9.5 Nagelkerke \(R^2\)

n_bin <- nrow(train_pima)

data.frame(
  Model = c("Model Penuh", "Model Reduksi"),
  `Nagelkerke R²` = c(
    nagelkerke_r2(fit_bin_full, null_bin, n_bin),
    nagelkerke_r2(fit_bin_red, null_bin, n_bin)
  ),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Nagelkerke R² kedua model")

Nagelkerke R² kedua model

Model	Nagelkerke R²
Model Penuh	0.459
Model Reduksi	0.442

Nilai Nagelkerke R² pada model penuh sebesar 0,459, sedangkan pada model reduksi sebesar 0,442. Hal ini menunjukkan bahwa model penuh mampu menjelaskan sekitar 45,9% variasi status diabetes, sementara model reduksi menjelaskan sekitar 44,2% variasi. Perbedaan nilai yang relatif kecil menunjukkan bahwa model reduksi tetap memiliki kemampuan penjelasan yang hampir sama dengan model penuh meskipun menggunakan lebih sedikit variabel

1.10 Odds Ratio Model Final (Model Reduksi)

or_bin_final <- broom::tidy(fit_bin_red) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_low = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error)
  ) %>%
  arrange(p.value) %>%
  transmute(
    Variabel = term,
    `Odds Ratio` = round(odds_ratio, 3),
    `IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
    `p-value` = signif(p.value, 3)
  )

or_bin_final %>%
  kable_rapi(caption = "Odds Ratio — Model Final")

Odds Ratio — Model Final

Variabel	Odds Ratio	IK 95%	p-value
glu	1.035	1.02 - 1.05	0.000
age	1.072	1.029 - 1.117	0.001
bmi	1.146	1.04 - 1.262	0.006
skin	0.965	0.922 - 1.011	0.135

Berdasarkan model final, variabel glukosa (glu), usia (age), dan BMI (bmi) berpengaruh signifikan terhadap status diabetes (p-value < 0,05). Variabel glukosa memiliki OR = 1,035 yang berarti setiap kenaikan 1 satuan glukosa meningkatkan odds diabetes sebesar 3,5%. Variabel usia memiliki OR = 1,072 sehingga setiap tambahan 1 tahun usia meningkatkan odds diabetes sebesar 7,2%, sedangkan BMI memiliki OR = 1,146 yang menunjukkan bahwa setiap kenaikan 1 satuan BMI meningkatkan odds diabetes sebesar 14,6%. Sementara itu, variabel skin tidak signifikan (p-value = 0,135), sehingga belum terdapat bukti yang cukup bahwa ketebalan lipatan kulit berpengaruh terhadap status diabetes pada model final.

1.11 Prediksi Probabilitas

prob_bin <- predict(fit_bin_red, newdata = test_pima, type = "response")
pred_bin <- ifelse(prob_bin >= 0.5, "Yes", "No") %>%
  factor(levels = c("No", "Yes"))

cm_bin <- table(Aktual = test_pima$type, Prediksi = pred_bin)
cm_bin

##       Prediksi
## Aktual No Yes
##    No  20   7
##    Yes  5   9

akurasi_bin <- akurasi(test_pima$type, pred_bin)

data.frame(
  Ukuran = "Akurasi",
  Nilai = round(akurasi_bin, 4)
) %>%
  kable_rapi(caption = "Akurasi prediksi model logistik biner")

Akurasi prediksi model logistik biner

Ukuran	Nilai
Akurasi	0.707

Nilai akurasi sebesar 70,7% menunjukkan bahwa model mampu memprediksi status diabetes dengan benar pada sekitar 7 dari 10 observasi data pengujian, sehingga model memiliki kemampuan klasifikasi yang cukup baik.

1.12 Kesimpulan

Berdasarkan hasil analisis regresi logistik biner menggunakan dataset Pima.tr, diperoleh bahwa model reduksi yang terdiri dari variabel glukosa (glu), BMI (bmi), usia (age), dan ketebalan lipatan kulit (skin) merupakan model yang lebih parsimonious dengan nilai AIC lebih rendah dibandingkan model penuh. Hasil pengujian menunjukkan bahwa variabel glukosa, BMI, dan usia berpengaruh signifikan terhadap status diabetes, sedangkan variabel skin tidak signifikan secara statistik. Uji asumsi menunjukkan tidak terdapat masalah multikolinearitas, model memiliki kecocokan yang baik berdasarkan uji Hosmer-Lemeshow, serta mampu menjelaskan sekitar 44,2% variasi status diabetes berdasarkan nilai Nagelkerke R². Selain itu, model menghasilkan tingkat akurasi sebesar 70,7% pada data pengujian, yang menunjukkan kemampuan klasifikasi yang cukup baik. Dengan demikian, kadar glukosa, BMI, dan usia merupakan faktor utama yang berkontribusi terhadap kemungkinan seseorang menderita diabetes pada dataset ini.

2 Regresi Logistik Multinomial

2.1 Deskripsi Data

Dataset yang digunakan adalah iris dari package datasets. Dataset ini merupakan data nyata yang memuat pengukuran morfologi bunga iris dari tiga spesies, yaitu setosa, versicolor, dan virginica. Variabel respons adalah Species, yang memiliki tiga kategori nominal. Karena kategori respons lebih dari dua dan tidak memiliki urutan alami, metode yang digunakan adalah regresi logistik multinomial.

Variabel dependen: Species

• setosa
• versicolor
• virginica

Variabel prediktor:

tibble(
  Variabel = names(iris),
  Keterangan = c(
    "Panjang sepal",
    "Lebar sepal",
    "Panjang petal",
    "Lebar petal",
    "Spesies bunga iris"
  )
) %>%
  kable_rapi(caption = "Keterangan variabel pada dataset iris")

Keterangan variabel pada dataset iris

Variabel	Keterangan
Sepal.Length	Panjang sepal
Sepal.Width	Lebar sepal
Petal.Length	Panjang petal
Petal.Width	Lebar petal
Species	Spesies bunga iris

2.2 Load Library dan Data

raw_iris <- iris
glimpse(raw_iris)

## Rows: 150
## Columns: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
## $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
## $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…
## $ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s…

Interpretasi output: Dataset iris terdiri dari 150 observasi dan 5 variabel. Variabel Species merupakan variabel kategorik nominal dengan tiga kelas.

2.3 Persiapan Data

2.3.1 Preprocessing dan Pemberian Label

iris_dat <- raw_iris %>%
  mutate(Species = factor(Species))

glimpse(iris_dat)

## Rows: 150
## Columns: 5
## $ Sepal.Length <dbl> 5.1, 4.9, 4.7, 4.6, 5.0, 5.4, 4.6, 5.0, 4.4, 4.9, 5.4, 4.…
## $ Sepal.Width  <dbl> 3.5, 3.0, 3.2, 3.1, 3.6, 3.9, 3.4, 3.4, 2.9, 3.1, 3.7, 3.…
## $ Petal.Length <dbl> 1.4, 1.4, 1.3, 1.5, 1.4, 1.7, 1.4, 1.5, 1.4, 1.5, 1.5, 1.…
## $ Petal.Width  <dbl> 0.2, 0.2, 0.2, 0.2, 0.2, 0.4, 0.3, 0.2, 0.2, 0.1, 0.2, 0.…
## $ Species      <fct> setosa, setosa, setosa, setosa, setosa, setosa, setosa, s…

2.3.2 Pengecekan Data

cat("Missing values per kolom:\n")

## Missing values per kolom:

colSums(is.na(iris_dat))

## Sepal.Length  Sepal.Width Petal.Length  Petal.Width      Species 
##            0            0            0            0            0

iris_dat %>%
  count(Species) %>%
  mutate(Proporsi = round(n / sum(n), 4)) %>%
  kable_rapi(caption = "Distribusi variabel dependen")

Distribusi variabel dependen

Species	n	Proporsi
setosa	50	0.333
versicolor	50	0.333
virginica	50	0.333

2.4 Statistik Deskriptif

2.4.1 Distribusi Variabel Prediktor

iris_dat %>%
  dplyr::select(Sepal.Length, Sepal.Width, Petal.Length, Petal.Width) %>%
  summary()

##   Sepal.Length    Sepal.Width     Petal.Length    Petal.Width   
##  Min.   :4.300   Min.   :2.000   Min.   :1.000   Min.   :0.100  
##  1st Qu.:5.100   1st Qu.:2.800   1st Qu.:1.600   1st Qu.:0.300  
##  Median :5.800   Median :3.000   Median :4.350   Median :1.300  
##  Mean   :5.843   Mean   :3.057   Mean   :3.758   Mean   :1.199  
##  3rd Qu.:6.400   3rd Qu.:3.300   3rd Qu.:5.100   3rd Qu.:1.800  
##  Max.   :7.900   Max.   :4.400   Max.   :6.900   Max.   :2.500

Berdasarkan statistik deskriptif dataset iris, rata-rata panjang sepal (Sepal.Length) adalah 5,843 cm dengan rentang 4,3–7,9 cm, sedangkan rata-rata lebar sepal (Sepal.Width) sebesar 3,057 cm dengan rentang 2,0–4,4 cm. Untuk bagian mahkota bunga, rata-rata panjang petal (Petal.Length) adalah 3,758 cm dengan rentang 1,0–6,9 cm, dan rata-rata lebar petal (Petal.Width) sebesar 1,199 cm dengan rentang 0,1–2,5 cm. Variabel Petal.Length dan Petal.Width memiliki rentang nilai yang relatif lebih lebar dibandingkan variabel sepal, yang mengindikasikan adanya perbedaan karakteristik yang cukup jelas antar spesies bunga iris.

2.4.2 Tabulasi Silang: Proporsi Species per Kategori

iris_cat <- iris_dat %>%
  mutate(
    PetalLength_cat = cut(Petal.Length, breaks = 3, labels = c("Pendek", "Sedang", "Panjang")),
    PetalWidth_cat  = cut(Petal.Width, breaks = 3, labels = c("Kecil", "Sedang", "Besar"))
  )

iris_cat %>%
  count(PetalLength_cat, Species) %>%
  group_by(PetalLength_cat) %>%
  mutate(Proporsi = round(n / sum(n), 3)) %>%
  kable_rapi(caption = "Proporsi spesies menurut kategori panjang petal")

Proporsi spesies menurut kategori panjang petal

PetalLength_cat	Species	n	Proporsi
Pendek	setosa	50	1.000
Sedang	versicolor	48	0.889
Sedang	virginica	6	0.111
Panjang	versicolor	2	0.043
Panjang	virginica	44	0.957

iris_cat %>%
  count(PetalWidth_cat, Species) %>%
  group_by(PetalWidth_cat) %>%
  mutate(Proporsi = round(n / sum(n), 3)) %>%
  kable_rapi(caption = "Proporsi spesies menurut kategori lebar petal")

Proporsi spesies menurut kategori lebar petal

PetalWidth_cat	Species	n	Proporsi
Kecil	setosa	50	1.000
Sedang	versicolor	49	0.907
Sedang	virginica	5	0.093
Besar	versicolor	1	0.022
Besar	virginica	45	0.978

2.5 Eksplorasi Visual

ggplot(iris_dat, aes(x = Species, fill = Species)) +
  geom_bar(alpha = 0.85) +
  scale_fill_manual(values = pal[1:3]) +
  labs(title = "Distribusi Spesies Iris",
       x = "Spesies", y = "Frekuensi", fill = "Spesies") +
  theme_analisis()

Grafik menunjukkan bahwa ketiga spesies bunga iris, yaitu setosa, versicolor, dan virginica, memiliki jumlah observasi yang sama, masing-masing sebanyak 50 data. Distribusi kelas yang seimbang ini menguntungkan dalam analisis regresi logistik multinomial karena model tidak akan cenderung lebih berpihak pada salah satu kategori spesies tertentu.

ggplot(iris_dat, aes(x = Petal.Length, y = Petal.Width, color = Species)) +
  geom_point(size = 2.5, alpha = 0.85) +
  scale_color_manual(values = pal[1:3]) +
  labs(title = "Sebaran Petal Length dan Petal Width",
       x = "Petal Length", y = "Petal Width", color = "Spesies") +
  theme_analisis()

Grafik menunjukkan adanya hubungan positif antara Petal Length dan Petal Width, di mana semakin panjang petal maka semakin lebar petalnya. Selain itu, terlihat bahwa ketiga spesies iris membentuk kelompok yang relatif terpisah, sehingga kedua variabel ini berpotensi menjadi prediktor penting dalam membedakan spesies bunga iris.

ggplot(iris_dat, aes(x = Species, y = Petal.Length, fill = Species)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = pal[1:3]) +
  labs(title = "Perbandingan Petal Length berdasarkan Species",
       x = "Spesies", y = "Petal Length", fill = "Spesies") +
  theme_analisis()

Boxplot menunjukkan bahwa Petal Length berbeda cukup jelas antar spesies. Spesies setosa memiliki panjang petal paling pendek, versicolor berada di tingkat sedang, dan virginica memiliki panjang petal paling besar. Perbedaan yang jelas ini menunjukkan bahwa Petal Length merupakan variabel yang penting untuk membedakan spesies bunga iris.

2.6 Pembagian Training dan Testing

set.seed(42)
train_id_iris <- stratified_split(iris_dat$Species, prop = 0.8, seed = 42)
train_iris <- iris_dat[train_id_iris, ]
test_iris  <- iris_dat[-train_id_iris, ]

bind_rows(
  train_iris %>% count(Species) %>% mutate(Data = "Training"),
  test_iris %>% count(Species) %>% mutate(Data = "Testing")
) %>%
  group_by(Data) %>%
  mutate(Proporsi = percent(n / sum(n), accuracy = 0.1)) %>%
  ungroup() %>%
  dplyr::select(Data, Species, Jumlah = n, Proporsi) %>%
  kable_rapi(caption = "Distribusi kelas pada training dan testing")

Distribusi kelas pada training dan testing

Data	Species	Jumlah	Proporsi
Training	setosa	40	33.3%
Training	versicolor	40	33.3%
Training	virginica	40	33.3%
Testing	setosa	10	33.3%
Testing	versicolor	10	33.3%
Testing	virginica	10	33.3%

2.7 Model Regresi Logistik Multinomial

2.7.1 Rumus Model Multinomial

Misalkan \(Y_i\) adalah spesies iris dengan kategori \(1,2,\ldots,J\). Salah satu kategori dijadikan referensi. Untuk kategori \(j\), model multinomial ditulis:

\[ \log \left(\frac{P(Y_i=j)}{P(Y_i=ref)}\right) = \beta_{0j} + \beta_{1j}X_{1i} + \cdots + \beta_{kj}X_{ki} \]

Ukuran efek pada regresi multinomial dinyatakan dengan Relative Risk Ratio (RRR):

\[ RRR = \exp(\hat{\beta}) \]

2.7.2 Model Penuh

fit_multi_full <- multinom(
  Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
  data = train_iris,
  trace = FALSE
)

null_multi <- multinom(Species ~ 1, data = train_iris, trace = FALSE)

ringkasan_multi <- data.frame(
  Keterangan = c("Jumlah observasi training", "LogLik model null",
                 "LogLik model penuh", "AIC"),
  Nilai = c(
    nrow(train_iris),
    round(as.numeric(logLik(null_multi)), 3),
    round(as.numeric(logLik(fit_multi_full)), 3),
    round(AIC(fit_multi_full), 3)
  )
)

ringkasan_multi %>%
  kable_rapi(caption = "Ringkasan kecocokan model penuh")

Ringkasan kecocokan model penuh

Keterangan	Nilai
Jumlah observasi training	120.000
LogLik model null	-131.833
LogLik model penuh	-0.195
AIC	20.391

Model regresi logistik multinomial diestimasi menggunakan 120 observasi data training. Nilai LogLik model penuh (-0,195) jauh lebih besar dibandingkan LogLik model null (-131,833), yang menunjukkan bahwa penambahan variabel prediktor memberikan peningkatan kecocokan model yang sangat baik dibandingkan model tanpa prediktor. Selain itu, nilai AIC sebesar 20,391 yang relatif kecil mengindikasikan bahwa model memiliki keseimbangan yang baik antara kecocokan dan kompleksitas model. Secara umum, hasil ini menunjukkan bahwa model mampu membedakan kategori spesies iris dengan sangat baik.

z_multi <- summary(fit_multi_full)$coefficients / summary(fit_multi_full)$standard.errors
p_multi <- 2 * (1 - pnorm(abs(z_multi)))

coef_multi <- broom::tidy(fit_multi_full) %>%
  mutate(
    RRR = exp(estimate),
    p_value = as.vector(t(p_multi))
  ) %>%
  transmute(
    Kelas = y.level,
    Variabel = term,
    `Relative Risk Ratio` = round(RRR, 3),
    `p-value` = signif(p_value, 3)
  )

coef_multi %>%
  kable_rapi(caption = "Relative Risk Ratio dan p-value — Model Penuh")

Relative Risk Ratio dan p-value — Model Penuh

Kelas	Variabel	Relative Risk Ratio	p-value
versicolor	(Intercept)	7.380668e+20	0.931
versicolor	Sepal.Length	1.300000e-02	0.983
versicolor	Sepal.Width	0.000000e+00	0.935
versicolor	Petal.Length	1.311609e+09	0.931
versicolor	Petal.Width	0.000000e+00	0.971
virginica	(Intercept)	0.000000e+00	0.934
virginica	Sepal.Length	0.000000e+00	0.908
virginica	Sepal.Width	0.000000e+00	0.825
virginica	Petal.Length	8.332876e+26	0.808
virginica	Petal.Width	2.554927e+24	0.924

Hasil model menunjukkan bahwa seluruh variabel memiliki p-value lebih besar dari 0,05, sehingga tidak terdapat bukti yang cukup bahwa masing-masing prediktor berpengaruh signifikan dalam membedakan spesies versicolor maupun virginica terhadap kategori referensi. Selain itu, nilai Relative Risk Ratio (RRR) yang sangat besar atau sangat kecil mengindikasikan adanya ketidakstabilan estimasi parameter, yang kemungkinan disebabkan oleh pemisahan data (separation) yang sangat baik pada dataset iris. Kondisi ini umum terjadi karena karakteristik spesies iris dapat dibedakan dengan sangat jelas berdasarkan ukuran petalnya.

2.8 Seleksi Variabel: Stepwise Backward (AIC)

cat("== Stepwise Backward Selection (AIC) ==\n")

## == Stepwise Backward Selection (AIC) ==

fit_multi_red <- stepAIC(fit_multi_full, direction = "backward", trace = TRUE)

## Start:  AIC=20.39
## Species ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width
## 
##                Df    AIC
## - Sepal.Length  2 17.156
## <none>            20.391
## - Sepal.Width   2 24.185
## - Petal.Width   2 28.405
## - Petal.Length  2 34.920
## 
## Step:  AIC=17.16
## Species ~ Sepal.Width + Petal.Length + Petal.Width
## 
##                Df    AIC
## <none>            17.156
## - Sepal.Width   2 23.234
## - Petal.Width   2 29.644
## - Petal.Length  2 33.106

cat("\n== Formula Model Reduksi (hasil stepAIC) ==\n")

## 
## == Formula Model Reduksi (hasil stepAIC) ==

print(formula(fit_multi_red))

## Species ~ Sepal.Width + Petal.Length + Petal.Width
## attr(,"variables")
## list(Species, Sepal.Width, Petal.Length, Petal.Width)
## attr(,"factors")
##              Sepal.Width Petal.Length Petal.Width
## Species                0            0           0
## Sepal.Width            1            0           0
## Petal.Length           0            1           0
## Petal.Width            0            0           1
## attr(,"term.labels")
## [1] "Sepal.Width"  "Petal.Length" "Petal.Width" 
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,"predvars")
## list(Species, Sepal.Width, Petal.Length, Petal.Width)
## attr(,"dataClasses")
##      Species  Sepal.Width Petal.Length  Petal.Width 
##     "factor"    "numeric"    "numeric"    "numeric"

data.frame(
  Keterangan = c("AIC", "LogLik"),
  `Model Penuh` = c(round(AIC(fit_multi_full), 3),
                    round(as.numeric(logLik(fit_multi_full)), 3)),
  `Model Reduksi` = c(round(AIC(fit_multi_red), 3),
                      round(as.numeric(logLik(fit_multi_red)), 3)),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Perbandingan ringkasan Model Penuh vs Reduksi")

Perbandingan ringkasan Model Penuh vs Reduksi

Keterangan	Model Penuh	Model Reduksi
AIC	20.391	17.156
LogLik	-0.195	-0.578

Hasil seleksi variabel menunjukkan bahwa model reduksi memiliki nilai AIC lebih kecil (17,156) dibandingkan model penuh (20,391), sehingga model reduksi lebih baik dalam menyeimbangkan kecocokan model dan kompleksitas parameter. Meskipun nilai LogLik model reduksi sedikit lebih kecil dibandingkan model penuh, perbedaannya relatif kecil. Oleh karena itu, model reduksi lebih disarankan karena lebih sederhana namun tetap mampu menjelaskan data dengan baik.

2.9 Uji Asumsi

2.9.1 Multikolinearitas — VIF

vif_multi_lm <- lm(as.numeric(Species) ~ Sepal.Length + Sepal.Width + Petal.Length + Petal.Width,
                   data = train_iris)
vif(vif_multi_lm)

## Sepal.Length  Sepal.Width Petal.Length  Petal.Width 
##     7.161789     2.052739    28.931854    14.639661

Nilai VIF untuk variabel Sepal.Length (7,16), Petal.Length (28,93), dan Petal.Width (14,64) melebihi batas umum 5, bahkan beberapa melebihi 10. Hal ini menunjukkan adanya multikolinearitas yang kuat antar variabel prediktor, terutama pada variabel Petal.Length dan Petal.Width. Sebaliknya, variabel Sepal.Width memiliki VIF sebesar 2,05 sehingga tidak menunjukkan masalah multikolinearitas. Dengan demikian, terdapat indikasi bahwa beberapa prediktor saling berkorelasi tinggi sehingga dapat menyebabkan ketidakstabilan estimasi parameter pada model regresi logistik multinomial.

2.9.2 Likelihood Ratio Test (LRT)

cat("== LRT: Model Penuh vs Null ==\n")

## == LRT: Model Penuh vs Null ==

lrtest(null_multi, fit_multi_full)

cat("\n== LRT: Model Reduksi vs Null ==\n")

## 
## == LRT: Model Reduksi vs Null ==

lrtest(null_multi, fit_multi_red)

cat("\n== LRT: Model Reduksi vs Model Penuh ==\n")

## 
## == LRT: Model Reduksi vs Model Penuh ==

lrtest(fit_multi_red, fit_multi_full)

Hasil Likelihood Ratio Test (LRT) menunjukkan bahwa baik model penuh maupun model reduksi secara signifikan lebih baik dibandingkan model null (p-value < 0,001). Selain itu, perbandingan antara model reduksi dan model penuh menghasilkan p-value = 0,682 (> 0,05), yang menunjukkan bahwa penghapusan variabel Sepal.Length tidak menurunkan performa model secara signifikan. Oleh karena itu, model reduksi dapat dipilih karena lebih sederhana namun tetap memiliki kemampuan yang sangat baik dalam membedakan spesies bunga iris.

2.9.3 Pseudo \(R^2\)

data.frame(
  Model = c("Model Penuh", "Model Reduksi"),
  `McFadden R²` = c(
    mcfadden_r2(fit_multi_full, null_multi),
    mcfadden_r2(fit_multi_red, null_multi)
  ),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Pseudo R² model multinomial")

Pseudo R² model multinomial

Model	McFadden R²
Model Penuh	0.999
Model Reduksi	0.996

Nilai McFadden R² pada model penuh sebesar 0,999 dan pada model reduksi sebesar 0,996. Nilai yang sangat mendekati 1 menunjukkan bahwa kedua model memiliki kemampuan yang sangat baik dalam menjelaskan variasi kategori spesies bunga iris. Perbedaan nilai yang sangat kecil mengindikasikan bahwa model reduksi tetap mampu mempertahankan performa yang hampir sama dengan model penuh meskipun menggunakan lebih sedikit variabel.

2.10 Relative Risk Ratio Model Final

z_multi_red <- summary(fit_multi_red)$coefficients / summary(fit_multi_red)$standard.errors
p_multi_red <- 2 * (1 - pnorm(abs(z_multi_red)))

broom::tidy(fit_multi_red) %>%
  mutate(
    RRR = exp(estimate),
    p_value = as.vector(t(p_multi_red))
  ) %>%
  transmute(
    Kelas = y.level,
    Variabel = term,
    `Relative Risk Ratio` = round(RRR, 3),
    `p-value` = signif(p_value, 3)
  ) %>%
  kable_rapi(caption = "Relative Risk Ratio — Model Final")

Relative Risk Ratio — Model Final

Kelas	Variabel	Relative Risk Ratio	p-value
versicolor	(Intercept)	1.512382e+26	0.278
versicolor	Sepal.Width	0.000000e+00	0.000
versicolor	Petal.Length	1.840219e+17	0.000
versicolor	Petal.Width	3.000000e-03	0.846
virginica	(Intercept)	0.000000e+00	0.067
virginica	Sepal.Width	0.000000e+00	0.000
virginica	Petal.Length	8.027058e+27	0.000
virginica	Petal.Width	3.079057e+31	0.014

Berdasarkan model final, variabel Sepal.Width dan Petal.Length berpengaruh signifikan dalam membedakan spesies bunga iris karena memiliki p-value < 0,05 pada kedua kategori. Selain itu, Petal.Width juga berpengaruh signifikan untuk membedakan spesies virginica terhadap kategori referensi (p-value = 0,014). Nilai Relative Risk Ratio (RRR) yang sangat besar atau sangat kecil menunjukkan bahwa variabel-variabel tersebut memiliki kemampuan yang sangat kuat dalam memisahkan kategori spesies. Hasil ini sejalan dengan karakteristik dataset iris yang memang memiliki perbedaan morfologi bunga yang sangat jelas antar spesies.

2.11 Prediksi Kelas

pred_multi <- predict(fit_multi_red, newdata = test_iris, type = "class")

cm_multi <- table(Aktual = test_iris$Species, Prediksi = pred_multi)
cm_multi

##             Prediksi
## Aktual       setosa versicolor virginica
##   setosa         10          0         0
##   versicolor      0          9         1
##   virginica       0          1         9

akurasi_multi <- akurasi(test_iris$Species, pred_multi)

data.frame(
  Ukuran = "Akurasi",
  Nilai = round(akurasi_multi, 4)
) %>%
  kable_rapi(caption = "Akurasi prediksi model logistik multinomial")

Akurasi prediksi model logistik multinomial

Ukuran	Nilai
Akurasi	0.933

Nilai akurasi sebesar 93,3% menunjukkan bahwa model mampu mengklasifikasikan spesies bunga iris dengan benar pada sekitar 93 dari 100 observasi data pengujian. Nilai ini mengindikasikan bahwa model regresi logistik multinomial memiliki kemampuan klasifikasi yang sangat baik dalam membedakan spesies setosa, versicolor, dan virginica.

2.12 Kesimpulan

Berdasarkan hasil analisis regresi logistik multinomial menggunakan dataset iris, diperoleh bahwa model mampu mengklasifikasikan spesies bunga iris dengan sangat baik. Hasil seleksi variabel menunjukkan bahwa model reduksi lebih disarankan karena memiliki nilai AIC yang lebih rendah dan performa yang hampir sama dengan model penuh. Variabel Sepal.Width, Petal.Length, dan Petal.Width merupakan prediktor utama yang berperan dalam membedakan spesies bunga iris. Nilai McFadden R² sebesar 0,996 menunjukkan bahwa model memiliki kemampuan yang sangat tinggi dalam menjelaskan variasi kategori spesies, sedangkan nilai akurasi sebesar 93,3% menunjukkan kemampuan klasifikasi yang sangat baik pada data pengujian. Dengan demikian, model regresi logistik multinomial yang diperoleh efektif digunakan untuk mengidentifikasi dan memprediksi spesies bunga iris berdasarkan karakteristik morfologinya..

3 Regresi Logistik Ordinal

3.1 Deskripsi Data

Dataset yang digunakan adalah housing dari package MASS. Dataset ini merupakan data survei kepuasan terhadap kondisi perumahan di Copenhagen. Variabel respons adalah Sat, yaitu tingkat kepuasan penghuni rumah dengan kategori Low, Medium, dan High. Karena kategori respons memiliki urutan alami dari rendah ke tinggi, metode yang tepat digunakan adalah regresi logistik ordinal.

Variabel dependen: Sat

• Low = kepuasan rendah
• Medium = kepuasan sedang
• High = kepuasan tinggi

Variabel prediktor:

data(housing, package = "MASS")

tibble(
  Variabel = names(housing),
  Keterangan = c(
    "Tingkat kepuasan penghuni",
    "Pengaruh/tingkat pengaruh manajemen",
    "Tipe perumahan",
    "Kontak dengan penghuni lain",
    "Frekuensi observasi pada kombinasi kategori"
  )
) %>%
  kable_rapi(caption = "Keterangan variabel pada dataset housing")

Keterangan variabel pada dataset housing

Variabel	Keterangan
Sat	Tingkat kepuasan penghuni
Infl	Pengaruh/tingkat pengaruh manajemen
Type	Tipe perumahan
Cont	Kontak dengan penghuni lain
Freq	Frekuensi observasi pada kombinasi kategori

3.2 Load Library dan Data

raw_housing <- housing
glimpse(raw_housing)

## Rows: 72
## Columns: 5
## $ Sat  <ord> Low, Medium, High, Low, Medium, High, Low, Medium, High, Low, Med…
## $ Infl <fct> Low, Low, Low, Medium, Medium, Medium, High, High, High, Low, Low…
## $ Type <fct> Tower, Tower, Tower, Tower, Tower, Tower, Tower, Tower, Tower, Ap…
## $ Cont <fct> Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, …
## $ Freq <int> 21, 21, 28, 34, 22, 36, 10, 11, 36, 61, 23, 17, 43, 35, 40, 26, 1…

Interpretasi output: Dataset housing tersedia dalam bentuk data berkelompok dengan kolom Freq sebagai frekuensi. Untuk analisis individual, data akan diperluas berdasarkan frekuensi tersebut.

3.3 Persiapan Data

3.3.1 Preprocessing dan Pemberian Label

housing_dat <- raw_housing %>%
  tidyr::uncount(weights = Freq) %>%
  mutate(
    Sat = ordered(Sat, levels = c("Low", "Medium", "High")),
    Infl = factor(Infl),
    Type = factor(Type),
    Cont = factor(Cont)
  )

glimpse(housing_dat)

## Rows: 1,681
## Columns: 4
## $ Sat  <ord> Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, …
## $ Infl <fct> Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, …
## $ Type <fct> Tower, Tower, Tower, Tower, Tower, Tower, Tower, Tower, Tower, To…
## $ Cont <fct> Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, Low, …

3.3.2 Pengecekan Data

cat("Missing values per kolom:\n")

## Missing values per kolom:

colSums(is.na(housing_dat))

##  Sat Infl Type Cont 
##    0    0    0    0

housing_dat %>%
  count(Sat) %>%
  mutate(Proporsi = round(n / sum(n), 4)) %>%
  kable_rapi(caption = "Distribusi variabel dependen")

Distribusi variabel dependen

Sat	n	Proporsi
Low	567	0.337
Medium	446	0.265
High	668	0.397

3.4 Statistik Deskriptif

3.4.1 Distribusi Variabel Prediktor

housing_dat %>%
  dplyr::select(Infl, Type, Cont) %>%
  summary()

##      Infl            Type       Cont    
##  Low   :627   Tower    :400   Low :713  
##  Medium:659   Apartment:765   High:968  
##  High  :395   Atrium   :239             
##               Terrace  :277

Berdasarkan statistik deskriptif dataset housing, variabel respon Sat (tingkat kepuasan penghuni) terdiri dari 627 kategori Low, 659 kategori Medium, dan 395 kategori High. Sementara itu, variabel Type menunjukkan bahwa jenis hunian Apartment merupakan yang paling banyak diamati (765 unit), diikuti Tower (400 unit), Terrace (277 unit), dan Atrium (239 unit). Untuk variabel Cont (kontak dengan penghuni lain), terdapat 713 observasi dengan tingkat kontak Low dan 968 observasi dengan tingkat kontak High. Secara umum, distribusi data menunjukkan variasi yang cukup baik pada setiap kategori sehingga layak digunakan untuk analisis regresi logistik ordinal.

3.4.2 Tabulasi Silang: Proporsi Kepuasan per Kategori

tab_ord <- function(var, label) {
  housing_dat %>%
    count(!!sym(var), Sat) %>%
    group_by(!!sym(var)) %>%
    mutate(Proporsi = round(n / sum(n), 3)) %>%
    kable_rapi(caption = paste("Proporsi tingkat kepuasan menurut", label))
}

tab_ord("Infl", "pengaruh manajemen")

Proporsi tingkat kepuasan menurut pengaruh manajemen

Infl	Sat	n	Proporsi
Low	Low	282	0.450
Low	Medium	170	0.271
Low	High	175	0.279
Medium	Low	206	0.313
Medium	Medium	189	0.287
Medium	High	264	0.401
High	Low	79	0.200
High	Medium	87	0.220
High	High	229	0.580

tab_ord("Type", "tipe perumahan")

Proporsi tingkat kepuasan menurut tipe perumahan

Type	Sat	n	Proporsi
Tower	Low	99	0.248
Tower	Medium	101	0.252
Tower	High	200	0.500
Apartment	Low	271	0.354
Apartment	Medium	192	0.251
Apartment	High	302	0.395
Atrium	Low	64	0.268
Atrium	Medium	79	0.331
Atrium	High	96	0.402
Terrace	Low	133	0.480
Terrace	Medium	74	0.267
Terrace	High	70	0.253

tab_ord("Cont", "kontak sosial")

Proporsi tingkat kepuasan menurut kontak sosial

Cont	Sat	n	Proporsi
Low	Low	262	0.367
Low	Medium	178	0.250
Low	High	273	0.383
High	Low	305	0.315
High	Medium	268	0.277
High	High	395	0.408

3.5 Eksplorasi Visual

ggplot(housing_dat, aes(x = Sat, fill = Sat)) +
  geom_bar(alpha = 0.85) +
  scale_fill_manual(values = pal[1:3]) +
  labs(title = "Distribusi Tingkat Kepuasan",
       x = "Tingkat Kepuasan", y = "Frekuensi", fill = "Kepuasan") +
  theme_analisis()

Grafik menunjukkan bahwa kategori kepuasan tinggi (High) memiliki frekuensi terbesar dibandingkan kategori Low dan Medium. Hal ini mengindikasikan bahwa sebagian besar penghuni merasa puas terhadap kondisi tempat tinggalnya. Selain itu, ketiga kategori kepuasan memiliki jumlah observasi yang cukup memadai sehingga analisis regresi logistik ordinal dapat dilakukan dengan baik.

ggplot(housing_dat, aes(x = Infl, fill = Sat)) +
  geom_bar(position = "fill", alpha = 0.85) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = pal[1:3]) +
  labs(title = "Proporsi Kepuasan berdasarkan Infl",
       x = "Infl", y = "Proporsi", fill = "Kepuasan") +
  theme_analisis()

menunjukkan bahwa tingkat kepuasan penghuni cenderung meningkat seiring meningkatnya tingkat pengaruh (Infl). Pada kategori Infl tinggi, proporsi kepuasan High terlihat lebih besar dibandingkan kategori lainnya, sedangkan pada Infl rendah proporsi kepuasan Low masih cukup dominan. Pola ini mengindikasikan bahwa variabel Infl berpotensi berhubungan dengan tingkat kepuasan penghuni dan dapat menjadi prediktor penting dalam model regresi logistik ordinal.

ggplot(housing_dat, aes(x = Type, fill = Sat)) +
  geom_bar(position = "fill", alpha = 0.85) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = pal[1:3]) +
  labs(title = "Proporsi Kepuasan berdasarkan Tipe Perumahan",
       x = "Tipe", y = "Proporsi", fill = "Kepuasan") +
  theme_analisis()

Grafik menunjukkan bahwa tingkat kepuasan penghuni berbeda pada setiap jenis hunian (Type). Hunian Tower memiliki proporsi kepuasan tinggi yang relatif paling besar, sedangkan Terrace memiliki proporsi kepuasan rendah yang lebih dominan dibandingkan jenis hunian lainnya. Perbedaan pola ini mengindikasikan bahwa jenis hunian berpotensi memengaruhi tingkat kepuasan penghuni dan layak dipertimbangkan dalam model regresi logistik ordinal.

3.6 Pembagian Training dan Testing

set.seed(42)
train_id_housing <- stratified_split(housing_dat$Sat, prop = 0.8, seed = 42)
train_housing <- housing_dat[train_id_housing, ]
test_housing  <- housing_dat[-train_id_housing, ]

bind_rows(
  train_housing %>% count(Sat) %>% mutate(Data = "Training"),
  test_housing %>% count(Sat) %>% mutate(Data = "Testing")
) %>%
  group_by(Data) %>%
  mutate(Proporsi = percent(n / sum(n), accuracy = 0.1)) %>%
  ungroup() %>%
  dplyr::select(Data, Sat, Jumlah = n, Proporsi) %>%
  kable_rapi(caption = "Distribusi kelas pada training dan testing")

Distribusi kelas pada training dan testing

Data	Sat	Jumlah	Proporsi
Training	Low	453	33.7%
Training	Medium	356	26.5%
Training	High	534	39.8%
Testing	Low	114	33.7%
Testing	Medium	90	26.6%
Testing	High	134	39.6%

3.7 Model Regresi Logistik Ordinal

3.7.1 Rumus Model Proportional Odds

Model proportional odds untuk respons ordinal dituliskan sebagai:

\[ \log \left(\frac{P(Y_i \leq j)}{P(Y_i > j)}\right) = \alpha_j - \beta_1X_{1i} - \cdots - \beta_kX_{ki} \]

dengan \(j = 1,\ldots,J-1\). Pada MASS::polr, tanda koefisien mengikuti bentuk \(\alpha_j - X\beta\). Oleh karena itu, interpretasi odds ratio sering memperhatikan arah model dan konteks kategori respons.

3.7.2 Model Penuh

fit_ord_full <- polr(
  Sat ~ Infl + Type + Cont,
  data = train_housing,
  Hess = TRUE,
  method = "logistic"
)

null_ord <- polr(Sat ~ 1, data = train_housing, Hess = TRUE, method = "logistic")

ringkasan_ord <- data.frame(
  Keterangan = c("Jumlah observasi training", "LogLik model null",
                 "LogLik model penuh", "AIC"),
  Nilai = c(
    nrow(train_housing),
    round(as.numeric(logLik(null_ord)), 3),
    round(as.numeric(logLik(fit_ord_full)), 3),
    round(AIC(fit_ord_full), 3)
  )
)

ringkasan_ord %>%
  kable_rapi(caption = "Ringkasan kecocokan model penuh")

Ringkasan kecocokan model penuh

Keterangan	Nilai
Jumlah observasi training	1343.000
LogLik model null	-1457.468
LogLik model penuh	-1385.638
AIC	2787.275

Model regresi logistik ordinal diestimasi menggunakan 1.343 observasi data training. Nilai LogLik model penuh (-1385,638) lebih besar dibandingkan LogLik model null (-1457,468), yang menunjukkan bahwa variabel prediktor mampu meningkatkan kecocokan model dibandingkan model tanpa prediktor. Selain itu, nilai AIC sebesar 2787,275 dapat digunakan sebagai acuan dalam pemilihan model, di mana nilai AIC yang lebih kecil menunjukkan model yang lebih baik. Secara umum, hasil ini menunjukkan bahwa model memiliki kemampuan yang cukup baik dalam menjelaskan variasi tingkat kepuasan penghuni..

ct_ord <- coef(summary(fit_ord_full))
p_ord <- pnorm(abs(ct_ord[, "t value"]), lower.tail = FALSE) * 2

data.frame(ct_ord) %>%
  rownames_to_column("Parameter") %>%
  mutate(`p-value` = signif(p_ord, 3)) %>%
  kable_rapi(caption = "Koefisien dan p-value — Model Penuh")

Koefisien dan p-value — Model Penuh

Parameter	Value	Std..Error	t.value	p-value
InflMedium	0.524	0.117	4.484	0.000
InflHigh	1.285	0.143	8.989	0.000
TypeApartment	-0.640	0.133	-4.812	0.000
TypeAtrium	-0.434	0.173	-2.504	0.012
TypeTerrace	-1.181	0.169	-6.982	0.000
ContHigh	0.420	0.107	3.936	0.000
Low|Medium	-0.542	0.138	-3.925	0.000
Medium|High	0.649	0.139	4.684	0.000

Hasil model menunjukkan bahwa variabel Infl, Type, dan Cont berpengaruh signifikan terhadap tingkat kepuasan penghuni karena memiliki p-value < 0,05. Koefisien positif pada InflMedium, InflHigh, dan ContHigh menunjukkan bahwa peningkatan tingkat pengaruh dan kontak dengan penghuni lain cenderung meningkatkan peluang berada pada kategori kepuasan yang lebih tinggi. Sebaliknya, koefisien negatif pada TypeApartment, TypeAtrium, dan TypeTerrace menunjukkan bahwa jenis hunian tersebut memiliki kecenderungan tingkat kepuasan yang lebih rendah dibandingkan kategori referensi (Tower). Secara keseluruhan, faktor pengaruh lingkungan, jenis hunian, dan tingkat kontak sosial memiliki peran penting dalam menentukan tingkat kepuasan penghuni.

3.8 Seleksi Variabel: Stepwise Backward (AIC)

cat("== Stepwise Backward Selection (AIC) ==\n")

## == Stepwise Backward Selection (AIC) ==

fit_ord_red <- stepAIC(fit_ord_full, direction = "backward", trace = TRUE)

## Start:  AIC=2787.28
## Sat ~ Infl + Type + Cont
## 
##        Df    AIC
## <none>    2787.3
## - Cont  1 2800.9
## - Type  3 2834.1
## - Infl  2 2868.1

cat("\n== Formula Model Reduksi (hasil stepAIC) ==\n")

## 
## == Formula Model Reduksi (hasil stepAIC) ==

print(formula(fit_ord_red))

## Sat ~ Infl + Type + Cont
## attr(,"variables")
## list(Sat, Infl, Type, Cont)
## attr(,"factors")
##      Infl Type Cont
## Sat     0    0    0
## Infl    1    0    0
## Type    0    1    0
## Cont    0    0    1
## attr(,"term.labels")
## [1] "Infl" "Type" "Cont"
## attr(,"order")
## [1] 1 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Sat, Infl, Type, Cont)
## attr(,"dataClasses")
##       Sat      Infl      Type      Cont 
## "ordered"  "factor"  "factor"  "factor"

data.frame(
  Keterangan = c("AIC", "LogLik"),
  `Model Penuh` = c(round(AIC(fit_ord_full), 3),
                    round(as.numeric(logLik(fit_ord_full)), 3)),
  `Model Reduksi` = c(round(AIC(fit_ord_red), 3),
                      round(as.numeric(logLik(fit_ord_red)), 3)),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Perbandingan ringkasan Model Penuh vs Reduksi")

Perbandingan ringkasan Model Penuh vs Reduksi

Keterangan	Model Penuh	Model Reduksi
AIC	2787.275	2787.275
LogLik	-1385.638	-1385.638

Hasil perbandingan menunjukkan bahwa model penuh dan model reduksi memiliki nilai AIC dan LogLik yang sama, yaitu AIC sebesar 2787,275 dan LogLik sebesar -1385,638. Hal ini menunjukkan bahwa pengurangan variabel tidak menyebabkan perubahan pada kecocokan model. Oleh karena itu, model reduksi lebih disarankan karena lebih sederhana namun tetap memiliki performa yang sama dengan model penuh.

3.9 Uji Asumsi

3.9.1 Multikolinearitas — VIF

vif_ord_lm <- lm(as.numeric(Sat) ~ Infl + Type + Cont, data = train_housing)
vif(vif_ord_lm)

##          GVIF Df GVIF^(1/(2*Df))
## Infl 1.016344  2        1.004061
## Type 1.031644  3        1.005206
## Cont 1.033510  1        1.016617

Seluruh nilai VIF mendekati 1 dan jauh di bawah 5, sehingga tidak terdapat masalah multikolinearitas pada model. Dengan demikian, variabel Infl, Type, dan Cont dapat digunakan bersama dalam model regresi logistik ordinal.

3.9.2 Likelihood Ratio Test (LRT)

cat("== LRT: Model Penuh vs Null ==\n")

## == LRT: Model Penuh vs Null ==

lrtest(null_ord, fit_ord_full)

cat("\n== LRT: Model Reduksi vs Null ==\n")

## 
## == LRT: Model Reduksi vs Null ==

lrtest(null_ord, fit_ord_red)

cat("\n== LRT: Model Reduksi vs Model Penuh ==\n")

## 
## == LRT: Model Reduksi vs Model Penuh ==

lrtest(fit_ord_red, fit_ord_full)

Hasil Likelihood Ratio Test (LRT) menunjukkan bahwa model dengan prediktor Infl, Type, dan Cont secara signifikan lebih baik dibandingkan model null (p-value < 0,001). Sementara itu, model reduksi dan model penuh memiliki nilai log-likelihood yang sama sehingga tidak terdapat perbedaan performa di antara keduanya. Dengan demikian, model yang digunakan sudah memadai untuk menjelaskan tingkat kepuasan penghuni.

3.9.3 Uji Parallel Lines

fit_clm <- clm(Sat ~ Infl + Type + Cont, data = train_housing, link = "logit")

nominal_test(fit_clm)

Output: Hasil uji parallel lines (nominal effects) menunjukkan bahwa seluruh variabel memiliki p-value lebih besar dari 0,05, yaitu Infl (0,8085), Type (0,2558), dan Cont (0,4228). Hal ini menunjukkan bahwa asumsi proportional odds (parallel lines) tidak dilanggar, sehingga model regresi logistik ordinal yang digunakan sudah sesuai untuk memodelkan tingkat kepuasan penghuni.

3.9.4 Pseudo \(R^2\)

data.frame(
  Model = c("Model Penuh", "Model Reduksi"),
  `McFadden R²` = c(
    mcfadden_r2(fit_ord_full, null_ord),
    mcfadden_r2(fit_ord_red, null_ord)
  ),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Pseudo R² model ordinal")

Pseudo R² model ordinal

Model	McFadden R²
Model Penuh	0.049
Model Reduksi	0.049

Nilai McFadden R² sebesar 0,049 pada model penuh maupun model reduksi menunjukkan bahwa model mampu menjelaskan sekitar 4,9% variasi tingkat kepuasan penghuni. Nilai yang relatif kecil ini mengindikasikan bahwa masih terdapat faktor lain di luar model yang memengaruhi tingkat kepuasan penghuni. Namun, kedua model memiliki kemampuan penjelasan yang sama karena menghasilkan nilai McFadden R² yang identik.

3.10 Odds Ratio Model Final

coef_ord <- coef(summary(fit_ord_red))
p_ord_red <- pnorm(abs(coef_ord[, "t value"]), lower.tail = FALSE) * 2

data.frame(coef_ord) %>%
  rownames_to_column("Parameter") %>%
  filter(!str_detect(Parameter, "\\|")) %>%
  mutate(
    `Odds Ratio` = round(exp(Value), 3),
    `p-value` = signif(p_ord_red[Parameter], 3)
  ) %>%
  dplyr::select(Parameter, `Odds Ratio`, `p-value`) %>%
  kable_rapi(caption = "Odds Ratio — Model Final")

Odds Ratio — Model Final

Parameter	Odds Ratio	p-value
InflMedium	1.688	0.000
InflHigh	3.615	0.000
TypeApartment	0.527	0.000
TypeAtrium	0.648	0.012
TypeTerrace	0.307	0.000
ContHigh	1.522	0.000

Berdasarkan model final, seluruh variabel memiliki p-value < 0,05 sehingga berpengaruh signifikan terhadap tingkat kepuasan penghuni. Variabel InflMedium dan InflHigh memiliki odds ratio lebih besar dari 1, yang menunjukkan bahwa semakin tinggi tingkat pengaruh (Infl), semakin besar peluang penghuni berada pada kategori kepuasan yang lebih tinggi. Variabel ContHigh juga meningkatkan peluang kepuasan sebesar 1,522 kali dibandingkan kategori referensi. Sebaliknya, variabel TypeApartment, TypeAtrium, dan TypeTerrace memiliki odds ratio kurang dari 1, yang menunjukkan bahwa jenis hunian tersebut cenderung memiliki peluang kepuasan yang lebih rendah dibandingkan jenis hunian referensi (Tower).

3.11 Prediksi Kategori

pred_ord <- predict(fit_ord_red, newdata = test_housing, type = "class")
pred_ord <- factor(as.character(pred_ord), levels = levels(test_housing$Sat), ordered = TRUE)

cm_ord <- table(Aktual = test_housing$Sat, Prediksi = pred_ord)
cm_ord

##         Prediksi
## Aktual   Low Medium High
##   Low     72      0   42
##   Medium  43      0   47
##   High    43      0   91

akurasi_ord <- akurasi(test_housing$Sat, pred_ord)

data.frame(
  Ukuran = "Akurasi",
  Nilai = round(akurasi_ord, 4)
) %>%
  kable_rapi(caption = "Akurasi prediksi model logistik ordinal")

Akurasi prediksi model logistik ordinal

Ukuran	Nilai
Akurasi	0.482

Nilai akurasi sebesar 48,2% menunjukkan bahwa model mampu mengklasifikasikan tingkat kepuasan penghuni dengan benar pada sekitar 48 dari 100 observasi data pengujian. Nilai ini tergolong rendah, sehingga kemampuan prediksi model dalam menentukan kategori kepuasan masih terbatas meskipun beberapa variabel prediktor terbukti signifikan secara statistik.

3.12 Kesimpulan

Berdasarkan hasil analisis regresi logistik ordinal menggunakan dataset housing, diperoleh bahwa variabel Infl, Type, dan Cont berpengaruh signifikan terhadap tingkat kepuasan penghuni. Hasil pengujian menunjukkan bahwa model memenuhi asumsi yang diperlukan, termasuk tidak adanya masalah multikolinearitas dan terpenuhinya asumsi proportional odds. Variabel dengan tingkat pengaruh dan kontak sosial yang lebih tinggi cenderung meningkatkan peluang kepuasan, sedangkan beberapa jenis hunian memiliki peluang kepuasan yang lebih rendah dibandingkan kategori referensi. Meskipun demikian, nilai McFadden R² sebesar 0,049 dan akurasi sebesar 48,2% menunjukkan bahwa kemampuan model dalam menjelaskan dan memprediksi tingkat kepuasan masih terbatas. Oleh karena itu, faktor-faktor lain di luar model kemungkinan juga berperan dalam menentukan tingkat kepuasan penghuni..

4 Regresi Poisson

4.1 Deskripsi Data

Dataset yang digunakan adalah quine dari package MASS. Dataset ini merupakan data nyata mengenai jumlah hari tidak masuk sekolah pada anak-anak di Australia. Variabel respons adalah Days, yaitu jumlah hari tidak masuk sekolah. Karena variabel respons berbentuk data cacah nonnegatif, regresi Poisson digunakan untuk memodelkan ekspektasi jumlah hari tidak masuk berdasarkan prediktor yang tersedia.

Variabel dependen: Days

Variabel prediktor:

data(quine, package = "MASS")

tibble(
  Variabel = names(quine),
  Keterangan = c(
    "Etnis siswa",
    "Jenis kelamin",
    "Kelompok usia",
    "Status belajar",
    "Jumlah hari tidak masuk sekolah"
  )
) %>%
  kable_rapi(caption = "Keterangan variabel pada dataset quine")

Keterangan variabel pada dataset quine

Variabel	Keterangan
Eth	Etnis siswa
Sex	Jenis kelamin
Age	Kelompok usia
Lrn	Status belajar
Days	Jumlah hari tidak masuk sekolah

4.2 Load Library dan Data

raw_quine <- quine
glimpse(raw_quine)

## Rows: 146
## Columns: 5
## $ Eth  <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A,…
## $ Sex  <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M,…
## $ Age  <fct> F0, F0, F0, F0, F0, F0, F0, F0, F1, F1, F1, F1, F1, F2, F2, F2, F…
## $ Lrn  <fct> SL, SL, SL, AL, AL, AL, AL, AL, SL, SL, SL, AL, AL, SL, SL, SL, S…
## $ Days <int> 2, 11, 14, 5, 5, 13, 20, 22, 6, 6, 15, 7, 14, 6, 32, 53, 57, 14, …

Dataset quine terdiri dari 146 observasi dan 5 variabel. Variabel Days merupakan variabel count sehingga sesuai dianalisis menggunakan regresi Poisson.

4.3 Persiapan Data

4.3.1 Preprocessing dan Pemberian Label

quine_dat <- raw_quine %>%
  mutate(
    Eth = factor(Eth),
    Sex = factor(Sex),
    Age = factor(Age),
    Lrn = factor(Lrn)
  )

glimpse(quine_dat)

## Rows: 146
## Columns: 5
## $ Eth  <fct> A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A, A,…
## $ Sex  <fct> M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M, M,…
## $ Age  <fct> F0, F0, F0, F0, F0, F0, F0, F0, F1, F1, F1, F1, F1, F2, F2, F2, F…
## $ Lrn  <fct> SL, SL, SL, AL, AL, AL, AL, AL, SL, SL, SL, AL, AL, SL, SL, SL, S…
## $ Days <int> 2, 11, 14, 5, 5, 13, 20, 22, 6, 6, 15, 7, 14, 6, 32, 53, 57, 14, …

4.3.2 Pengecekan Data

cat("Missing values per kolom:\n")

## Missing values per kolom:

colSums(is.na(quine_dat))

##  Eth  Sex  Age  Lrn Days 
##    0    0    0    0    0

summary(quine_dat$Days)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##    0.00    5.00   11.00   16.46   22.75   81.00

data.frame(
  Rata_rata = mean(quine_dat$Days),
  Varians = var(quine_dat$Days),
  Minimum = min(quine_dat$Days),
  Maksimum = max(quine_dat$Days)
) %>%
  kable_rapi(caption = "Ringkasan variabel respons Days")

Ringkasan variabel respons Days

Rata_rata	Varians	Minimum	Maksimum
16.459	264.167	0	81

4.4 Statistik Deskriptif

4.4.1 Distribusi Variabel Prediktor

quine_dat %>%
  dplyr::select(Eth, Sex, Age, Lrn) %>%
  summary()

##  Eth    Sex    Age     Lrn    
##  A:69   F:80   F0:27   AL:83  
##  N:77   M:66   F1:46   SL:63  
##                F2:40          
##                F3:33

Dataset quine terdiri dari 146 siswa dengan distribusi yang relatif seimbang berdasarkan etnis (A = 69, N = 77) dan jenis kelamin (F = 80, M = 66). Variabel usia terbagi ke dalam empat kelompok (F0–F3), sedangkan tingkat pembelajaran terdiri dari kategori AL (83 siswa) dan SL (63 siswa). Distribusi data yang cukup merata ini mendukung analisis regresi Poisson untuk memodelkan jumlah hari ketidakhadiran siswa.

4.4.2 Hubungan Respon dan Prediktor

mean_days <- function(var, label) {
  quine_dat %>%
    group_by(!!sym(var)) %>%
    summarise(
      n = n(),
      Rata_rata_Days = round(mean(Days), 3),
      Median_Days = median(Days),
      .groups = "drop"
    ) %>%
    kable_rapi(caption = paste("Rata-rata Days menurut", label))
}

mean_days("Eth", "etnis")

Rata-rata Days menurut etnis

Eth	n	Rata_rata_Days	Median_Days
A	69	21.232	15
N	77	12.182	7

mean_days("Sex", "jenis kelamin")

Rata-rata Days menurut jenis kelamin

Sex	n	Rata_rata_Days	Median_Days
F	80	15.225	10
M	66	17.955	14

mean_days("Age", "kelompok usia")

Rata-rata Days menurut kelompok usia

Age	n	Rata_rata_Days	Median_Days
F0	27	14.852	11
F1	46	11.152	6
F2	40	21.050	14
F3	33	19.606	18

mean_days("Lrn", "status belajar")

Rata-rata Days menurut status belajar

Lrn	n	Rata_rata_Days	Median_Days
AL	83	15.819	12
SL	63	17.302	10

Siswa dengan kategori Slow Learner (SL) memiliki rata-rata jumlah hari tidak masuk sekolah sebesar 17,302 hari, lebih tinggi dibandingkan siswa Average Learner (AL) yang memiliki rata-rata 15,819 hari. Hal ini menunjukkan bahwa siswa yang mengalami kesulitan belajar cenderung memiliki tingkat ketidakhadiran yang lebih tinggi dibandingkan siswa dengan kemampuan belajar rata-rata.

4.5 Eksplorasi Visual

ggplot(quine_dat, aes(x = Days)) +
  geom_histogram(bins = 25, fill = pal[1], alpha = 0.85) +
  labs(title = "Distribusi Jumlah Hari Tidak Masuk Sekolah",
       x = "Days", y = "Frekuensi") +
  theme_analisis()

Histogram menunjukkan bahwa sebagian besar siswa memiliki jumlah hari tidak masuk sekolah yang relatif rendah, dengan frekuensi tertinggi berada pada rentang 0–10 hari. Distribusi data terlihat miring ke kanan (right-skewed), di mana hanya sedikit siswa yang memiliki jumlah ketidakhadiran sangat tinggi hingga lebih dari 50 hari. Pola ini sesuai dengan karakteristik data cacah (count data) sehingga regresi Poisson menjadi metode yang sesuai untuk memodelkan jumlah hari ketidakhadiran siswa.

ggplot(quine_dat, aes(x = Age, y = Days, fill = Age)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = pal[1:4]) +
  labs(title = "Jumlah Hari Tidak Masuk berdasarkan Kelompok Usia",
       x = "Kelompok Usia", y = "Days", fill = "Usia") +
  theme_analisis()

Boxplot menunjukkan bahwa jumlah hari tidak masuk sekolah cenderung meningkat pada kelompok usia yang lebih tinggi. Kelompok F2 dan F3 memiliki median serta variasi jumlah ketidakhadiran yang lebih besar dibandingkan F0 dan F1. Selain itu, terdapat beberapa pencilan (outlier) pada setiap kelompok usia, yang menunjukkan adanya siswa dengan tingkat ketidakhadiran yang jauh lebih tinggi dibandingkan siswa lainnya. Pola ini mengindikasikan bahwa usia berpotensi memengaruhi jumlah hari tidak masuk sekolah.

ggplot(quine_dat, aes(x = Lrn, y = Days, fill = Lrn)) +
  geom_boxplot(alpha = 0.85) +
  scale_fill_manual(values = pal[3:4]) +
  labs(title = "Jumlah Hari Tidak Masuk berdasarkan Status Belajar",
       x = "Status Belajar", y = "Days", fill = "Status") +
  theme_analisis()

Boxplot menunjukkan bahwa siswa Slow Learner (SL) cenderung memiliki jumlah hari tidak masuk sekolah yang sedikit lebih tinggi dibandingkan siswa Average Learner (AL). Selain itu, kelompok SL memiliki lebih banyak pencilan (outlier) dengan tingkat ketidakhadiran yang sangat tinggi. Hal ini mengindikasikan bahwa status belajar berpotensi memengaruhi jumlah hari ketidakhadiran siswa.

4.6 Pembagian Training dan Testing

set.seed(42)
train_id_quine <- sample(seq_len(nrow(quine_dat)), size = floor(0.8 * nrow(quine_dat)))
train_quine <- quine_dat[train_id_quine, ]
test_quine  <- quine_dat[-train_id_quine, ]

data.frame(
  Data = c("Training", "Testing"),
  Jumlah = c(nrow(train_quine), nrow(test_quine)),
  Rata_rata_Days = c(mean(train_quine$Days), mean(test_quine$Days))
) %>%
  kable_rapi(caption = "Ringkasan pembagian data training dan testing")

Ringkasan pembagian data training dan testing

Data	Jumlah	Rata_rata_Days
Training	116	16.405
Testing	30	16.667

4.7 Model Regresi Poisson

4.7.1 Rumus Model Poisson

Misalkan \(Y_i\) adalah jumlah hari tidak masuk sekolah siswa ke-\(i\), dengan:

\[ Y_i \sim \text{Poisson}(\mu_i) \]

Model regresi Poisson menggunakan link log:

\[ \log(\mu_i) = \beta_0 + \beta_1X_{1i} + \cdots + \beta_kX_{ki} \]

Sehingga:

\[ \mu_i = \exp(\beta_0 + \beta_1X_{1i} + \cdots + \beta_kX_{ki}) \]

Ukuran efek pada regresi Poisson disebut Incidence Rate Ratio (IRR):

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

4.7.2 Model Penuh

fit_pois_full <- glm(
  Days ~ Eth + Sex + Age + Lrn,
  data = train_quine,
  family = poisson(link = "log")
)

null_pois <- glm(Days ~ 1, data = train_quine, family = poisson(link = "log"))

ringkasan_pois <- data.frame(
  Keterangan = c("Jumlah observasi training", "Null deviance",
                 "Residual deviance", "Derajat bebas residual", "AIC"),
  Nilai = c(
    nobs(fit_pois_full),
    round(fit_pois_full$null.deviance, 3),
    round(fit_pois_full$deviance, 3),
    fit_pois_full$df.residual,
    round(AIC(fit_pois_full), 3)
  )
)

ringkasan_pois %>%
  kable_rapi(caption = "Ringkasan kecocokan model penuh")

Ringkasan kecocokan model penuh

Keterangan	Nilai
Jumlah observasi training	116.000
Null deviance	1632.167
Residual deviance	1326.191
Derajat bebas residual	109.000
AIC	1802.464

Model regresi Poisson diestimasi menggunakan 116 observasi data training. Nilai residual deviance sebesar 1326,191 lebih kecil dibandingkan null deviance sebesar 1632,167, yang menunjukkan bahwa variabel prediktor mampu meningkatkan kecocokan model dibandingkan model tanpa prediktor. Selain itu, nilai AIC sebesar 1802,464 digunakan sebagai ukuran pemilihan model, di mana nilai yang lebih kecil menunjukkan model yang lebih baik. Secara umum, model mampu menjelaskan variasi jumlah hari tidak masuk sekolah lebih baik dibandingkan model null.

irr_pois_full <- broom::tidy(fit_pois_full) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    IRR = exp(estimate),
    ci_low = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error)
  ) %>%
  arrange(p.value) %>%
  transmute(
    Variabel = term,
    `Incidence Rate Ratio` = round(IRR, 3),
    `IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
    `p-value` = signif(p.value, 3)
  )

irr_pois_full %>%
  kable_rapi(caption = "IRR dan p-value — Model Penuh")

IRR dan p-value — Model Penuh

Variabel	Incidence Rate Ratio	IK 95%	p-value
EthN	0.667	0.609 - 0.731	0.000
LrnSL	1.558	1.388 - 1.748	0.000
AgeF3	1.721	1.492 - 1.985	0.000
AgeF1	0.628	0.538 - 0.733	0.000
SexM	1.190	1.084 - 1.307	0.000
AgeF2	1.201	1.051 - 1.372	0.007

Berdasarkan model Poisson, seluruh variabel yang ditampilkan berpengaruh signifikan terhadap jumlah hari tidak masuk sekolah (p-value < 0,05). Siswa dengan kategori SL (Slow Learner) memiliki tingkat ketidakhadiran 1,558 kali lebih tinggi dibandingkan kategori referensi. Selain itu, siswa pada kelompok usia F3 memiliki tingkat ketidakhadiran 1,721 kali lebih tinggi dibandingkan kelompok referensi, sedangkan siswa laki-laki (SexM) memiliki tingkat ketidakhadiran 1,190 kali lebih tinggi dibandingkan perempuan. Sebaliknya, siswa dengan etnis N memiliki tingkat ketidakhadiran yang lebih rendah (IRR = 0,667) dibandingkan etnis referensi. Hasil ini menunjukkan bahwa etnis, status belajar, usia, dan jenis kelamin berpengaruh terhadap jumlah hari ketidakhadiran siswa.

4.8 Seleksi Variabel: Stepwise Backward (AIC)

cat("== Stepwise Backward Selection (AIC) ==\n")

## == Stepwise Backward Selection (AIC) ==

fit_pois_red <- stepAIC(fit_pois_full, direction = "backward", trace = TRUE)

## Start:  AIC=1802.46
## Days ~ Eth + Sex + Age + Lrn
## 
##        Df Deviance    AIC
## <none>      1326.2 1802.5
## - Sex   1   1339.6 1813.9
## - Lrn   1   1384.5 1858.7
## - Eth   1   1402.5 1876.7
## - Age   3   1508.5 1978.7

cat("\n== Formula Model Reduksi (hasil stepAIC) ==\n")

## 
## == Formula Model Reduksi (hasil stepAIC) ==

print(formula(fit_pois_red))

## Days ~ Eth + Sex + Age + Lrn

data.frame(
  Keterangan = c("Null deviance", "Residual deviance",
                 "Derajat bebas residual", "AIC"),
  `Model Penuh` = c(
    round(fit_pois_full$null.deviance, 3),
    round(fit_pois_full$deviance, 3),
    fit_pois_full$df.residual,
    round(AIC(fit_pois_full), 3)
  ),
  `Model Reduksi` = c(
    round(fit_pois_red$null.deviance, 3),
    round(fit_pois_red$deviance, 3),
    fit_pois_red$df.residual,
    round(AIC(fit_pois_red), 3)
  ),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Perbandingan ringkasan Model Penuh vs Reduksi")

Perbandingan ringkasan Model Penuh vs Reduksi

Keterangan	Model Penuh	Model Reduksi
Null deviance	1632.167	1632.167
Residual deviance	1326.191	1326.191
Derajat bebas residual	109.000	109.000
AIC	1802.464	1802.464

Hasil perbandingan menunjukkan bahwa model penuh dan model reduksi memiliki nilai yang sama untuk seluruh ukuran evaluasi, yaitu null deviance, residual deviance, dan AIC. Hal ini menunjukkan bahwa pengurangan variabel tidak memengaruhi performa model. Oleh karena itu, model reduksi dapat dipilih karena lebih sederhana namun tetap memberikan kecocokan yang sama dengan model penuh.

4.9 Uji Asumsi

4.9.1 Multikolinearitas — VIF

cat("== VIF Model Penuh ==\n")

## == VIF Model Penuh ==

vif(fit_pois_full)

##         GVIF Df GVIF^(1/(2*Df))
## Eth 1.014027  1        1.006989
## Sex 1.069431  1        1.034133
## Age 1.695608  3        1.091996
## Lrn 1.614461  1        1.270615

cat("\n== VIF Model Reduksi ==\n")

## 
## == VIF Model Reduksi ==

vif(fit_pois_red)

##         GVIF Df GVIF^(1/(2*Df))
## Eth 1.014027  1        1.006989
## Sex 1.069431  1        1.034133
## Age 1.695608  3        1.091996
## Lrn 1.614461  1        1.270615

Hasil pengujian multikolinearitas menunjukkan bahwa seluruh nilai VIF pada model penuh maupun model reduksi berada di bawah 5, bahkan mendekati 1. Hal ini menunjukkan bahwa tidak terdapat masalah multikolinearitas antar variabel Eth, Sex, Age, dan Lrn, sehingga seluruh prediktor dapat digunakan bersama dalam model regresi Poisson.

4.9.2 Overdispersion

cat("== Dispersion Test Model Penuh ==\n")

## == Dispersion Test Model Penuh ==

dispersiontest(fit_pois_full)

## 
##  Overdispersion test
## 
## data:  fit_pois_full
## z = 5.2123, p-value = 9.323e-08
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   11.70069

cat("\n== Dispersion Test Model Reduksi ==\n")

## 
## == Dispersion Test Model Reduksi ==

dispersiontest(fit_pois_red)

## 
##  Overdispersion test
## 
## data:  fit_pois_red
## z = 5.2123, p-value = 9.323e-08
## alternative hypothesis: true dispersion is greater than 1
## sample estimates:
## dispersion 
##   11.70069

dispersion_ratio <- sum(residuals(fit_pois_full, type = "pearson")^2) / fit_pois_red$df.residual

data.frame(
  Keterangan = "Pearson dispersion ratio model penuh",
  Nilai = round(dispersion_ratio, 3)
) %>%
  kable_rapi(caption = "Rasio dispersi")

Rasio dispersi

Keterangan	Nilai
Pearson dispersion ratio model penuh	12.452

Nilai Pearson Dispersion Ratio sebesar 12,452 jauh lebih besar dari 1, yang menunjukkan adanya overdispersion pada data. Hal ini berarti varians jumlah hari tidak masuk sekolah lebih besar daripada nilai rata-ratanya, sehingga asumsi dasar regresi Poisson tidak sepenuhnya terpenuhi.

4.9.3 Likelihood Ratio Test (LRT)

cat("== LRT: Model Penuh vs Null ==\n")

## == LRT: Model Penuh vs Null ==

lrtest(null_pois, fit_pois_full)

cat("\n== LRT: Model Reduksi vs Null ==\n")

## 
## == LRT: Model Reduksi vs Null ==

lrtest(null_pois, fit_pois_red)

cat("\n== LRT: Model Reduksi vs Model Penuh ==\n")

## 
## == LRT: Model Reduksi vs Model Penuh ==

lrtest(fit_pois_red, fit_pois_full)

Hasil Likelihood Ratio Test (LRT) menunjukkan bahwa model dengan prediktor Eth, Sex, Age, dan Lrn secara signifikan lebih baik dibandingkan model null (p-value < 0,001). Sementara itu, perbandingan model reduksi dan model penuh menghasilkan p-value = 1, yang menunjukkan bahwa kedua model memiliki performa yang identik. Dengan demikian, model yang digunakan sudah memadai untuk menjelaskan variasi jumlah hari tidak masuk sekolah

4.9.4 Pseudo \(R^2\)

data.frame(
  Model = c("Model Penuh", "Model Reduksi"),
  `McFadden R²` = c(
    mcfadden_r2(fit_pois_full, null_pois),
    mcfadden_r2(fit_pois_red, null_pois)
  ),
  check.names = FALSE
) %>%
  kable_rapi(caption = "Pseudo R² model Poisson")

Pseudo R² model Poisson

Model	McFadden R²
Model Penuh	0.146
Model Reduksi	0.146

Nilai McFadden R² sebesar 0,146 menunjukkan bahwa model mampu menjelaskan sekitar 14,6% variasi jumlah hari tidak masuk sekolah. Nilai ini mengindikasikan bahwa model memiliki kemampuan penjelasan yang moderat, namun masih terdapat faktor lain di luar model yang turut memengaruhi tingkat ketidakhadiran siswa. Karena model penuh dan model reduksi memiliki nilai yang sama, keduanya memberikan kemampuan penjelasan yang identik

4.10 Incidence Rate Ratio (IRR)

irr_pois_final <- broom::tidy(fit_pois_red) %>%
  filter(term != "(Intercept)") %>%
  mutate(
    IRR = exp(estimate),
    ci_low = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error)
  ) %>%
  arrange(p.value) %>%
  transmute(
    Variabel = term,
    `Incidence Rate Ratio` = round(IRR, 3),
    `IK 95%` = paste0(round(ci_low, 3), " - ", round(ci_high, 3)),
    `p-value` = signif(p.value, 3)
  )

irr_pois_final %>%
  kable_rapi(caption = "Incidence Rate Ratio — Model Final")

Incidence Rate Ratio — Model Final

Variabel	Incidence Rate Ratio	IK 95%	p-value
EthN	0.667	0.609 - 0.731	0.000
LrnSL	1.558	1.388 - 1.748	0.000
AgeF3	1.721	1.492 - 1.985	0.000
AgeF1	0.628	0.538 - 0.733	0.000
SexM	1.190	1.084 - 1.307	0.000
AgeF2	1.201	1.051 - 1.372	0.007

Berdasarkan model final, seluruh variabel memiliki p-value < 0,05 sehingga berpengaruh signifikan terhadap jumlah hari tidak masuk sekolah. Siswa dengan kategori SL (Slow Learner) memiliki tingkat ketidakhadiran 1,558 kali lebih tinggi dibandingkan AL (Average Learner). Selain itu, siswa pada kelompok usia F3 memiliki tingkat ketidakhadiran 1,721 kali lebih tinggi dibandingkan kelompok referensi, sedangkan siswa laki-laki (SexM) memiliki tingkat ketidakhadiran 1,190 kali lebih tinggi dibandingkan perempuan. Sebaliknya, siswa dengan etnis N memiliki tingkat ketidakhadiran yang lebih rendah (IRR = 0,667) dibandingkan etnis referensi. Dengan demikian, etnis, status belajar, usia, dan jenis kelamin berpengaruh signifikan terhadap jumlah hari ketidakhadiran siswa.

4.11 Prediksi Jumlah Kejadian

pred_pois <- predict(fit_pois_red, newdata = test_quine, type = "response")

eval_pois <- data.frame(
  MAE = mean(abs(test_quine$Days - pred_pois)),
  RMSE = sqrt(mean((test_quine$Days - pred_pois)^2))
)

eval_pois %>%
  kable_rapi(caption = "Evaluasi prediksi model Poisson")

Evaluasi prediksi model Poisson

MAE	RMSE
11.609	16.623

data.frame(
  Aktual = test_quine$Days,
  Prediksi = round(pred_pois, 2)
) %>%
  head(10) %>%
  kable_rapi(caption = "Contoh hasil prediksi jumlah hari tidak masuk")

Contoh hasil prediksi jumlah hari tidak masuk

	Aktual	Prediksi
7	20	16.52
11	15	16.17
12	7	10.38
23	43	19.84
28	28	28.42
37	5	13.58
39	6	13.58
45	53	13.58
51	19	8.72
52	8	25.96

abel menunjukkan hasil prediksi jumlah hari tidak masuk sekolah pada beberapa observasi data pengujian. Secara umum, model mampu mengikuti pola ketidakhadiran siswa, namun masih terdapat selisih antara nilai aktual dan nilai prediksi pada beberapa observasi. Misalnya, pada observasi ke-28 nilai prediksi (28,42 hari) sangat dekat dengan nilai aktual (28 hari), sedangkan pada observasi lain seperti ke-23 dan ke-45 terdapat perbedaan yang cukup besar. Hal ini menunjukkan bahwa model Poisson mampu menangkap pola umum data, tetapi belum sepenuhnya akurat dalam memprediksi jumlah ketidakhadiran untuk setiap siswa secara indiv

4.12 Kesimpulan

Berdasarkan hasil analisis regresi Poisson menggunakan dataset quine, diperoleh bahwa variabel Eth, Sex, Age, dan Lrn berpengaruh signifikan terhadap jumlah hari tidak masuk sekolah. Hasil model menunjukkan bahwa siswa dengan kategori Slow Learner (SL), siswa laki-laki, serta siswa pada kelompok usia yang lebih tinggi cenderung memiliki tingkat ketidakhadiran yang lebih besar. Pengujian asumsi menunjukkan tidak terdapat masalah multikolinearitas, namun terdeteksi adanya overdispersion pada data. Nilai McFadden R² sebesar 14,6% menunjukkan bahwa model mampu menjelaskan sebagian variasi jumlah hari tidak masuk sekolah, meskipun masih terdapat faktor lain di luar model yang memengaruhinya. Secara keseluruhan, model regresi Poisson mampu mengidentifikasi faktor-faktor yang berhubungan dengan tingkat ketidakhadiran siswa dan memberikan gambaran umum yang cukup baik terhadap pola data yang diamati.