logit--probit

# ============================================================
#   LOGIT DAN PROBIT — TEORI & VISUALISASI LENGKAP
#   Siap salin dan jalankan di R / RStudio
# ============================================================
# Paket yang dibutuhkan:
#   install.packages(c("ggplot2", "dplyr", "gridExtra",
#                      "margins", "pROC", "tidyr"))
# ============================================================

library(ggplot2)

Warning: package 'ggplot2' was built under R version 4.5.3

library(dplyr)

Warning: package 'dplyr' was built under R version 4.5.3


Attaching package: 'dplyr'

The following objects are masked from 'package:stats':

    filter, lag

The following objects are masked from 'package:base':

    intersect, setdiff, setequal, union

library(gridExtra)

Warning: package 'gridExtra' was built under R version 4.5.3


Attaching package: 'gridExtra'

The following object is masked from 'package:dplyr':

    combine

library(margins)

Warning: package 'margins' was built under R version 4.5.3

library(pROC)

Warning: package 'pROC' was built under R version 4.5.3

Type 'citation("pROC")' for a citation.


Attaching package: 'pROC'

The following objects are masked from 'package:stats':

    cov, smooth, var

library(tidyr)

Warning: package 'tidyr' was built under R version 4.5.3

# ─────────────────────────────────────────────
# 0. TEMA GRAFIK SERAGAM
# ─────────────────────────────────────────────
tema <- theme_minimal(base_size = 13) +
  theme(
    plot.title    = element_text(face = "bold", size = 14),
    plot.subtitle = element_text(color = "grey50", size = 11),
    axis.title    = element_text(size = 11),
    legend.position = "bottom",
    panel.grid.minor = element_blank()
  )


# ============================================================
# BAGIAN 1: FUNGSI LINK — LOGIT vs PROBIT
# ============================================================

z <- seq(-5, 5, length.out = 300)

df_link <- data.frame(
  z      = rep(z, 2),
  P      = c(plogis(z),         # fungsi logistik
             pnorm(z)),          # CDF normal standar
  Model  = rep(c("Logit (Logistik)", "Probit (Normal CDF)"), each = 300)
)

p1 <- ggplot(df_link, aes(x = z, y = P, color = Model, linetype = Model)) +
  geom_line(linewidth = 1.2) +
  geom_hline(yintercept = c(0, 1), linetype = "dotted", color = "grey70") +
  geom_vline(xintercept = 0, linetype = "dotted", color = "grey70") +
  scale_color_manual(values = c("Logit (Logistik)" = "#5046A8",
                                "Probit (Normal CDF)" = "#0D6B50")) +
  scale_linetype_manual(values = c("Logit (Logistik)" = "solid",
                                   "Probit (Normal CDF)" = "dashed")) +
  labs(
    title    = "Fungsi Link: Logit vs Probit",
    subtitle = "Kedua kurva sangat mirip di tengah; Logit memiliki ekor yang sedikit lebih tebal",
    x        = "Indeks linear  z = β₀ + β₁X₁ + …",
    y        = "P(Y = 1)",
    color    = NULL, linetype = NULL
  ) +
  tema

print(p1)

# ============================================================
# BAGIAN 2: PERBANDINGAN EKOR (TAIL BEHAVIOR)
# ============================================================

df_tail <- data.frame(
  z     = rep(z, 2),
  pdf   = c(dlogis(z),   # PDF distribusi logistik
            dnorm(z)),   # PDF normal standar
  Model = rep(c("Logistik (Logit)", "Normal (Probit)"), each = 300)
)

p2 <- ggplot(df_tail, aes(x = z, y = pdf, color = Model, linetype = Model)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c("Logistik (Logit)" = "#5046A8",
                                "Normal (Probit)"  = "#0D6B50")) +
  scale_linetype_manual(values = c("Logistik (Logit)" = "solid",
                                   "Normal (Probit)"  = "dashed")) +
  labs(
    title    = "PDF Distribusi Error Asumsi",
    subtitle = "Distribusi logistik memiliki varians π²/3 ≈ 3.29; normal standar = 1",
    x = "z", y = "Densitas",
    color = NULL, linetype = NULL
  ) +
  tema

print(p2)

# ============================================================
# BAGIAN 3: SIMULASI DATA & ESTIMASI MODEL
# ============================================================

set.seed(42)
n  <- 500
X1 <- rnorm(n, mean = 0, sd = 1)       # variabel kontinu
X2 <- rbinom(n, 1, 0.5)                # variabel dummy (0/1)

# Probabilitas sejati menggunakan model logistik
z_true <- -0.5 + 1.2 * X1 + 0.8 * X2
p_true <- plogis(z_true)

Y  <- rbinom(n, 1, p_true)             # variabel dependen biner

df <- data.frame(Y = Y, X1 = X1, X2 = factor(X2, labels = c("Grup A", "Grup B")))

cat("\n===== RINGKASAN DATA =====\n")


===== RINGKASAN DATA =====

cat("N obs:", n, "\n")

N obs: 500

cat("Y = 1:", sum(Y), " | Y = 0:", sum(Y == 0), "\n")

Y = 1: 255  | Y = 0: 245

# ── Estimasi Model ──────────────────────────────────────────
model_logit  <- glm(Y ~ X1 + X2, data = df, family = binomial(link = "logit"))
model_probit <- glm(Y ~ X1 + X2, data = df, family = binomial(link = "probit"))

cat("\n===== RINGKASAN MODEL LOGIT =====\n")


===== RINGKASAN MODEL LOGIT =====

print(summary(model_logit))


Call:
glm(formula = Y ~ X1 + X2, family = binomial(link = "logit"), 
    data = df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -0.1249     0.1394  -0.896   0.3701    
X1            1.0865     0.1255   8.659   <2e-16 ***
X2Grup B      0.4243     0.1997   2.125   0.0336 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 692.95  on 499  degrees of freedom
Residual deviance: 589.20  on 497  degrees of freedom
AIC: 595.2

Number of Fisher Scoring iterations: 4

cat("\n===== RINGKASAN MODEL PROBIT =====\n")


===== RINGKASAN MODEL PROBIT =====

print(summary(model_probit))


Call:
glm(formula = Y ~ X1 + X2, family = binomial(link = "probit"), 
    data = df)

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept) -0.07728    0.08407  -0.919   0.3580    
X1           0.66083    0.07132   9.266   <2e-16 ***
X2Grup B     0.26227    0.12040   2.178   0.0294 *  
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 692.95  on 499  degrees of freedom
Residual deviance: 588.54  on 497  degrees of freedom
AIC: 594.54

Number of Fisher Scoring iterations: 4

# ── Odds Ratio (khusus Logit) ────────────────────────────────
cat("\n===== ODDS RATIO (LOGIT) =====\n")


===== ODDS RATIO (LOGIT) =====

OR <- exp(cbind(OR = coef(model_logit), confint(model_logit)))

Waiting for profiling to be done...

print(round(OR, 4))

                OR  2.5 % 97.5 %
(Intercept) 0.8825 0.6707 1.1595
X1          2.9638 2.3367 3.8242
X2Grup B    1.5286 1.0351 2.2666

# ── Perbandingan Koefisien ───────────────────────────────────
cat("\n===== PERBANDINGAN KOEFISIEN =====\n")


===== PERBANDINGAN KOEFISIEN =====

koef <- data.frame(
  Parameter  = names(coef(model_logit)),
  Logit      = round(coef(model_logit), 4),
  Probit     = round(coef(model_probit), 4),
  Rasio_L_P  = round(coef(model_logit) / coef(model_probit), 3)
)
print(koef)

              Parameter   Logit  Probit Rasio_L_P
(Intercept) (Intercept) -0.1249 -0.0773     1.617
X1                   X1  1.0865  0.6608     1.644
X2Grup B       X2Grup B  0.4243  0.2623     1.618

cat("\nCatatan: β_logit ≈ π/√3 ≈ 1.814 × β_probit (konversi teoritis)\n")


Catatan: β_logit ≈ π/√3 ≈ 1.814 × β_probit (konversi teoritis)

# ============================================================
# BAGIAN 4: GRAFIK KOEFISIEN (COEFFICIENT PLOT)
# ============================================================

df_koef <- data.frame(
  term  = rep(names(coef(model_logit))[-1], 2),   # hapus intercept
  est   = c(coef(model_logit)[-1], coef(model_probit)[-1]),
  lower = c(confint(model_logit)[-1, 1], confint(model_probit)[-1, 1]),
  upper = c(confint(model_logit)[-1, 2], confint(model_probit)[-1, 2]),
  Model = rep(c("Logit", "Probit"), each = 2)
)

Waiting for profiling to be done...

Waiting for profiling to be done...
Waiting for profiling to be done...
Waiting for profiling to be done...

p3 <- ggplot(df_koef, aes(x = est, y = term, color = Model, xmin = lower, xmax = upper)) +
  geom_vline(xintercept = 0, linetype = "dashed", color = "grey60") +
  geom_errorbarh(height = 0.2, linewidth = 0.8, position = position_dodge(0.5)) +
  geom_point(size = 3.5, position = position_dodge(0.5)) +
  scale_color_manual(values = c(Logit = "#5046A8", Probit = "#0D6B50")) +
  labs(
    title    = "Estimasi Koefisien + Interval Kepercayaan 95%",
    subtitle = "Skala berbeda karena distribusi error berbeda (Logit ≈ 1.8 × Probit)",
    x = "Nilai Koefisien", y = NULL, color = NULL
  ) +
  tema

Warning: `geom_errorbarh()` was deprecated in ggplot2 4.0.0.
ℹ Please use the `orientation` argument of `geom_errorbar()` instead.

print(p3)

`height` was translated to `width`.

# ============================================================
# BAGIAN 5: MARGINAL EFFECTS
# ============================================================

# Average Marginal Effects (AME)
ame_logit  <- margins(model_logit)
ame_probit <- margins(model_probit)

cat("\n===== AVERAGE MARGINAL EFFECTS — LOGIT =====\n")


===== AVERAGE MARGINAL EFFECTS — LOGIT =====

print(summary(ame_logit))

   factor    AME     SE       z      p  lower  upper
       X1 0.2199 0.0174 12.6080 0.0000 0.1857 0.2540
 X2Grup B 0.0861 0.0402  2.1432 0.0321 0.0074 0.1649

cat("\n===== AVERAGE MARGINAL EFFECTS — PROBIT =====\n")


===== AVERAGE MARGINAL EFFECTS — PROBIT =====

print(summary(ame_probit))

   factor    AME     SE       z      p  lower  upper
       X1 0.2208 0.0171 12.9478 0.0000 0.1874 0.2542
 X2Grup B 0.0879 0.0401  2.1922 0.0284 0.0093 0.1665

# Visualisasi: Marginal Effect X1 di berbagai nilai X1
x1_seq <- seq(-3, 3, length.out = 200)
x2_fix <- 0   # X2 dikunci pada nilai referensi

z_logit  <- coef(model_logit)["(Intercept)"]  + coef(model_logit)["X1"]  * x1_seq + coef(model_logit)["X2Grup B"]  * x2_fix
z_probit <- coef(model_probit)["(Intercept)"] + coef(model_probit)["X1"] * x1_seq + coef(model_probit)["X2Grup B"] * x2_fix

me_logit  <- coef(model_logit)["X1"]  * dlogis(z_logit)
me_probit <- coef(model_probit)["X1"] * dnorm(z_probit)

df_me <- data.frame(
  X1    = rep(x1_seq, 2),
  ME    = c(me_logit, me_probit),
  Model = rep(c("Logit", "Probit"), each = 200)
)

p4 <- ggplot(df_me, aes(x = X1, y = ME, color = Model, linetype = Model)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c(Logit = "#5046A8", Probit = "#0D6B50")) +
  scale_linetype_manual(values = c(Logit = "solid", Probit = "dashed")) +
  labs(
    title    = "Marginal Effect X1 terhadap P(Y=1)",
    subtitle = "Efek tidak konstan — terbesar di tengah (P ≈ 0.5), mengecil di ekor",
    x = "Nilai X1", y = "∂P(Y=1)/∂X1",
    color = NULL, linetype = NULL
  ) +
  tema

print(p4)

# ============================================================
# BAGIAN 6: PREDIKSI & PROBABILITAS FITTED
# ============================================================

df$prob_logit  <- fitted(model_logit)
df$prob_probit <- fitted(model_probit)

# Scatter: prediksi logit vs probit
p5 <- ggplot(df, aes(x = prob_logit, y = prob_probit)) +
  geom_point(alpha = 0.35, color = "#5046A8", size = 1.5) +
  geom_abline(slope = 1, intercept = 0, linetype = "dashed", color = "grey50") +
  labs(
    title    = "Prediksi Probabilitas: Logit vs Probit",
    subtitle = "Titik hampir sejajar dengan garis 45° — prediksi sangat mirip",
    x = "P̂ dari Logit", y = "P̂ dari Probit"
  ) +
  tema

print(p5)

# Distribusi probabilitas fitted, dibedakan Y=0 dan Y=1
df_long <- df %>%
  pivot_longer(cols = c(prob_logit, prob_probit),
               names_to  = "Model",
               values_to = "Probabilitas") %>%
  mutate(Model = recode(Model,
                        prob_logit  = "Logit",
                        prob_probit = "Probit"),
         Y = factor(Y, labels = c("Y = 0", "Y = 1")))

p6 <- ggplot(df_long, aes(x = Probabilitas, fill = Y)) +
  geom_histogram(bins = 30, alpha = 0.7, position = "identity") +
  facet_wrap(~ Model) +
  scale_fill_manual(values = c("Y = 0" = "#C8C3F0", "Y = 1" = "#5046A8")) +
  labs(
    title    = "Distribusi Probabilitas Fitted",
    subtitle = "Kedua model memisahkan Y=0 dan Y=1 secara serupa",
    x = "Probabilitas Prediksi", y = "Frekuensi", fill = NULL
  ) +
  tema

print(p6)

# ============================================================
# BAGIAN 7: KURVA ROC & AUC
# ============================================================

roc_logit  <- roc(df$Y, df$prob_logit,  quiet = TRUE)
roc_probit <- roc(df$Y, df$prob_probit, quiet = TRUE)

cat("\n===== AUC =====\n")


===== AUC =====

cat("Logit  AUC:", round(auc(roc_logit),  4), "\n")

Logit  AUC: 0.7455

cat("Probit AUC:", round(auc(roc_probit), 4), "\n")

Probit AUC: 0.7456

df_roc <- rbind(
  data.frame(
    FPR   = 1 - roc_logit$specificities,
    TPR   = roc_logit$sensitivities,
    Model = paste0("Logit (AUC = ", round(auc(roc_logit), 3), ")")
  ),
  data.frame(
    FPR   = 1 - roc_probit$specificities,
    TPR   = roc_probit$sensitivities,
    Model = paste0("Probit (AUC = ", round(auc(roc_probit), 3), ")")
  )
)

p7 <- ggplot(df_roc, aes(x = FPR, y = TPR, color = Model, linetype = Model)) +
  geom_line(linewidth = 1.2) +
  geom_abline(slope = 1, intercept = 0, linetype = "dotted", color = "grey60") +
  scale_color_manual(
    name = NULL,
    values = setNames(
      c("#5046A8", "#0D6B50"),
      c(paste0("Logit (AUC = ", round(auc(roc_logit), 3), ")"),
        paste0("Probit (AUC = ", round(auc(roc_probit), 3), ")"))
    )
  ) +
  scale_linetype_manual(
    name = NULL,
    values = setNames(
      c("solid", "dashed"),
      c(paste0("Logit (AUC = ", round(auc(roc_logit), 3), ")"),
        paste0("Probit (AUC = ", round(auc(roc_probit), 3), ")"))
    )
  ) +
  labs(
    title = "Kurva ROC",
    subtitle = "Area under the curve (AUC) mengukur kemampuan diskriminasi model",
    x = "False Positive Rate (1 - Spesifisitas)",
    y = "True Positive Rate (Sensitivitas)",
    color = NULL,
    linetype = NULL
  )
  tema

<theme> List of 144
 $ line                            : <ggplot2::element_line>
  ..@ colour       : chr "black"
  ..@ linewidth    : num 0.591
  ..@ linetype     : num 1
  ..@ lineend      : chr "butt"
  ..@ linejoin     : chr "round"
  ..@ arrow        : logi FALSE
  ..@ arrow.fill   : chr "black"
  ..@ inherit.blank: logi TRUE
 $ rect                            : <ggplot2::element_rect>
  ..@ fill         : chr "white"
  ..@ colour       : chr "black"
  ..@ linewidth    : num 0.591
  ..@ linetype     : num 1
  ..@ linejoin     : chr "round"
  ..@ inherit.blank: logi TRUE
 $ text                            : <ggplot2::element_text>
  ..@ family       : chr ""
  ..@ face         : chr "plain"
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : chr "black"
  ..@ size         : num 13
  ..@ hjust        : num 0.5
  ..@ vjust        : num 0.5
  ..@ angle        : num 0
  ..@ lineheight   : num 0.9
  ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 0
  ..@ debug        : logi FALSE
  ..@ inherit.blank: logi TRUE
 $ title                           : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : NULL
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ point                           : <ggplot2::element_point>
  ..@ colour       : chr "black"
  ..@ shape        : num 19
  ..@ size         : num 1.77
  ..@ fill         : chr "white"
  ..@ stroke       : num 0.591
  ..@ inherit.blank: logi TRUE
 $ polygon                         : <ggplot2::element_polygon>
  ..@ fill         : chr "white"
  ..@ colour       : chr "black"
  ..@ linewidth    : num 0.591
  ..@ linetype     : num 1
  ..@ linejoin     : chr "round"
  ..@ inherit.blank: logi TRUE
 $ geom                            : <ggplot2::element_geom>
  ..@ ink        : chr "black"
  ..@ paper      : chr "white"
  ..@ accent     : chr "#3366FF"
  ..@ linewidth  : num 0.591
  ..@ borderwidth: num 0.591
  ..@ linetype   : int 1
  ..@ bordertype : int 1
  ..@ family     : chr ""
  ..@ fontsize   : num 4.57
  ..@ pointsize  : num 1.77
  ..@ pointshape : num 19
  ..@ colour     : NULL
  ..@ fill       : NULL
 $ spacing                         : 'simpleUnit' num 6.5points
  ..- attr(*, "unit")= int 8
 $ margins                         : <ggplot2::margin> num [1:4] 6.5 6.5 6.5 6.5
 $ aspect.ratio                    : NULL
 $ axis.title                      : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : num 11
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : NULL
  ..@ debug        : NULL
  ..@ inherit.blank: logi FALSE
 $ axis.title.x                    : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : num 1
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 3.25 0 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.title.x.top                : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : num 0
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 0 3.25 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.title.x.bottom             : NULL
 $ axis.title.y                    : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : num 1
  ..@ angle        : num 90
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 3.25 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.title.y.left               : NULL
 $ axis.title.y.right              : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : num 1
  ..@ angle        : num -90
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 3.25
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text                       : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : chr "#4D4D4DFF"
  ..@ size         : 'rel' num 0.8
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : NULL
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.x                     : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : num 1
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 2.6 0 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.x.top                 : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 0 5.85 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.x.bottom              : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 5.85 0 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.y                     : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : num 1
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 2.6 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.y.left                : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 5.85 0 0
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.y.right               : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 0 0 5.85
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.text.theta                 : NULL
 $ axis.text.r                     : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : num 0.5
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : <ggplot2::margin> num [1:4] 0 2.6 0 2.6
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ axis.ticks                      : <ggplot2::element_blank>
 $ axis.ticks.x                    : NULL
 $ axis.ticks.x.top                : NULL
 $ axis.ticks.x.bottom             : NULL
 $ axis.ticks.y                    : NULL
 $ axis.ticks.y.left               : NULL
 $ axis.ticks.y.right              : NULL
 $ axis.ticks.theta                : NULL
 $ axis.ticks.r                    : NULL
 $ axis.minor.ticks.x.top          : NULL
 $ axis.minor.ticks.x.bottom       : NULL
 $ axis.minor.ticks.y.left         : NULL
 $ axis.minor.ticks.y.right        : NULL
 $ axis.minor.ticks.theta          : NULL
 $ axis.minor.ticks.r              : NULL
 $ axis.ticks.length               : 'rel' num 0.5
 $ axis.ticks.length.x             : NULL
 $ axis.ticks.length.x.top         : NULL
 $ axis.ticks.length.x.bottom      : NULL
 $ axis.ticks.length.y             : NULL
 $ axis.ticks.length.y.left        : NULL
 $ axis.ticks.length.y.right       : NULL
 $ axis.ticks.length.theta         : NULL
 $ axis.ticks.length.r             : NULL
 $ axis.minor.ticks.length         : 'rel' num 0.75
 $ axis.minor.ticks.length.x       : NULL
 $ axis.minor.ticks.length.x.top   : NULL
 $ axis.minor.ticks.length.x.bottom: NULL
 $ axis.minor.ticks.length.y       : NULL
 $ axis.minor.ticks.length.y.left  : NULL
 $ axis.minor.ticks.length.y.right : NULL
 $ axis.minor.ticks.length.theta   : NULL
 $ axis.minor.ticks.length.r       : NULL
 $ axis.line                       : <ggplot2::element_blank>
 $ axis.line.x                     : NULL
 $ axis.line.x.top                 : NULL
 $ axis.line.x.bottom              : NULL
 $ axis.line.y                     : NULL
 $ axis.line.y.left                : NULL
 $ axis.line.y.right               : NULL
 $ axis.line.theta                 : NULL
 $ axis.line.r                     : NULL
 $ legend.background               : <ggplot2::element_blank>
 $ legend.margin                   : NULL
 $ legend.spacing                  : 'rel' num 2
 $ legend.spacing.x                : NULL
 $ legend.spacing.y                : NULL
 $ legend.key                      : <ggplot2::element_blank>
 $ legend.key.size                 : 'simpleUnit' num 1.2lines
  ..- attr(*, "unit")= int 3
 $ legend.key.height               : NULL
 $ legend.key.width                : NULL
 $ legend.key.spacing              : NULL
 $ legend.key.spacing.x            : NULL
 $ legend.key.spacing.y            : NULL
 $ legend.key.justification        : NULL
 $ legend.frame                    : NULL
 $ legend.ticks                    : NULL
 $ legend.ticks.length             : 'rel' num 0.2
 $ legend.axis.line                : NULL
 $ legend.text                     : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : 'rel' num 0.8
  ..@ hjust        : NULL
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : NULL
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ legend.text.position            : NULL
 $ legend.title                    : <ggplot2::element_text>
  ..@ family       : NULL
  ..@ face         : NULL
  ..@ italic       : chr NA
  ..@ fontweight   : num NA
  ..@ fontwidth    : num NA
  ..@ colour       : NULL
  ..@ size         : NULL
  ..@ hjust        : num 0
  ..@ vjust        : NULL
  ..@ angle        : NULL
  ..@ lineheight   : NULL
  ..@ margin       : NULL
  ..@ debug        : NULL
  ..@ inherit.blank: logi TRUE
 $ legend.title.position           : NULL
 $ legend.position                 : chr "bottom"
 $ legend.position.inside          : NULL
 $ legend.direction                : NULL
 $ legend.byrow                    : NULL
 $ legend.justification            : chr "center"
 $ legend.justification.top        : NULL
 $ legend.justification.bottom     : NULL
 $ legend.justification.left       : NULL
 $ legend.justification.right      : NULL
 $ legend.justification.inside     : NULL
  [list output truncated]
 @ complete: logi TRUE
 @ validate: logi TRUE

print(p7)

# ============================================================
# BAGIAN 8: PSEUDO R² DAN GOODNESS-OF-FIT
# ============================================================

# McFadden Pseudo-R²
mcfadden <- function(model) {
  ll_null <- logLik(update(model, . ~ 1))
  ll_full <- logLik(model)
  1 - as.numeric(ll_full / ll_null)
}

cat("\n===== GOODNESS-OF-FIT =====\n")


===== GOODNESS-OF-FIT =====

gof <- data.frame(
  Ukuran         = c("Log-Likelihood", "AIC", "BIC", "McFadden Pseudo-R²"),
  Logit          = c(
    round(as.numeric(logLik(model_logit)),  3),
    round(AIC(model_logit),  3),
    round(BIC(model_logit),  3),
    round(mcfadden(model_logit),  4)
  ),
  Probit         = c(
    round(as.numeric(logLik(model_probit)), 3),
    round(AIC(model_probit), 3),
    round(BIC(model_probit), 3),
    round(mcfadden(model_probit), 4)
  )
)
print(gof)

              Ukuran     Logit    Probit
1     Log-Likelihood -294.6010 -294.2700
2                AIC  595.2020  594.5410
3                BIC  607.8460  607.1840
4 McFadden Pseudo-R²    0.1497    0.1507

# ============================================================
# BAGIAN 9: PANEL RINGKASAN (4 GRAFIK SEKALIGUS)
# ============================================================

grid.arrange(p1, p2, p4, p7,
             ncol = 2,
             top  = "Logit & Probit — Ringkasan Visual")

# ============================================================
# BAGIAN 10: PREDIKSI UNTUK OBSERVASI BARU
# ============================================================

data_baru <- data.frame(
  X1 = c(-2, -1, 0, 1, 2),
  X2 = factor(c("Grup A", "Grup A", "Grup B", "Grup B", "Grup A"),
              levels = c("Grup A", "Grup B"))
)

pred_baru <- data_baru
pred_baru$P_logit  <- round(predict(model_logit,  newdata = data_baru, type = "response"), 4)
pred_baru$P_probit <- round(predict(model_probit, newdata = data_baru, type = "response"), 4)
pred_baru$Selisih  <- round(abs(pred_baru$P_logit - pred_baru$P_probit), 4)

cat("\n===== PREDIKSI DATA BARU =====\n")


===== PREDIKSI DATA BARU =====

print(pred_baru)

  X1     X2 P_logit P_probit Selisih
1 -2 Grup A  0.0913   0.0809  0.0104
2 -1 Grup A  0.2295   0.2302  0.0007
3  0 Grup B  0.5743   0.5734  0.0009
4  1 Grup B  0.7999   0.8012  0.0013
5  2 Grup A  0.8857   0.8933  0.0076

cat("\n✓ Selesai! Semua grafik dan output telah ditampilkan.\n")


✓ Selesai! Semua grafik dan output telah ditampilkan.

cat("  Gunakan Ctrl+Z (undo) di RStudio jika ingin kembali ke grafik sebelumnya.\n")

  Gunakan Ctrl+Z (undo) di RStudio jika ingin kembali ke grafik sebelumnya.