Data Respon Biner
## Sumber Data: Agresti (2015), hlm: 188, "Dose–Response Study"
dosis <- c(1.691, 1.724, 1.755, 1.784, 1.811, 1.837, 1.861, 1.884)
mati <- c(6, 13, 18, 28, 52, 53, 61, 60)
n <- c(59, 60, 62, 56, 63, 59, 62, 60)
hidup <- n - mati
data1 <- matrix(append(mati,hidup),ncol=2)
data1
## [,1] [,2]
## [1,] 6 53
## [2,] 13 47
## [3,] 18 44
## [4,] 28 28
## [5,] 52 11
## [6,] 53 6
## [7,] 61 1
## [8,] 60 0
Model Logit
model.logit <- glm(data1 ~ dosis, family=binomial(link=logit))
summary(model.logit)
##
## Call:
## glm(formula = data1 ~ dosis, family = binomial(link = logit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.740 5.182 -11.72 <2e-16 ***
## dosis 34.286 2.913 11.77 <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: 284.202 on 7 degrees of freedom
## Residual deviance: 11.116 on 6 degrees of freedom
## AIC: 41.314
##
## Number of Fisher Scoring iterations: 4
Model Probit
model.probit <- glm(data1 ~ dosis, family=binomial(link=probit))
summary(model.probit)
##
## Call:
## glm(formula = data1 ~ dosis, family = binomial(link = probit))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -34.956 2.649 -13.20 <2e-16 ***
## dosis 19.741 1.488 13.27 <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: 284.202 on 7 degrees of freedom
## Residual deviance: 9.987 on 6 degrees of freedom
## AIC: 40.185
##
## Number of Fisher Scoring iterations: 4
Model C-log-log
model.cloglog <- glm(data1 ~ dosis, family=binomial(link=cloglog))
summary(model.cloglog)
##
## Call:
## glm(formula = data1 ~ dosis, family = binomial(link = cloglog))
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -39.522 3.236 -12.21 <2e-16 ***
## dosis 22.015 1.797 12.25 <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: 284.2024 on 7 degrees of freedom
## Residual deviance: 3.5143 on 6 degrees of freedom
## AIC: 33.712
##
## Number of Fisher Scoring iterations: 4
Nilai Prediksi
## Prediksi pi(i), peluang(Y=1)
fitted.values(model.logit)
## 1 2 3 4 5 6 7
## 0.05937747 0.16366723 0.36162283 0.60490961 0.79440490 0.90405532 0.95546748
## 8
## 0.97925643
fitted.values(model.probit)
## 1 2 3 4 5 6 7 8
## 0.0577367 0.1781060 0.3780390 0.6032833 0.7866532 0.9045852 0.9626183 0.9873227
fitted.values(model.cloglog)
## 1 2 3 4 5 6 7
## 0.09582195 0.18802653 0.33777217 0.54177644 0.75683967 0.91843509 0.98575181
## 8
## 0.99913561
## Prediksi Yi
n*fitted.values(model.logit)
## 1 2 3 4 5 6 7 8
## 3.503271 9.820034 22.420616 33.874938 50.047508 53.339264 59.238984 58.755386
n*fitted.values(model.probit)
## 1 2 3 4 5 6 7 8
## 3.406465 10.686363 23.438420 33.783862 49.559153 53.370529 59.682337 59.239363
n*fitted.values(model.cloglog)
## 1 2 3 4 5 6 7 8
## 5.653495 11.281592 20.941875 30.339480 47.680899 54.187670 61.116612 59.948137
Kurva: Peluang atau proporsi aktual berupa titik
ggplot() +
## Garis halus untuk model
geom_smooth(data = df_model,
aes(x = X_dosis, y = Nilai, color = Jenis, linetype = Jenis),
method = "loess", se = FALSE, linewidth = 0.9) +
## Titik untuk aktual
geom_point(data = df_aktual,
aes(x = X_dosis, y = Nilai),
size = 3, color = "black") +
## Warna dan tipe garis manual
scale_color_manual(values = c(
"P_logit" = "red",
"P_probit" = "green",
"P_cloglog" = "blue"
)) +
scale_linetype_manual(values = c(
"P_logit" = "solid",
"P_probit" = "solid",
"P_cloglog" = "solid"
)) +
labs(
title = "Hubungan X_dosis dengan Proporsi(Y=1): Aktual & Prediksi",
x = "X_dosis",
y = "Proporsi(Y=1), Aktual dan Prediksi",
color = "Prediksi P:",
linetype = "Prediksi P:"
) +
## Pengaturan tema grafik
theme_classic(base_size = 13) +
theme(
legend.position = "bottom",
legend.text = element_text(size = 10),
axis.text = element_text(color = "black"),
panel.grid.major = element_line(color = "gray80", size = 0.4),
panel.grid.minor = element_line(color = "gray90", size = 0.1),
plot.margin = margin(10, 20, 10, 10)
)
