1 Persiapan Data

1.1 Deskripsi Data

Data yang digunakan adalah General Social Survey (GSS) 2024, yaitu data survei nasional Amerika Serikat yang mengumpulkan informasi tentang karakteristik demografis, perilaku sosial, dan opini masyarakat.

Variabel	Tipe	Peran	Keterangan
`educ`	Kategorik/Ordinal	Prediktor (X)	Tingkat pendidikan tertinggi
`gunlaw`	Biner	Respons (Y1)	Sikap terhadap izin kepemilikan senjata
`marital`	Nominal	Respons (Y2)	Status pernikahan
`happy`	Ordinal	Respons (Y3)	Tingkat kebahagiaan
`tvhours`	Count	Respons (Y4)	Jam menonton TV per hari

1.2 Import dan Cleaning Data

# Import data dari file Excel
gss_raw <- read_excel("GSS.xlsx", sheet = "Data")

# Konversi educ ke numerik (jumlah tahun sekolah)
educ_map <- c(
  "No formal schooling" = 0,
  "1st grade" = 1, "2nd grade" = 2, "3rd grade" = 3, "4th grade" = 4,
  "5th grade" = 5, "6th grade" = 6, "7th grade" = 7, "8th grade" = 8,
  "9th grade" = 9, "10th grade" = 10, "11th grade" = 11, "12th grade" = 12,
  "1 year of college" = 13, "2 years of college" = 14, "3 years of college" = 15,
  "4 years of college" = 16, "5 years of college" = 17, "6 years of college" = 18,
  "7 years of college" = 19, "8 or more years of college" = 20
)

# Bersihkan tvhours (menghapus 0 hours dan 24 hours karne bukan numerik)
gss <- gss_raw %>%
  mutate(
    educ_num = educ_map[educ],
    
    tvhours_num = as.numeric(gsub("[^0-9]", "", tvhours)),
    
    gunlaw = factor(gunlaw, levels = c("OPPOSE", "FAVOR")),
    
    marital = factor(marital,
      levels = c("Married", "Widowed", "Divorced", "Separated", "Never married")),
    
    happy = factor(happy,
      levels = c("Not too happy", "Pretty happy", "Very happy"),
      ordered = TRUE)
  ) %>%
  filter(!is.na(educ_num), !is.na(gunlaw), !is.na(marital),
         !is.na(happy), !is.na(tvhours_num))

gss[, c("educ","educ_num","gunlaw","marital","happy","tvhours_num")] |> head(8) |>
  kable(caption = "Contoh 8 baris pertama data GSS (setelah cleaning)",
        col.names = c("educ", "educ_num", "gunlaw", "marital", "happy", "tvhours")) %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Contoh 8 baris pertama data GSS (setelah cleaning)
educ	educ_num	gunlaw	marital	happy	tvhours
12th grade	12	FAVOR	Never married	Pretty happy	2
12th grade	12	FAVOR	Never married	Pretty happy	24
12th grade	12	FAVOR	Divorced	Very happy	0
7 years of college	19	OPPOSE	Married	Very happy	2
1 year of college	13	FAVOR	Married	Very happy	6
8th grade	8	FAVOR	Divorced	Pretty happy	1
1 year of college	13	FAVOR	Divorced	Pretty happy	4
4 years of college	16	FAVOR	Married	Pretty happy	12

2 Regresi Logistik Biner

2.1 Definisi dan Data

Regresi Logistik Biner digunakan ketika variabel respons bersifat dikotomis (dua kategori).

Variabel dalam analisis ini:

Y gunlaw — Sikap terhadap izin kepemilikan senjata: FAVOR vs OPPOSE
X educ_num — Lama pendidikan (0–20 tahun)

# Distribusi variabel respons
gss %>%
  count(gunlaw) %>%
  mutate(pct = n / sum(n)) %>%
  ggplot(aes(x = gunlaw, y = n, fill = gunlaw)) +
  geom_col(width = 0.55, color = "white") +
  geom_text(aes(label = paste0(n, "\n(", percent(pct, 1), ")")),
            vjust = -0.3, fontface = "bold", size = 4) +
  scale_fill_manual(values = c("OPPOSE" = "#e74c3c", "FAVOR" = "#2980b9")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title = "Distribusi Sikap terhadap Izin Kepemilikan Senjata",
       x = NULL, y = "Frekuensi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", color = "#1a3a5c"))

2.2 Model dan Interpretasi

# Membuat variabel biner (1 = FAVOR, 0 = OPPOSE)
gss <- gss %>% mutate(favor_gun = ifelse(gunlaw == "FAVOR", 1, 0))

# Fitting model
fit_biner <- glm(favor_gun ~ educ_num, data = gss, family = binomial(link = "logit"))
summary(fit_biner)

## 
## Call:
## glm(formula = favor_gun ~ educ_num, family = binomial(link = "logit"), 
##     data = gss)
## 
## Coefficients:
##             Estimate Std. Error z value Pr(>|z|)  
## (Intercept)  0.16670    0.34947   0.477   0.6334  
## educ_num     0.04618    0.02423   1.906   0.0567 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1290.4  on 1048  degrees of freedom
## Residual deviance: 1286.8  on 1047  degrees of freedom
## AIC: 1290.8
## 
## Number of Fisher Scoring iterations: 4

# Ringkasan koefisien
tidy_biner <- broom::tidy(fit_biner) %>%
  mutate(
    odds_ratio = exp(estimate),
    ci_low  = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error),
    p_value = signif(p.value, 3)
  )

tidy_biner %>%
  transmute(
    Parameter    = term,
    Koefisien    = round(estimate, 4),
    `Std. Error` = round(std.error, 4),
    `Odds Ratio` = round(odds_ratio, 4),
    `95% CI`     = paste0("[", round(ci_low, 3), ", ", round(ci_high, 3), "]"),
    `p-value`    = p_value
  ) %>%
  kable(caption = "Hasil Regresi Logistik Biner: gunlaw ~ educ_num") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Hasil Regresi Logistik Biner: gunlaw ~ educ_num
Parameter	Koefisien	Std. Error	Odds Ratio	95% CI	p-value
(Intercept)	0.1667	0.3495	1.1814	[0.596, 2.344]	0.6330
educ_num	0.0462	0.0242	1.0473	[0.999, 1.098]	0.0567

Hasil Regresi Logistik Biner

Model: logit(P(FAVOR)) = 0.1667 + 0.0462 × educ_num

Berdasarkan output model:

Intercept (β₀=0.1667): Log-odds seseorang mendukung izin kepemilikan senjata saat educ_num = 0adalah 0.1667
Koefisien educ_num (β₁=0.0462): Setiap penambahan 1 tahun pendidikan, log-odds mendukung izin senjata berubah sebesar 0.0462
Odds Ratio (OR=1.0473): Setiap penambahan 1 tahun pendidikan, odds mendukung FAVOR naik 4.73%. Karena OR > 1, makasemakin berpendidikan semakin cenderung mendukung izin kepemilikan senjata
p-value (0.0567): Pada taraf signifikansi 5%, pendidikan tidak berpengaruh signifikan terhadap sikap responden karena p > 0,05. Interval kepercayaan 95% odds ratio, yaitu [0,999; 1,098] yang mencakup nilai 1 menunjukkan bahwa tidak terdapat perbedaan antar tingkat pendidikan.

2.3 Visualisasi dan Interpretasi

# Kurva probabilitas prediksi
grid_biner <- data.frame(educ_num = seq(0, 20, length.out = 200))
grid_biner$prob <- predict(fit_biner, newdata = grid_biner, type = "response")

ggplot() +
  geom_point(data = gss, aes(x = educ_num, y = favor_gun),
             alpha = 0.08, color = "#7f8c8d", shape = 16, size = 2) +
  geom_line(data = grid_biner, aes(x = educ_num, y = prob),
            color = "#2980b9", linewidth = 1.4) +
  scale_y_continuous(labels = percent_format()) +
  labs(title = "Probabilitas Prediksi: FAVOR Izin Kepemilikan Senjata",
       x = "Lama Pendidikan (tahun)", y = "P(FAVOR)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a3a5c"))

3 Regresi Logistik Multinomial

3.1 Definisi dan Data

Regresi Logistik Multinomial adalah model regresi untuk variabel respons yang memiliki tiga kategori atau lebih dan tidak memiliki urutan (nominal).

Variabel dalam analisis ini:

Y marital — Status pernikahan: Married / Widowed / Divorced / Separated / Never married
X educ_num — Lama pendidikan (0–20 tahun)
Kategori referensi: Married

gss %>%
  count(marital) %>%
  mutate(pct = n / sum(n)) %>%
  ggplot(aes(x = reorder(marital, n), y = n, fill = marital)) +
  geom_col(width = 0.65, color = "white") +
  geom_text(aes(label = paste0(n, " (", percent(pct, 1), ")")),
            hjust = -0.08, size = 3.5, fontface = "bold") +
  coord_flip() +
  scale_fill_manual(values = c("#2980b9","#27ae60","#e67e22","#8e44ad","#c0392b")) +
  scale_x_discrete(labels = function(x) tools::toTitleCase(tolower(x))) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.25))) +
  labs(title = "Distribusi Status Pernikahan (marital)",
       x = NULL, y = "Frekuensi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", color = "#1a3a5c"))

3.2 Model dan Interpretasi

# Set kategori referensi = Married
gss$marital <- relevel(gss$marital, ref = "Married")

# Fitting model multinomial
fit_multi <- nnet::multinom(marital ~ educ_num, data = gss, trace = FALSE)

# Koefisien dan RRR
multi_sum <- summary(fit_multi)

coef_m <- as.data.frame(multi_sum$coefficients)
se_m   <- as.data.frame(multi_sum$standard.errors)

result_multi <- coef_m %>%
  tibble::rownames_to_column("kategori") %>%
  pivot_longer(cols = -kategori, names_to = "variabel", values_to = "estimate") %>%
  left_join(
    se_m %>%
      tibble::rownames_to_column("kategori") %>%
      pivot_longer(cols = -kategori, names_to = "variabel", values_to = "std_error"),
    by = c("kategori", "variabel")
  ) %>%
  mutate(
    z_value = estimate / std_error,
    p_value = 2 * (1 - pnorm(abs(z_value))),
    RRR = exp(estimate),
    CI   = paste0("[", round(exp(estimate - 1.96*std_error), 3),
                  ", ", round(exp(estimate + 1.96*std_error), 3), "]")
  )

result_multi %>%
  filter(variabel != "(Intercept)") %>%
  transmute(
    `Kategori (vs Married)` = kategori,
    Koefisien = round(estimate, 4),
    `Std. Error` = round(std_error, 4),
    `z-value` = round(z_value, 3),
    RRR = round(RRR, 4),
    `95% CI` = CI,
    `p-value` = signif(p_value, 3)
  ) %>%
  kable(caption = "Koefisien Regresi Logistik Multinomial: marital ~ educ_num (Referensi: Married)") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Koefisien Regresi Logistik Multinomial: marital ~ educ_num (Referensi: Married)
Kategori (vs Married)	Koefisien	Std. Error	z-value	RRR	95% CI	p-value
Widowed	-0.1327	0.0452	-2.935	0.8757	[0.801, 0.957]	0.00334
Divorced	-0.0577	0.0337	-1.714	0.9439	[0.884, 1.008]	0.08660
Separated	-0.0881	0.0627	-1.404	0.9157	[0.81, 1.035]	0.16000
Never married	-0.0855	0.0266	-3.208	0.9181	[0.871, 0.967]	0.00134

Hasil Regresi Logistik Multinomial

Model: Regresi Logistik Multinomial, kategori referensi = Married

log(P(Widowed)/P(Married)) = β₀ + (-0.1327) × educ_num

Widowed (RRR = 0.8757, p = 0.003): Setiap tambah 1 tahun pendidikan, odds menjadi janda/duda dibanding menikah turun 12.43%. Signifikan (p < 0.05). CI 95% [0.801, 0.957] tidak mencakup angka 1, mengonfirmasi hasil signifikan.

log(P(Divorced)/P(Married)) = β₀ + (-0.0577) × educ_num

Divorced (RRR = 0.9439, p = 0.087): Setiap tambah 1 tahun pendidikan, odds menjadi bercerai dibanding menikah turun 5.61%. Tidak signifikan (p > 0.05). CI 95% [0.884, 1.008] mencakup angka 1, mengonfirmasi hasil tidak signifikan.

log(P(Separated)/P(Married)) = β₀ + (-0.0881) × educ_num

Separated (RRR = 0.9157, p = 0.160): Setiap tambah 1 tahun pendidikan, odds menjadi pisah dibanding menikah turun 8.43%. Tidak signifikan (p > 0.05). CI 95% [0.810, 1.035] mencakup angka 1, mengonfirmasi hasil tidak signifikan.

log(P(Never married)/P(Married)) = β₀ + (-0.0855) × educ_num

Never married (RRR = 0.9181, p = 0.001): Setiap tambah 1 tahun pendidikan, odds belum pernah menikah dibanding menikah turun 8.19%. Signifikan (p < 0.05). CI 95% [0.871, 0.967] tidak mencakup angka 1, mengonfirmasi hasil signifikan.

Kesimpulan: Semakin tinggi tingkat pendidikan, seseorang cenderung lebih besar kemungkinannya berstatus Married dibanding status lainnya. Pengaruh pendidikan hanya signifikan terhadap kategori Widowed dan Never married sementara Divorced dan Separated tidak menunjukkan pengaruh yang signifikan secara statistik.

3.3 Visualisasi dan Interpretasi

# Probabilitas prediksi sepanjang educ_num
grid_m <- data.frame(educ_num = seq(0, 20, length.out = 200))
prob_m  <- predict(fit_multi, newdata = grid_m, type = "probs")

plot_m <- grid_m %>%
  bind_cols(as.data.frame(prob_m)) %>%
  pivot_longer(cols = -educ_num, names_to = "Status", values_to = "Probabilitas")

ggplot(plot_m, aes(x = educ_num, y = Probabilitas, color = Status)) +
  geom_line(linewidth = 1.2) +
  scale_color_manual(values = c("#2980b9","#27ae60","#e67e22","#8e44ad","#c0392b")) +
  scale_y_continuous(labels = percent_format()) +
  labs(title = "Probabilitas Prediksi Status Pernikahan berdasarkan Pendidikan",
       x = "Lama Pendidikan (tahun)", y = "Probabilitas", color = "Status") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a3a5c"),
        legend.position = "right")

Interpretasi Grafik:

Married (hijau): Probabilitas berstatus menikah terus naik seiring meningkatnya pendidikan, dari sekitar 17% (0 tahun) hingga lebih dari 52% (20 tahun).
Never married (oranye): Probabilitas berstatus belum pernah menikah terus turun seiring pendidikan meningkat, dari sekitar 46% (0 tahun) menjadi sekitar 26% (20 tahun). Garis Married dan Never married saling bersilangan di sekitar 12 tahun pendidikan.
Widowed (merah): Probabilitas berstatus janda/duda turun seiring pendidikan, dari sekitar 19% (0 tahun) menjadi sekitar 5% (20 tahun).
Divorced (biru): Probabilitas berstatus bercerai relatif stabil di angka 15%-16% sepanjang rentang pendidikan.
Separated (ungu): Probabilitas berstatus pisah sedikit turun dan tetap rendah di sekitar 4%-6% sepanjang rentang pendidikan.

Kesimpulan grafik: Semakin tinggi pendidikan, seseorang semakin besar kemungkinannya berstatus Married dan semakin kecil kemungkinannya berstatus Never married maupun Widowed. Hal ini konsisten dengan hasil tabel koefisien yang menunjukkan RRR < 1 untuk semua kategori dibanding Married.

4 Regresi Logistik Ordinal

4.1 Definisi dan Data

Regresi Logistik Ordinal (Proportional Odds Model) digunakan ketika variabel respons memiliki tiga kategori atau lebih dengan urutan (ordinal).

Variabel dalam analisis ini:

Y happy — Tingkat kebahagiaan: Not too happy < Pretty happy < Very happy
X educ_num — Lama pendidikan (0–20 tahun)

gss %>%
  count(happy) %>%
  mutate(pct = n / sum(n)) %>%
  ggplot(aes(x = happy, y = n, fill = happy)) +
  geom_col(width = 0.55, color = "white") +
  geom_text(aes(label = paste0(n, "\n(", percent(pct, 1), ")")),
            vjust = -0.3, fontface = "bold", size = 3.8) +
  scale_fill_manual(values = c("#e74c3c", "#f39c12", "#27ae60")) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.15))) +
  labs(title = "Distribusi Tingkat Kebahagiaan (happy)",
       x = NULL, y = "Frekuensi") +
  theme_minimal(base_size = 12) +
  theme(legend.position = "none",
        plot.title = element_text(face = "bold", color = "#1a3a5c"))

4.2 Model dan Interpretasi

# Fitting model ordinal (proportional odds)
fit_ord <- MASS::polr(
  happy ~ educ_num,
  data   = gss,
  method = "logistic",
  Hess   = TRUE
)

# Ringkasan hasil
ord_coef <- coef(summary(fit_ord))

result_ord <- as.data.frame(ord_coef) %>%
  tibble::rownames_to_column("parameter") %>%
  rename(estimate = Value, std_error = `Std. Error`, t_value = `t value`) %>%
  mutate(
    p_value = 2 * (1 - pnorm(abs(t_value))),
    jenis   = ifelse(grepl("\\|", parameter), "Cutpoint (τ)", "Koefisien"),
    OR      = ifelse(jenis == "Koefisien", exp(estimate), NA_real_),
    `95% CI` = ifelse(
      jenis == "Koefisien",
      paste0("[", round(exp(estimate - 1.96*std_error), 3),
             ", ", round(exp(estimate + 1.96*std_error), 3), "]"),
      "—"
    )
  )

result_ord %>%
  transmute(
    Parameter = parameter,
    Jenis = jenis,
    Koefisien = round(estimate, 4),
    `Std. Error` = round(std_error, 4),
    `t-value` = round(t_value, 3),
    OR = round(OR, 4),
    `95% CI`,
    `p-value` = signif(p_value, 3)
  ) %>%
  kable(caption = "Hasil Regresi Logistik Ordinal: happy ~ educ_num") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE) %>%
  row_spec(which(result_ord$jenis == "Cutpoint (τ)"), background = "#f8f9fa")

Hasil Regresi Logistik Ordinal: happy ~ educ_num
Parameter	Jenis	Koefisien	Std. Error	t-value	OR	95% CI	p-value
educ_num	Koefisien	0.0396	0.0225	1.761	1.0404	[0.996, 1.087]	7.83e-02
Not too happy\|Pretty happy	Cutpoint (τ)	-1.0113	0.3290	-3.074	NA	—	2.11e-03
Pretty happy\|Very happy	Cutpoint (τ)	1.7691	0.3330	5.312	NA	—	1.00e-07

Hasil Regresi Logistik Ordinal

Model: log(P(Y ≤ j) / P(Y > j)) = τⱼ − (0.0396) × educ_num

educ_num (OR = 1.0404, p = 0.0783): Setiap tambah 1 tahun pendidikan, odds berada di kategori kebahagiaan lebih tinggi dibanding kategori di bawahnya naik 4.04%. Tidak signifikan (p > 0.05). CI 95% [0.996, 1.087] mencakup angka 1, mengonfirmasi hasil tidak signifikan.

log(P(Y ≤ "Not too happy") / P(Y > "Not too happy")) = −1.0113

Cutpoint τ₁ = −1.0113 (p = 0.002): Titik pisah antara kategoriNot too happy dan Pretty happy. Log-odds kumulatif bernilai negatif, artinya sebagian besar responden berada di atas kategori terendah. Signifikan (p < 0.05).

log(P(Y ≤ "Pretty happy") / P(Y > "Pretty happy")) = 1.7691

Cutpoint τ₂ = 1.7691 (p < 0.001): Titik pisah antara kategori Pretty happy dan Very happy. Log-odds kumulatif bernilai positif, artinya mayoritas responden berada di bawah kategori tertinggi. Sangat signifikan (p < 0.05).

Kesimpulan: Arah hubungan pendidikan terhadap kebahagiaan bersifat positif — semakin tinggi pendidikan, semakin besar kemungkinan seseorang berada di kategori kebahagiaan yang lebih tinggi (OR = 1.0404). Namun pengaruh ini tidak signifikan secara statistik (p = 0.0783 > 0.05), sehingga tidak cukup bukti untuk menyimpulkan bahwa pendidikan berpengaruh nyata terhadap tingkat kebahagiaan secara ordinal dalam data ini.

# Probabilitas prediksi ordinal
grid_ord  <- data.frame(educ_num = seq(0, 20, length.out = 200))
prob_ord  <- predict(fit_ord, newdata = grid_ord, type = "probs")

plot_ord <- grid_ord %>%
  bind_cols(as.data.frame(prob_ord)) %>%
  pivot_longer(cols = -educ_num, names_to = "Kebahagiaan", values_to = "Probabilitas") %>%
  mutate(Kebahagiaan = factor(Kebahagiaan,
    levels = c("Not too happy", "Pretty happy", "Very happy")))

ggplot(plot_ord, aes(x = educ_num, y = Probabilitas, color = Kebahagiaan)) +
  geom_line(linewidth = 1.3) +
  scale_color_manual(values = c("#e74c3c", "#f39c12", "#27ae60")) +
  scale_y_continuous(labels = percent_format()) +
  labs(title = "Probabilitas Prediksi Tingkat Kebahagiaan berdasarkan Pendidikan",
       x = "Lama Pendidikan (tahun)", y = "Probabilitas", color = "Kebahagiaan") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a3a5c"),
        legend.position = "right")

4.3 Interpretasi

Interpretasi Grafik

Not too happy (merah): Probabilitas berada di kategori Not too happy menurun seiring bertambahnya lama pendidikan, yaitu dari sekitar 27% pada 0 tahun pendidikan turun menjadi sekitar 15% pada 20 tahun pendidikan. Ini menunjukkan bahwa pendidikan yang lebih tinggi mengurangi kemungkinan seseorang berada di kategori kebahagiaan terendah.
Pretty happy (oranye): Probabilitas berada di kategori Pretty happy relatif stabil di kisaran 59–60% sepanjang rentang pendidikan 0–20 tahun, dengan sedikit penurunan di ujung kanan. Kategori ini mendominasi di seluruh level pendidikan, artinya sebagian besar responden cenderung berada di kategori tengah terlepas dari tingkat pendidikannya.
Very happy (hijau): Probabilitas berada di kategori Very happy meningkat seiring bertambahnya lama pendidikan — dari sekitar 15% pada 0 tahun pendidikan naik menjadi sekitar 28% pada 20 tahun pendidikan. Ini menunjukkan bahwa pendidikan yang lebih tinggi meningkatkan kemungkinan seseorang berada di kategori kebahagiaan tertinggi.
Titik persilangan (~10 tahun): Kurva Not too happy dan Very happy berpotongan di sekitar 10 tahun pendidikan. Di bawah 10 tahun, probabilitas Not too happy lebih besar dari Very happy; di atas 10 tahun, kondisinya berbalik — Very happy melampaui Not too happy.

Kesimpulan: Grafik memperlihatkan pola yang konsisten dengan koefisien positif model ordinal (OR = 1.0404), semakin lama pendidikan, distribusi probabilitas bergeser dari kategori kebahagiaan rendah menuju kategori kebahagiaan tinggi. Meski pengaruh pendidikan tidak signifikan secara statistik (p = 0.0783), arah tren dalam grafik ini tetap menunjukkan hubungan positif yang substantif antara pendidikan dan tingkat kebahagiaan.

5 Regresi Poisson

5.1 Definisi dan Data

Regresi Poisson digunakan untuk memodelkan variabel respons berupa data cacah (count) atau bilangan bulat non-negatif yang mengikuti distribusi Poisson.

Variabel dalam analisis ini:

Y tvhours — Jumlah jam menonton TV per hari (0, 1, 2, …, 24)
X educ_num — Lama pendidikan (0–20 tahun)

gss %>%
  count(tvhours_num) %>%
  filter(tvhours_num <= 15) %>%
  ggplot(aes(x = tvhours_num, y = n)) +
  geom_col(fill = "#2980b9", width = 0.75, alpha = 0.88, color = "white") +
  geom_text(aes(label = n), vjust = -0.4, size = 3, fontface = "bold") +
  scale_x_continuous(breaks = 0:15) +
  scale_y_continuous(expand = expansion(mult = c(0, 0.12))) +
  labs(title = "Distribusi Jam Menonton TV per Hari (tvhours)",
       subtitle = "Data cacah cocok untuk Regresi Poisson",
       x = "Jam/hari", y = "Frekuensi") +
  theme_minimal(base_size = 12) +
  theme(plot.title    = element_text(face = "bold", color = "#1a3a5c"),
        plot.subtitle = element_text(color = "#7f8c8d"))

5.2 Model dan Interpretasi

# Fitting Regresi Poisson
fit_pois <- glm(
  tvhours_num ~ educ_num,
  data   = gss,
  family = poisson(link = "log")
)

# Ringkasan koefisien
tidy_pois <- broom::tidy(fit_pois) %>%
  mutate(
    IRR     = exp(estimate),
    ci_low  = exp(estimate - 1.96 * std.error),
    ci_high = exp(estimate + 1.96 * std.error)
  )

tidy_pois %>%
  transmute(
    Parameter = term,
    Koefisien = round(estimate, 4),
    `Std. Error` = round(std.error, 4),
    `z-value` = round(statistic, 3),
    IRR = round(IRR, 4),
    `95% CI` = paste0("[", round(ci_low, 3), ", ", round(ci_high, 3), "]"),
    `p-value` = signif(p.value, 3)
  ) %>%
  kable(caption = "Hasil Regresi Poisson: tvhours ~ educ_num") %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed"), full_width = FALSE)

Hasil Regresi Poisson: tvhours ~ educ_num
Parameter	Koefisien	Std. Error	z-value	IRR	95% CI	p-value
(Intercept)	2.0474	0.0857	23.896	7.7481	[6.55, 9.165]	0
educ_num	-0.0634	0.0061	-10.426	0.9385	[0.927, 0.95]	0

# Informasi model
cat("Null Deviance    :", round(fit_pois$null.deviance, 2), 
    "(df =", fit_pois$df.null, ")\n")

## Null Deviance    : 2599.01 (df = 1048 )

cat("Residual Deviance:", round(fit_pois$deviance, 2),
    "(df =", fit_pois$df.residual, ")\n")

## Residual Deviance: 2492.53 (df = 1047 )

cat("AIC              :", round(AIC(fit_pois), 2), "\n")

## AIC              : 5258.72

Hasil Regresi Poisson: tvhours ~ educ_num

Model: Regresi Poisson, memodelkan log rata-rata jam menonton TV sebagai fungsi pendidikan:

log(μ) = 2.0474 + (−0.0634) × educ_num

Intercept (IRR = 7.7481, p < 0.001): Ketika lama pendidikan = 0 tahun, rata-rata jam menonton TV per hari diprediksi sebesar 7.75 jam. Sangat signifikan (p < 0.05). CI 95% [6.55, 9.165] tidak mencakup angka 0, mengonfirmasi hasil signifikan.
educ_num (IRR = 0.9385, p < 0.001): Setiap tambah 1 tahun pendidikan, rata-rata jam menonton TV turun 6.15% (IRR = 0.9385). Signifikan (p < 0.05). CI 95% [0.927, 0.950] tidak mencakup angka 1, mengonfirmasi hasil signifikan.
Goodness of Fit: Null deviance = 2599.01 turun menjadi residual deviance = 2492.53 setelah memasukkan educ_num, dengan selisih 106.48 pada 1 df. Penurunan ini menunjukkan bahwa pendidikan mampu menjelaskan sebagian variasi dalam jam menonton TV. AIC model = 5258.72 dapat digunakan sebagai pembanding jika ada model alternatif.

Kesimpulan: Tingkat pendidikan berpengaruh negatif dan signifikan terhadap jumlah jam menonton TV (p < 0.001). Semakin tinggi pendidikan seseorang, semakin sedikit waktu yang dihabiskan untuk menonton TV per harinya. Setiap tambahan 1 tahun pendidikan secara rata-rata mengurangi jam menonton TV sebesar 6.15%.

# Rata-rata prediksi
grid_pois <- data.frame(educ_num = seq(0, 20, length.out = 200))
grid_pois$mu <- predict(fit_pois, newdata = grid_pois, type = "response")

ggplot() +
  geom_jitter(data = gss %>% filter(tvhours_num <= 20),
              aes(x = educ_num, y = tvhours_num),
              alpha = 0.12, color = "#7f8c8d", height = 0.25, width = 0.3, size = 1.5) +
  geom_line(data = grid_pois, aes(x = educ_num, y = mu),
            color = "#2980b9", linewidth = 1.5) +
  labs(title = "Rata-rata Prediksi Jam Menonton TV berdasarkan Pendidikan",
       x = "Lama Pendidikan (tahun)", y = "Jam TV per Hari (prediksi rata-rata)") +
  theme_minimal(base_size = 12) +
  theme(plot.title = element_text(face = "bold", color = "#1a3a5c"))

Interpretasi Grafik

Garis prediksi biru menunjukkan pola menurun eksponensial, yaitu dari sekitar 7–8 jam/hari pada 0 tahun pendidikan hingga 2–3 jam/hari pada 20 tahun pendidikan, konsisten dengan koefisien negatif model (IRR = 0.9385). Titik-titik observasi tersebar cukup lebar di atas garis prediksi, terutama pada rentang 10–20 tahun pendidikan dengan beberapa nilai ekstrem hingga 20 jam/hari, mengindikasikan kemungkinan overdispersi sehingga perlu dipertimbangkan penggunaan model Negative Binomial sebagai alternatif.

6 Kesimpulan

Kesimpulan Keempat Model Regresi

Keempat model regresi di atas secara konsisten menunjukkan bahwa tingkat pendidikan memiliki peran yang beragam terhadap berbagai aspek kehidupan sosial responden GSS 2024.

Regresi Logistik Biner (gunlaw): Pendidikan berpengaruh positif terhadap dukungan izin kepemilikan senjata (OR = 1.0473), namun tidak signifikan (p = 0.057). Arah hubungan menunjukkan bahwa semakin berpendidikan, seseorang sedikit lebih cenderung mendukung regulasi senjata.
Regresi Logistik Multinomial (marital): Pendidikan berpengaruh negatif terhadap semua status pernikahan non-Married. Pengaruh signifikan hanya ditemukan pada kategori Widowed (RRR = 0.876, p = 0.003) dan Never married (RRR = 0.918, p = 0.001), yang menunjukkan bahwa semakin tinggi pendidikan, maka seseorang semakin cenderung berstatus menikah.
Regresi Logistik Ordinal (happy): Pendidikan berpengaruh positif terhadap tingkat kebahagiaan (OR = 1.040), namun tidak signifikan (p = 0.078). Arah tren tetap menunjukkan bahwa semakin berpendidikan, maka seseorang sedikit lebih mungkin berada di kategori kebahagiaan yang lebih tinggi.
Regresi Poisson (tvhours): Pendidikan berpengaruh negatif dan signifikan terhadap jam menonton TV (IRR = 0.939, p < 0.001). Setiap tambahan 1 tahun pendidikan, rata-rata jam menonton TV berkurang 6.15%. Ini merupakan hubungan yang paling kuat di antara keempat model.

Secara keseluruhan, pendidikan paling nyata pengaruhnya terhadap perilaku menonton TV dan status pernikahan, sementara pengaruhnya terhadap sikap kepemilikan senjata dan tingkat kebahagiaan belum cukup terbukti secara statistik dalam data ini.

Analisis Regresi Logistik

Analisis Regresi Logistik biner, Multinomial, Ordinal, dan Poisson

Fathiya Syakira

01 June 2026

1 Persiapan Data

1.1 Deskripsi Data

1.2 Import dan Cleaning Data

2 Regresi Logistik Biner

2.1 Definisi dan Data

2.2 Model dan Interpretasi

2.3 Visualisasi dan Interpretasi

3 Regresi Logistik Multinomial

3.1 Definisi dan Data

3.2 Model dan Interpretasi

3.3 Visualisasi dan Interpretasi

4 Regresi Logistik Ordinal

4.1 Definisi dan Data

4.2 Model dan Interpretasi

4.3 Interpretasi

5 Regresi Poisson

5.1 Definisi dan Data

5.2 Model dan Interpretasi

6 Kesimpulan