# Memanggil dataset esoph yang sudah otomatis ada di base R
data(esoph)
# Melihat struktur dataset
str(esoph)
## 'data.frame': 88 obs. of 5 variables:
## $ agegp : Ord.factor w/ 6 levels "25-34"<"35-44"<..: 1 1 1 1 1 1 1 1 1 1 ...
## $ alcgp : Ord.factor w/ 4 levels "0-39g/day"<"40-79"<..: 1 1 1 1 2 2 2 2 3 3 ...
## $ tobgp : Ord.factor w/ 4 levels "0-9g/day"<"10-19"<..: 1 2 3 4 1 2 3 4 1 2 ...
## $ ncases : num 0 0 0 0 0 0 0 0 0 0 ...
## $ ncontrols: num 40 10 6 5 27 7 4 7 2 1 ...
head(esoph)
## agegp alcgp tobgp ncases ncontrols
## 1 25-34 0-39g/day 0-9g/day 0 40
## 2 25-34 0-39g/day 10-19 0 10
## 3 25-34 0-39g/day 20-29 0 6
## 4 25-34 0-39g/day 30+ 0 5
## 5 25-34 40-79 0-9g/day 0 27
## 6 25-34 40-79 10-19 0 7
#level masing-masing variabel
cat("Level agegp (Kelompok Usia):\n")
## Level agegp (Kelompok Usia):
print(levels(esoph$agegp))
## [1] "25-34" "35-44" "45-54" "55-64" "65-74" "75+"
cat("\nLevel alcgp (Kelompok Konsumsi Alkohol):\n")
##
## Level alcgp (Kelompok Konsumsi Alkohol):
print(levels(esoph$alcgp))
## [1] "0-39g/day" "40-79" "80-119" "120+"
cat("\nLevel tobgp (Kelompok Konsumsi Tembakau):\n")
##
## Level tobgp (Kelompok Konsumsi Tembakau):
print(levels(esoph$tobgp))
## [1] "0-9g/day" "10-19" "20-29" "30+"
# Membuat proporsi langsung di dalam esoph
esoph$Total_People <- esoph$ncases + esoph$ncontrols
esoph$Prop_Cases <- esoph$ncases / esoph$Total_People
# Pemodelan dengan Logit
esoph$agegp <- factor(esoph$agegp, ordered = TRUE)
esoph$alcgp <- factor(esoph$alcgp, ordered = TRUE)
esoph$tobgp <- factor(esoph$tobgp, ordered = TRUE)
options(contrasts = c("contr.treatment", "contr.treatment"))
model_esoph1 <- glm(
Prop_Cases ~ agegp + alcgp + tobgp,
family = binomial(link = "logit"), data = esoph, weights = Total_People
)
summary(model_esoph1)
##
## Call:
## glm(formula = Prop_Cases ~ agegp + alcgp + tobgp, family = binomial(link = "logit"),
## data = esoph, weights = Total_People)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.8954 1.0859 -6.350 2.16e-10 ***
## agegp35-44 1.9809 1.1041 1.794 0.072786 .
## agegp45-54 3.7763 1.0680 3.536 0.000407 ***
## agegp55-64 4.3352 1.0650 4.070 4.69e-05 ***
## agegp65-74 4.8964 1.0764 4.549 5.39e-06 ***
## agegp75+ 4.8265 1.1213 4.304 1.67e-05 ***
## alcgp40-79 1.4346 0.2501 5.737 9.63e-09 ***
## alcgp80-119 1.9807 0.2848 6.956 3.51e-12 ***
## alcgp120+ 3.6029 0.3850 9.357 < 2e-16 ***
## tobgp10-19 0.4381 0.2283 1.919 0.055039 .
## tobgp20-29 0.5126 0.2730 1.878 0.060398 .
## tobgp30+ 1.6410 0.3441 4.769 1.85e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 367.953 on 87 degrees of freedom
## Residual deviance: 82.337 on 76 degrees of freedom
## AIC: 221.39
##
## Number of Fisher Scoring iterations: 6
model_esoph2 <- glm(
cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp,
data = esoph,
family = binomial(link = "logit")
)
summary(model_esoph2)
##
## Call:
## glm(formula = cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp,
## family = binomial(link = "logit"), data = esoph)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -6.8954 1.0859 -6.350 2.16e-10 ***
## agegp35-44 1.9809 1.1041 1.794 0.072786 .
## agegp45-54 3.7763 1.0680 3.536 0.000407 ***
## agegp55-64 4.3352 1.0650 4.070 4.69e-05 ***
## agegp65-74 4.8964 1.0764 4.549 5.39e-06 ***
## agegp75+ 4.8265 1.1213 4.304 1.67e-05 ***
## alcgp40-79 1.4346 0.2501 5.737 9.63e-09 ***
## alcgp80-119 1.9807 0.2848 6.956 3.51e-12 ***
## alcgp120+ 3.6029 0.3850 9.357 < 2e-16 ***
## tobgp10-19 0.4381 0.2283 1.919 0.055039 .
## tobgp20-29 0.5126 0.2730 1.878 0.060398 .
## tobgp30+ 1.6410 0.3441 4.769 1.85e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 367.953 on 87 degrees of freedom
## Residual deviance: 82.337 on 76 degrees of freedom
## AIC: 221.39
##
## Number of Fisher Scoring iterations: 6
coef(model_esoph1) ; coef(model_esoph2)
## (Intercept) agegp35-44 agegp45-54 agegp55-64 agegp65-74 agegp75+
## -6.8954152 1.9808846 3.7762865 4.3351817 4.8964059 4.8265420
## alcgp40-79 alcgp80-119 alcgp120+ tobgp10-19 tobgp20-29 tobgp30+
## 1.4346287 1.9807173 3.6028688 0.4380525 0.5126181 1.6409973
## (Intercept) agegp35-44 agegp45-54 agegp55-64 agegp65-74 agegp75+
## -6.8954152 1.9808846 3.7762865 4.3351817 4.8964059 4.8265420
## alcgp40-79 alcgp80-119 alcgp120+ tobgp10-19 tobgp20-29 tobgp30+
## 1.4346287 1.9807173 3.6028688 0.4380525 0.5126181 1.6409973
cat("odds ratio adalah: ")
## odds ratio adalah:
oddsratio <- round(exp(coef(model_esoph1)), 4)
oddsratio
## (Intercept) agegp35-44 agegp45-54 agegp55-64 agegp65-74 agegp75+
## 0.0010 7.2492 43.6536 76.3388 133.8080 124.7787
## alcgp40-79 alcgp80-119 alcgp120+ tobgp10-19 tobgp20-29 tobgp30+
## 4.1981 7.2479 36.7034 1.5497 1.6697 5.1603
options(contrasts = c("contr.poly", "contr.poly"))
# Simultan test
# Model null (tanpa prediktor)
model_null <- glm(
cbind(ncases, ncontrols) ~ 1,
data = esoph,
family = binomial(link = "logit")
)
model_full <- model_esoph2
# Likelihood Ratio Test
LRTest <- anova(model_null, model_full, test = "Chisq")
LRTest
## Analysis of Deviance Table
##
## Model 1: cbind(ncases, ncontrols) ~ 1
## Model 2: cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp
## Resid. Df Resid. Dev Df Deviance Pr(>Chi)
## 1 87 367.95
## 2 76 82.34 11 285.62 < 2.2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# Ambil p-value otomatis
p_LRTest <- LRTest$`Pr(>Chi)`[2]
# Keputusan
if (p_LRTest < 0.05) {
cat("\n Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan. \n")
} else {
cat("\n Gagal Tolak H_0 karena nilai p-value > alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, s tidak ada pengaruh simultan dari seluruh variabel prediktor \n")
}
##
## Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan.
#LR Test cara manual
# Ambil nilai log-likelihood
LL_null <- as.numeric(logLik(model_null))
LL_full <- as.numeric(logLik(model_full))
# Hitung Likelihood Ratio Statistic
LR_Test1 <- -2 * (LL_null - LL_full)
# Hitung derajat bebas (df)
df <- attr(logLik(model_full), "df") - attr(logLik(model_null), "df")
# Hitung p-value
p_LRTest1 <- 1 - pchisq(LR_Test1, df)
# Tampilkan hasil
cat("Pengujian LR Test secara manual \n")
## Pengujian LR Test secara manual
cat("LogLik Null =", LL_null, "\n")
## LogLik Null = -241.5042
cat("LogLik Full =", LL_full, "\n")
## LogLik Full = -98.6959
cat("LR Statistic =", LR_Test1, "\n")
## LR Statistic = 285.6166
cat("df =", df, "\n")
## df = 11
cat("p-value =", p_LRTest1, "\n")
## p-value = 0
# Keputusan
if (p_LRTest1 < 0.05) {
cat("\n Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan. \n")
} else {
cat("\n Gagal Tolak H_0 karena nilai p-value > alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, s tidak ada pengaruh simultan dari seluruh variabel prediktor \n")
}
##
## Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan.
#Uji Parsial
ringkasan <- summary(model_esoph2)
wald <- ringkasan$coefficients
# Buatkan Data Frame Uji Parsial Wald Test
wald_df <- as.data.frame(wald)
# Hitung Odds Ratio
wald_df$Odds_Ratio <- exp(wald_df$Estimate)
# Merapihkan angka di belakang koma
wald_df$Estimate <- round(wald_df$Estimate, 4)
wald_df$Std_Error <- round(wald_df$`Std. Error`, 4)
wald_df$z_value <- round(wald_df$`z value`, 4)
wald_df$p_value <- round(wald_df$`Pr(>|z|)`, 4)
wald_df$Odds_Ratio <- round(wald_df$Odds_Ratio, 4)
# Tambah kolom keputusan otomatis
wald_df$Decision <- ifelse(wald_df$`Pr(>|z|)` < 0.05,
"Signifikan",
"Tidak signifikan")
# Rename kolom agar lebih rapi
wald_df <- wald_df[, c("Estimate","Std_Error","z_value","p_value","Odds_Ratio","Decision")]
# Cetak tabel
cat("\n Hasil Pengujian Parsial dengan Wald Test \n")
##
## Hasil Pengujian Parsial dengan Wald Test
print(wald_df)
## Estimate Std_Error z_value p_value Odds_Ratio Decision
## (Intercept) -6.8954 1.0859 -6.3497 0.0000 0.0010 Signifikan
## agegp35-44 1.9809 1.1041 1.7942 0.0728 7.2492 Tidak signifikan
## agegp45-54 3.7763 1.0680 3.5357 0.0004 43.6536 Signifikan
## agegp55-64 4.3352 1.0650 4.0704 0.0000 76.3388 Signifikan
## agegp65-74 4.8964 1.0764 4.5490 0.0000 133.8080 Signifikan
## agegp75+ 4.8265 1.1213 4.3044 0.0000 124.7787 Signifikan
## alcgp40-79 1.4346 0.2501 5.7371 0.0000 4.1981 Signifikan
## alcgp80-119 1.9807 0.2848 6.9557 0.0000 7.2479 Signifikan
## alcgp120+ 3.6029 0.3850 9.3572 0.0000 36.7034 Signifikan
## tobgp10-19 0.4381 0.2283 1.9186 0.0550 1.5497 Tidak signifikan
## tobgp20-29 0.5126 0.2730 1.8779 0.0604 1.6697 Tidak signifikan
## tobgp30+ 1.6410 0.3441 4.7688 0.0000 5.1603 Signifikan
#Uji Devians
# Ambil deviance dan df
dev_val <- model_esoph2$deviance
df_val <- model_esoph2$df.residual
# Hitung p-value
p_val_dev <- 1 - pchisq(dev_val, df_val)
dev_val_r <- round(dev_val, 4)
df_val_r <- round(df_val, 4)
p_val_r <- round(p_val_dev, 4)
# Keputusan
decision <- ifelse(p_val_dev > 0.05,
"Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit",
"Gagal Tolak H_0, artinya dengan tingkat signifikansi 0.05 tidak terdapat cukup bukti unutk menyatakan bahwa Model Fit")
# Cetak hasil
cat("\n Goodness of Fit Test dengan Uji Devians \n")
##
## Goodness of Fit Test dengan Uji Devians
cat("Residual Deviance : ", dev_val_r, "\n")
## Residual Deviance : 82.3369
cat("Residual DF : ", df_val_r, "\n")
## Residual DF : 76
cat("p-value : ", p_val_r, "\n")
## p-value : 0.2898
cat("Keputusan : ", decision, "\n")
## Keputusan : Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit
#Model Probit
options(contrasts = c("contr.treatment", "contr.treatment"))
model_probit <- glm(
cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp,
family = binomial(link = "probit"),
data = esoph
)
options(contrasts = c("contr.poly", "contr.poly"))
# Ringkasan model
summary(model_probit)
##
## Call:
## glm(formula = cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp,
## family = binomial(link = "probit"), data = esoph)
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.7990 0.5251 -7.235 4.65e-13 ***
## agegp35-44 1.0343 0.5300 1.952 0.050994 .
## agegp45-54 1.9677 0.5146 3.824 0.000131 ***
## agegp55-64 2.3020 0.5125 4.492 7.06e-06 ***
## agegp65-74 2.6295 0.5207 5.050 4.41e-07 ***
## agegp75+ 2.5850 0.5515 4.688 2.76e-06 ***
## alcgp40-79 0.8110 0.1362 5.953 2.64e-09 ***
## alcgp80-119 1.1259 0.1596 7.056 1.72e-12 ***
## alcgp120+ 2.0762 0.2112 9.833 < 2e-16 ***
## tobgp10-19 0.2935 0.1302 2.254 0.024203 *
## tobgp20-29 0.3146 0.1574 1.999 0.045610 *
## tobgp30+ 0.9348 0.1966 4.754 1.99e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 367.953 on 87 degrees of freedom
## Residual deviance: 80.562 on 76 degrees of freedom
## AIC: 219.62
##
## Number of Fisher Scoring iterations: 6
library(margins)
# Menghitung efek marginal
margins_probit <- margins(model_probit)
# Menampilkan hasil efek marginal
summary(margins_probit)
## factor AME SE z p lower upper
## agegp35-44 0.1077 0.0403 2.6755 0.0075 0.0288 0.1866
## agegp45-54 0.3224 0.0412 7.8304 0.0000 0.2417 0.4030
## agegp55-64 0.4208 0.0394 10.6922 0.0000 0.3437 0.4979
## agegp65-74 0.5209 0.0445 11.6981 0.0000 0.4336 0.6082
## agegp75+ 0.5073 0.0729 6.9570 0.0000 0.3644 0.6502
## alcgp120+ 0.5153 0.0483 10.6680 0.0000 0.4207 0.6100
## alcgp40-79 0.1686 0.0270 6.2345 0.0000 0.1156 0.2215
## alcgp80-119 0.2536 0.0370 6.8551 0.0000 0.1811 0.3261
## tobgp10-19 0.0656 0.0295 2.2283 0.0259 0.0079 0.1233
## tobgp20-29 0.0706 0.0362 1.9505 0.0511 -0.0003 0.1415
## tobgp30+ 0.2240 0.0477 4.6943 0.0000 0.1305 0.3176
#Uji Simultan
# Misal kita punya dua model:
# model_full: model dengan semua prediktor
# model_reduced: model dengan prediktor dikurangi
model_full <- glm(cbind(ncases, ncontrols) ~ agegp + alcgp + tobgp, family = binomial(link="probit"), data=esoph)
model_reduced <- glm(cbind(ncases, ncontrols) ~ 1, family = binomial(link="probit"), data=esoph)
# Hitung deviance
deviance_full <- deviance(model_full)
deviance_reduced <- deviance(model_reduced)
# Hitung statistik uji G (likelihood ratio)
G_statistic <- deviance_reduced - deviance_full
# Hitung degrees of freedom
df <- df.residual(model_reduced) - df.residual(model_full)
# Hitung p-value
p_value <- pchisq(G_statistic, df, lower.tail = FALSE)
# Tampilkan hasil
cat("G-Statistic:", G_statistic, "\n")
## G-Statistic: 287.3911
cat("Degrees of Freedom:", df, "\n")
## Degrees of Freedom: 11
cat("P-value:", p_value, "\n")
## P-value: 3.955173e-55
# Keputusan
if (p_value < 0.05) {
cat("\n Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan. \n")
} else {
cat("\n Gagal Tolak H_0 karena nilai p-value > alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, s tidak ada pengaruh simultan dari seluruh variabel prediktor \n")
}
##
## Tolak H_0 karena nilai p-value < alpha = 0.05. Artinya, dengan tingkat signifikansi 0.05, seluruh variabel prediktornya berpengaruh secara simultan.
#Uji Parsial
# Fungsi uji Wald otomatis untuk semua koefisien
wald_test_all <- function(model) {
coefs <- coef(summary(model))
wald_results <- data.frame(
Variable = rownames(coefs),
Estimate = coefs[, "Estimate"],
Std_Error = coefs[, "Std. Error"],
Wald_Statistic = NA,
P_value = NA
)
for (i in 1:nrow(wald_results)) {
wald_stat <- (wald_results$Estimate[i] / wald_results$Std_Error[i])^2
p_val <- pchisq(wald_stat, df = 1, lower.tail = FALSE)
wald_results$Wald_Statistic[i] <- wald_stat
wald_results$P_value[i] <- round(p_val, 4)
}
return(wald_results)
}
# Hitung uji Wald untuk semua koefisien
wald_results <- wald_test_all(model_probit)
# Tampilkan hasil
print(wald_results)
## Variable Estimate Std_Error Wald_Statistic P_value
## (Intercept) (Intercept) -3.7990494 0.5250758 52.348690 0.0000
## agegp35-44 agegp35-44 1.0342740 0.5299819 3.808458 0.0510
## agegp45-54 agegp45-54 1.9677446 0.5145526 14.624398 0.0001
## agegp55-64 agegp55-64 2.3020214 0.5125009 20.175740 0.0000
## agegp65-74 agegp65-74 2.6295301 0.5206645 25.505885 0.0000
## agegp75+ agegp75+ 2.5850246 0.5514662 21.973113 0.0000
## alcgp40-79 alcgp40-79 0.8109709 0.1362335 35.435823 0.0000
## alcgp80-119 alcgp80-119 1.1259036 0.1595717 49.784117 0.0000
## alcgp120+ alcgp120+ 2.0761633 0.2111521 96.679019 0.0000
## tobgp10-19 tobgp10-19 0.2935003 0.1302198 5.079991 0.0242
## tobgp20-29 tobgp20-29 0.3146132 0.1573861 3.995957 0.0456
## tobgp30+ tobgp30+ 0.9347740 0.1966232 22.601833 0.0000
library(caret)
## Warning: package 'caret' was built under R version 4.4.1
## Loading required package: ggplot2
## Loading required package: lattice
# Model proporsional sudah dibuat: model_probit
# Buat variabel Y biner dari dataset asli
esoph$Y_binary <- ifelse(esoph$ncases > 0, 1, 0)
# 2a. Prediksi Probabilitas (proporsi kasus)
prob_probit <- predict(model_probit, type = "response")
# 2b. Konversi probabilitas ke kelas biner (threshold 0.1)
pred_class_probit <- ifelse(prob_probit >= 0.1, 1, 0)
# 2c. Ubah kelas prediksi menjadi faktor
pred_class_probit_factor <- factor(pred_class_probit, levels = c(0, 1))
# 2d. Outcome asli sebagai faktor
Outcome_factor <- factor(esoph$Y_binary, levels = c(0, 1))
# 2e. Buat Confusion Matrix
conf_matrix_probit <- confusionMatrix(pred_class_probit_factor, Outcome_factor)
# 2f. Tampilkan hasil
print(conf_matrix_probit)
## Confusion Matrix and Statistics
##
## Reference
## Prediction 0 1
## 0 22 6
## 1 7 53
##
## Accuracy : 0.8523
## 95% CI : (0.7606, 0.9189)
## No Information Rate : 0.6705
## P-Value [Acc > NIR] : 9.252e-05
##
## Kappa : 0.6627
##
## Mcnemar's Test P-Value : 1
##
## Sensitivity : 0.7586
## Specificity : 0.8983
## Pos Pred Value : 0.7857
## Neg Pred Value : 0.8833
## Prevalence : 0.3295
## Detection Rate : 0.2500
## Detection Prevalence : 0.3182
## Balanced Accuracy : 0.8285
##
## 'Positive' Class : 0
##
# Model proporsional sudah dibuat: model_probit
# deviance residual dan null deviance sudah tersedia di model
# 1. Hitung deviance residual
deviance_model <- deviance(model_probit) # Deviance model
df_model <- df.residual(model_probit) # Degrees of freedom residual
# 2. Null deviance (model hanya dengan intercept)
null_deviance <- model_probit$null.deviance
df_null <- model_probit$df.null
# 3. Tampilkan deviance dan df
cat("Null Deviance:", null_deviance, "on", df_null, "degrees of freedom\n")
## Null Deviance: 367.9535 on 87 degrees of freedom
cat("Residual Deviance:", deviance_model, "on", df_model, "degrees of freedom\n")
## Residual Deviance: 80.56233 on 76 degrees of freedom
# 4. Uji deviance (chi-square) untuk goodness-of-fit
p_valuedeviance <- round(pchisq(deviance_model, df = df_model, lower.tail = FALSE),4)
cat("P-value untuk uji deviance:", p_valuedeviance, "\n")
## P-value untuk uji deviance: 0.3384
ifelse( p_valuedeviance > 0.05,
"Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit",
"Gagal Tolak H_0, artinya dengan tingkat signifikansi 0.05 tidak terdapat cukup bukti unutk menyatakan bahwa Model Fit")
## [1] "Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit"
library(ResourceSelection)
## ResourceSelection 0.3-6 2023-06-27
# Model probit biner
# Outcome asli sebagai vektor
y <- esoph$Y_binary
# Prediksi probabilitas dari model
model_probit_bin <- glm(Y_binary ~ agegp + alcgp + tobgp, family = binomial(link="probit"), data=esoph)
## Warning: glm.fit: fitted probabilities numerically 0 or 1 occurred
prob <- predict(model_probit_bin, type = "response")
# Lakukan Hosmer-Lemeshow test
hl_test <- hoslem.test(y, prob, g = 10) # g = jumlah grup decile (default 10)
# Tampilkan hasil
hl_test
##
## Hosmer and Lemeshow goodness of fit (GOF) test
##
## data: y, prob
## X-squared = 2.0782, df = 8, p-value = 0.9785
ifelse(hl_test$p.value > 0.05,
"Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit",
"Gagal Tolak H_0, artinya dengan tingkat signifikansi 0.05 tidak terdapat cukup bukti unutk menyatakan bahwa Model Fit")
## [1] "Tolak H_0, artinya dengan tingkat signifikansi 0.05 terdapat cukup bukti untuk menyatakan bahwa Model Fit"