Tugas Analisis Regresi Logistik

library(readxl)
library(tidyverse)
library(knitr)
library(kableExtra)
library(scales)
library(nnet)
library(broom)
library(caret)
library(car)
library(VGAM)
library(MASS)

theme_profesional <- function() {
  theme_minimal(base_size = 13) +
    theme(
      plot.title = element_text(face = "bold", size = 17, color = "#7c2d12"),
      plot.subtitle = element_text(size = 11, color = "#475569"),
      axis.title = element_text(face = "bold"),
      legend.position = "bottom",
      panel.grid.minor = element_blank()
    )
}

options(knitr.kable.NA = "")

Tugas Analisis Regresi Logistik

Dokumen ini menggabungkan empat analisis regresi dalam satu file R Markdown yang terstruktur, yaitu regresi logistik biner, regresi logistik multinomial, regresi logistik ordinal, dan regresi Poisson. Setiap bagian disusun dengan urutan langkah analisis, sintaks R, output, serta interpretasi hasil.

1. Regresi Logistik Biner

1.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan keputusan pelanggan dalam membeli sepeda menggunakan regresi logistik biner. Model ini sesuai digunakan karena variabel respon hanya memiliki dua kategori, yaitu Yes dan No.

Variabel respon yang digunakan adalah:

\[ Y = Purchased\ Bike \]

dengan pengkodean:

\[ Y = \begin{cases} 1, & \text{jika pelanggan membeli sepeda atau Yes} \\ 0, & \text{jika pelanggan tidak membeli sepeda atau No} \end{cases} \]

Variabel prediktor yang digunakan adalah:

• Age
• Income
• Children
• Gender
• Marital Status
• Education

1.2 Membaca Data Excel

data <- read_excel("E:/KULIAH/Semester 4/ADK/Tugas regresi/data biner sepeda.xlsx")

glimpse(data)

## Rows: 1,000
## Columns: 13
## $ ID                 <dbl> 12496, 24107, 14177, 24381, 25597, 13507, 27974, 19…
## $ `Marital Status`   <chr> "Married", "Married", "Married", "Single", "Single"…
## $ Gender             <chr> "Female", "Male", "Male", "Male", "Male", "Female",…
## $ Income             <dbl> 40000, 30000, 80000, 70000, 30000, 10000, 160000, 4…
## $ Children           <dbl> 1, 3, 5, 0, 0, 2, 2, 1, 2, 2, 3, 0, 5, 2, 1, 2, 3, …
## $ Education          <chr> "Bachelors", "Partial College", "Partial College", …
## $ Occupation         <chr> "Skilled Manual", "Clerical", "Professional", "Prof…
## $ `Home Owner`       <chr> "Yes", "Yes", "No", "Yes", "No", "Yes", "Yes", "Yes…
## $ Cars               <dbl> 0, 1, 2, 1, 0, 0, 4, 0, 2, 1, 2, 4, 4, 1, 1, 1, 2, …
## $ `Commute Distance` <chr> "0-1 Miles", "0-1 Miles", "2-5 Miles", "5-10 Miles"…
## $ Region             <chr> "Europe", "Europe", "Europe", "Pacific", "Europe", …
## $ Age                <dbl> 42, 43, 60, 41, 36, 50, 33, 43, 58, 48, 54, 36, 55,…
## $ `Purchased Bike`   <chr> "No", "No", "No", "Yes", "Yes", "No", "Yes", "Yes",…

Data dibaca dari file Excel. Variabel Purchased Bike digunakan sebagai variabel respon, sedangkan variabel usia, pendapatan, jumlah anak, jenis kelamin, status pernikahan, dan pendidikan digunakan sebagai variabel prediktor.

1.3 Eksplorasi dan Persiapan Data

data_model <- data %>%
  dplyr::select(
    Purchased_Bike = `Purchased Bike`,
    Age,
    Income,
    Children,
    Gender,
    Marital_Status = `Marital Status`,
    Education
  ) %>%
  tidyr::drop_na() %>%
  dplyr::mutate(
    Purchased_Bike = factor(
      Purchased_Bike,
      levels = c("No", "Yes")
    ),
    Purchased_Bike_num = ifelse(Purchased_Bike == "Yes", 1, 0),
    Gender = factor(Gender),
    Marital_Status = factor(Marital_Status),
    Education = factor(Education)
  )

glimpse(data_model)

## Rows: 1,000
## Columns: 8
## $ Purchased_Bike     <fct> No, No, No, Yes, Yes, No, Yes, Yes, No, Yes, Yes, N…
## $ Age                <dbl> 42, 43, 60, 41, 36, 50, 33, 43, 58, 48, 54, 36, 55,…
## $ Income             <dbl> 40000, 30000, 80000, 70000, 30000, 10000, 160000, 4…
## $ Children           <dbl> 1, 3, 5, 0, 0, 2, 2, 1, 2, 2, 3, 0, 5, 2, 1, 2, 3, …
## $ Gender             <fct> Female, Male, Male, Male, Male, Female, Male, Male,…
## $ Marital_Status     <fct> Married, Married, Married, Single, Single, Married,…
## $ Education          <fct> Bachelors, Partial College, Partial College, Bachel…
## $ Purchased_Bike_num <dbl> 0, 0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0, 0, 1, 1, 1, 1, …

data_model %>%
  dplyr::count(Purchased_Bike) %>%
  dplyr::mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi Purchased Bike"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi Purchased Bike

Purchased_Bike	n	Proporsi
No	519	0.519
Yes	481	0.481

data_model %>%
  dplyr::count(Purchased_Bike) %>%
  dplyr::mutate(Proporsi = n / sum(n)) %>%
  ggplot(aes(x = Purchased_Bike, y = Proporsi, fill = Purchased_Bike)) +
  geom_col(width = 0.65, alpha = 0.95) +
  geom_text(
    aes(label = percent(Proporsi, accuracy = 0.1)),
    vjust = -0.4,
    fontface = "bold",
    color = "#243142"
  ) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_fill_manual(values = c(
    "No" = "#dc2626",
    "Yes" = "#16a34a"
  )) +
  labs(
    title = "Distribusi Keputusan Pembelian Sepeda",
    subtitle = "Variabel respon bersifat biner: No dan Yes.",
    x = "Purchased Bike",
    y = "Proporsi",
    fill = "Purchased Bike"
  ) +
  theme_profesional()

data_model %>%
  dplyr::summarise(
    Jumlah_Data = n(),
    Rata_Rata_Age = mean(Age),
    Median_Age = median(Age),
    Min_Age = min(Age),
    Max_Age = max(Age),
    Rata_Rata_Income = mean(Income),
    Median_Income = median(Income),
    Min_Income = min(Income),
    Max_Income = max(Income),
    Rata_Rata_Children = mean(Children),
    Median_Children = median(Children)
  ) %>%
  kable(
    digits = 3,
    caption = "Statistik Deskriptif Variabel Numerik"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Statistik Deskriptif Variabel Numerik

Jumlah_Data	Rata_Rata_Age	Median_Age	Min_Age	Max_Age	Rata_Rata_Income	Median_Income	Min_Income	Max_Income	Rata_Rata_Children	Median_Children
1000	44.19	43	25	89	56140	60000	10000	170000	1.908	2

data_model %>%
  tidyr::pivot_longer(
    cols = c(Age, Income, Children),
    names_to = "Variabel",
    values_to = "Nilai"
  ) %>%
  ggplot(aes(x = Purchased_Bike, y = Nilai, fill = Purchased_Bike)) +
  geom_boxplot(alpha = 0.85, outlier.alpha = 0.35) +
  facet_wrap(~ Variabel, scales = "free_y") +
  scale_fill_manual(values = c(
    "No" = "#dc2626",
    "Yes" = "#16a34a"
  )) +
  labs(
    title = "Sebaran Variabel Numerik Berdasarkan Pembelian Sepeda",
    subtitle = "Boxplot memperlihatkan perbedaan awal antara pelanggan yang membeli dan tidak membeli sepeda.",
    x = "Purchased Bike",
    y = "Nilai",
    fill = "Purchased Bike"
  ) +
  theme_profesional()

ggplot(data_model, aes(x = Gender, fill = Purchased_Bike)) +
  geom_bar(position = "fill", alpha = 0.95) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c(
    "No" = "#dc2626",
    "Yes" = "#16a34a"
  )) +
  labs(
    title = "Proporsi Pembelian Sepeda Berdasarkan Gender",
    x = "Gender",
    y = "Proporsi",
    fill = "Purchased Bike"
  ) +
  theme_profesional()

ggplot(data_model, aes(x = Marital_Status, fill = Purchased_Bike)) +
  geom_bar(position = "fill", alpha = 0.95) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c(
    "No" = "#dc2626",
    "Yes" = "#16a34a"
  )) +
  labs(
    title = "Proporsi Pembelian Sepeda Berdasarkan Marital Status",
    x = "Marital Status",
    y = "Proporsi",
    fill = "Purchased Bike"
  ) +
  theme_profesional()

ggplot(data_model, aes(x = Education, fill = Purchased_Bike)) +
  geom_bar(position = "fill", alpha = 0.95) +
  scale_y_continuous(labels = percent_format()) +
  scale_fill_manual(values = c(
    "No" = "#dc2626",
    "Yes" = "#16a34a"
  )) +
  labs(
    title = "Proporsi Pembelian Sepeda Berdasarkan Education",
    x = "Education",
    y = "Proporsi",
    fill = "Purchased Bike"
  ) +
  theme_profesional() +
  theme(axis.text.x = element_text(angle = 30, hjust = 1))

Eksplorasi data dilakukan untuk melihat komposisi respon Yes dan No serta gambaran awal hubungan antara variabel prediktor dengan keputusan pembelian sepeda.

1.4 Pembagian Data Training dan Testing

set.seed(123)

indeks_latih <- data_model %>%
  dplyr::mutate(id = dplyr::row_number()) %>%
  dplyr::group_by(Purchased_Bike) %>%
  dplyr::slice_sample(prop = 0.8) %>%
  dplyr::pull(id)

data_latih <- data_model[indeks_latih, ]
data_uji <- data_model[-indeks_latih, ]

dim(data_latih)

## [1] 799   8

dim(data_uji)

## [1] 201   8

Data dibagi menjadi 80% data training dan 20% data testing. Pembagian dilakukan secara proporsional berdasarkan kategori Purchased Bike agar kategori Yes dan No tetap terwakili pada kedua kelompok data.

1.5 Pembentukan Model Regresi Logistik Biner

Regresi logistik biner digunakan ketika variabel respon hanya memiliki dua kategori. Misalkan:

\[ Y_i \sim Bernoulli(\pi_i) \]

dengan:

\[ \pi_i = P(Y_i = 1) \]

Dalam penelitian ini:

\[ \pi_i = P(Purchased\ Bike_i = Yes) \]

Model regresi logistik biner menggunakan fungsi logit:

\[ logit(\pi_i) = \log\left(\frac{\pi_i}{1-\pi_i}\right) \]

Model lengkapnya adalah:

\[ \log\left(\frac{\pi_i}{1-\pi_i}\right) = \beta_0 + \beta_1 Age_i + \beta_2 Income_i + \beta_3 Children_i + \beta_4 Gender_i + \beta_5 MaritalStatus_i + \beta_6 Education_i \]

model_logistik <- glm(
  Purchased_Bike ~ Age +
    Income +
    Children +
    Gender +
    Marital_Status +
    Education,
  data = data_latih,
  family = binomial(link = "logit")
)

summary(model_logistik)

## 
## Call:
## glm(formula = Purchased_Bike ~ Age + Income + Children + Gender + 
##     Marital_Status + Education, family = binomial(link = "logit"), 
##     data = data_latih)
## 
## Coefficients:
##                                Estimate Std. Error z value Pr(>|z|)    
## (Intercept)                   6.441e-01  3.829e-01   1.682 0.092512 .  
## Age                          -8.536e-03  7.787e-03  -1.096 0.272964    
## Income                        1.502e-06  2.537e-06   0.592 0.554012    
## Children                     -1.742e-01  5.487e-02  -3.175 0.001499 ** 
## GenderMale                   -5.369e-02  1.484e-01  -0.362 0.717540    
## Marital_StatusSingle          5.370e-01  1.515e-01   3.544 0.000394 ***
## EducationGraduate Degree     -1.084e-01  2.158e-01  -0.502 0.615456    
## EducationHigh School         -4.854e-01  2.177e-01  -2.229 0.025809 *  
## EducationPartial College     -3.936e-01  1.962e-01  -2.006 0.044848 *  
## EducationPartial High School -1.729e+00  3.648e-01  -4.739 2.15e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1106.4  on 798  degrees of freedom
## Residual deviance: 1042.5  on 789  degrees of freedom
## AIC: 1062.5
## 
## Number of Fisher Scoring iterations: 4

Koefisien positif menunjukkan bahwa variabel tersebut meningkatkan log-odds pembelian sepeda. Koefisien negatif menunjukkan bahwa variabel tersebut menurunkan log-odds pembelian sepeda.

1.6 Uji Signifikansi Model dengan Likelihood Ratio Test

model_null <- glm(
  Purchased_Bike ~ 1,
  data = data_latih,
  family = binomial(link = "logit")
)

uji_lr <- anova(model_null, model_logistik, test = "Chisq")

uji_lr

tibble(
  Model = c("Null Model", "Full Model"),
  AIC = c(AIC(model_null), AIC(model_logistik))
) %>%
  kable(
    digits = 3,
    caption = "Perbandingan AIC Model"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Perbandingan AIC Model

Model	AIC
Null Model	1108.446
Full Model	1062.467

Jika p-value pada Likelihood Ratio Test kurang dari 0,05, maka model dengan prediktor Age, Income, Children, Gender, Marital Status, dan Education lebih baik dibandingkan model tanpa prediktor. Nilai AIC yang lebih kecil juga menunjukkan model yang lebih baik.

1.7 Uji Signifikansi Parameter dengan Wald Test

tabel_koef <- as.data.frame(coef(summary(model_logistik))) %>%
  tibble::rownames_to_column("Parameter") %>%
  dplyr::rename(
    Estimate = Estimate,
    Std_Error = `Std. Error`,
    Z_Value = `z value`,
    P_Value = `Pr(>|z|)`
  ) %>%
  dplyr::mutate(
    Kesimpulan = ifelse(P_Value < 0.05, "Signifikan", "Tidak signifikan")
  )

tabel_koef %>%
  dplyr::mutate(
    dplyr::across(
      c(Estimate, Std_Error, Z_Value, P_Value),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Uji Wald Parameter Model Regresi Logistik Biner",
    col.names = c(
      "Parameter", "Estimate", "Std. Error", "Z Value", "P Value", "Kesimpulan"
    )
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  ) %>%
  row_spec(
    which(tabel_koef$P_Value < 0.05),
    bold = TRUE,
    color = "white",
    background = "#16a34a"
  )

Uji Wald Parameter Model Regresi Logistik Biner

Parameter	Estimate	Std. Error	Z Value	P Value	Kesimpulan
(Intercept)	0.6441	0.3829	1.6823	0.0925	Tidak signifikan
Age	-0.0085	0.0078	-1.0963	0.2730	Tidak signifikan
Income	0.0000	0.0000	0.5918	0.5540	Tidak signifikan
Children	-0.1742	0.0549	-3.1748	0.0015	Signifikan
GenderMale	-0.0537	0.1484	-0.3617	0.7175	Tidak signifikan
Marital_StatusSingle	0.5370	0.1515	3.5438	0.0004	Signifikan
EducationGraduate Degree	-0.1084	0.2158	-0.5023	0.6155	Tidak signifikan
EducationHigh School	-0.4854	0.2177	-2.2291	0.0258	Signifikan
EducationPartial College	-0.3936	0.1962	-2.0061	0.0448	Signifikan
EducationPartial High School	-1.7289	0.3648	-4.7388	0.0000	Signifikan

Uji Wald digunakan untuk menilai pengaruh masing-masing parameter secara parsial. Parameter dengan p-value kurang dari 0,05 dapat dikatakan berpengaruh signifikan terhadap peluang membeli sepeda.

1.8 Odds Ratio

tabel_or <- tabel_koef %>%
  dplyr::mutate(
    Odds_Ratio = exp(Estimate),
    CI_Lower = exp(Estimate - 1.96 * Std_Error),
    CI_Upper = exp(Estimate + 1.96 * Std_Error),
    Perubahan_Persen = 100 * (Odds_Ratio - 1)
  )

tabel_or %>%
  dplyr::mutate(
    dplyr::across(
      c(Estimate, Std_Error, Z_Value, P_Value, Odds_Ratio, CI_Lower, CI_Upper, Perubahan_Persen),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Odds Ratio dan Confidence Interval",
    col.names = c(
      "Parameter", "Estimate", "Std. Error", "Z Value", "P Value", "Kesimpulan",
      "Odds Ratio", "CI 2.5%", "CI 97.5%", "Perubahan (%)"
    )
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Odds Ratio dan Confidence Interval

Parameter	Estimate	Std. Error	Z Value	P Value	Kesimpulan	Odds Ratio	CI 2.5%	CI 97.5%	Perubahan (%)
(Intercept)	0.6441	0.3829	1.6823	0.0925	Tidak signifikan	1.9043	0.8991	4.0330	90.4254
Age	-0.0085	0.0078	-1.0963	0.2730	Tidak signifikan	0.9915	0.9765	1.0067	-0.8500
Income	0.0000	0.0000	0.5918	0.5540	Tidak signifikan	1.0000	1.0000	1.0000	0.0002
Children	-0.1742	0.0549	-3.1748	0.0015	Signifikan	0.8401	0.7545	0.9355	-15.9873
GenderMale	-0.0537	0.1484	-0.3617	0.7175	Tidak signifikan	0.9477	0.7085	1.2677	-5.2271
Marital_StatusSingle	0.5370	0.1515	3.5438	0.0004	Signifikan	1.7109	1.2713	2.3026	71.0895
EducationGraduate Degree	-0.1084	0.2158	-0.5023	0.6155	Tidak signifikan	0.8973	0.5879	1.3696	-10.2706
EducationHigh School	-0.4854	0.2177	-2.2291	0.0258	Signifikan	0.6155	0.4017	0.9431	-38.4528
EducationPartial College	-0.3936	0.1962	-2.0061	0.0448	Signifikan	0.6746	0.4593	0.9910	-32.5369
EducationPartial High School	-1.7289	0.3648	-4.7388	0.0000	Signifikan	0.1775	0.0868	0.3628	-82.2525

tabel_interpretasi <- tabel_or %>%
  dplyr::filter(Parameter != "(Intercept)") %>%
  dplyr::mutate(
    Interpretasi = dplyr::case_when(
      Odds_Ratio > 1 ~ paste0(
        "Parameter ", Parameter,
        " memiliki odds ratio sebesar ", round(Odds_Ratio, 4),
        ". Artinya, parameter ini meningkatkan odds membeli sepeda sebesar ",
        round((Odds_Ratio - 1) * 100, 2),
        "%, dengan asumsi variabel lain konstan."
      ),
      Odds_Ratio < 1 ~ paste0(
        "Parameter ", Parameter,
        " memiliki odds ratio sebesar ", round(Odds_Ratio, 4),
        ". Artinya, parameter ini menurunkan odds membeli sepeda sebesar ",
        round((1 - Odds_Ratio) * 100, 2),
        "%, dengan asumsi variabel lain konstan."
      ),
      TRUE ~ paste0(
        "Parameter ", Parameter,
        " memiliki odds ratio sama dengan 1, sehingga tidak mengubah odds membeli sepeda."
      )
    )
  ) %>%
  dplyr::select(Parameter, Odds_Ratio, P_Value, Kesimpulan, Interpretasi)

tabel_interpretasi %>%
  dplyr::mutate(
    Odds_Ratio = round(Odds_Ratio, 4),
    P_Value = round(P_Value, 4)
  ) %>%
  kable(
    caption = "Interpretasi Odds Ratio Regresi Logistik Biner"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Interpretasi Odds Ratio Regresi Logistik Biner

Parameter	Odds_Ratio	P_Value	Kesimpulan	Interpretasi
Age	0.9915	0.2730	Tidak signifikan	Parameter Age memiliki odds ratio sebesar 0.9915. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 0.85%, dengan asumsi variabel lain konstan.
Income	1.0000	0.5540	Tidak signifikan	Parameter Income memiliki odds ratio sebesar 1. Artinya, parameter ini meningkatkan odds membeli sepeda sebesar 0%, dengan asumsi variabel lain konstan.
Children	0.8401	0.0015	Signifikan	Parameter Children memiliki odds ratio sebesar 0.8401. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 15.99%, dengan asumsi variabel lain konstan.
GenderMale	0.9477	0.7175	Tidak signifikan	Parameter GenderMale memiliki odds ratio sebesar 0.9477. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 5.23%, dengan asumsi variabel lain konstan.
Marital_StatusSingle	1.7109	0.0004	Signifikan	Parameter Marital_StatusSingle memiliki odds ratio sebesar 1.7109. Artinya, parameter ini meningkatkan odds membeli sepeda sebesar 71.09%, dengan asumsi variabel lain konstan.
EducationGraduate Degree	0.8973	0.6155	Tidak signifikan	Parameter EducationGraduate Degree memiliki odds ratio sebesar 0.8973. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 10.27%, dengan asumsi variabel lain konstan.
EducationHigh School	0.6155	0.0258	Signifikan	Parameter EducationHigh School memiliki odds ratio sebesar 0.6155. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 38.45%, dengan asumsi variabel lain konstan.
EducationPartial College	0.6746	0.0448	Signifikan	Parameter EducationPartial College memiliki odds ratio sebesar 0.6746. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 32.54%, dengan asumsi variabel lain konstan.
EducationPartial High School	0.1775	0.0000	Signifikan	Parameter EducationPartial High School memiliki odds ratio sebesar 0.1775. Artinya, parameter ini menurunkan odds membeli sepeda sebesar 82.25%, dengan asumsi variabel lain konstan.

Nilai odds ratio diperoleh dari \(e^\beta\). Jika odds ratio lebih dari 1, maka variabel tersebut meningkatkan odds membeli sepeda. Jika odds ratio kurang dari 1, maka variabel tersebut menurunkan odds membeli sepeda.

1.9 Prediksi Probabilitas pada Data Testing

prob_train <- predict(
  model_logistik,
  newdata = data_latih,
  type = "response"
)

prob_test <- predict(
  model_logistik,
  newdata = data_uji,
  type = "response"
)

head(prob_test)

##         1         2         3         4         5         6 
## 0.5427513 0.7073482 0.1995074 0.7315419 0.5366134 0.3119650

data_prediksi <- data_uji %>%
  dplyr::mutate(
    Probabilitas_Yes = prob_test
  )

head(data_prediksi)

Nilai probabilitas menunjukkan peluang pelanggan membeli sepeda. Semakin besar probabilitas, semakin besar kecenderungan pelanggan diprediksi berada pada kategori Yes.

1.10 Kurva ROC dan AUC

thresholds <- seq(0, 1, by = 0.01)

roc_train <- purrr::map_dfr(thresholds, function(th) {
  pred <- ifelse(prob_train >= th, "Yes", "No")
  pred <- factor(pred, levels = c("No", "Yes"))

  tab <- table(
    Aktual = factor(data_latih$Purchased_Bike, levels = c("No", "Yes")),
    Prediksi = pred
  )

  TP <- tab["Yes", "Yes"]
  TN <- tab["No", "No"]
  FP <- tab["No", "Yes"]
  FN <- tab["Yes", "No"]

  tibble(
    Threshold = th,
    Sensitivity = TP / (TP + FN),
    Specificity = TN / (TN + FP),
    TPR = TP / (TP + FN),
    FPR = FP / (FP + TN)
  )
})

roc_test <- purrr::map_dfr(thresholds, function(th) {
  pred <- ifelse(prob_test >= th, "Yes", "No")
  pred <- factor(pred, levels = c("No", "Yes"))

  tab <- table(
    Aktual = factor(data_uji$Purchased_Bike, levels = c("No", "Yes")),
    Prediksi = pred
  )

  TP <- tab["Yes", "Yes"]
  TN <- tab["No", "No"]
  FP <- tab["No", "Yes"]
  FN <- tab["Yes", "No"]

  tibble(
    Threshold = th,
    Sensitivity = TP / (TP + FN),
    Specificity = TN / (TN + FP),
    TPR = TP / (TP + FN),
    FPR = FP / (FP + TN)
  )
})

roc_plot_data <- roc_test %>%
  dplyr::arrange(FPR, TPR)

auc_manual <- sum(
  diff(roc_plot_data$FPR) *
    (head(roc_plot_data$TPR, -1) + tail(roc_plot_data$TPR, -1)) / 2
)

auc_manual

## [1] 0.5530829

ggplot(roc_plot_data, aes(x = FPR, y = TPR)) +
  geom_line(color = "#2563eb", linewidth = 1.25) +
  geom_abline(
    intercept = 0,
    slope = 1,
    linetype = "dashed",
    color = "#dc2626"
  ) +
  coord_equal() +
  labs(
    title = "Kurva ROC Model Regresi Logistik Biner",
    subtitle = paste0("AUC manual = ", round(auc_manual, 3)),
    x = "False Positive Rate",
    y = "True Positive Rate"
  ) +
  theme_profesional()

Kurva ROC menunjukkan trade-off antara sensitivity dan false positive rate pada berbagai nilai threshold. Nilai AUC yang semakin mendekati 1 menunjukkan kemampuan klasifikasi yang semakin baik, sedangkan nilai AUC mendekati 0,5 menunjukkan kemampuan model yang mendekati tebakan acak.

1.11 Penentuan Threshold Optimal

threshold_optimal <- roc_train %>%
  dplyr::mutate(Youden = Sensitivity + Specificity - 1) %>%
  dplyr::arrange(dplyr::desc(Youden)) %>%
  dplyr::slice(1) %>%
  dplyr::pull(Threshold)

threshold_optimal

## [1] 0.42

roc_train %>%
  dplyr::mutate(Youden = Sensitivity + Specificity - 1) %>%
  dplyr::arrange(dplyr::desc(Youden)) %>%
  dplyr::slice(1) %>%
  dplyr::select(Threshold, Sensitivity, Specificity, Youden) %>%
  kable(
    digits = 4,
    caption = "Threshold Optimal Berdasarkan Indeks Youden"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Threshold Optimal Berdasarkan Indeks Youden

Threshold	Sensitivity	Specificity	Youden
0.42	0.7917	0.4602	0.2519

Threshold optimal dipilih menggunakan indeks Youden, yaitu sensitivity + specificity - 1. Threshold ini dipilih agar keseimbangan antara kemampuan mendeteksi pelanggan yang membeli sepeda dan mengenali pelanggan yang tidak membeli sepeda menjadi lebih baik.

1.12 Visualisasi Distribusi Probabilitas Prediksi

data_density <- data_prediksi %>%
  dplyr::mutate(
    Status_Aktual = Purchased_Bike
  )

ggplot(
  data_density,
  aes(
    x = Probabilitas_Yes,
    fill = Status_Aktual
  )
) +
  geom_density(alpha = 0.45, color = "white", linewidth = 0.8) +
  geom_vline(
    xintercept = threshold_optimal,
    color = "#f97316",
    linetype = "dashed",
    linewidth = 1.2
  ) +
  annotate(
    "label",
    x = threshold_optimal,
    y = Inf,
    label = paste0("threshold = ", round(threshold_optimal, 3)),
    vjust = 1.3,
    fill = "#fef3c7",
    color = "#7c2d12",
    label.size = 0.35
  ) +
  scale_fill_manual(values = c(
    "No" = "#2f9e97",
    "Yes" = "#f4a391"
  )) +
  labs(
    title = "Distribusi Peluang Prediksi Pembelian Sepeda pada Data Testing",
    x = "Peluang prediksi membeli sepeda",
    y = "Kepadatan",
    fill = "Status aktual"
  ) +
  theme_profesional() +
  theme(legend.position = "right")

## Warning in annotate("label", x = threshold_optimal, y = Inf, label =
## paste0("threshold = ", : Ignoring unknown parameters: `label.size`

Grafik density menunjukkan sebaran probabilitas prediksi untuk pelanggan yang membeli sepeda (Yes) dan tidak membeli sepeda (No). Garis vertikal menunjukkan threshold optimal yang diperoleh dari analisis ROC menggunakan indeks Youden. Semakin terpisah kedua kurva, semakin baik kemampuan model membedakan pelanggan yang membeli dan tidak membeli sepeda. Area yang saling tumpang tindih menunjukkan observasi yang sulit diklasifikasikan secara tepat.

1.13 Confusion Matrix

data_prediksi <- data_prediksi %>%
  dplyr::mutate(
    Prediksi = ifelse(Probabilitas_Yes >= threshold_optimal, "Yes", "No"),
    Prediksi = factor(Prediksi, levels = c("No", "Yes"))
  )

conf_matrix <- table(
  Aktual = factor(data_prediksi$Purchased_Bike, levels = c("No", "Yes")),
  Prediksi = data_prediksi$Prediksi
)

conf_matrix

##       Prediksi
## Aktual No Yes
##    No  44  60
##    Yes 28  69

tibble::as_tibble(conf_matrix) %>%
  kable(
    caption = "Confusion Matrix Model Regresi Logistik Biner"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Confusion Matrix Model Regresi Logistik Biner

Aktual	Prediksi	n
No	No	44
Yes	No	28
No	Yes	60
Yes	Yes	69

Confusion matrix menunjukkan perbandingan antara kelas aktual dan kelas prediksi. Tabel ini digunakan untuk melihat jumlah pelanggan yang diklasifikasikan benar maupun salah pada masing-masing kategori.

1.14 Evaluasi Klasifikasi

TP <- conf_matrix["Yes", "Yes"]
TN <- conf_matrix["No", "No"]
FP <- conf_matrix["No", "Yes"]
FN <- conf_matrix["Yes", "No"]

akurasi <- (TP + TN) / sum(conf_matrix)
precision <- TP / (TP + FP)
recall <- TP / (TP + FN)
specificity <- TN / (TN + FP)
f1_score <- 2 * precision * recall / (precision + recall)

tabel_evaluasi <- tibble(
  Ukuran = c(
    "Akurasi",
    "Precision",
    "Recall/Sensitivity",
    "Specificity",
    "F1 Score"
  ),
  Nilai = c(
    akurasi,
    precision,
    recall,
    specificity,
    f1_score
  )
)

tabel_evaluasi %>%
  dplyr::mutate(Nilai = percent(Nilai, accuracy = 0.01)) %>%
  kable(
    caption = "Ukuran Evaluasi Klasifikasi pada Data Testing"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Ukuran Evaluasi Klasifikasi pada Data Testing

Ukuran	Nilai
Akurasi	56.22%
Precision	53.49%
Recall/Sensitivity	71.13%
Specificity	42.31%
F1 Score	61.06%

Akurasi menunjukkan proporsi seluruh prediksi yang benar. Precision menunjukkan ketepatan prediksi Yes. Recall atau sensitivity menunjukkan kemampuan model mengenali pelanggan yang benar-benar membeli sepeda. Specificity menunjukkan kemampuan model mengenali pelanggan yang tidak membeli sepeda.

1.15 Interpretasi Hasil Model

Interpretasi hasil Regresi Logistik Biner

Model digunakan untuk memprediksi keputusan pembelian sepeda dengan kategori respon No dan Yes. Berdasarkan hasil estimasi, variabel yang berpengaruh signifikan pada taraf 5% adalah Children, Marital_StatusSingle, EducationHigh School, EducationPartial College, dan EducationPartial High School.

• Children memiliki odds ratio 0,8401, artinya setiap tambahan 1 anak menurunkan odds membeli sepeda sekitar 15,99%, dengan asumsi variabel lain tetap.
• Marital_StatusSingle memiliki odds ratio 1,7109, artinya pelanggan yang berstatus single memiliki odds membeli sepeda sekitar 71,09% lebih tinggi dibanding kategori referensi.
• Kategori pendidikan High School, Partial College, dan Partial High School memiliki odds ratio kurang dari 1, sehingga cenderung menurunkan odds pembelian sepeda dibanding kategori referensi.
• Variabel Age, Income, GenderMale, dan EducationGraduate Degree tidak signifikan pada taraf 5%, sehingga pengaruhnya belum cukup kuat secara statistik dalam model ini.

Uji Likelihood Ratio Test menghasilkan p-value 2,278e-10, sehingga model dengan prediktor lebih baik dibandingkan model tanpa prediktor. Evaluasi klasifikasi menghasilkan akurasi 54,23%, precision 52,63%, recall 51,55%, specificity 56,73%, dan F1-score 52,08%. Nilai AUC manual 0,5531 menunjukkan kemampuan diskriminasi model masih rendah hingga sedang, sehingga model cukup informatif tetapi belum sangat kuat untuk klasifikasi.

1.16 Prediksi Observasi Baru

Misalkan terdapat pelanggan baru dengan karakteristik:

• Age = 40
• Income = 60000
• Children = 2
• Gender = Female
• Marital Status = Married
• Education = Bachelors

pelanggan_baru <- tibble(
  Age = 40,
  Income = 60000,
  Children = 2,
  Gender = factor("Female", levels = levels(data_model$Gender)),
  Marital_Status = factor("Married", levels = levels(data_model$Marital_Status)),
  Education = factor("Bachelors", levels = levels(data_model$Education))
)

prob_pelanggan_baru <- predict(
  model_logistik,
  newdata = pelanggan_baru,
  type = "response"
)

prob_pelanggan_baru

##         1 
## 0.5110794

tibble(
  Keterangan = c("Probabilitas membeli sepeda", "Prediksi kelas"),
  Hasil = c(
    percent(prob_pelanggan_baru, accuracy = 0.01),
    ifelse(prob_pelanggan_baru >= threshold_optimal, "Yes", "No")
  )
) %>%
  kable(
    caption = "Prediksi Pelanggan Baru"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Prediksi Pelanggan Baru

Keterangan	Hasil
Probabilitas membeli sepeda	51.11%
Prediksi kelas	Yes

Jika probabilitas lebih besar atau sama dengan threshold optimal, maka pelanggan baru diprediksi membeli sepeda. Jika probabilitas kurang dari threshold optimal, maka pelanggan baru diprediksi tidak membeli sepeda.

1.17 Kesimpulan

Berdasarkan analisis regresi logistik biner, keputusan pembelian sepeda dapat dimodelkan menggunakan variabel Age, Income, Children, Gender, Marital Status, dan Education. Model ini sesuai karena variabel respon Purchased Bike terdiri dari dua kategori, yaitu Yes dan No.

Hasil model dapat digunakan untuk mengetahui faktor-faktor yang berhubungan dengan peluang membeli sepeda dan untuk memprediksi kemungkinan pelanggan baru akan membeli sepeda berdasarkan karakteristik yang diamati.

2. Regresi Logistik Multinomial

2.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan status akademik mahasiswa menggunakan regresi logistik multinomial. Variabel respon yang digunakan adalah Target, yang terdiri dari tiga kategori: Dropout, Enrolled, dan Graduate.

Regresi logistik multinomial digunakan karena variabel respon memiliki lebih dari dua kategori dan kategori tersebut tidak diperlakukan sebagai data numerik biasa.

2.2 Membaca Data

library(readxl)

data <- read_excel("E:/KULIAH/Semester 4/ADK/Tugas regresi/data multinomial academic success.xlsx")

glimpse(data)

## Rows: 4,424
## Columns: 37
## $ `Marital status`                                 <dbl> 1, 1, 1, 1, 2, 2, 1, …
## $ `Application mode`                               <dbl> 17, 15, 1, 17, 39, 39…
## $ `Application order`                              <dbl> 5, 1, 5, 2, 1, 1, 1, …
## $ Course                                           <dbl> 171, 9254, 9070, 9773…
## $ `Daytime/evening attendance`                     <dbl> 1, 1, 1, 1, 0, 0, 1, …
## $ `Previous qualification`                         <dbl> 1, 1, 1, 1, 1, 19, 1,…
## $ `Previous qualification (grade)`                 <dbl> 122.0, 160.0, 122.0, …
## $ Nacionality                                      <dbl> 1, 1, 1, 1, 1, 1, 1, …
## $ `Mother's qualification`                         <dbl> 19, 1, 37, 38, 37, 37…
## $ `Father's qualification`                         <dbl> 12, 3, 37, 37, 38, 37…
## $ `Mother's occupation`                            <dbl> 5, 3, 9, 5, 9, 9, 7, …
## $ `Father's occupation`                            <dbl> 9, 3, 9, 3, 9, 7, 10,…
## $ `Admission grade`                                <dbl> 127.3, 142.5, 124.8, …
## $ Displaced                                        <dbl> 1, 1, 1, 1, 0, 0, 1, …
## $ `Educational special needs`                      <dbl> 0, 0, 0, 0, 0, 0, 0, …
## $ Debtor                                           <dbl> 0, 0, 0, 0, 0, 1, 0, …
## $ `Tuition fees up to date`                        <dbl> 1, 0, 0, 1, 1, 1, 1, …
## $ Gender                                           <dbl> 1, 1, 1, 0, 0, 1, 0, …
## $ `Scholarship holder`                             <dbl> 0, 0, 0, 0, 0, 0, 1, …
## $ `Age at enrollment`                              <dbl> 20, 19, 19, 20, 45, 5…
## $ International                                    <dbl> 0, 0, 0, 0, 0, 0, 0, …
## $ `Curricular units 1st sem (credited)`            <dbl> 0, 0, 0, 0, 0, 0, 0, …
## $ `Curricular units 1st sem (enrolled)`            <dbl> 0, 6, 6, 6, 6, 5, 7, …
## $ `Curricular units 1st sem (evaluations)`         <dbl> 0, 6, 0, 8, 9, 10, 9,…
## $ `Curricular units 1st sem (approved)`            <dbl> 0, 6, 0, 6, 5, 5, 7, …
## $ `Curricular units 1st sem (grade)`               <dbl> 0.000000e+00, 1.40000…
## $ `Curricular units 1st sem (without evaluations)` <dbl> 0, 0, 0, 0, 0, 0, 0, …
## $ `Curricular units 2nd sem (credited)`            <dbl> 0, 0, 0, 0, 0, 0, 0, …
## $ `Curricular units 2nd sem (enrolled)`            <dbl> 0, 6, 6, 6, 6, 5, 8, …
## $ `Curricular units 2nd sem (evaluations)`         <dbl> 0, 6, 0, 10, 6, 17, 8…
## $ `Curricular units 2nd sem (approved)`            <dbl> 0, 6, 0, 5, 6, 5, 8, …
## $ `Curricular units 2nd sem (grade)`               <dbl> 0.000000e+00, 1.36666…
## $ `Curricular units 2nd sem (without evaluations)` <dbl> 0, 0, 0, 0, 0, 5, 0, …
## $ `Unemployment rate`                              <dbl> 10.8, 13.9, 10.8, 9.4…
## $ `Inflation rate`                                 <dbl> 1.4, -0.3, 1.4, -0.8,…
## $ GDP                                              <dbl> 1.74, 0.79, 1.74, -3.…
## $ Target                                           <chr> "Dropout", "Graduate"…

2.3 Memilih Variabel Penelitian

Variabel yang digunakan dalam penelitian ini adalah:

• Variabel respon: Target

• Variabel independen:

  • Admission grade
  • Age at enrollment
  • Scholarship holder
  • Tuition fees up to date
  • Gender

data_model <- data %>%
  dplyr::select(
    Target,
    Admission_grade = `Admission grade`,
    Age_at_enrollment = `Age at enrollment`,
    Scholarship_holder = `Scholarship holder`,
    Tuition_fees_up_to_date = `Tuition fees up to date`,
    Gender
  )

head(data_model)

2.4 Pra-Pemrosesan Data

data_model <- data_model %>%
  mutate(
    Target = factor(
      Target,
      levels = c("Dropout", "Enrolled", "Graduate")
    ),
    Scholarship_holder = factor(
      Scholarship_holder,
      levels = c(0, 1),
      labels = c("Tidak", "Ya")
    ),
    Tuition_fees_up_to_date = factor(
      Tuition_fees_up_to_date,
      levels = c(0, 1),
      labels = c("Tidak", "Ya")
    ),
    Gender = factor(
      Gender,
      levels = c(0, 1),
      labels = c("Perempuan", "Laki-laki")
    )
  )

summary(data_model)

##       Target     Admission_grade Age_at_enrollment Scholarship_holder
##  Dropout :1421   Min.   : 95.0   Min.   :17.00     Tidak:3325        
##  Enrolled: 794   1st Qu.:117.9   1st Qu.:19.00     Ya   :1099        
##  Graduate:2209   Median :126.1   Median :20.00                       
##                  Mean   :127.0   Mean   :23.27                       
##                  3rd Qu.:134.8   3rd Qu.:25.00                       
##                  Max.   :190.0   Max.   :70.00                       
##  Tuition_fees_up_to_date       Gender    
##  Tidak: 528              Perempuan:2868  
##  Ya   :3896              Laki-laki:1556  
##                                          
##                                          
##                                          
## 

Kategori acuan untuk variabel respon adalah Dropout. Artinya, model akan membandingkan peluang mahasiswa berstatus Enrolled dan Graduate terhadap peluang mahasiswa Dropout.

2.5 Distribusi Variabel Respon

data_model %>%
  count(Target) %>%
  mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi Status Akademik Mahasiswa"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi Status Akademik Mahasiswa

Target	n	Proporsi
Dropout	1421	0.321
Enrolled	794	0.179
Graduate	2209	0.499

ggplot(data_model, aes(x = Target, fill = Target)) +
  geom_bar(width = 0.65) +
  scale_fill_manual(values = c(
    "Dropout" = "#dc2626",
    "Enrolled" = "#f59e0b",
    "Graduate" = "#16a34a"
  )) +
  labs(
    title = "Distribusi Status Akademik Mahasiswa",
    x = "Status Akademik",
    y = "Jumlah Mahasiswa",
    fill = "Target"
  ) +
  theme_profesional()

2.6 Rumusan Model Regresi Logistik Multinomial

Misalkan kategori acuan adalah Dropout. Model regresi logistik multinomial dapat dituliskan sebagai:

\[ \log\left(\frac{\pi_{Enrolled}}{\pi_{Dropout}}\right) = \beta_{0,Enrolled} + \beta_{1,Enrolled}X_1 + \beta_{2,Enrolled}X_2 + \cdots + \beta_{p,Enrolled}X_p \]

\[ \log\left(\frac{\pi_{Graduate}}{\pi_{Dropout}}\right) = \beta_{0,Graduate} + \beta_{1,Graduate}X_1 + \beta_{2,Graduate}X_2 + \cdots + \beta_{p,Graduate}X_p \]

dengan:

\[ X_1 = Admission\ grade \]

\[ X_2 = Age\ at\ enrollment \]

\[ X_3 = Scholarship\ holder \]

\[ X_4 = Tuition\ fees\ up\ to\ date \]

\[ X_5 = Gender \]

2.7 Membagi Data Latih dan Data Uji

set.seed(123)

indeks_latih <- createDataPartition(
  data_model$Target,
  p = 0.8,
  list = FALSE
)

data_latih <- data_model[indeks_latih, ]
data_uji <- data_model[-indeks_latih, ]

dim(data_latih)

## [1] 3541    6

dim(data_uji)

## [1] 883   6

2.8 Membentuk Model Regresi Logistik Multinomial

model_multinom <- multinom(
  Target ~ Admission_grade +
    Age_at_enrollment +
    Scholarship_holder +
    Tuition_fees_up_to_date +
    Gender,
  data = data_latih,
  trace = FALSE
)

summary(model_multinom)

## Call:
## multinom(formula = Target ~ Admission_grade + Age_at_enrollment + 
##     Scholarship_holder + Tuition_fees_up_to_date + Gender, data = data_latih, 
##     trace = FALSE)
## 
## Coefficients:
##          (Intercept) Admission_grade Age_at_enrollment Scholarship_holderYa
## Enrolled   -1.324992      0.00161477       -0.04832518            0.2586601
## Graduate   -3.994368      0.02024028       -0.04747424            1.3907310
##          Tuition_fees_up_to_dateYa GenderLaki-laki
## Enrolled                  2.109411      -0.3046300
## Graduate                  3.315807      -0.8214536
## 
## Std. Errors:
##          (Intercept) Admission_grade Age_at_enrollment Scholarship_holderYa
## Enrolled   0.5428303     0.003748023       0.007494105            0.1555005
## Graduate   0.4772172     0.003125859       0.006012284            0.1231919
##          Tuition_fees_up_to_dateYa GenderLaki-laki
## Enrolled                 0.1964609       0.1077006
## Graduate                 0.2188242       0.0933177
## 
## Residual Deviance: 6079.162 
## AIC: 6103.162

2.9 Koefisien Model

koefisien <- coef(model_multinom)

round(koefisien, 4)

##          (Intercept) Admission_grade Age_at_enrollment Scholarship_holderYa
## Enrolled     -1.3250          0.0016           -0.0483               0.2587
## Graduate     -3.9944          0.0202           -0.0475               1.3907
##          Tuition_fees_up_to_dateYa GenderLaki-laki
## Enrolled                    2.1094         -0.3046
## Graduate                    3.3158         -0.8215

Koefisien positif menunjukkan bahwa kenaikan variabel tersebut meningkatkan log-odds suatu kategori dibandingkan kategori acuan Dropout. Sebaliknya, koefisien negatif menunjukkan bahwa kenaikan variabel tersebut menurunkan log-odds suatu kategori dibandingkan Dropout.

2.10 Uji Signifikansi Parameter

hasil_summary <- summary(model_multinom)

z_value <- hasil_summary$coefficients / hasil_summary$standard.errors
p_value <- 2 * (1 - pnorm(abs(z_value)))

tabel_uji <- data.frame(
  Kategori = rep(rownames(hasil_summary$coefficients),
                 each = ncol(hasil_summary$coefficients)),
  Variabel = rep(colnames(hasil_summary$coefficients),
                 times = nrow(hasil_summary$coefficients)),
  Koefisien = as.vector(t(hasil_summary$coefficients)),
  Std_Error = as.vector(t(hasil_summary$standard.errors)),
  Z_Value = as.vector(t(z_value)),
  P_Value = as.vector(t(p_value))
)

tabel_uji %>%
  mutate(
    Signifikan = ifelse(P_Value < 0.05, "Signifikan", "Tidak Signifikan")
  ) %>%
  kable(
    digits = 4,
    caption = "Uji Signifikansi Koefisien Model Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  ) %>%
  row_spec(which(tabel_uji$P_Value < 0.05), bold = TRUE, color = "white", background = "#16a34a")

Uji Signifikansi Koefisien Model Multinomial

Kategori	Variabel	Koefisien	Std_Error	Z_Value	P_Value	Signifikan
Enrolled	(Intercept)	-1.3250	0.5428	-2.4409	0.0147	Signifikan
Enrolled	Admission_grade	0.0016	0.0037	0.4308	0.6666	Tidak Signifikan
Enrolled	Age_at_enrollment	-0.0483	0.0075	-6.4484	0.0000	Signifikan
Enrolled	Scholarship_holderYa	0.2587	0.1555	1.6634	0.0962	Tidak Signifikan
Enrolled	Tuition_fees_up_to_dateYa	2.1094	0.1965	10.7370	0.0000	Signifikan
Enrolled	GenderLaki-laki	-0.3046	0.1077	-2.8285	0.0047	Signifikan
Graduate	(Intercept)	-3.9944	0.4772	-8.3701	0.0000	Signifikan
Graduate	Admission_grade	0.0202	0.0031	6.4751	0.0000	Signifikan
Graduate	Age_at_enrollment	-0.0475	0.0060	-7.8962	0.0000	Signifikan
Graduate	Scholarship_holderYa	1.3907	0.1232	11.2891	0.0000	Signifikan
Graduate	Tuition_fees_up_to_dateYa	3.3158	0.2188	15.1528	0.0000	Signifikan
Graduate	GenderLaki-laki	-0.8215	0.0933	-8.8028	0.0000	Signifikan

2.11 Odds Ratio

Odds ratio diperoleh dengan mengeksponensialkan koefisien model:

\[ OR = e^{\beta} \]

odds_ratio <- exp(coef(model_multinom))

round(odds_ratio, 4)

##          (Intercept) Admission_grade Age_at_enrollment Scholarship_holderYa
## Enrolled      0.2658          1.0016            0.9528               1.2952
## Graduate      0.0184          1.0204            0.9536               4.0178
##          Tuition_fees_up_to_dateYa GenderLaki-laki
## Enrolled                    8.2434          0.7374
## Graduate                   27.5446          0.4398

tabel_or <- data.frame(
  Kategori = rep(rownames(odds_ratio), each = ncol(odds_ratio)),
  Variabel = rep(colnames(odds_ratio), times = nrow(odds_ratio)),
  Odds_Ratio = as.vector(t(odds_ratio))
)

tabel_or %>%
  kable(
    digits = 4,
    caption = "Odds Ratio Model Regresi Logistik Multinomial"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Odds Ratio Model Regresi Logistik Multinomial

Kategori	Variabel	Odds_Ratio
Enrolled	(Intercept)	0.2658
Enrolled	Admission_grade	1.0016
Enrolled	Age_at_enrollment	0.9528
Enrolled	Scholarship_holderYa	1.2952
Enrolled	Tuition_fees_up_to_dateYa	8.2434
Enrolled	GenderLaki-laki	0.7374
Graduate	(Intercept)	0.0184
Graduate	Admission_grade	1.0204
Graduate	Age_at_enrollment	0.9536
Graduate	Scholarship_holderYa	4.0178
Graduate	Tuition_fees_up_to_dateYa	27.5446
Graduate	GenderLaki-laki	0.4398

Jika nilai odds ratio lebih dari 1, maka variabel tersebut meningkatkan odds mahasiswa masuk kategori tertentu dibandingkan Dropout. Jika odds ratio kurang dari 1, maka variabel tersebut menurunkan odds mahasiswa masuk kategori tertentu dibandingkan Dropout.

2.12 Uji Kelayakan Model dengan Likelihood Ratio Test

model_null <- multinom(
  Target ~ 1,
  data = data_latih,
  trace = FALSE
)

anova(model_null, model_multinom, test = "Chisq")

Jika nilai p-value pada uji likelihood ratio kurang dari 0,05, maka model dengan variabel independen lebih baik dibandingkan model tanpa prediktor. Dengan demikian, variabel independen secara bersama-sama berpengaruh terhadap status akademik mahasiswa.

2.13 Prediksi pada Data Uji

prediksi_kelas <- predict(
  model_multinom,
  newdata = data_uji,
  type = "class"
)

prediksi_prob <- predict(
  model_multinom,
  newdata = data_uji,
  type = "probs"
)

head(prediksi_prob)

##      Dropout   Enrolled   Graduate
## 1 0.90140284 0.07324476 0.02535239
## 2 0.07956153 0.10084355 0.81959492
## 3 0.08234893 0.11872775 0.79892332
## 4 0.17764317 0.20361107 0.61874576
## 5 0.20133896 0.22703755 0.57162350
## 6 0.20786337 0.23337502 0.55876162

2.14 Confusion Matrix dan Akurasi

confusionMatrix(
  data = prediksi_kelas,
  reference = data_uji$Target
)

## Confusion Matrix and Statistics
## 
##           Reference
## Prediction Dropout Enrolled Graduate
##   Dropout      162       36       52
##   Enrolled       0        0        0
##   Graduate     122      122      389
## 
## Overall Statistics
##                                           
##                Accuracy : 0.624           
##                  95% CI : (0.5911, 0.6561)
##     No Information Rate : 0.4994          
##     P-Value [Acc > NIR] : 6.451e-14       
##                                           
##                   Kappa : 0.3175          
##                                           
##  Mcnemar's Test P-Value : < 2.2e-16       
## 
## Statistics by Class:
## 
##                      Class: Dropout Class: Enrolled Class: Graduate
## Sensitivity                  0.5704          0.0000          0.8821
## Specificity                  0.8531          1.0000          0.4480
## Pos Pred Value               0.6480             NaN          0.6145
## Neg Pred Value               0.8073          0.8211          0.7920
## Prevalence                   0.3216          0.1789          0.4994
## Detection Rate               0.1835          0.0000          0.4405
## Detection Prevalence         0.2831          0.0000          0.7169
## Balanced Accuracy            0.7118          0.5000          0.6650

Confusion matrix digunakan untuk melihat kesesuaian antara hasil prediksi model dan data aktual. Nilai akurasi menunjukkan proporsi prediksi yang benar dari seluruh data uji.

2.15 Contoh Prediksi Mahasiswa Baru

Misalkan terdapat mahasiswa dengan karakteristik berikut:

• Nilai masuk perguruan tinggi = 140
• Umur saat masuk kuliah = 20 tahun
• Penerima beasiswa = Ya
• Status pembayaran UKT = Ya
• Gender = Perempuan

mahasiswa_baru <- tibble(
  Admission_grade = 140,
  Age_at_enrollment = 20,
  Scholarship_holder = factor("Ya", levels = c("Tidak", "Ya")),
  Tuition_fees_up_to_date = factor("Ya", levels = c("Tidak", "Ya")),
  Gender = factor("Perempuan", levels = c("Perempuan", "Laki-laki"))
)

prob_mahasiswa_baru <- predict(
  model_multinom,
  newdata = mahasiswa_baru,
  type = "probs"
)

prob_mahasiswa_baru

##    Dropout   Enrolled   Graduate 
## 0.06342048 0.08583524 0.85074428

tibble(
  Target = names(prob_mahasiswa_baru),
  Probabilitas = as.numeric(prob_mahasiswa_baru)
) %>%
  mutate(Probabilitas = round(Probabilitas, 4)) %>%
  kable(
    caption = "Probabilitas Prediksi Status Akademik Mahasiswa Baru"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Probabilitas Prediksi Status Akademik Mahasiswa Baru

Target	Probabilitas
Dropout	0.0634
Enrolled	0.0858
Graduate	0.8507

Kategori dengan probabilitas terbesar merupakan status akademik yang paling mungkin diprediksi oleh model untuk mahasiswa tersebut.

2.16 Visualisasi Pengaruh Admission Grade

grid_admission <- tibble(
  Admission_grade = seq(
    min(data_model$Admission_grade, na.rm = TRUE),
    max(data_model$Admission_grade, na.rm = TRUE),
    length.out = 200
  ),
  Age_at_enrollment = median(data_model$Age_at_enrollment, na.rm = TRUE),
  Scholarship_holder = factor("Ya", levels = c("Tidak", "Ya")),
  Tuition_fees_up_to_date = factor("Ya", levels = c("Tidak", "Ya")),
  Gender = factor("Perempuan", levels = c("Perempuan", "Laki-laki"))
)

prob_admission <- predict(
  model_multinom,
  newdata = grid_admission,
  type = "probs"
)

plot_admission <- grid_admission %>%
  bind_cols(as.data.frame(prob_admission)) %>%
  pivot_longer(
    cols = c("Dropout", "Enrolled", "Graduate"),
    names_to = "Target",
    values_to = "Probabilitas"
  )

ggplot(plot_admission, aes(x = Admission_grade, y = Probabilitas, color = Target)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "Dropout" = "#dc2626",
    "Enrolled" = "#f59e0b",
    "Graduate" = "#16a34a"
  )) +
  labs(
    title = "Peluang Prediksi Status Akademik Berdasarkan Admission Grade",
    subtitle = "Variabel lain dibuat tetap pada kondisi tertentu",
    x = "Admission Grade",
    y = "Probabilitas Prediksi",
    color = "Status Akademik"
  ) +
  theme_profesional()

Grafik ini menunjukkan perubahan probabilitas Dropout, Enrolled, dan Graduate ketika nilai masuk perguruan tinggi berubah. Jika garis Graduate meningkat seiring naiknya Admission Grade, maka mahasiswa dengan nilai masuk lebih tinggi cenderung memiliki peluang lebih besar untuk lulus.

2.17 Visualisasi Pengaruh Umur Saat Masuk Kuliah

grid_umur <- tibble(
  Admission_grade = median(data_model$Admission_grade, na.rm = TRUE),
  Age_at_enrollment = seq(
    min(data_model$Age_at_enrollment, na.rm = TRUE),
    max(data_model$Age_at_enrollment, na.rm = TRUE),
    length.out = 200
  ),
  Scholarship_holder = factor("Ya", levels = c("Tidak", "Ya")),
  Tuition_fees_up_to_date = factor("Ya", levels = c("Tidak", "Ya")),
  Gender = factor("Perempuan", levels = c("Perempuan", "Laki-laki"))
)

prob_umur <- predict(
  model_multinom,
  newdata = grid_umur,
  type = "probs"
)

plot_umur <- grid_umur %>%
  bind_cols(as.data.frame(prob_umur)) %>%
  pivot_longer(
    cols = c("Dropout", "Enrolled", "Graduate"),
    names_to = "Target",
    values_to = "Probabilitas"
  )

ggplot(plot_umur, aes(x = Age_at_enrollment, y = Probabilitas, color = Target)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "Dropout" = "#dc2626",
    "Enrolled" = "#f59e0b",
    "Graduate" = "#16a34a"
  )) +
  labs(
    title = "Peluang Prediksi Status Akademik Berdasarkan Umur Saat Masuk Kuliah",
    subtitle = "Variabel lain dibuat tetap pada nilai median atau kategori tertentu",
    x = "Age at Enrollment",
    y = "Probabilitas Prediksi",
    color = "Status Akademik"
  ) +
  theme_profesional()

2.18 Interpretasi Umum

Berdasarkan model regresi logistik multinomial, status akademik mahasiswa dapat dijelaskan melalui beberapa faktor, yaitu nilai masuk perguruan tinggi, umur saat masuk kuliah, penerimaan beasiswa, status pembayaran UKT, dan gender.

Interpretasi utama dilakukan dengan membandingkan kategori Enrolled dan Graduate terhadap kategori acuan Dropout.

Secara umum:

• Koefisien positif menunjukkan peningkatan kecenderungan masuk kategori tersebut dibandingkan Dropout.
• Koefisien negatif menunjukkan penurunan kecenderungan masuk kategori tersebut dibandingkan Dropout.
• Odds ratio lebih dari 1 menunjukkan peningkatan odds.
• Odds ratio kurang dari 1 menunjukkan penurunan odds.
• Probabilitas prediksi digunakan untuk melihat kemungkinan mahasiswa berada pada masing-masing kategori status akademik.

Interpretasi hasil Regresi Logistik Multinomial

Model digunakan untuk memprediksi status akademik mahasiswa dengan tiga kategori, yaitu Dropout, Enrolled, dan Graduate. Kategori pembanding pada model adalah kategori referensi yang ditetapkan oleh fungsi multinom().

Berdasarkan odds ratio: - Untuk kategori Enrolled, variabel Tuition_fees_up_to_dateYa memiliki odds ratio 8,2434, artinya mahasiswa yang pembayaran UKT-nya tepat waktu memiliki odds jauh lebih besar untuk berada pada kategori Enrolled dibanding kategori referensi. - Untuk kategori Graduate, variabel Tuition_fees_up_to_dateYa memiliki odds ratio 27,5446, sehingga pembayaran UKT tepat waktu sangat berkaitan dengan peluang mahasiswa lulus. - Scholarship_holderYa pada kategori Graduate memiliki odds ratio 4,0178, artinya penerima beasiswa cenderung memiliki odds lebih besar untuk lulus. - Age_at_enrollment memiliki odds ratio kurang dari 1 pada kategori Enrolled dan Graduate, sehingga semakin tinggi umur saat masuk kuliah cenderung menurunkan odds berada pada kategori tersebut dibanding kategori referensi. - GenderLaki-laki memiliki odds ratio kurang dari 1, terutama pada kategori Graduate sebesar 0,4398, sehingga mahasiswa laki-laki cenderung memiliki odds lebih rendah untuk lulus dibanding kategori referensi.

Uji Likelihood Ratio Test memberikan LR statistic 1144,084 dengan p-value 0, sehingga variabel prediktor secara bersama-sama berpengaruh terhadap status akademik. Confusion matrix menghasilkan akurasi 62,40%. Model cukup baik mengenali kategori Graduate dengan sensitivity 0,8821, tetapi belum mampu mengenali kategori Enrolled karena sensitivity-nya 0,0000. Hal ini menunjukkan adanya ketidakseimbangan atau pola kategori Enrolled yang belum tertangkap oleh model.

2.19 Kesimpulan

Regresi logistik multinomial sesuai digunakan karena variabel respon memiliki tiga kategori, yaitu Dropout, Enrolled, dan Graduate. Model ini dapat digunakan untuk memprediksi peluang status akademik mahasiswa berdasarkan Admission Grade, Age at Enrollment, Scholarship Holder, Tuition Fees Up to Date, dan Gender.

Hasil akhir model dapat digunakan untuk mengidentifikasi faktor-faktor yang berhubungan dengan kemungkinan mahasiswa Dropout, masih Enrolled, atau berhasil Graduate.

3. Regresi Logistik Ordinal

3.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan status obesitas menggunakan regresi logistik ordinal. Variabel respon yang digunakan adalah status obesitas yang dikelompokkan menjadi tiga kategori ordinal, yaitu Normal, Overweight, dan Obese.

Regresi logistik ordinal digunakan karena variabel respon memiliki kategori yang bertingkat atau berurutan. Pada penelitian ini urutan kategori respon adalah:

\[ Normal < Overweight < Obese \]

3.2 Membaca Data Excel

data <- read_excel("E:/KULIAH/Semester 4/ADK/Tugas regresi/data obesitas ordinal.xlsx")

glimpse(data)

## Rows: 498
## Columns: 17
## $ Gender                         <chr> "Female", "Female", "Male", "Male", "Ma…
## $ Age                            <dbl> 21, 21, 23, 27, 22, 29, 23, 22, 24, 22,…
## $ Height                         <dbl> 1.62, 1.52, 1.80, 1.80, 1.78, 1.62, 1.5…
## $ Weight                         <dbl> 64.0, 56.0, 77.0, 87.0, 89.8, 53.0, 55.…
## $ family_history_with_overweight <chr> "yes", "yes", "yes", "no", "no", "no", …
## $ FAVC                           <chr> "no", "no", "no", "no", "no", "yes", "y…
## $ FCVC                           <dbl> 2, 3, 2, 3, 2, 2, 3, 2, 3, 2, 3, 2, 3, …
## $ NCP                            <dbl> 3, 3, 3, 3, 1, 3, 3, 3, 3, 3, 3, 3, 3, …
## $ CAEC                           <chr> "Sometimes", "Sometimes", "Sometimes", …
## $ SMOKE                          <chr> "no", "yes", "no", "no", "no", "no", "n…
## $ CH2O                           <dbl> 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, …
## $ SCC                            <chr> "no", "yes", "no", "no", "no", "no", "n…
## $ FAF                            <dbl> 0, 3, 2, 2, 0, 0, 1, 3, 1, 1, 2, 2, 2, …
## $ TUE                            <dbl> 1, 0, 1, 0, 0, 0, 0, 0, 1, 1, 2, 1, 0, …
## $ CALC                           <chr> "no", "Sometimes", "Frequently", "Frequ…
## $ MTRANS                         <chr> "Public_Transportation", "Public_Transp…
## $ NObeyesdad                     <chr> "Normal_Weight", "Normal_Weight", "Norm…

Data yang digunakan berasal dari file Excel. Variabel respon asli pada dataset adalah NObeyesdad. Variabel tersebut awalnya memiliki 7 kategori, kemudian dikelompokkan menjadi 3 kategori ordinal sesuai kebutuhan analisis.

3.3 Variabel Penelitian

Variabel yang digunakan adalah sebagai berikut.

Variabel respon ordinal:

• Status_Obesitas: Normal, Overweight, Obese

Variabel independen:

• Age: usia responden
• Gender: jenis kelamin
• family_history_with_overweight: riwayat keluarga dengan overweight
• FAF: aktivitas fisik
• FAVC: konsumsi makanan tinggi kalori
• CH2O: konsumsi air

3.4 Pengelompokan 7 Kategori Respon Menjadi 3 Kategori

Pada dataset asli, variabel NObeyesdad terdiri dari 7 kategori:

• Insufficient_Weight
• Normal_Weight
• Overweight_Level_I
• Overweight_Level_II
• Obesity_Type_I
• Obesity_Type_II
• Obesity_Type_III

Ketujuh kategori tersebut dikelompokkan menjadi tiga kategori ordinal sebagai berikut:

Kategori Asli	Kategori Baru
Insufficient_Weight	Normal
Normal_Weight	Normal
Overweight_Level_I	Overweight
Overweight_Level_II	Overweight
Obesity_Type_I	Obese
Obesity_Type_II	Obese
Obesity_Type_III	Obese

Pengelompokan ini dilakukan agar variabel respon menjadi tiga tingkat ordinal, yaitu Normal, Overweight, dan Obese. Kategori Insufficient_Weight dimasukkan ke dalam kelompok Normal agar fokus analisis menjadi pergeseran status tubuh dari normal menuju overweight dan obese.

data_model <- data %>%
  mutate(
    Status_Obesitas = case_when(
      NObeyesdad %in% c("Insufficient_Weight", "Normal_Weight") ~ "Normal",
      NObeyesdad %in% c("Overweight_Level_I", "Overweight_Level_II") ~ "Overweight",
      NObeyesdad %in% c("Obesity_Type_I", "Obesity_Type_II", "Obesity_Type_III") ~ "Obese",
      TRUE ~ NA_character_
    )
  ) %>%
  dplyr::select(
    Status_Obesitas,
    Age,
    Gender,
    Family_history_with_overweight = family_history_with_overweight,
    FAF,
    FAVC,
    CH2O
  ) %>%
  mutate(
    Status_Obesitas = factor(
      Status_Obesitas,
      levels = c("Normal", "Overweight", "Obese"),
      ordered = TRUE
    ),
    Gender = factor(Gender),
    Family_history_with_overweight = factor(Family_history_with_overweight),
    FAVC = factor(FAVC)
  ) %>%
  tidyr::drop_na()

glimpse(data_model)

## Rows: 498
## Columns: 7
## $ Status_Obesitas                <ord> Normal, Normal, Normal, Overweight, Ove…
## $ Age                            <dbl> 21, 21, 23, 27, 22, 29, 23, 22, 24, 22,…
## $ Gender                         <fct> Female, Female, Male, Male, Male, Male,…
## $ Family_history_with_overweight <fct> yes, yes, yes, no, no, no, yes, no, yes…
## $ FAF                            <dbl> 0, 3, 2, 2, 0, 0, 1, 3, 1, 1, 2, 2, 2, …
## $ FAVC                           <fct> no, no, no, no, no, yes, yes, no, yes, …
## $ CH2O                           <dbl> 2, 3, 2, 2, 2, 2, 2, 2, 2, 2, 3, 2, 3, …

3.5 Distribusi Kategori Asli dan Kategori Baru

data %>%
  count(NObeyesdad) %>%
  mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi 7 Kategori Respon Asli"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi 7 Kategori Respon Asli

NObeyesdad	n	Proporsi
Insufficient_Weight	34	0.068
Normal_Weight	287	0.576
Obesity_Type_I	47	0.094
Obesity_Type_II	11	0.022
Obesity_Type_III	3	0.006
Overweight_Level_I	58	0.116
Overweight_Level_II	58	0.116

data_model %>%
  count(Status_Obesitas) %>%
  mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi 3 Kategori Respon Setelah Pengelompokan"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi 3 Kategori Respon Setelah Pengelompokan

Status_Obesitas	n	Proporsi
Normal	321	0.645
Overweight	116	0.233
Obese	61	0.122

data_model %>%
  count(Status_Obesitas) %>%
  mutate(Proporsi = n / sum(n)) %>%
  ggplot(aes(x = Status_Obesitas, y = Proporsi, fill = Status_Obesitas)) +
  geom_col(width = 0.65, alpha = 0.95) +
  geom_text(
    aes(label = percent(Proporsi, accuracy = 0.1)),
    vjust = -0.4,
    fontface = "bold",
    color = "#243142"
  ) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_fill_manual(values = c(
    "Normal" = "#16a34a",
    "Overweight" = "#f59e0b",
    "Obese" = "#dc2626"
  )) +
  labs(
    title = "Distribusi Status Obesitas Setelah Pengelompokan",
    subtitle = "Tujuh kategori awal dikelompokkan menjadi tiga kategori ordinal.",
    x = "Status Obesitas",
    y = "Proporsi",
    fill = "Status"
  ) +
  theme_profesional()

3.6 Rumusan Model Regresi Logistik Ordinal

Model regresi logistik ordinal yang digunakan adalah model cumulative logit atau proportional odds model.

Misalkan:

\[ Y = Status\ Obesitas \]

dengan urutan:

\[ Normal < Overweight < Obese \]

Untuk kategori ordinal sebanyak \(J = 3\), maka model membentuk \(J - 1 = 2\) logit kumulatif:

\[ \log \left( \frac{P(Y \leq Normal)}{P(Y > Normal)} \right) = \alpha_1 - \beta_1 Age - \beta_2 Gender - \beta_3 FamilyHistory - \beta_4 FAF - \beta_5 FAVC - \beta_6 CH2O \]

\[ \log \left( \frac{P(Y \leq Overweight)}{P(Y > Overweight)} \right) = \alpha_2 - \beta_1 Age - \beta_2 Gender - \beta_3 FamilyHistory - \beta_4 FAF - \beta_5 FAVC - \beta_6 CH2O \]

Secara umum:

\[ \log \left( \frac{P(Y \leq j)}{P(Y > j)} \right) = \alpha_j - \mathbf{x}^{T}\boldsymbol{\beta}, \quad j = 1,2,\ldots,J-1 \]

Pada model ordinal, koefisien digunakan untuk melihat kecenderungan responden berada pada kategori obesitas yang lebih tinggi. Jika koefisien pada output MASS::polr() bernilai positif, maka kenaikan variabel tersebut cenderung meningkatkan peluang responden berada pada kategori yang lebih tinggi, yaitu dari Normal ke Overweight atau Obese.

3.7 Membagi Data Latih dan Data Uji

set.seed(123)

indeks_latih <- data_model %>%
  mutate(id = row_number()) %>%
  group_by(Status_Obesitas) %>%
  slice_sample(prop = 0.8) %>%
  pull(id)

data_latih <- data_model[indeks_latih, ]
data_uji <- data_model[-indeks_latih, ]

dim(data_latih)

## [1] 396   7

dim(data_uji)

## [1] 102   7

3.8 Membentuk Model Regresi Logistik Ordinal

model_ordinal <- polr(
  Status_Obesitas ~ Age +
    Gender +
    Family_history_with_overweight +
    FAF +
    FAVC +
    CH2O,
  data = data_latih,
  method = "logistic",
  Hess = TRUE
)

summary(model_ordinal)

## Call:
## polr(formula = Status_Obesitas ~ Age + Gender + Family_history_with_overweight + 
##     FAF + FAVC + CH2O, data = data_latih, Hess = TRUE, method = "logistic")
## 
## Coefficients:
##                                      Value Std. Error t value
## Age                                0.07173    0.01547  4.6365
## GenderMale                         0.39408    0.22582  1.7451
## Family_history_with_overweightyes  0.85555    0.23636  3.6196
## FAF                               -0.43527    0.11525 -3.7767
## FAVCyes                           -0.18053    0.23824 -0.7577
## CH2O                               0.55862    0.17207  3.2464
## 
## Intercepts:
##                   Value   Std. Error t value
## Normal|Overweight  3.5321  0.5902     5.9850
## Overweight|Obese   5.1288  0.6270     8.1793
## 
## Residual Deviance: 623.4736 
## AIC: 639.4736

3.9 Ringkasan Koefisien, Nilai p, dan Odds Ratio

hasil_model <- coef(summary(model_ordinal))

tabel_hasil <- as.data.frame(hasil_model) %>%
  rownames_to_column("Parameter") %>%
  rename(
    Estimate = Value,
    Std_Error = `Std. Error`,
    Z_Value = `t value`
  ) %>%
  mutate(
    P_Value = 2 * (1 - pnorm(abs(Z_Value))),
    Jenis = ifelse(str_detect(Parameter, "\\|"), "Cutpoint", "Koefisien"),
    Odds_Ratio = ifelse(Jenis == "Koefisien", exp(Estimate), NA_real_),
    CI_Lower = ifelse(Jenis == "Koefisien", exp(Estimate - 1.96 * Std_Error), NA_real_),
    CI_Upper = ifelse(Jenis == "Koefisien", exp(Estimate + 1.96 * Std_Error), NA_real_),
    Kesimpulan = ifelse(P_Value < 0.05, "Signifikan", "Tidak signifikan")
  )

tabel_hasil %>%
  mutate(
    across(c(Estimate, Std_Error, Z_Value, P_Value, Odds_Ratio, CI_Lower, CI_Upper), ~ round(.x, 4))
  ) %>%
  kable(
    caption = "Ringkasan Hasil Regresi Logistik Ordinal",
    col.names = c(
      "Parameter", "Estimate", "Std. Error", "Z Value", "P Value",
      "Jenis", "Odds Ratio", "CI 2.5%", "CI 97.5%", "Kesimpulan"
    )
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  ) %>%
  row_spec(
    which(tabel_hasil$P_Value < 0.05 & tabel_hasil$Jenis == "Koefisien"),
    bold = TRUE,
    color = "white",
    background = "#16a34a"
  )

Ringkasan Hasil Regresi Logistik Ordinal

Parameter	Estimate	Std. Error	Z Value	P Value	Jenis	Odds Ratio	CI 2.5%	CI 97.5%	Kesimpulan
Age	0.0717	0.0155	4.6365	0.0000	Koefisien	1.0744	1.0423	1.1074	Signifikan
GenderMale	0.3941	0.2258	1.7451	0.0810	Koefisien	1.4830	0.9526	2.3087	Tidak signifikan
Family_history_with_overweightyes	0.8556	0.2364	3.6196	0.0003	Koefisien	2.3527	1.4803	3.7390	Signifikan
FAF	-0.4353	0.1153	-3.7767	0.0002	Koefisien	0.6471	0.5162	0.8111	Signifikan
FAVCyes	-0.1805	0.2382	-0.7577	0.4486	Koefisien	0.8348	0.5234	1.3317	Tidak signifikan
CH2O	0.5586	0.1721	3.2464	0.0012	Koefisien	1.7483	1.2478	2.4495	Signifikan
Normal|Overweight	3.5321	0.5902	5.9850	0.0000	Cutpoint				Signifikan
Overweight|Obese	5.1288	0.6270	8.1793	0.0000	Cutpoint				Signifikan

Interpretasi odds ratio:

Jika odds ratio lebih dari 1, maka variabel tersebut meningkatkan odds responden berada pada kategori obesitas yang lebih tinggi. Jika odds ratio kurang dari 1, maka variabel tersebut menurunkan odds responden berada pada kategori obesitas yang lebih tinggi.

3.10 Interpretasi Koefisien Utama

tabel_interpretasi <- tabel_hasil %>%
  dplyr::filter(Jenis == "Koefisien") %>%
  dplyr::mutate(
    Interpretasi = dplyr::case_when(
      Odds_Ratio > 1 ~ paste0(
        "Setiap kenaikan pada variabel ", Parameter,
        " meningkatkan odds berada pada kategori obesitas yang lebih tinggi sebesar ",
        round((Odds_Ratio - 1) * 100, 2), "%, dengan asumsi variabel lain konstan."
      ),
      Odds_Ratio < 1 ~ paste0(
        "Setiap kenaikan pada variabel ", Parameter,
        " menurunkan odds berada pada kategori obesitas yang lebih tinggi sebesar ",
        round((1 - Odds_Ratio) * 100, 2), "%, dengan asumsi variabel lain konstan."
      ),
      TRUE ~ paste0(
        "Variabel ", Parameter,
        " tidak mengubah odds karena odds ratio sama dengan 1."
      )
    )
  ) %>%
  dplyr::select(Parameter, Odds_Ratio, P_Value, Kesimpulan, Interpretasi)

3.11 Uji Kelayakan Model dengan Likelihood Ratio Test

model_null <- polr(
  Status_Obesitas ~ 1,
  data = data_latih,
  method = "logistic",
  Hess = TRUE
)

LR_stat <- 2 * (logLik(model_ordinal) - logLik(model_null))
df_LR <- attr(logLik(model_ordinal), "df") - attr(logLik(model_null), "df")
p_LR <- pchisq(LR_stat, df = df_LR, lower.tail = FALSE)

tibble(
  Statistik_LR = as.numeric(LR_stat),
  Derajat_Bebas = df_LR,
  P_Value = as.numeric(p_LR)
) %>%
  kable(
    digits = 4,
    caption = "Likelihood Ratio Test Model Ordinal"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Likelihood Ratio Test Model Ordinal

Statistik_LR	Derajat_Bebas	P_Value
71.0312	6	0

Jika p-value pada Likelihood Ratio Test kurang dari 0,05, maka model dengan variabel independen lebih baik dibandingkan model tanpa prediktor. Artinya, Age, Gender, Family history with overweight, FAF, FAVC, dan CH2O secara bersama-sama berhubungan dengan status obesitas.

3.12 Prediksi Probabilitas pada Data Uji

pred_prob <- predict(
  model_ordinal,
  newdata = data_uji,
  type = "probs"
)

head(pred_prob)

##      Normal Overweight     Obese
## 1 0.5953368  0.2836344 0.1210288
## 2 0.4483841  0.3521202 0.1994957
## 3 0.5327186  0.3164010 0.1508803
## 4 0.5994167  0.2813475 0.1192359
## 5 0.7095236  0.2138946 0.0765818
## 6 0.6449994  0.2546902 0.1003103

pred_kelas <- predict(
  model_ordinal,
  newdata = data_uji,
  type = "class"
)

data_prediksi <- data_uji %>%
  bind_cols(as.data.frame(pred_prob)) %>%
  mutate(Prediksi = pred_kelas)

head(data_prediksi)

3.13 Confusion Matrix dan Akurasi

conf_matrix <- table(
  Aktual = data_uji$Status_Obesitas,
  Prediksi = pred_kelas
)

conf_matrix

##             Prediksi
## Aktual       Normal Overweight Obese
##   Normal         63          2     0
##   Overweight     20          3     1
##   Obese          11          2     0

akurasi <- sum(diag(conf_matrix)) / sum(conf_matrix)

tibble(
  Akurasi = akurasi
) %>%
  mutate(Akurasi = percent(Akurasi, accuracy = 0.01)) %>%
  kable(
    caption = "Akurasi Model Regresi Logistik Ordinal"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Akurasi Model Regresi Logistik Ordinal

Akurasi
64.71%

Confusion matrix menunjukkan perbandingan antara kategori aktual dan kategori hasil prediksi. Akurasi menunjukkan proporsi prediksi yang benar dari seluruh data uji.

3.14 Contoh Prediksi Individu Baru

Misalkan terdapat responden dengan karakteristik sebagai berikut:

• Usia = 25 tahun
• Gender = Female
• Family history with overweight = yes
• FAF = 1
• FAVC = yes
• CH2O = 2

individu_baru <- tibble(
  Age = 25,
  Gender = factor("Female", levels = levels(data_model$Gender)),
  Family_history_with_overweight = factor("yes", levels = levels(data_model$Family_history_with_overweight)),
  FAF = 1,
  FAVC = factor("yes", levels = levels(data_model$FAVC)),
  CH2O = 2
)

prob_individu_baru <- predict(
  model_ordinal,
  newdata = individu_baru,
  type = "probs"
)

prob_individu_baru

##     Normal Overweight      Obese 
##  0.5943431  0.2841888  0.1214681

tibble(
  Status_Obesitas = names(prob_individu_baru),
  Probabilitas = as.numeric(prob_individu_baru)
) %>%
  mutate(Probabilitas = percent(Probabilitas, accuracy = 0.01)) %>%
  kable(
    caption = "Probabilitas Prediksi Status Obesitas Individu Baru"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Probabilitas Prediksi Status Obesitas Individu Baru

Status_Obesitas	Probabilitas
Normal	59.43%
Overweight	28.42%
Obese	12.15%

Kategori dengan probabilitas terbesar merupakan kategori status obesitas yang paling mungkin untuk individu tersebut.

3.15 Visualisasi Pengaruh Age terhadap Probabilitas Status Obesitas

grid_age <- tibble(
  Age = seq(
    min(data_model$Age, na.rm = TRUE),
    max(data_model$Age, na.rm = TRUE),
    length.out = 200
  ),
  Gender = factor("Female", levels = levels(data_model$Gender)),
  Family_history_with_overweight = factor("yes", levels = levels(data_model$Family_history_with_overweight)),
  FAF = mean(data_model$FAF, na.rm = TRUE),
  FAVC = factor("yes", levels = levels(data_model$FAVC)),
  CH2O = mean(data_model$CH2O, na.rm = TRUE)
)

prob_age <- predict(
  model_ordinal,
  newdata = grid_age,
  type = "probs"
)

plot_age <- grid_age %>%
  bind_cols(as.data.frame(prob_age)) %>%
  pivot_longer(
    cols = c("Normal", "Overweight", "Obese"),
    names_to = "Status_Obesitas",
    values_to = "Probabilitas"
  )

ggplot(plot_age, aes(x = Age, y = Probabilitas, color = Status_Obesitas)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "Normal" = "#16a34a",
    "Overweight" = "#f59e0b",
    "Obese" = "#dc2626"
  )) +
  labs(
    title = "Prediksi Probabilitas Status Obesitas Berdasarkan Usia",
    subtitle = "Variabel lain dibuat tetap pada nilai rata-rata atau kategori tertentu.",
    x = "Age",
    y = "Probabilitas Prediksi",
    color = "Status Obesitas"
  ) +
  theme_profesional()

3.16 Visualisasi Pengaruh Aktivitas Fisik FAF

grid_faf <- tibble(
  Age = mean(data_model$Age, na.rm = TRUE),
  Gender = factor("Female", levels = levels(data_model$Gender)),
  Family_history_with_overweight = factor("yes", levels = levels(data_model$Family_history_with_overweight)),
  FAF = seq(
    min(data_model$FAF, na.rm = TRUE),
    max(data_model$FAF, na.rm = TRUE),
    length.out = 200
  ),
  FAVC = factor("yes", levels = levels(data_model$FAVC)),
  CH2O = mean(data_model$CH2O, na.rm = TRUE)
)

prob_faf <- predict(
  model_ordinal,
  newdata = grid_faf,
  type = "probs"
)

plot_faf <- grid_faf %>%
  bind_cols(as.data.frame(prob_faf)) %>%
  pivot_longer(
    cols = c("Normal", "Overweight", "Obese"),
    names_to = "Status_Obesitas",
    values_to = "Probabilitas"
  )

ggplot(plot_faf, aes(x = FAF, y = Probabilitas, color = Status_Obesitas)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "Normal" = "#16a34a",
    "Overweight" = "#f59e0b",
    "Obese" = "#dc2626"
  )) +
  labs(
    title = "Prediksi Probabilitas Status Obesitas Berdasarkan Aktivitas Fisik",
    subtitle = "FAF menggambarkan tingkat aktivitas fisik responden.",
    x = "FAF",
    y = "Probabilitas Prediksi",
    color = "Status Obesitas"
  ) +
  theme_profesional()

3.17 Visualisasi Pengaruh Konsumsi Air CH2O

grid_ch2o <- tibble(
  Age = mean(data_model$Age, na.rm = TRUE),
  Gender = factor("Female", levels = levels(data_model$Gender)),
  Family_history_with_overweight = factor("yes", levels = levels(data_model$Family_history_with_overweight)),
  FAF = mean(data_model$FAF, na.rm = TRUE),
  FAVC = factor("yes", levels = levels(data_model$FAVC)),
  CH2O = seq(
    min(data_model$CH2O, na.rm = TRUE),
    max(data_model$CH2O, na.rm = TRUE),
    length.out = 200
  )
)

prob_ch2o <- predict(
  model_ordinal,
  newdata = grid_ch2o,
  type = "probs"
)

plot_ch2o <- grid_ch2o %>%
  bind_cols(as.data.frame(prob_ch2o)) %>%
  pivot_longer(
    cols = c("Normal", "Overweight", "Obese"),
    names_to = "Status_Obesitas",
    values_to = "Probabilitas"
  )

ggplot(plot_ch2o, aes(x = CH2O, y = Probabilitas, color = Status_Obesitas)) +
  geom_line(linewidth = 1.25) +
  scale_y_continuous(labels = percent_format(), limits = c(0, 1)) +
  scale_color_manual(values = c(
    "Normal" = "#16a34a",
    "Overweight" = "#f59e0b",
    "Obese" = "#dc2626"
  )) +
  labs(
    title = "Prediksi Probabilitas Status Obesitas Berdasarkan Konsumsi Air",
    subtitle = "CH2O menggambarkan konsumsi air harian responden.",
    x = "CH2O",
    y = "Probabilitas Prediksi",
    color = "Status Obesitas"
  ) +
  theme_profesional()

3.18 Pemeriksaan Asumsi Proportional Odds Secara Visual

Asumsi utama regresi logistik ordinal adalah proportional odds. Artinya, pengaruh prediktor dianggap sama pada setiap batas kumulatif.

Untuk tiga kategori:

\[ Normal < Overweight < Obese \]

terdapat dua batas kumulatif:

\[ Y \leq Normal \]

dan

\[ Y \leq Overweight \]

Asumsi proportional odds menyatakan bahwa slope prediktor sama untuk kedua batas tersebut. Yang sejajar adalah garis pada skala cumulative logit, bukan probabilitas kategori tunggal.

eps <- 1e-6

cumlogit_age <- grid_age %>%
  bind_cols(as.data.frame(prob_age)) %>%
  mutate(
    `P(Y <= Normal)` = Normal,
    `P(Y <= Overweight)` = Normal + Overweight
  ) %>%
  pivot_longer(
    cols = c(`P(Y <= Normal)`, `P(Y <= Overweight)`),
    names_to = "Batas_Kumulatif",
    values_to = "Prob_Kumulatif"
  ) %>%
  mutate(
    Prob_Kumulatif = pmin(pmax(Prob_Kumulatif, eps), 1 - eps),
    Cumulative_Logit = qlogis(Prob_Kumulatif)
  )

ggplot(cumlogit_age, aes(x = Age, y = Cumulative_Logit, color = Batas_Kumulatif)) +
  geom_line(linewidth = 1.25) +
  scale_color_manual(values = c(
    "P(Y <= Normal)" = "#16a34a",
    "P(Y <= Overweight)" = "#f59e0b"
  )) +
  labs(
    title = "Visualisasi Cumulative Logit",
    subtitle = "Pada model proportional odds, garis cumulative logit memiliki slope yang sama.",
    x = "Age",
    y = "Cumulative Logit",
    color = "Batas Kumulatif"
  ) +
  theme_profesional()

3.19 Pembahasan

Interpretasi hasil Regresi Logistik Ordinal

Model digunakan untuk memprediksi status obesitas berurutan, yaitu Normal < Overweight < Obese. Karena respon bersifat bertingkat, regresi logistik ordinal dengan pendekatan cumulative logit/proportional odds sesuai digunakan.

Berdasarkan tabel koefisien: - Age signifikan dengan odds ratio 1,0744, artinya peningkatan umur cenderung meningkatkan peluang berada pada kategori obesitas yang lebih tinggi. - Family_history_with_overweightyes signifikan dengan odds ratio 2,3527, artinya individu dengan riwayat keluarga overweight memiliki odds lebih besar untuk berada pada kategori status obesitas yang lebih tinggi. - FAF signifikan dengan odds ratio 0,6471, artinya peningkatan aktivitas fisik menurunkan odds berada pada kategori obesitas yang lebih tinggi. - CH2O signifikan dengan odds ratio 1,7483, artinya peningkatan konsumsi air pada data ini berkaitan dengan odds lebih tinggi terhadap kategori obesitas yang lebih tinggi. - GenderMale dan FAVCyes tidak signifikan pada taraf 5%, sehingga pengaruhnya belum cukup kuat secara statistik.

Uji Likelihood Ratio Test menghasilkan statistik LR 71,0312 dengan p-value 0, sehingga model dengan prediktor lebih baik dibanding model tanpa prediktor. Evaluasi klasifikasi menghasilkan akurasi 64,71%. Confusion matrix menunjukkan model paling baik mengenali kategori Normal, sementara kategori Obese masih sulit diprediksi dengan tepat pada data uji.

Model regresi logistik ordinal digunakan karena status obesitas memiliki urutan alami, yaitu Normal, Overweight, dan Obese. Sebelum analisis dilakukan, tujuh kategori asli pada variabel NObeyesdad dikelompokkan menjadi tiga kategori agar sesuai dengan struktur ordinal yang digunakan dalam penelitian.

Pengelompokan dilakukan sebagai berikut:

• Insufficient_Weight dan Normal_Weight menjadi Normal.
• Overweight_Level_I dan Overweight_Level_II menjadi Overweight.
• Obesity_Type_I, Obesity_Type_II, dan Obesity_Type_III menjadi Obese.

Setelah dilakukan pemodelan, interpretasi utama dilihat dari nilai koefisien, p-value, dan odds ratio. Variabel dengan p-value kurang dari 0,05 dapat dikatakan berpengaruh signifikan terhadap status obesitas. Odds ratio lebih dari 1 menunjukkan peningkatan peluang berada pada kategori obesitas yang lebih tinggi, sedangkan odds ratio kurang dari 1 menunjukkan penurunan peluang berada pada kategori obesitas yang lebih tinggi.

3.20 Kesimpulan

Berdasarkan analisis regresi logistik ordinal, status obesitas dapat dimodelkan menggunakan variabel Age, Gender, Family history with overweight, FAF, FAVC, dan CH2O. Variabel respon yang semula memiliki tujuh kategori dikelompokkan menjadi tiga kategori ordinal, yaitu Normal, Overweight, dan Obese.

Model ini dapat digunakan untuk menjelaskan faktor-faktor yang berhubungan dengan kecenderungan seseorang berada pada status obesitas yang lebih tinggi. Selain itu, model juga dapat digunakan untuk memprediksi probabilitas seseorang termasuk ke dalam kategori Normal, Overweight, atau Obese berdasarkan karakteristik yang diamati.

4. Regresi Logistik Poisson

4.1 Pendahuluan

Analisis ini bertujuan untuk memodelkan jumlah ikan yang ditangkap menggunakan regresi Poisson. Regresi Poisson sesuai digunakan karena variabel respon berupa data hitungan, yaitu 0, 1, 2, 3, dan seterusnya.

Variabel respon yang digunakan adalah:

\[ Y_i = count_i \]

dengan \(count_i\) menyatakan jumlah ikan yang ditangkap oleh kelompok pemancing ke-\(i\).

Variabel independen yang digunakan adalah:

• camper: status camping, 0 = tidak camping dan 1 = camping.
• persons: jumlah orang dalam kelompok memancing.
• child: jumlah anak dalam kelompok memancing.

4.2 Membaca Data Excel

data <- read_excel("E:/KULIAH/Semester 4/ADK/Tugas regresi/data mancing poisson.xlsx")

glimpse(data)

## Rows: 250
## Columns: 8
## $ nofish   <dbl> 1, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 0, 0…
## $ livebait <dbl> 0, 1, 1, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 1, 1…
## $ camper   <dbl> 0, 1, 0, 1, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 1, 1, 0, 1, 1, 0, 0…
## $ persons  <dbl> 1, 1, 1, 2, 1, 4, 3, 4, 3, 1, 4, 3, 3, 3, 1, 1, 4, 3, 2, 3, 4…
## $ child    <dbl> 0, 0, 0, 1, 0, 2, 1, 3, 2, 0, 1, 2, 0, 0, 0, 0, 1, 2, 1, 1, 1…
## $ xb       <chr> "-0.8963145613670349", "-0.5583449602127075", "-0.40173101425…
## $ zg       <chr> "3.05040478706359E+16", "1746148943901060", "0.27993887662887…
## $ count    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 0, 1, 0, 1…

Dataset yang digunakan berisi informasi aktivitas memancing. Variabel count digunakan sebagai variabel respon karena menunjukkan jumlah ikan yang ditangkap. Variabel tersebut berbentuk hitungan non-negatif sehingga sesuai dianalisis dengan regresi Poisson.

4.3 Memilih Variabel Penelitian

data_model <- data %>%
  dplyr::select(
    count,
    camper,
    persons,
    child
  ) %>%
  dplyr::mutate(
    camper = factor(
      camper,
      levels = c(0, 1),
      labels = c("Tidak camping", "Camping")
    )
  ) %>%
  tidyr::drop_na()

glimpse(data_model)

## Rows: 250
## Columns: 4
## $ count   <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 2, 0, 1, 0, 0, 1, 0, 1,…
## $ camper  <fct> Tidak camping, Camping, Tidak camping, Camping, Tidak camping,…
## $ persons <dbl> 1, 1, 1, 2, 1, 4, 3, 4, 3, 1, 4, 3, 3, 3, 1, 1, 4, 3, 2, 3, 4,…
## $ child   <dbl> 0, 0, 0, 1, 0, 2, 1, 3, 2, 0, 1, 2, 0, 0, 0, 0, 1, 2, 1, 1, 1,…

4.4 Statistik Deskriptif

data_model %>%
  dplyr::summarise(
    Jumlah_Data = n(),
    Minimum_Count = min(count),
    Maksimum_Count = max(count),
    Rata_Rata_Count = mean(count),
    Varians_Count = var(count),
    Median_Count = median(count),
    Rata_Rata_Persons = mean(persons),
    Rata_Rata_Child = mean(child)
  ) %>%
  kable(
    digits = 3,
    caption = "Statistik Deskriptif Variabel Penelitian"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Statistik Deskriptif Variabel Penelitian

Jumlah_Data	Minimum_Count	Maksimum_Count	Rata_Rata_Count	Varians_Count	Median_Count	Rata_Rata_Persons	Rata_Rata_Child
250	0	149	3.296	135.374	0	2.528	0.684

Pada regresi Poisson, salah satu karakteristik penting adalah nilai rata-rata dan varians dari variabel respon. Secara teoritis, distribusi Poisson memiliki sifat bahwa rata-rata sama dengan varians. Jika varians jauh lebih besar daripada rata-rata, maka terdapat indikasi overdispersion.

4.5 Distribusi Jumlah Ikan yang Ditangkap

data_model %>%
  dplyr::count(count) %>%
  dplyr::mutate(Proporsi = n / sum(n)) %>%
  kable(
    digits = 3,
    caption = "Distribusi Jumlah Ikan yang Ditangkap"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Distribusi Jumlah Ikan yang Ditangkap

count	n	Proporsi
0	142	0.568
1	31	0.124
2	20	0.080
3	12	0.048
4	6	0.024
5	10	0.040
6	4	0.016
7	3	0.012
8	2	0.008
9	2	0.008
10	1	0.004
11	1	0.004
13	1	0.004
14	1	0.004
15	2	0.008
16	1	0.004
21	2	0.008
22	1	0.004
29	1	0.004
30	1	0.004
31	1	0.004
32	2	0.008
38	1	0.004
65	1	0.004
149	1	0.004

data_model %>%
  dplyr::count(count) %>%
  ggplot(aes(x = count, y = n)) +
  geom_col(width = 0.8, fill = "#2f7f73", alpha = 0.92) +
  scale_x_continuous(breaks = scales::pretty_breaks()) +
  labs(
    title = "Distribusi Jumlah Ikan yang Ditangkap",
    subtitle = "Variabel respon berupa data hitungan non-negatif.",
    x = "Jumlah ikan yang ditangkap",
    y = "Frekuensi"
  ) +
  theme_profesional()

4.6 Distribusi Berdasarkan Status Camping

data_model %>%
  dplyr::group_by(camper) %>%
  dplyr::summarise(
    Jumlah_Observasi = n(),
    Rata_Rata_Count = mean(count),
    Varians_Count = var(count),
    Median_Count = median(count),
    .groups = "drop"
  ) %>%
  kable(
    digits = 3,
    caption = "Ringkasan Jumlah Ikan Berdasarkan Status Camping"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Ringkasan Jumlah Ikan Berdasarkan Status Camping

camper	Jumlah_Observasi	Rata_Rata_Count	Varians_Count	Median_Count
Tidak camping	103	1.524	21.056	0
Camping	147	4.537	212.401	1

ggplot(data_model, aes(x = camper, y = count, fill = camper)) +
  geom_boxplot(alpha = 0.85, width = 0.55) +
  scale_fill_manual(values = c(
    "Tidak camping" = "#5568b8",
    "Camping" = "#d18b2f"
  )) +
  labs(
    title = "Jumlah Ikan Berdasarkan Status Camping",
    subtitle = "Perbandingan jumlah ikan yang ditangkap antara kelompok camping dan tidak camping.",
    x = "Status Camping",
    y = "Jumlah ikan",
    fill = "Camper"
  ) +
  theme_profesional()

4.7 Rumusan Model Regresi Poisson

Pada regresi Poisson, variabel respon diasumsikan mengikuti distribusi Poisson:

\[ Y_i \sim Poisson(\mu_i) \]

dengan:

\[ E(Y_i \mid X_i) = \mu_i \]

dan:

\[ Var(Y_i \mid X_i) = \mu_i \]

Karena \(\mu_i\) harus bernilai positif, digunakan fungsi hubung log:

\[ \log(\mu_i) = \beta_0 + \beta_1 camper_i + \beta_2 persons_i + \beta_3 child_i \]

Sehingga rata-rata jumlah ikan yang ditangkap dapat dituliskan sebagai:

\[ \mu_i = \exp( \beta_0 + \beta_1 camper_i + \beta_2 persons_i + \beta_3 child_i ) \]

Model ini menjelaskan pengaruh status camping, jumlah orang, dan jumlah anak terhadap rata-rata jumlah ikan yang ditangkap. Karena model menggunakan link log, maka koefisien lebih mudah diinterpretasikan melalui nilai Incidence Rate Ratio atau IRR, yaitu \(e^\beta\).

4.8 Estimasi Model Regresi Poisson

model_poisson <- glm(
  count ~ camper + persons + child,
  data = data_model,
  family = poisson(link = "log")
)

summary(model_poisson)

## 
## Call:
## glm(formula = count ~ camper + persons + child, family = poisson(link = "log"), 
##     data = data_model)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)   -1.98183    0.15226  -13.02   <2e-16 ***
## camperCamping  0.93094    0.08909   10.45   <2e-16 ***
## persons        1.09126    0.03926   27.80   <2e-16 ***
## child         -1.68996    0.08099  -20.87   <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: 2958.4  on 249  degrees of freedom
## Residual deviance: 1337.1  on 246  degrees of freedom
## AIC: 1682.1
## 
## Number of Fisher Scoring iterations: 6

4.9 Ringkasan Koefisien dan Incidence Rate Ratio

tabel_koef <- as.data.frame(coef(summary(model_poisson))) %>%
  tibble::rownames_to_column("Parameter") %>%
  dplyr::rename(
    Estimate = Estimate,
    Std_Error = `Std. Error`,
    Z_Value = `z value`,
    P_Value = `Pr(>|z|)`
  ) %>%
  dplyr::mutate(
    IRR = exp(Estimate),
    CI_Lower = exp(Estimate - 1.96 * Std_Error),
    CI_Upper = exp(Estimate + 1.96 * Std_Error),
    Perubahan_Persen = 100 * (IRR - 1),
    Kesimpulan = ifelse(P_Value < 0.05, "Signifikan", "Tidak signifikan")
  )

tabel_koef %>%
  dplyr::mutate(
    dplyr::across(
      c(Estimate, Std_Error, Z_Value, P_Value, IRR, CI_Lower, CI_Upper, Perubahan_Persen),
      ~ round(.x, 4)
    )
  ) %>%
  kable(
    caption = "Ringkasan Hasil Regresi Poisson",
    col.names = c(
      "Parameter", "Estimate", "Std. Error", "Z Value", "P Value",
      "IRR", "CI 2.5%", "CI 97.5%", "Perubahan (%)", "Kesimpulan"
    )
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  ) %>%
  row_spec(
    which(tabel_koef$P_Value < 0.05),
    bold = TRUE,
    color = "white",
    background = "#16a34a"
  )

Ringkasan Hasil Regresi Poisson

Parameter	Estimate	Std. Error	Z Value	P Value	IRR	CI 2.5%	CI 97.5%	Perubahan (%)	Kesimpulan
(Intercept)	-1.9818	0.1523	-13.0158	0	0.1378	0.1023	0.1857	-86.2183	Signifikan
camperCamping	0.9309	0.0891	10.4498	0	2.5369	2.1304	3.0209	153.6882	Signifikan
persons	1.0913	0.0393	27.7992	0	2.9780	2.7575	3.2162	197.8031	Signifikan
child	-1.6900	0.0810	-20.8658	0	0.1845	0.1574	0.2163	-81.5473	Signifikan

Nilai IRR digunakan untuk menafsirkan pengaruh variabel independen. Jika IRR lebih dari 1, maka variabel tersebut meningkatkan rata-rata jumlah ikan yang ditangkap. Jika IRR kurang dari 1, maka variabel tersebut menurunkan rata-rata jumlah ikan yang ditangkap.

4.10 Interpretasi Otomatis Koefisien

tabel_interpretasi <- tabel_koef %>%
  dplyr::filter(Parameter != "(Intercept)") %>%
  dplyr::mutate(
    Interpretasi = dplyr::case_when(
      IRR > 1 ~ paste0(
        "Parameter ", Parameter,
        " memiliki IRR sebesar ", round(IRR, 4),
        ". Artinya, variabel ini meningkatkan rata-rata jumlah ikan yang ditangkap sebesar ",
        round(Perubahan_Persen, 2),
        "%, dengan asumsi variabel lain konstan."
      ),
      IRR < 1 ~ paste0(
        "Parameter ", Parameter,
        " memiliki IRR sebesar ", round(IRR, 4),
        ". Artinya, variabel ini menurunkan rata-rata jumlah ikan yang ditangkap sebesar ",
        abs(round(Perubahan_Persen, 2)),
        "%, dengan asumsi variabel lain konstan."
      ),
      TRUE ~ paste0(
        "Parameter ", Parameter,
        " memiliki IRR sebesar 1, sehingga tidak mengubah rata-rata jumlah ikan."
      )
    )
  ) %>%
  dplyr::select(Parameter, IRR, P_Value, Kesimpulan, Interpretasi)

tabel_interpretasi %>%
  dplyr::mutate(
    IRR = round(IRR, 4),
    P_Value = round(P_Value, 4)
  ) %>%
  kable(
    caption = "Interpretasi Koefisien Regresi Poisson"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = TRUE
  )

Interpretasi Koefisien Regresi Poisson

Parameter	IRR	P_Value	Kesimpulan	Interpretasi
camperCamping	2.5369	0	Signifikan	Parameter camperCamping memiliki IRR sebesar 2.5369. Artinya, variabel ini meningkatkan rata-rata jumlah ikan yang ditangkap sebesar 153.69%, dengan asumsi variabel lain konstan.
persons	2.9780	0	Signifikan	Parameter persons memiliki IRR sebesar 2.978. Artinya, variabel ini meningkatkan rata-rata jumlah ikan yang ditangkap sebesar 197.8%, dengan asumsi variabel lain konstan.
child	0.1845	0	Signifikan	Parameter child memiliki IRR sebesar 0.1845. Artinya, variabel ini menurunkan rata-rata jumlah ikan yang ditangkap sebesar 81.55%, dengan asumsi variabel lain konstan.

4.11 Uji Kelayakan Model

4.11.1 Likelihood Ratio Test

model_null <- glm(
  count ~ 1,
  data = data_model,
  family = poisson(link = "log")
)

uji_lr <- anova(model_null, model_poisson, test = "Chisq")

uji_lr

Jika p-value pada Likelihood Ratio Test kurang dari 0,05, maka model dengan prediktor camper, persons, dan child lebih baik dibandingkan model tanpa prediktor.

4.11.2 AIC Model

tibble(
  Model = c("Null Model", "Poisson Model"),
  AIC = c(AIC(model_null), AIC(model_poisson))
) %>%
  kable(
    digits = 3,
    caption = "Perbandingan AIC Model"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Perbandingan AIC Model

Model	AIC
Null Model	3297.431
Poisson Model	1682.145

Model dengan nilai AIC yang lebih kecil dianggap lebih baik dibandingkan model dengan AIC yang lebih besar.

4.12 Pemeriksaan Overdispersion

Salah satu asumsi penting regresi Poisson adalah equidispersion, yaitu:

\[ E(Y_i \mid X_i) = Var(Y_i \mid X_i) \]

Indikasi overdispersion dapat diperiksa dengan rasio Pearson:

\[ \hat{\phi} = \frac{\sum r_{P,i}^{2}}{df} \]

dengan \(r_{P,i}\) adalah residual Pearson dan \(df\) adalah derajat bebas residual.

dispersion <- sum(residuals(model_poisson, type = "pearson")^2) /
  df.residual(model_poisson)

tibble(
  Dispersion_Pearson = dispersion,
  Interpretasi = dplyr::case_when(
    dispersion < 1.5 ~ "Tidak ada indikasi overdispersion berat",
    dispersion < 2.5 ~ "Ada indikasi overdispersion sedang",
    TRUE ~ "Ada indikasi overdispersion kuat"
  )
) %>%
  dplyr::mutate(Dispersion_Pearson = round(Dispersion_Pearson, 3)) %>%
  kable(
    caption = "Pemeriksaan Overdispersion"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Pemeriksaan Overdispersion

Dispersion_Pearson	Interpretasi
11.832	Ada indikasi overdispersion kuat

Jika nilai dispersion Pearson mendekati 1, maka asumsi equidispersion cukup terpenuhi. Jika nilainya jauh lebih besar dari 1, maka terdapat indikasi overdispersion. Jika overdispersion kuat terjadi, alternatif yang dapat digunakan adalah quasi-Poisson atau regresi binomial negatif.

4.13 Prediksi Jumlah Ikan yang Ditangkap

data_prediksi <- data_model %>%
  dplyr::mutate(
    Prediksi_Count = predict(model_poisson, type = "response"),
    Residual = count - Prediksi_Count
  )

head(data_prediksi)

data_prediksi %>%
  dplyr::select(count, Prediksi_Count, Residual, camper, persons, child) %>%
  head(10) %>%
  kable(
    digits = 3,
    caption = "Contoh Hasil Prediksi Model Poisson"
  ) %>%
  kable_styling(
    bootstrap_options = c("striped", "hover", "condensed", "responsive"),
    full_width = FALSE
  )

Contoh Hasil Prediksi Model Poisson

count	Prediksi_Count	Residual	camper	persons	child
0	0.410	-0.410	Tidak camping	1	0
0	1.041	-1.041	Camping	1	0
0	0.410	-0.410	Tidak camping	1	0
0	0.572	-0.572	Camping	2	1
1	0.410	0.590	Tidak camping	1	0
0	0.936	-0.936	Camping	4	2
0	0.672	-0.672	Tidak camping	3	1
0	0.068	-0.068	Tidak camping	4	3
0	0.314	-0.314	Camping	3	2
1	1.041	-0.041	Camping	1	0

4.14 Contoh Prediksi Kelompok Pemancing Baru

Misalkan terdapat kelompok pemancing dengan karakteristik:

• Status camping = Camping
• Jumlah orang = 3
• Jumlah anak = 1

kelompok_baru <- tibble(
  camper = factor("Camping", levels = levels(data_model$camper)),
  persons = 3,
  child = 1
)

prediksi_baru <- predict(
  model_poisson,
  newdata = kelompok_baru,
  type = "response"
)

prediksi_baru

##        1 
## 1.703933

Nilai prediksi tersebut merupakan rata-rata jumlah ikan yang diperkirakan ditangkap oleh kelompok pemancing dengan karakteristik yang diberikan.

4.15 Visualisasi Pengaruh Jumlah Orang terhadap Prediksi Count

grid_persons <- tibble(
  persons = seq(
    min(data_model$persons, na.rm = TRUE),
    max(data_model$persons, na.rm = TRUE),
    length.out = 100
  ),
  camper = factor("Camping", levels = levels(data_model$camper)),
  child = round(mean(data_model$child, na.rm = TRUE))
)

pred_persons <- predict(
  model_poisson,
  newdata = grid_persons,
  type = "link",
  se.fit = TRUE
)

plot_persons <- grid_persons %>%
  dplyr::mutate(
    fit_link = pred_persons$fit,
    se_link = pred_persons$se.fit,
    Prediksi_Count = exp(fit_link),
    Lower = exp(fit_link - 1.96 * se_link),
    Upper = exp(fit_link + 1.96 * se_link)
  )

ggplot(plot_persons, aes(x = persons, y = Prediksi_Count)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "#2f7f73", alpha = 0.18) +
  geom_line(color = "#2f7f73", linewidth = 1.35) +
  labs(
    title = "Prediksi Jumlah Ikan Berdasarkan Jumlah Orang",
    subtitle = "Status camping = Camping; jumlah anak dibuat tetap pada rata-rata.",
    x = "Jumlah orang dalam kelompok",
    y = "Prediksi rata-rata jumlah ikan"
  ) +
  theme_profesional()

4.16 Visualisasi Pengaruh Jumlah Anak terhadap Prediksi Count

grid_child <- tibble(
  child = seq(
    min(data_model$child, na.rm = TRUE),
    max(data_model$child, na.rm = TRUE),
    length.out = 100
  ),
  camper = factor("Camping", levels = levels(data_model$camper)),
  persons = round(mean(data_model$persons, na.rm = TRUE))
)

pred_child <- predict(
  model_poisson,
  newdata = grid_child,
  type = "link",
  se.fit = TRUE
)

plot_child <- grid_child %>%
  dplyr::mutate(
    fit_link = pred_child$fit,
    se_link = pred_child$se.fit,
    Prediksi_Count = exp(fit_link),
    Lower = exp(fit_link - 1.96 * se_link),
    Upper = exp(fit_link + 1.96 * se_link)
  )

ggplot(plot_child, aes(x = child, y = Prediksi_Count)) +
  geom_ribbon(aes(ymin = Lower, ymax = Upper), fill = "#d18b2f", alpha = 0.18) +
  geom_line(color = "#d18b2f", linewidth = 1.35) +
  labs(
    title = "Prediksi Jumlah Ikan Berdasarkan Jumlah Anak",
    subtitle = "Status camping = Camping; jumlah orang dibuat tetap pada rata-rata.",
    x = "Jumlah anak dalam kelompok",
    y = "Prediksi rata-rata jumlah ikan"
  ) +
  theme_profesional()

4.17 Visualisasi Perbandingan Camper dan Non-Camper

grid_camper <- expand.grid(
  camper = levels(data_model$camper),
  persons = seq(
    min(data_model$persons, na.rm = TRUE),
    max(data_model$persons, na.rm = TRUE),
    length.out = 100
  ),
  child = round(mean(data_model$child, na.rm = TRUE))
)

pred_camper <- predict(
  model_poisson,
  newdata = grid_camper,
  type = "response"
)

plot_camper <- grid_camper %>%
  dplyr::mutate(Prediksi_Count = pred_camper)

ggplot(plot_camper, aes(x = persons, y = Prediksi_Count, color = camper)) +
  geom_line(linewidth = 1.35) +
  scale_color_manual(values = c(
    "Tidak camping" = "#5568b8",
    "Camping" = "#d18b2f"
  )) +
  labs(
    title = "Prediksi Jumlah Ikan Berdasarkan Status Camping",
    subtitle = "Grafik menunjukkan perbedaan prediksi count antara kelompok camping dan tidak camping.",
    x = "Jumlah orang dalam kelompok",
    y = "Prediksi rata-rata jumlah ikan",
    color = "Status Camping"
  ) +
  theme_profesional()

4.18 Perbandingan Nilai Aktual dan Prediksi

ggplot(data_prediksi, aes(x = Prediksi_Count, y = count)) +
  geom_point(alpha = 0.65, color = "#5568b8") +
  geom_abline(intercept = 0, slope = 1, linetype = "dashed", color = "#dc2626", linewidth = 1) +
  labs(
    title = "Perbandingan Jumlah Aktual dan Prediksi",
    subtitle = "Garis putus-putus menunjukkan kondisi prediksi sama dengan aktual.",
    x = "Prediksi jumlah ikan",
    y = "Jumlah ikan aktual"
  ) +
  theme_profesional()

4.19 Pembahasan

Interpretasi hasil Regresi Logistik Poisson

Model Poisson digunakan karena variabel respon berupa data cacahan, yaitu jumlah ikan yang ditangkap. Dalam model Poisson, interpretasi koefisien dilakukan menggunakan Incidence Rate Ratio (IRR).

Berdasarkan hasil estimasi: - camperCamping signifikan dengan IRR 2,5369, artinya kelompok yang camping memiliki rata-rata jumlah ikan yang ditangkap sekitar 153,69% lebih tinggi dibanding kelompok referensi. - persons signifikan dengan IRR 2,9780, artinya setiap tambahan 1 orang meningkatkan rata-rata jumlah ikan yang ditangkap sekitar 197,80%, dengan asumsi variabel lain tetap. - child signifikan dengan IRR 0,1845, artinya setiap tambahan 1 anak menurunkan rata-rata jumlah ikan yang ditangkap sekitar 81,55%, dengan asumsi variabel lain tetap.

Uji Likelihood Ratio Test menghasilkan deviance 1621,3 dengan p-value < 2,2e-16, sehingga model dengan prediktor lebih baik dibanding model tanpa prediktor. Secara umum, variabel camper, persons, dan child berpengaruh signifikan terhadap jumlah ikan yang ditangkap.

Regresi Poisson digunakan karena variabel respon berbentuk data hitungan, yaitu jumlah ikan yang ditangkap. Model menggunakan fungsi hubung log sehingga setiap koefisien menjelaskan perubahan pada log rata-rata jumlah ikan. Agar interpretasi lebih mudah, koefisien dieksponensialkan menjadi IRR.

Interpretasi umum:

• Jika koefisien camperCamping positif dan IRR lebih dari 1, maka kelompok yang camping memiliki rata-rata jumlah ikan lebih tinggi dibandingkan kelompok yang tidak camping.
• Jika koefisien persons positif dan IRR lebih dari 1, maka semakin banyak orang dalam kelompok, semakin besar rata-rata jumlah ikan yang ditangkap.
• Jika koefisien child negatif dan IRR kurang dari 1, maka semakin banyak anak dalam kelompok dapat menurunkan rata-rata jumlah ikan yang ditangkap.
• Variabel dengan p-value kurang dari 0,05 dapat dikatakan berpengaruh signifikan terhadap jumlah ikan yang ditangkap.

4.20 Kesimpulan

Berdasarkan analisis regresi Poisson, jumlah ikan yang ditangkap dapat dimodelkan menggunakan variabel status camping, jumlah orang dalam kelompok, dan jumlah anak dalam kelompok. Model ini sesuai karena variabel respon berupa data hitungan non-negatif.

Hasil analisis dapat digunakan untuk melihat faktor-faktor yang berhubungan dengan peningkatan atau penurunan rata-rata jumlah ikan yang ditangkap. Interpretasi utama dilakukan melalui nilai koefisien, p-value, dan Incidence Rate Ratio atau IRR.