## ============================================================
## 2. Design Matrix per Strata (tanpa Pemisah)
## ============================================================
X_jalan <- model.matrix(~ X1 + X2 + X3 + X4 + X5, data = data_jalan)
X_inap <- model.matrix(~ X1 + X2 + X3 + X4 + X5, data = data_inap)
## Fungsi loglik sama:
loglik_geo <- function(beta, X, Y){
eta <- as.vector(X %*% beta)
h <- 1 / (1 + exp(-eta))
ll <- sum(log(h) + (Y - 1) * log(1 - h))
return(-ll)
}
## Fungsi inverse aman sama:
if (!"MASS" %in% rownames(installed.packages())) install.packages("MASS")
library(MASS)
safe_inverse <- function(H){
tryCatch(
solve(H),
error = function(e){
message("Hessian singular → gunakan generalized inverse (ginv).")
MASS::ginv(H)
}
)
}
## ============================================================
## 3. MODEL KHUSUS RAWAT JALAN (Pemisah = 0)
## ============================================================
start_jalan <- rep(0, ncol(X_jalan))
fit_jalan <- optim(par = start_jalan,
fn = loglik_geo,
X = X_jalan,
Y = Y_jalan,
hessian = TRUE)
coef_jalan <- fit_jalan$par
names(coef_jalan) <- colnames(X_jalan)
vcov_jalan <- safe_inverse(fit_jalan$hessian)
## Hessian singular → gunakan generalized inverse (ginv).
se_jalan <- sqrt(diag(vcov_jalan))
z_jalan <- coef_jalan / se_jalan
p_jalan <- 2 * pnorm(-abs(z_jalan))
HR_jalan <- exp(coef_jalan)
lowerCI_jalan <- exp(coef_jalan - 1.96 * se_jalan)
upperCI_jalan <- exp(coef_jalan + 1.96 * se_jalan)
hasil_jalan <- data.frame(
Parameter = names(coef_jalan),
Estimate = coef_jalan,
Std.Error = se_jalan,
z.value = z_jalan,
p.value = p_jalan,
HR = HR_jalan,
CI.lower = lowerCI_jalan,
CI.upper = upperCI_jalan
)
cat("\n======================================\n")
##
## ======================================
cat("MODEL KHUSUS RAWAT JALAN (Pemisah=0)\n")
## MODEL KHUSUS RAWAT JALAN (Pemisah=0)
cat("======================================\n")
## ======================================
print(hasil_jalan, digits = 4)
## Parameter Estimate Std.Error z.value p.value HR CI.lower
## (Intercept) (Intercept) 3.27042 1.885e+00 1.735e+00 0.08280 26.32238 0.65387
## X1 X1 0.01004 3.329e-02 3.016e-01 0.76299 1.01009 0.94628
## X2 X2 0.88133 1.903e-16 4.632e+15 0.00000 2.41412 2.41412
## X31 X31 0.21820 7.329e-01 2.977e-01 0.76592 1.24383 0.29573
## X41 X41 -2.34502 1.073e+00 -2.185e+00 0.02891 0.09585 0.01169
## X51 X51 -0.91961 7.417e-01 -1.240e+00 0.21505 0.39868 0.09316
## CI.upper
## (Intercept) 1059.6347
## X1 1.0782
## X2 2.4141
## X31 5.2316
## X41 0.7857
## X51 1.7061
LL_jalan <- -fit_jalan$value
k_jalan <- length(coef_jalan)
n_jalan <- nrow(data_jalan)
AIC_jalan <- -2*LL_jalan + 2*k_jalan
BIC_jalan <- -2*LL_jalan + k_jalan*log(n_jalan)
cat("\n=== VALIDASI MODEL RAWAT JALAN ===\n")
##
## === VALIDASI MODEL RAWAT JALAN ===
cat("Log-likelihood:", LL_jalan, "\n")
## Log-likelihood: -43.44246
cat("AIC:", AIC_jalan, "\n")
## AIC: 98.88491
cat("BIC:", BIC_jalan, "\n")
## BIC: 114.5756
## ============================================================
## 4. MODEL KHUSUS RAWAT INAP (Pemisah = 1)
## ============================================================
start_inap <- rep(0, ncol(X_inap))
fit_inap <- optim(par = start_inap,
fn = loglik_geo,
X = X_inap,
Y = Y_inap,
hessian = TRUE)
coef_inap <- fit_inap$par
names(coef_inap) <- colnames(X_inap)
vcov_inap <- safe_inverse(fit_inap$hessian)
se_inap <- sqrt(diag(vcov_inap))
z_inap <- coef_inap / se_inap
p_inap <- 2 * pnorm(-abs(z_inap))
HR_inap <- exp(coef_inap)
lowerCI_inap <- exp(coef_inap - 1.96 * se_inap)
upperCI_inap <- exp(coef_inap + 1.96 * se_inap)
hasil_inap <- data.frame(
Parameter = names(coef_inap),
Estimate = coef_inap,
Std.Error = se_inap,
z.value = z_inap,
p.value = p_inap,
HR = HR_inap,
CI.lower = lowerCI_inap,
CI.upper = upperCI_inap
)
cat("\n======================================\n")
##
## ======================================
cat("MODEL KHUSUS RAWAT INAP (Pemisah=1)\n")
## MODEL KHUSUS RAWAT INAP (Pemisah=1)
cat("======================================\n")
## ======================================
print(hasil_inap, digits = 4)
## Parameter Estimate Std.Error z.value p.value HR CI.lower
## (Intercept) (Intercept) 0.9998 0.77292 1.2935 1.958e-01 2.7177 0.5974
## X1 X1 -0.0110 0.01343 -0.8191 4.127e-01 0.9891 0.9634
## X2 X2 -0.4060 0.08050 -5.0429 4.585e-07 0.6663 0.5691
## X31 X31 -0.2871 0.33682 -0.8522 3.941e-01 0.7505 0.3878
## X41 X41 -0.6127 0.27255 -2.2481 2.457e-02 0.5419 0.3176
## X51 X51 0.2548 0.34850 0.7312 4.646e-01 1.2902 0.6517
## CI.upper
## (Intercept) 12.3632
## X1 1.0154
## X2 0.7802
## X31 1.4522
## X41 0.9245
## X51 2.5546
LL_inap <- -fit_inap$value
k_inap <- length(coef_inap)
n_inap <- nrow(data_inap)
AIC_inap <- -2*LL_inap + 2*k_inap
BIC_inap <- -2*LL_inap + k_inap*log(n_inap)
cat("\n=== VALIDASI MODEL RAWAT INAP ===\n")
##
## === VALIDASI MODEL RAWAT INAP ===
cat("Log-likelihood:", LL_inap, "\n")
## Log-likelihood: -212.9369
cat("AIC:", AIC_inap, "\n")
## AIC: 437.8737
cat("BIC:", BIC_inap, "\n")
## BIC: 453.322