Estimation of the frequency of neurocysticercosis among people with headache
Before talking about Bayesian binomial model, let’s review some important concepts.
The binomial distribution describes the probability of observing a certain number of successes in a fixed number of independent trials, each with the same probability of success. The binomial distribution is frequently used to model the number of successes in a sample drawn with replacement (that is important) from a population.
In simpler terms, it’s used when you have an experiment that:
Let’s walk through a simple example. Suppose we want to estimate the prevalence of Disease X—that is, the proportion of people in Montréal, Canada, in the year 2025 who have the disease at a specific point in time. The true prevalence in the entire population (often referred to as the parameter) is unknown. Since it is not practical to test every individual living in Montréal, we instead collect data from a sample of the population and use that sample to estimate the true prevalence.
Let’s assume we randomly select 1,000 individuals from the population and test them. This situation fits perfectly within a binomial distribution framework because:
n, is fixed at 1,000.p, of having the diseaseIf 30 people test positive, then we have observed y = 30 “Sick” outcomes. The full range of possible outcomes (from 0 to 1,000 positive results), along with their corresponding probabilities, forms a binomial distribution.
For more information about the binomial distribution:https://www.itl.nist.gov/div898/handbook/eda/section3/eda366i.htm
So, what is the prevalence estimate for the above example?
Answer:
## So, we have -
y <- 30 # Number of people with disease
n <- 1000 # Number of people tested
## Prevalence of X disease in Montréal
Prevalence <- y/n
Prevalence[1] 0.03## If we want to express in percentage
paste0("Prevalence of X disease: ",Prevalence*100,"%")[1] "Prevalence of X disease: 3%"So, we can say that the prevalence of disease X in the population of Montréal is 3%.
However, this interpretation is incomplete. Because, a single proportion estimate does not account for uncertainty — such as random sampling variation, test accuracy, or selection bias. In practice, we want to quantify this uncertainty by constructing a confidence interval (CI) around our estimate. This uncertainty reflects how confident we are that the true prevalence in the population lies within a certain range.
There are several ways to calculate confidence intervals (CIs), such as the Wald approximation, Wilson method, and Clopper–Pearson method.
Let’s calculate CI for our prevalence estimate,
binom.test(y, n, conf.level = 0.95)
Exact binomial test
data: y and n
number of successes = 30, number of trials = 1000, p-value < 2.2e-16
alternative hypothesis: true probability of success is not equal to 0.5
95 percent confidence interval:
0.02033049 0.04255140
sample estimates:
probability of success
0.03 Now, we can say that “The prevalence of X disease in the population of Montréal is 3% (95% CI: 2% – 4%)”.
What does this mean?
This result means the following:
If we were to repeat this same study many times — each time randomly selecting 1000 individuals and calculating a 95% confidence interval — then approximately 95% of those intervals would contain the true population prevalence, while 5% would not.
In other words, the method we used to construct the interval has a long-run success rate of 95% in capturing the true value.
However, for our single study, we cannot say with certainty whether the true prevalence lies inside or outside this particular interval — we only know that the procedure has a 95% chance of producing a correct interval over many repetitions.
😔
That’s not quite what we usually want to know, is it?
In practice, we often wish to make a direct probability statement about the true prevalence, such as:
“We are 95% confident that the true prevalence lies between 2% and 4%.”
Unfortunately, this statement is not technically correct for our previous analysis.
Why?. Because the true prevalence is fixed — it’s not something that changes from one study to another. The uncertainty lies in our sample, NOT in the parameter itself.
Of note, the analysis we performed so far represents the frequentist approach to estimating prevalence for binomial data.
There is another powerful framework called the Bayesian approach, which allows us to directly express uncertainty about the parameter itself — for example, by saying that the probability of the prevalence being between 2% and 4% is 95%.
This interpretation aligns much better with how researchers naturally think about uncertainty.
Bayesian statistics provide a way to update our beliefs about an unknown parameter (such as disease prevalence) based on observed data.
It starts with a prior distribution, which represents what we believe about the parameter before seeing the data (for example, previous studies might suggest the prevalence is around 5%).
Where can we get prior information? - From previous studies, reasonable guesses from experts, or we can simply say that “we don’t know.”
Does the prior influence the final results? - Yes — it matters when the data are weak, and not much otherwise.
Then, we collect new data (e.g., 30 positive cases out of 1000 tests), and use Bayes’ theorem to combine this prior information with the data’s likelihood.
The result is a posterior distribution — a new, updated belief about the parameter after considering the data.
Mathematically:
\[ \text{Posterior} \propto \text{Likelihood} \times \text{Prior} \]
From the posterior, we can calculate a credible interval, which gives the range of values within which the true parameter lies with a given probability . Unlike frequentist confidence intervals, Bayesian credible intervals directly express uncertainty about the parameter itself. So, in Bayesian statistics, we assume that “true prevalence” is not a fixed value (*****That is very important****). It has its own distribution.
For a bit more about Bayesian statistics visit: https://jalal2025.quarto.pub/my-web/Bayesian_statistics_course.html
The Bayesian binomial model is a statistical approach used to estimate the proportion or prevalence of a binary outcome, such as the presence or absence of a disease. It combines observed data with prior knowledge through Bayes’ theorem to generate a posterior distribution for the parameter of interest.
In the context of disease prevalence:
y be the number of cases observed in a sample of size n.y ~ Binomial(n, p), where p is the true prevalence.p to reflect previous knowledge or uncertainty (commonly Beta or uniform).p integrates the prior and the observed data, providing a probabilistic estimate of prevalence with credible intervals.Advantages of the Bayesian approach include:
Applications: Estimating disease prevalence, infection rates, or any scenario with binary outcomes in epidemiology, clinical trials, or public health research.
Let’s calculate again the prevalence of X disease.
First, we will assume that we not know what is the true prevalence of X disease in the population of Montréal. But, as are talking about proportion, so, it will be between 0 and 100. So,
Beta(1,1) is the approximate representation of our prior information. People usually say such information as “non-informative or vague prior”. Notably, non-informative prior in one scale can be very informative prior in another scale.## Prior: Uniform prior (Beta(1,1))
alpha_prior <- 1
beta_prior <- 1
## Data
positive_cases <- 30
total_tests <- 1000
## Posterior parameters (Beta distribution)
alpha_posterior <- alpha_prior + positive_cases
beta_posterior <- beta_prior + total_tests - positive_cases
## Calculate posterior mean (point estimate of prevalence)
posterior_mean <- round((alpha_posterior / (alpha_posterior + beta_posterior))*100,3)
## Calculate posterior median
posterior_median <- round((qbeta(0.5, alpha_posterior, beta_posterior))*100,3)
## Calculate posterior mode
posterior_mode <- round(((alpha_posterior - 1) / (alpha_posterior + beta_posterior - 2))*100,3)
## Calculate 95% credible interval
credible_interval <- round((qbeta(c(0.025, 0.975), alpha_posterior, beta_posterior))*100,3)
## Results table
results <- data.frame(
Statistic = c("Mean", "Median", "Mode", "95% CrI (lower)", "95% CrI (upper)"),
Value = c(posterior_mean, posterior_median, posterior_mode, credible_interval[1], credible_interval[2])
)
results Statistic Value
1 Mean 3.094
2 Median 3.063
3 Mode 3.000
4 95% CrI (lower) 2.114
5 95% CrI (upper) 4.251We can plot the posterior distribution:
# Use unscaled (proportion) values for plotting
posterior_mean_prop <- posterior_mean / 100
posterior_median_prop <- posterior_median / 100
posterior_mode_prop <- posterior_mode / 100
credible_interval_prop <- credible_interval / 100
# Plot posterior distribution
curve(dbeta(x, alpha_posterior, beta_posterior),
from = 0.02, to = 0.045,
main = "Posterior Distribution of Prevalence\nwith Central Tendency Measures",
xlab = "Prevalence (proportion)", ylab = "Density",
col = "blue", lwd = 2)
# Add vertical lines
abline(v = posterior_mean_prop, col = "red", lty = 1, lwd = 2)
abline(v = posterior_median_prop, col = "darkgreen", lty = 2, lwd = 2)
abline(v = posterior_mode_prop, col = "purple", lty = 3, lwd = 2)
# Add credible interval
abline(v = credible_interval_prop, col = "gray", lty = 3)
# Legend
legend("topright",
legend = c("Posterior", "Mean", "Median", "Mode", "95% CrI"),
col = c("blue", "red", "darkgreen", "purple", "gray"),
lty = c(1, 1, 2, 3, 3), lwd = 2)
Now, this analysis solve our fundamental problem. We can say that
Based on the Bayesian binomial model with a uniform prior, the estimated prevalence of X disease in Montréal is 3.09% (95% CrI: 2.11%–4.25%), meaning there is a 95% probability that the true prevalence lies within this range.
However, one can now want to know that -
“What is the probability that the true prevalence of X disease is 3% or less?”
Unfortunately, we cannot get the answer of such question form our previous analyses. Fortunately, there is way to answer this question by doing Monte Carlo simulation.
## MONTE CARLO APPROACH (rbeta sampling)
set.seed(123) # for reproducible results
n_sim <- 10000
posterior_samples <- rbeta(n_sim,
alpha_posterior,
beta_posterior)
# post_mode <- density(posterior_samples)$x[which.max(density(posterior_samples)$y)]
## Results table
results <- data.frame(
Statistic = c("Mean",
"Median",
"95% CrI (lower)",
"95% CrI (upper)"),
Value = c(mean(posterior_samples)*100,
median(posterior_samples)*100,
quantile(posterior_samples, c(0.025))*100,
quantile(posterior_samples, c(0.975))*100
))
results Statistic Value
1 Mean 3.095348
2 Median 3.061198
3 95% CrI (lower) 2.097944
4 95% CrI (upper) 4.239486Plot the results
# Monte Carlo histogram
hist(posterior_samples, breaks = 50, prob = TRUE,
main = NULL,
xlab = "Prevalence", ylab = "Density",
xlim = c(0.02, 0.045), col = "lightblue")
lines(density(posterior_samples), col = "red", lwd = 2)
legend("topright", legend = c("Sample density", "Histogram"),
col = c("red", "lightblue"), lwd = 2, fill = c(NA, "lightblue"))
Figure: Posterior distribution of prevalence.
Now, we can calculate “What is the probability that the true prevalence of X disease is 3% or less?”
## Probability that prevalence <= 3%
prob_leq_3 <- mean(posterior_samples <= 0.03)*100
prob_leq_3[1] 45.21We can interpret this as:
“Based on the observed data and our uniform prior, there is approximately 45% probability that the true prevalence is 3% or less.”
So far, we have learnt how Bayesian binomial model works. And we have calculated prevalence of disease X in one population. But in real world we may require to calculate the prevalence of same disease (or multiple disease) in different populations.
Let’s see, how we can do such analysis!!
For this purpose,
We will calculate prevalence of a disease, called “neurocysticercosis”, in different population using a real world data set.
We will present the results in a nice table.
We have recently completed a brief systematic review of research articles reporting the frequency of neurocysticercosis (NCC, a parasitic disease of human central nervous system caused by the larval stage of zoonotic tapeworm, Taenia solium) among people with headache (PWH). In that review, we identified 10 studies which met the inclusion criteria (data on NCC among PWH, and NCC diagnosis was based on brain imaging). The data set is given in the “Data” folder. The data is in Microsoft Excel (sheet = PWH). There are only five variables in this the data set:
id: Unique ID for each study
first_author: First author
study_setting: Study setting (“Clinic_or_hospital” OR “Community”)
n: Number of people with headache manifestation who underwent brain imaging.
y: Number of NCC-cases
First, import data (users are advised to keep all file within same folder, otherwise you may require to set working directory manually using setwd() function)
tdat
library(here)
tdat <- read_excel(here("Data", "NCC_among_PWH.xlsx"))
tdat# A tibble: 10 × 5
id first_author study_setting n y
<dbl> <chr> <chr> <dbl> <dbl>
1 500001 Pradhan et al Clinic_or_hospital 166 72
2 500002 Valtis et al Clinic_or_hospital 924 119
3 500003 deOliveiraTaveira et al Clinic_or_hospital 108 12
4 500004 DelBrutto et al Clinic_or_hospital 1017 48
5 500005 Sanchez et al Clinic_or_hospital 16 12
6 500006 Cruz et al Community 69 21
7 500007 Stelzle et al Clinic_or_hospital 17 0
8 500008 Moyano et al Community 3 2
9 500009 Valenca et al Clinic_or_hospital 78 3
10 500010 Unpublished Community 131 10## Working data set
df <- tdatLet’s check what are the data we have,
glimpse(df)Rows: 10
Columns: 5
$ id <dbl> 500001, 500002, 500003, 500004, 500005, 500006, 500007, …
$ first_author <chr> "Pradhan et al", "Valtis et al", "deOliveiraTaveira et a…
$ study_setting <chr> "Clinic_or_hospital", "Clinic_or_hospital", "Clinic_or_h…
$ n <dbl> 166, 924, 108, 1017, 16, 69, 17, 3, 78, 131
$ y <dbl> 72, 119, 12, 48, 12, 21, 0, 2, 3, 10head(df)# A tibble: 6 × 5
id first_author study_setting n y
<dbl> <chr> <chr> <dbl> <dbl>
1 500001 Pradhan et al Clinic_or_hospital 166 72
2 500002 Valtis et al Clinic_or_hospital 924 119
3 500003 deOliveiraTaveira et al Clinic_or_hospital 108 12
4 500004 DelBrutto et al Clinic_or_hospital 1017 48
5 500005 Sanchez et al Clinic_or_hospital 16 12
6 500006 Cruz et al Community 69 21summary(df) id first_author study_setting n
Min. :5e+05 Length:10 Length:10 Min. : 3.0
1st Qu.:5e+05 Class :character Class :character 1st Qu.: 30.0
Median :5e+05 Mode :character Mode :character Median : 93.0
Mean :5e+05 Mean : 252.9
3rd Qu.:5e+05 3rd Qu.: 157.2
Max. :5e+05 Max. :1017.0
y
Min. : 0.00
1st Qu.: 4.75
Median : 12.00
Mean : 29.90
3rd Qu.: 41.25
Max. :119.00 We see the inputs for variable “study_setting” are a bit incosistant. Let’s first check what are the inputs
unique(df$study_setting)[1] "Clinic_or_hospital" "Community" We see there are two types of inputs, `
"Clinic_or_hospital""Community" We will make Clinic_or_hospital a bit readable
df <- df %>%
mutate(study_setting = case_when(
study_setting == "Clinic_or_hospital" ~ "Health center",
TRUE ~ study_setting
))Now, we want to calculate the proportion of NCC among PWH for each row of the data set and express it as a percentage along with the associated 95% credible intervals. To complete the task,
Stan model. Then, we will compile the mode (it may take few minutes).The model assumes that the number of NCC cases \(y_i\) in study \(i\) follows a binomial distribution with study-specific sample size \(n_i\) and probability parameter \(p_i\):
\[ \begin{aligned} y_i &\sim \text{Binomial}(n_i, p_i) \\ p_i &\sim \text{Beta}(1, 1) \end{aligned} \]
A uniform \(\text{Beta}(1, 1)\) prior was specified for \(p_i\), representing an uninformative prior distribution that is equivalent to a uniform distribution over the interval \((0, 1)\).
model_code <- "
data{
int<lower=0> N_studies;
int<lower=0> n[N_studies];
int<lower=0> y[N_studies];
int<lower=0> alpha_prior;
int<lower=0> beta_prior;
}
parameters {
vector <lower=0, upper=1> [N_studies]p;
}
model{
p ~ beta(alpha_prior,beta_prior);
for (i in 1:N_studies){
y[i] ~ binomial(n[i],p[i]);
}
}
generated quantities{
vector[N_studies] y_rep;
for (i in 1:N_studies){
y_rep[i] = binomial_rng (n[i],p[i]);
}
}
"
compiled_model <- stan_model(model_code = model_code)df <- df %>%
mutate(st_id = paste0("St_",id))
df1 <- df %>%
select(st_id,n,y)
model_data <- list(
N_studies = nrow(df),
n = df$n,
y = df$y,
alpha_prior = 1,
beta_prior =1
)## Fit the model
fit <- sampling(
compiled_model,
data = model_data,
iter = 10000,
warmup = 5000,
chains = 4,
seed = 123
)It is better to save the results as .rds file. In that case, we will not need to run the model, again we can just, check the “Results”
library(here)
saveRDS(fit, file = here("Outputs", "Bayesian_model_results.rds"))First, reload the model and print the results
library(here)
x <- readRDS(here("Outputs","Bayesian_model_results.rds"))
print(x, pars = "p", probs = c(0.025,0.50,0.975) )Inference for Stan model: anon_model.
4 chains, each with iter=10000; warmup=5000; thin=1;
post-warmup draws per chain=5000, total post-warmup draws=20000.
mean se_mean sd 2.5% 50% 97.5% n_eff Rhat
p[1] 0.43 0 0.04 0.36 0.43 0.51 38436 1
p[2] 0.13 0 0.01 0.11 0.13 0.15 37357 1
p[3] 0.12 0 0.03 0.07 0.12 0.18 35724 1
p[4] 0.05 0 0.01 0.04 0.05 0.06 36390 1
p[5] 0.72 0 0.10 0.50 0.73 0.90 38862 1
p[6] 0.31 0 0.05 0.21 0.31 0.42 39721 1
p[7] 0.05 0 0.05 0.00 0.04 0.19 31598 1
p[8] 0.60 0 0.20 0.20 0.61 0.93 39946 1
p[9] 0.05 0 0.02 0.01 0.05 0.11 39272 1
p[10] 0.08 0 0.02 0.04 0.08 0.14 39839 1
Samples were drawn using NUTS(diag_e) at Thu Oct 23 00:47:51 2025.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at
convergence, Rhat=1).Now, check the convergence
How to check model convergence ?
In Bayesian statistics, convergence refers to two main concepts: the posterior distribution converging to the “true” distribution as data grows, and the MCMC sampling algorithm converging to its stationary distribution.
Source: https://ourcodingclub.github.io/tutorials/brms/
Bulk_ESS and Tail_ESS > 1000Rhat = 1We can easily check (1) and (2) from previous chunk of code. For different plots there is an amazing function launch_shinystan() from R-package {shinystan}.
launch_shinystan(x)So far, we have calculated the prevalence proportion for multiple studies using Bayesian binomial model.
Now, we want to present the results in a nice table.
# Extract posterior samples
posterior_samples <- extract(x)
# Get summary statistics for theta parameters
theta_samples <- posterior_samples$p
colnames(theta_samples) <- df1$st_id
# Calculate posterior summaries
posterior_summaries <- data.frame(
st_id = df1$st_id,
y = df1$y,
n = df1$n,
Posterior_Mean = apply(theta_samples, 2, mean),
Posterior_Median = apply(theta_samples, 2, median),
CI_2.5 = apply(theta_samples, 2, quantile, probs = 0.025),
CI_97.5 = apply(theta_samples, 2, quantile, probs = 0.975),
SD = apply(theta_samples, 2, sd)
)
df2 <- df %>%
left_join(posterior_summaries, by = "st_id") %>%
mutate(Prop = paste0(
sprintf("%.1f",Posterior_Median*100)," (",
sprintf("%.1f",CI_2.5*100)," - ",
sprintf("%.1f",CI_97.5*100),")")) %>%
select(id,first_author,study_setting,Prop)
table <- df2 %>%
arrange(id) %>%
gt() %>%
cols_label(
id = "Study ID",
first_author = "Reference",
study_setting = "Setting",
Prop = md("Bayesian estimate<br>%(95% CrI)")
) %>%
tab_style(
style = cell_text(weight = "bold"),
locations = cells_column_labels()
) %>%
tab_style(
style = cell_text(align = "center"),
locations = cells_column_labels()
) %>%
tab_style(
style = cell_text(align = "center"),
locations = cells_body()
) %>%
fmt_number(
columns = c(id),
decimals = 0
) %>%
tab_style(style = list(cell_fill(color = "lightgray")),
locations = cells_body(rows = seq(1, nrow(df2), 2))) %>%
tab_options(
table.font.size = px(12),
table.font.names = "Times New Roman"
)
table| Study ID | Reference | Setting | Bayesian estimate %(95% CrI) |
|---|---|---|---|
| 500,001 | Pradhan et al | Health center | 43.4 (35.9 - 51.1) |
| 500,002 | Valtis et al | Health center | 12.9 (10.9 - 15.1) |
| 500,003 | deOliveiraTaveira et al | Health center | 11.6 (6.6 - 18.4) |
| 500,004 | DelBrutto et al | Health center | 4.8 (3.6 - 6.2) |
| 500,005 | Sanchez et al | Health center | 73.1 (50.3 - 89.7) |
| 500,006 | Cruz et al | Community | 30.8 (20.9 - 42.1) |
| 500,007 | Stelzle et al | Health center | 3.8 (0.1 - 18.9) |
| 500,008 | Moyano et al | Community | 61.2 (19.6 - 93.1) |
| 500,009 | Valenca et al | Health center | 4.6 (1.4 - 10.8) |
| 500,010 | Unpublished | Community | 8.1 (4.3 - 13.6) |
Make the table a bit pretty
# Create a prettier table
table <- df2 %>%
arrange(id) %>%
gt() %>%
# Column labels with better formatting
cols_label(
id = md("**Study ID**"),
first_author = md("**Reference**"),
study_setting = md("**Study Setting**"),
Prop = md("**Bayesian Estimate**<br>**% (95% CrI)**")
) %>%
# Header styling
tab_header(
title = md("**Neurocysticercosis Prevalence Among People with Headaches**"),
subtitle = "Bayesian Beta-Binomial Model Estimates by Study"
) %>%
# Column alignment
cols_align(
align = "center",
columns = c(study_setting, Prop)
) %>%
cols_align(
align = "left",
columns = c(first_author)
) %>%
# Zebra striping with subtle colors
tab_style(
style = list(
cell_fill(color = "#f8f9fa"),
cell_text(weight = "normal")
),
locations = cells_body(rows = seq(1, nrow(df2), 2))
) %>%
# Header style
tab_style(
style = list(
cell_fill(color = "#2c3e50"),
cell_text(color = "white", weight = "bold", size = px(14))
),
locations = cells_column_labels()
) %>%
# Title styling
tab_style(
style = list(
cell_text(weight = "bold", size = px(16), color = "#2c3e50")
),
locations = cells_title(groups = "title")
) %>%
# Border styling
tab_style(
style = cell_borders(
sides = c("bottom", "top"),
color = "#dee2e6",
weight = px(1)
),
locations = cells_body()
) %>%
# Table options
tab_options(
table.font.names = "Times New Roman",
table.font.size = px(12),
heading.align = "center",
column_labels.border.top.color = "white",
column_labels.border.bottom.color = "#2c3e50",
column_labels.border.bottom.width = px(2),
table_body.border.bottom.color = "white",
data_row.padding = px(6)
)
table_html <- table %>% as_raw_html()
table| Neurocysticercosis Prevalence Among People with Headaches | |||
|---|---|---|---|
| Bayesian Beta-Binomial Model Estimates by Study | |||
| Study ID | Reference | Study Setting | Bayesian Estimate % (95% CrI) |
| 500001 | Pradhan et al | Health center | 43.4 (35.9 - 51.1) |
| 500002 | Valtis et al | Health center | 12.9 (10.9 - 15.1) |
| 500003 | deOliveiraTaveira et al | Health center | 11.6 (6.6 - 18.4) |
| 500004 | DelBrutto et al | Health center | 4.8 (3.6 - 6.2) |
| 500005 | Sanchez et al | Health center | 73.1 (50.3 - 89.7) |
| 500006 | Cruz et al | Community | 30.8 (20.9 - 42.1) |
| 500007 | Stelzle et al | Health center | 3.8 (0.1 - 18.9) |
| 500008 | Moyano et al | Community | 61.2 (19.6 - 93.1) |
| 500009 | Valenca et al | Health center | 4.6 (1.4 - 10.8) |
| 500010 | Unpublished | Community | 8.1 (4.3 - 13.6) |