EVALUASI DAMPAK UKURAN SAMPEL DAN MULTIKOLINEARITAS TERHADAP ESTIMASI PARAMETER KEMENANGAN RRQ HOSHI MELALUI PENDEKATAN REGRESI LOGISTIK

library(MASS)

## Warning: package 'MASS' was built under R version 4.5.2

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following object is masked from 'package:MASS':
## 
##     select

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(ggplot2)
library(tidyr)

## Warning: package 'tidyr' was built under R version 4.5.2

# 1. Konfigurasi Awal
set.seed(040)
n_seq <- c(35, 100, 500) 
rho_seq <- c(0.0, 0.85) 
replikasi <- 100
beta_true <- c(-3.0, 1.2, 0.8) # Intercept, Gold (b1), Objective (b2)

hasil_simulasi <- data.frame()

# 2. Loop Simulasi
for (n in n_seq) {
  for (rho in rho_seq) {
    
    m_bias_b1 <- m_se_b1 <- m_rse_b1 <- numeric(replikasi)
    m_bias_b2 <- m_se_b2 <- m_rse_b2 <- numeric(replikasi)
    
    for (i in 1:replikasi) {
      Sigma <- matrix(c(1, rho, rho, 1), 2, 2)
      X <- mvrnorm(n, mu = c(2, 2), Sigma = Sigma)
      X1 <- X[, 1]; X2 <- X[, 2]
      
      log_odds <- beta_true[1] + beta_true[2]*X1 + beta_true[3]*X2
      prob <- 1 / (1 + exp(-log_odds))
      Y <- rbinom(n, 1, prob)
      
      data_sim <- data.frame(Y, X1, X2)
      
      model <- tryCatch({
        glm(Y ~ X1 + X2, data = data_sim, family = binomial)
      }, warning = function(w) NULL, error = function(e) NULL)
      
      if (!is.null(model) && model$converged) {
        coef_est <- coef(model)
        se_est <- summary(model)$coefficients[, "Std. Error"]
        
        m_bias_b1[i] <- coef_est[2] - beta_true[2]
        m_se_b1[i]   <- se_est[2]
        
        m_bias_b2[i] <- coef_est[3] - beta_true[3]
        m_se_b2[i]   <- se_est[3]
      } else {
        m_bias_b1[i] <- m_se_b1[i] <- NA
        m_bias_b2[i] <- m_se_b2[i] <- NA
      }
    }
    
    # Simpan hasil Gold
    hasil_simulasi <- rbind(hasil_simulasi, data.frame(
      n = n,
      rho = as.factor(rho),
      Variable = "Gold",
      Mean_Bias = mean(m_bias_b1, na.rm = TRUE),
      Mean_SE = mean(m_se_b1, na.rm = TRUE)
    ))
    
    # Simpan hasil Objective
    hasil_simulasi <- rbind(hasil_simulasi, data.frame(
      n = n,
      rho = as.factor(rho),
      Variable = "Objective",
      Mean_Bias = mean(m_bias_b2, na.rm = TRUE),
      Mean_SE = mean(m_se_b2, na.rm = TRUE)
    ))
  }
}

# 3. Visualisasi

# Grafik 1: Tren Bias 
plot_bias_combined <- ggplot(hasil_simulasi, aes(x = n, y = Mean_Bias, color = rho, group = rho)) +
  geom_line(linewidth = 1.2) +
  geom_point(size = 3) +
  geom_hline(yintercept = 0, linetype = "dashed", color = "red") +
  facet_wrap(~ Variable) + 
  labs(x = "Ukuran Sampel (n)",
       y = "Mean Bias",
       color = "Korelasi (Rho)") +
  theme_minimal() +
  theme(legend.position = "bottom")

# Grafik 2: Tren Standard Error 
plot_se_combined <- ggplot(hasil_simulasi, aes(x = as.factor(n), y = Mean_SE, fill = rho)) +
  geom_bar(stat = "identity", position = "dodge") +
  facet_wrap(~ Variable) +
  labs(x = "Ukuran Sampel (n)",
       y = "Mean Standard Error",
       fill = "Korelasi (Rho)") +
  theme_minimal() +
  theme(legend.position = "bottom")

# 4. Cetak Tabel Hasil dan plot Visualisasi
print(hasil_simulasi)

##      n  rho  Variable     Mean_Bias   Mean_SE
## 1   35    0      Gold  2.552582e-01 0.6777105
## 2   35    0 Objective  1.768878e-01 0.5900855
## 3   35 0.85      Gold  1.692321e-01 1.1068881
## 4   35 0.85 Objective  1.729370e-01 1.0545515
## 5  100    0      Gold  1.006459e-01 0.3364469
## 6  100    0 Objective  7.423939e-02 0.2982835
## 7  100 0.85      Gold  1.240199e-01 0.5711522
## 8  100 0.85 Objective -1.014450e-02 0.5452383
## 9  500    0      Gold -6.871417e-05 0.1390296
## 10 500    0 Objective  1.088872e-02 0.1261068
## 11 500 0.85      Gold  3.379039e-03 0.2440246
## 12 500 0.85 Objective  2.144445e-02 0.2354799

print(plot_bias_combined)

print(plot_se_combined)