1 Persiapan Data

1.1 Load Library dan Data

library(car)      
library(ggplot2)

# Membaca data
data <- read.csv("C:/Users/Asus/Downloads/bahagia.csv")

# Melihat struktur awal
str(data)

## 'data.frame':    185331 obs. of  4 variables:
##  $ kepuasan: chr  "Sangat tidak puas" "Sangat tidak puas" "Sangat tidak puas" "Sangat tidak puas" ...
##  $ bahagia : chr  "Bahagia" "Bahagia" "Bahagia" "Bahagia" ...
##  $ usia    : int  59 29 39 16 12 3 1 29 30 6 ...
##  $ gender  : chr  "Laki-laki" "Perempuan" "Perempuan" "Perempuan" ...

head(data)

##            kepuasan bahagia usia    gender
## 1 Sangat tidak puas Bahagia   59 Laki-laki
## 2 Sangat tidak puas Bahagia   29 Perempuan
## 3 Sangat tidak puas Bahagia   39 Perempuan
## 4 Sangat tidak puas Bahagia   16 Perempuan
## 5 Sangat tidak puas Bahagia   12 Perempuan
## 6 Sangat tidak puas Bahagia    3 Laki-laki

1.2 Pra-pemrosesan Data

Sesuai soal:

Variabel respons: bahagia2 — biner: "Bahagia" jika individu bahagia, "Tidak" jika tidak bahagia.
Kategori asli "Bahagia" dan "Sangat bahagia" → Bahagia;
kategori "Tidak bahagia" dan "Sangat tidak bahagia" → Tidak.
Filter usia: 20–40 tahun.
Variabel penjelas: usia (kuantitatif) dan kepuasan (kategorik ordinal, 5 level).

# Membuat variabel respons biner bahagia2
data$bahagia2 <- ifelse(
  data$bahagia %in% c("Bahagia", "Sangat bahagia"), "Bahagia", "Tidak"
)

# Filter usia 20-40 tahun
df <- subset(data, usia >= 20 & usia <= 40)

# Relevel kepuasan sesuai urutan ordinal
df$kepuasan <- factor(
  df$kepuasan,
  levels = c("Sangat tidak puas", "Tidak puas", "Agak puas", "Puas", "Sangat puas")
)

# Variabel respons numerik: 1 = Bahagia, 0 = Tidak
df$y <- ifelse(df$bahagia2 == "Bahagia", 1, 0)

# Ringkasan data setelah filter
cat("Jumlah observasi:", nrow(df), "\n")

## Jumlah observasi: 78863

table(df$kepuasan)

## 
## Sangat tidak puas        Tidak puas         Agak puas              Puas 
##              1262              9755             33318             31194 
##       Sangat puas 
##              3334

table(df$bahagia2)

## 
## Bahagia   Tidak 
##   71847    7016

2 (a) Estimasi Model Regresi Logistik

2.1 Model

Model regresi logistik berganda dengan variabel penjelas usia (\(x\)) dan kepuasan (variabel kategorik dengan kategori referensi Sangat tidak puas) adalah:

\[ \text{logit}\,\hat{P}(Y=1) = \hat{\beta}_0 + \hat{\beta}_1 x + \hat{\beta}_2 c_2 + \hat{\beta}_3 c_3 + \hat{\beta}_4 c_4 + \hat{\beta}_5 c_5 \]

di mana variabel indikator kepuasan didefinisikan sebagai:

Indikator	Kategori
\(c_2 = 1\)	Tidak puas
\(c_3 = 1\)	Agak puas
\(c_4 = 1\)	Puas
\(c_5 = 1\)	Sangat puas

dan \(c_j = 0\) untuk lainnya. Kategori referensi adalah Sangat tidak puas.

# Estimasi model regresi logistik berganda
fit <- glm(y ~ usia + kepuasan, family = binomial, data = df)
summary(fit)

## 
## Call:
## glm(formula = y ~ usia + kepuasan, family = binomial, data = df)
## 
## Coefficients:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)         0.326651   0.089211   3.662 0.000251 ***
## usia                0.002409   0.002299   1.048 0.294749    
## kepuasanTidak puas  0.307056   0.061332   5.006 5.54e-07 ***
## kepuasanAgak puas   2.307192   0.061726  37.378  < 2e-16 ***
## kepuasanPuas        2.940344   0.065326  45.011  < 2e-16 ***
## kepuasanSangat puas 2.781718   0.105494  26.368  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47339  on 78862  degrees of freedom
## Residual deviance: 40095  on 78857  degrees of freedom
## AIC: 40107
## 
## Number of Fisher Scoring iterations: 6

2.1.1 Interpretasi Output

Dari hasil estimasi diperoleh koefisien model. Parameter \(\hat{\beta}_j\) pada variabel kepuasan menyatakan pengaruh kategori kepuasan tersebut terhadap log-odds individu bahagia, dibandingkan dengan kategori referensi (Sangat tidak puas), dengan usia dianggap konstan.

3 (b) Persamaan Regresi Logistik untuk Setiap Kategori Kepuasan

3.1 Pengambilan Koefisien

# Ekstrak koefisien
b0  <- coef(fit)["(Intercept)"]
b1  <- coef(fit)["usia"]
b_tp  <- coef(fit)["kepuasanTidak puas"]
b_ap  <- coef(fit)["kepuasanAgak puas"]
b_p   <- coef(fit)["kepuasanPuas"]
b_sp  <- coef(fit)["kepuasanSangat puas"]

cat("Intercept (β0)        :", round(b0, 4), "\n")

## Intercept (β0)        : 0.3267

cat("Usia (β1)             :", round(b1, 4), "\n")

## Usia (β1)             : 0.0024

cat("Tidak puas (β2)       :", round(b_tp, 4), "\n")

## Tidak puas (β2)       : 0.3071

cat("Agak puas (β3)        :", round(b_ap, 4), "\n")

## Agak puas (β3)        : 2.3072

cat("Puas (β4)             :", round(b_p, 4), "\n")

## Puas (β4)             : 2.9403

cat("Sangat puas (β5)      :", round(b_sp, 4), "\n")

## Sangat puas (β5)      : 2.7817

3.2 Persamaan per Kategori

Persamaan logistik dugaan untuk masing-masing kategori kepuasan diperoleh dengan memasukkan nilai indikator yang sesuai. Karena hanya satu \(c_j = 1\) pada tiap kategori, intersepnya bergeser sebesar \(\hat{\beta}_j\).

# Intersep masing-masing kategori kepuasan
int_stp <- b0               # Sangat tidak puas: c2=c3=c4=c5=0
int_tp  <- b0 + b_tp        # Tidak puas
int_ap  <- b0 + b_ap        # Agak puas
int_p   <- b0 + b_p         # Puas
int_sp  <- b0 + b_sp        # Sangat puas

cat("=== Persamaan Regresi Logistik Dugaan ===\n\n")

## === Persamaan Regresi Logistik Dugaan ===

cat(sprintf("Sangat tidak puas : logit P(Y=1) = %.4f + %.4f * usia\n", int_stp, b1))

## Sangat tidak puas : logit P(Y=1) = 0.3267 + 0.0024 * usia

cat(sprintf("Tidak puas        : logit P(Y=1) = %.4f + %.4f * usia\n", int_tp,  b1))

## Tidak puas        : logit P(Y=1) = 0.6337 + 0.0024 * usia

cat(sprintf("Agak puas         : logit P(Y=1) = %.4f + %.4f * usia\n", int_ap,  b1))

## Agak puas         : logit P(Y=1) = 2.6338 + 0.0024 * usia

cat(sprintf("Puas              : logit P(Y=1) = %.4f + %.4f * usia\n", int_p,   b1))

## Puas              : logit P(Y=1) = 3.2670 + 0.0024 * usia

cat(sprintf("Sangat puas       : logit P(Y=1) = %.4f + %.4f * usia\n", int_sp,  b1))

## Sangat puas       : logit P(Y=1) = 3.1084 + 0.0024 * usia

3.2.1 Interpretasi

Setiap persamaan di atas memiliki kemiringan yang sama (\(\hat{\beta}_1\) untuk usia) karena model tidak mengandung interaksi — artinya pengaruh usia terhadap log-odds bahagia identik di semua kategori kepuasan. Yang berbeda hanyalah intersep, yang mencerminkan perbedaan log-odds dasar antar kategori kepuasan pada nilai usia berapa pun.

4 (c) Plot Peluang Dugaan vs Usia untuk Semua Kategori Kepuasan

4.1 Formula Peluang

Peluang dugaan diperoleh dari transformasi logistik (sigmoid):

\[ \hat{P}(Y=1 \mid x, c_j) = \frac{e^{\hat{\beta}_0^{(j)} + \hat{\beta}_1 x}}{1 + e^{\hat{\beta}_0^{(j)} + \hat{\beta}_1 x}} \]

di mana \(\hat{\beta}_0^{(j)}\) adalah intersep untuk kategori kepuasan ke-\(j\).

# Membuat data prediksi
x_seq <- seq(20, 40, length.out = 200)

sigmoid <- function(eta) exp(eta) / (1 + exp(eta))

pred_df <- data.frame(
  usia      = rep(x_seq, 5),
  kepuasan  = rep(c("Sangat tidak puas","Tidak puas","Agak puas","Puas","Sangat puas"),
                  each = 200),
  pi_hat    = c(
    sigmoid(int_stp + b1 * x_seq),
    sigmoid(int_tp  + b1 * x_seq),
    sigmoid(int_ap  + b1 * x_seq),
    sigmoid(int_p   + b1 * x_seq),
    sigmoid(int_sp  + b1 * x_seq)
  )
)

pred_df$kepuasan <- factor(
  pred_df$kepuasan,
  levels = c("Sangat tidak puas","Tidak puas","Agak puas","Puas","Sangat puas")
)

ggplot(pred_df, aes(x = usia, y = pi_hat, color = kepuasan, linetype = kepuasan)) +
  geom_line(size = 1) +
  labs(
    title    = "Peluang Dugaan Individu Bahagia berdasarkan Usia",
    subtitle = "untuk Setiap Kategori Kepuasan",
    x        = "Usia (tahun)",
    y        = expression(hat(pi)),
    color    = "Kepuasan",
    linetype = "Kepuasan"
  ) +
  theme_bw() +
  theme(legend.position = "right") +
  scale_color_manual(values = c("red","orange","green","blue","purple"))

4.1.1 Interpretasi

Karena model tidak mengandung interaksi, kelima kurva paralel secara horizontal (tidak bersilangan). Individu dengan kepuasan Sangat puas memiliki peluang dugaan tertinggi untuk bahagia pada semua nilai usia, sementara individu dengan kepuasan Sangat tidak puas memiliki peluang terendah.

5 (d) Rasio Odds: “Sangat Puas” vs “Sangat Tidak Puas”

5.1 Formula Rasio Odds

Rasio odds antara dua kategori diperoleh dari selisih koefisien kedua kategori tersebut. Karena Sangat tidak puas adalah referensi (\(\hat{\beta} = 0\)), rasio odds-nya adalah:

\[ \text{OR}\left(\text{Sangat puas vs Sangat tidak puas}\right) = e^{\hat{\beta}_5 - 0} = e^{\hat{\beta}_5} \]

Lebih umum, rasio odds antara dua kelompok \(A\) dan \(B\) adalah:

\[ \text{OR}(A \text{ vs } B) = \frac{\text{odds}_A}{\text{odds}_B} = e^{\hat{\beta}_A - \hat{\beta}_B} \]

OR_sp_vs_stp <- exp(b_sp)  # Sangat puas vs Sangat tidak puas (referensi)

cat("Koefisien Sangat puas (β5)       :", round(b_sp, 4), "\n")

## Koefisien Sangat puas (β5)       : 2.7817

cat("Rasio Odds (Sangat puas / Sangat tidak puas):", round(OR_sp_vs_stp, 4), "\n")

## Rasio Odds (Sangat puas / Sangat tidak puas): 16.1467

# Confidence interval untuk OR
ci <- exp(confint(fit))
cat("\n95% CI untuk OR Sangat puas:\n")

## 
## 95% CI untuk OR Sangat puas:

print(ci["kepuasanSangat puas", ])

##    2.5 %   97.5 % 
## 13.17007 19.92058

5.1.1 Interpretasi

Pada usia yang sama, odds individu dengan kepuasan “Sangat puas” untuk bahagia adalah \(e^{\hat{\beta}_5}\) kali lebih besar dibandingkan individu dengan kepuasan “Sangat tidak puas”. Nilai OR > 1 menunjukkan bahwa tingkat kepuasan yang lebih tinggi berkaitan positif dengan kebahagiaan.

6 (e) Uji Signifikansi Model dengan Variabel Kategorik Kepuasan

6.1 Hipotesis

\[ H_0: \beta_2 = \beta_3 = \beta_4 = \beta_5 = 0 \quad \text{(kepuasan tidak berpengaruh signifikan)} \] \[ H_1: \exists\, \beta_c \neq 0, \quad c = 2, 3, 4, 5 \]

6.2 Statistik Uji Likelihood-Ratio (LR)

Statistik uji:

\[ -2 \log \Lambda = -2(\log L_0 - \log L_1) = D_0 - D_1 \sim \chi^2_{db} \]

di mana \(db = p_1 - p_0\) (selisih banyaknya parameter), dengan \(p_1 = 6\) (model lengkap) dan \(p_0 = 2\) (model tanpa kepuasan), sehingga \(db = 4\).

Kriteria keputusan: tolak \(H_0\) jika \(-2\log\Lambda > \chi^2_{0.05; 4} = 9.488\).

# Model tanpa kepuasan (model reduced)
fit.H0 <- glm(y ~ usia,           family = binomial, data = df)
# Model lengkap (dengan kepuasan)
fit.H1 <- glm(y ~ usia + kepuasan, family = binomial, data = df)

# Uji LR via anova
anova(fit.H0, fit.H1, test = "Chisq")

## Analysis of Deviance Table
## 
## Model 1: y ~ usia
## Model 2: y ~ usia + kepuasan
##   Resid. Df Resid. Dev Df Deviance  Pr(>Chi)    
## 1     78861      47337                          
## 2     78857      40095  4   7241.8 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Uji Type II menggunakan Anova() dari package car
Anova(fit.H1)

## Analysis of Deviance Table (Type II tests)
## 
## Response: y
##          LR Chisq Df Pr(>Chisq)    
## usia          1.1  1     0.2948    
## kepuasan   7241.8  4     <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

# Nilai kritis Chi-square df=4, alpha=0.05
chi_kritis <- qchisq(0.95, df = 4)
LR_stat    <- fit.H0$deviance - fit.H1$deviance
p_val      <- 1 - pchisq(LR_stat, df = 4)

cat("Statistik LR (-2 log Λ) :", round(LR_stat, 4), "\n")

## Statistik LR (-2 log Λ) : 7241.751

cat("Nilai kritis χ²(0.05;4) :", round(chi_kritis, 4), "\n")

## Nilai kritis χ²(0.05;4) : 9.4877

cat("p-value                 :", format(p_val, scientific = TRUE), "\n")

## p-value                 : 0e+00

6.2.1 Interpretasi

Karena \(-2\log\Lambda > \chi^2_{0.05;4}\) dan \(p\text{-value} < 0.05\), maka \(H_0\) ditolak. Terdapat cukup bukti untuk menyatakan bahwa variabel kepuasan berpengaruh signifikan terhadap peluang individu bahagia, setelah mengontrol usia.

7 (f) Variabel Biner “Puas”

Variabel biner puas didefinisikan sebagai:

\[ \text{puas} = \begin{cases} 1 & \text{jika kepuasan} \in \{\text{Puas, Sangat puas}\} \\ 0 & \text{jika kepuasan} \in \{\text{Agak puas, Tidak puas, Sangat tidak puas}\} \end{cases} \]

df$puas <- ifelse(df$kepuasan %in% c("Puas", "Sangat puas"), 1, 0)

# Distribusi variabel puas
cat("Distribusi variabel puas:\n")

## Distribusi variabel puas:

table(df$puas, dnn = "puas (1=Puas/Sangat puas, 0=lainnya)")

## puas (1=Puas/Sangat puas, 0=lainnya)
##     0     1 
## 44335 34528

cat("\nCross-tabulation kepuasan dan puas:\n")

## 
## Cross-tabulation kepuasan dan puas:

table(df$kepuasan, df$puas)

##                    
##                         0     1
##   Sangat tidak puas  1262     0
##   Tidak puas         9755     0
##   Agak puas         33318     0
##   Puas                  0 31194
##   Sangat puas           0  3334

8 (g) Model Interaksi: Usia, Puas, dan Usia × Puas

8.1 Model

Model dengan interaksi antara usia dan puas adalah:

\[ \text{logit}\,\hat{P}(Y=1) = \hat{\beta}_0 + \hat{\beta}_1 \cdot \text{usia} + \hat{\beta}_2 \cdot \text{puas} + \hat{\beta}_3 \cdot (\text{usia} \times \text{puas}) \]

fit.int   <- glm(y ~ usia + puas + usia:puas, family = binomial, data = df)
summary(fit.int)

## 
## Call:
## glm(formula = y ~ usia + puas + usia:puas, family = binomial, 
##     data = df)
## 
## Coefficients:
##               Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  1.7682118  0.0741658  23.841   <2e-16 ***
## usia         0.0040582  0.0024159   1.680    0.093 .  
## puas         1.4506999  0.1719624   8.436   <2e-16 ***
## usia:puas   -0.0005949  0.0056139  -0.106    0.916    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 47339  on 78862  degrees of freedom
## Residual deviance: 44880  on 78859  degrees of freedom
## AIC: 44888
## 
## Number of Fisher Scoring iterations: 6

8.2 Persamaan Dugaan per Kelompok

g0 <- coef(fit.int)["(Intercept)"]
g1 <- coef(fit.int)["usia"]
g2 <- coef(fit.int)["puas"]
g3 <- coef(fit.int)["usia:puas"]

cat("=== Koefisien Model Interaksi ===\n")

## === Koefisien Model Interaksi ===

cat(sprintf("Intercept (γ0)    : %.4f\n", g0))

## Intercept (γ0)    : 1.7682

cat(sprintf("Usia (γ1)         : %.4f\n", g1))

## Usia (γ1)         : 0.0041

cat(sprintf("Puas (γ2)         : %.4f\n", g2))

## Puas (γ2)         : 1.4507

cat(sprintf("Usia × Puas (γ3)  : %.4f\n", g3))

## Usia × Puas (γ3)  : -0.0006

cat("\n=== Persamaan Dugaan ===\n\n")

## 
## === Persamaan Dugaan ===

# Kelompok puas = 0
cat(sprintf("Kelompok TIDAK puas (puas=0):\n"))

## Kelompok TIDAK puas (puas=0):

cat(sprintf("  logit P(Y=1) = %.4f + %.4f * usia\n\n", g0, g1))

##   logit P(Y=1) = 1.7682 + 0.0041 * usia

# Kelompok puas = 1
int_puas <- g0 + g2
slope_puas <- g1 + g3
cat(sprintf("Kelompok PUAS (puas=1):\n"))

## Kelompok PUAS (puas=1):

cat(sprintf("  logit P(Y=1) = (%.4f + %.4f) + (%.4f + %.4f) * usia\n", g0, g2, g1, g3))

##   logit P(Y=1) = (1.7682 + 1.4507) + (0.0041 + -0.0006) * usia

cat(sprintf("               = %.4f + %.4f * usia\n\n", int_puas, slope_puas))

##                = 3.2189 + 0.0035 * usia

# Titik persilangan
x_cross <- (g0 - (g0 + g2)) / ((g1 + g3) - g1)
cat(sprintf("Titik persilangan dua kurva: usia = %.2f tahun\n", x_cross))

## Titik persilangan dua kurva: usia = 2438.65 tahun

8.2.1 Interpretasi

Pada model interaksi, pengaruh usia terhadap log-odds bahagia berbeda antara kelompok puas dan tidak puas:

Kelompok tidak puas (puas = 0): slope usia = \(\hat{\gamma}_1\).
Kelompok puas (puas = 1): slope usia = \(\hat{\gamma}_1 + \hat{\gamma}_3\).

Nilai \(\hat{\gamma}_3\) (koefisien interaksi) menunjukkan seberapa besar perbedaan pengaruh usia antara dua kelompok. Karena kedua kurva memiliki slope berbeda, mereka akan bersilangan pada nilai usia tertentu.

9 (h) Uji LR: Model Interaksi vs Model Tanpa Interaksi

9.1 Hipotesis

\[ H_0: \beta_3 = 0 \quad \text{(tidak ada interaksi antara usia dan puas)} \] \[ H_1: \beta_3 \neq 0 \quad \text{(ada interaksi)} \]

9.2 Statistik Uji

\[ -2\log\Lambda = D_{\text{tanpa interaksi}} - D_{\text{dengan interaksi}} \sim \chi^2_{1} \]

Kriteria keputusan: tolak \(H_0\) jika \(-2\log\Lambda > \chi^2_{0.05;1} = 3.841\).

fit.noint <- glm(y ~ usia + puas,            family = binomial, data = df)
fit.int   <- glm(y ~ usia + puas + usia:puas, family = binomial, data = df)

anova(fit.noint, fit.int, test = "LRT")

## Analysis of Deviance Table
## 
## Model 1: y ~ usia + puas
## Model 2: y ~ usia + puas + usia:puas
##   Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1     78860      44880                     
## 2     78859      44880  1 0.011229   0.9156

LR_int   <- fit.noint$deviance - fit.int$deviance
chi_1    <- qchisq(0.95, df = 1)
p_int    <- 1 - pchisq(LR_int, df = 1)

cat("Devians model tanpa interaksi :", round(fit.noint$deviance, 4), "\n")

## Devians model tanpa interaksi : 44879.84

cat("Devians model dengan interaksi:", round(fit.int$deviance, 4), "\n")

## Devians model dengan interaksi: 44879.83

cat("Beda devians (-2 log Λ)       :", round(LR_int, 4), "\n")

## Beda devians (-2 log Λ)       : 0.0112

cat("Nilai kritis χ²(0.05; 1)      :", round(chi_1, 4), "\n")

## Nilai kritis χ²(0.05; 1)      : 3.8415

cat("p-value                       :", round(p_int, 4), "\n")

## p-value                       : 0.9156

9.2.1 Kesimpulan Uji

if (LR_int > chi_1) {
  cat("Keputusan: TOLAK H0\n")
  cat("Kesimpulan: Ada bukti signifikan adanya interaksi antara usia dan kepuasan (puas).\n")
  cat("           Model DENGAN interaksi lebih baik digunakan.\n")
} else {
  cat("Keputusan: GAGAL TOLAK H0\n")
  cat("Kesimpulan: Tidak cukup bukti adanya interaksi antara usia dan kepuasan (puas).\n")
  cat("           Model TANPA interaksi lebih baik digunakan (lebih sederhana).\n")
}

## Keputusan: GAGAL TOLAK H0
## Kesimpulan: Tidak cukup bukti adanya interaksi antara usia dan kepuasan (puas).
##            Model TANPA interaksi lebih baik digunakan (lebih sederhana).

9.2.2 Interpretasi Akhir

Jika \(p\text{-value} > 0.05\), maka interaksi antara usia dan puas tidak signifikan secara statistik. Model yang lebih sederhana (tanpa interaksi) lebih direkomendasikan karena memiliki interpretasi yang lebih mudah dengan kecocokan data yang sebanding. Sebaliknya, jika \(p\text{-value} \leq 0.05\), maka pengaruh usia terhadap kebahagiaan berbeda secara signifikan antara individu yang puas dan tidak puas.

Analisis Data Kategorik — Program Studi Statistika, Universitas Sebelas Maret

Regresi Logistik Berganda

Analisis Data Kategorik — Universitas Sebelas Maret

Adyannisa

18 May 2026