# ==== 2. Fit Bayesian model for Sensitivity ====
fit_Se <- brm(
formula = TP | trials(TP + FN) ~ 1 + (1 | id),
data = ebpp,
family = binomial(),
prior = c(prior(normal(0, 5), class = "Intercept"),
prior(exponential(1), class = sd)),
chains = 4,
iter = 4000,
warmup = 1000,
cores = 4,
seed = 2025
)
## Compiling Stan program...
## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.0 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Trying to compile a simple C file
## Running "C:/PROGRA~1/R/R-45~1.0/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc -I"C:/PROGRA~1/R/R-45~1.0/include" -DNDEBUG -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
## from C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
## from C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.0/etc/x64/Makeconf:289: foo.o] Error 1
## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.0 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Start sampling
fit_Se
## Family: binomial
## Links: mu = logit
## Formula: TP | trials(TP + FN) ~ 1 + (1 | id)
## Data: ebpp (Number of observations: 42)
## Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.38 0.18 1.07 1.76 1.00 3012 4529
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.12 0.23 1.67 2.57 1.00 2118 3519
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# ==== 3. Fit Bayesian model for Specificity ====
fit_Sp <- brm(
formula = TN | trials(TN + FP) ~ 1 + (1 | id),
data = ebpp,
family = binomial(),
prior = c(prior(normal(0, 5), class = "Intercept"),
prior(exponential(1), class = sd)),
chains = 4,
iter = 4000,
warmup = 1000,
cores = 4,
seed = 2025
)
## Compiling Stan program...
## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.0 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Trying to compile a simple C file
## Running "C:/PROGRA~1/R/R-45~1.0/bin/x64/Rcmd.exe" SHLIB foo.c
## using C compiler: 'gcc.exe (GCC) 14.2.0'
## gcc -I"C:/PROGRA~1/R/R-45~1.0/include" -DNDEBUG -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/Rcpp/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/unsupported" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/BH/include" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/src/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/" -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppParallel/include/" -DRCPP_PARALLEL_USE_TBB=1 -I"C:/Users/Dell/AppData/Local/R/win-library/4.5/rstan/include" -DEIGEN_NO_DEBUG -DBOOST_DISABLE_ASSERTS -DBOOST_PENDING_INTEGER_LOG2_HPP -DSTAN_THREADS -DUSE_STANC3 -DSTRICT_R_HEADERS -DBOOST_PHOENIX_NO_VARIADIC_EXPRESSION -D_HAS_AUTO_PTR_ETC=0 -include "C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp" -std=c++1y -I"C:/rtools45/x86_64-w64-mingw32.static.posix/include" -O2 -Wall -std=gnu2x -mfpmath=sse -msse2 -mstackrealign -c foo.c -o foo.o
## cc1.exe: warning: command-line option '-std=c++14' is valid for C++/ObjC++ but not for C
## In file included from C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Core:19,
## from C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/Dense:1,
## from C:/Users/Dell/AppData/Local/R/win-library/4.5/StanHeaders/include/stan/math/prim/fun/Eigen.hpp:22,
## from <command-line>:
## C:/Users/Dell/AppData/Local/R/win-library/4.5/RcppEigen/include/Eigen/src/Core/util/Macros.h:679:10: fatal error: cmath: No such file or directory
## 679 | #include <cmath>
## | ^~~~~~~
## compilation terminated.
## make: *** [C:/PROGRA~1/R/R-45~1.0/etc/x64/Makeconf:289: foo.o] Error 1
## WARNING: Rtools is required to build R packages, but is not currently installed.
##
## Please download and install the appropriate version of Rtools for 4.5.0 from
## https://cran.r-project.org/bin/windows/Rtools/.
## Start sampling
fit_Sp
## Family: binomial
## Links: mu = logit
## Formula: TN | trials(TN + FP) ~ 1 + (1 | id)
## Data: ebpp (Number of observations: 42)
## Draws: 4 chains, each with iter = 4000; warmup = 1000; thin = 1;
## total post-warmup draws = 12000
##
## Multilevel Hyperparameters:
## ~id (Number of levels: 42)
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept) 1.23 0.15 0.97 1.58 1.00 1810 3079
##
## Regression Coefficients:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 1.71 0.19 1.33 2.10 1.01 1018 2231
##
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
# ==== 4. Posterior summaries ====
plogis(fixef(fit_Se)["Intercept", "Estimate"]) # Trung bình độ nhạy
## [1] 0.8923745
plogis(fixef(fit_Sp)["Intercept", "Estimate"]) # Trung bình độ đặc hiệu
## [1] 0.8473061
# ==== 5. Posterior distributions ====
draws_Se <- plogis(as_draws_df(fit_Se)$b_Intercept)
draws_Sp <- plogis(as_draws_df(fit_Sp)$b_Intercept)
p1 <- mcmc_areas(as_draws_df(fit_Se), pars = "b_Intercept") +
ggtitle("Posterior logit(Sensitivity)")
p2 <- mcmc_areas(as_draws_df(fit_Sp), pars = "b_Intercept") +
ggtitle("Posterior logit(Specificity)")
grid.arrange(p1, p2, ncol = 2)
# ==== 6. Histograms of posterior draws ====
hist(draws_Se, breaks = 40, col = "skyblue",
main = "Posterior of Sensitivity",
xlab = "Sensitivity")
hist(draws_Sp, breaks = 40, col = "salmon",
main = "Posterior of Specificity",
xlab = "Specificity")
# ==== 7. Forest plot - Sensitivity ====
mu_se <- fixef(fit_Se)["Intercept", "Estimate"]
post_Se <- as.data.frame(ranef(fit_Se)$id[, , "Intercept"]) %>%
rename(Estimate = Estimate, Q2.5 = Q2.5, Q97.5 = Q97.5)
id_names <- rownames(ranef(fit_Se)$id)
forest_data_Se <- data.frame(id = id_names,
logit_mean = mu_se + post_Se$Estimate,
logit_lower = mu_se + post_Se$Q2.5,
logit_upper = mu_se + post_Se$Q97.5) %>%
mutate(Se_mean = plogis(logit_mean),
Se_lower = plogis(logit_lower),
Se_upper = plogis(logit_upper)) %>%
left_join(id_study, by = "id") %>%
mutate(Study = ifelse(is.na(Study), id, Study))
ggplot(forest_data_Se, aes(x = Se_mean, y = reorder(Study, Se_mean))) +
geom_point() +
geom_errorbarh(aes(xmin = Se_lower, xmax = Se_upper), height = 0.2) +
geom_vline(xintercept = plogis(mu_se), linetype = "dashed", color = "red") +
labs(x = "Độ nhạy hậu định", y = "Nghiên cứu", title = "Forest plot - Sensitivity") +
theme_minimal()
# ==== 8. Forest plot - Specificity ====
mu_sp <- fixef(fit_Sp)["Intercept", "Estimate"]
post_Sp <- as.data.frame(ranef(fit_Sp)$id[, , "Intercept"]) %>%
rename(Estimate = Estimate, Q2.5 = Q2.5, Q97.5 = Q97.5)
id_names <- rownames(ranef(fit_Sp)$id)
forest_data_Sp <- data.frame(id = id_names,
logit_mean = mu_sp + post_Sp$Estimate,
logit_lower = mu_sp + post_Sp$Q2.5,
logit_upper = mu_sp + post_Sp$Q97.5) %>%
mutate(Sp_mean = plogis(logit_mean),
Sp_lower = plogis(logit_lower),
Sp_upper = plogis(logit_upper)) %>%
left_join(id_study, by = "id") %>%
mutate(Study = ifelse(is.na(Study), id, Study))
ggplot(forest_data_Sp, aes(x = Sp_mean, y = reorder(Study, Sp_mean))) +
geom_point() +
geom_errorbarh(aes(xmin = Sp_lower, xmax = Sp_upper), height = 0.2) +
geom_vline(xintercept = plogis(mu_sp), linetype = "dashed", color = "blue") +
labs(x = "Độ đặc hiệu hậu định",
y = "Nghiên cứu",
title = "Forest plot - Specificity") +
theme_minimal()
## het ##