Authors: Razieh Bidhendi Yarandi, Mohammad Ali Mansournia, Hojjat Zeraati, Kazem Mohammad
Categories: Methodology, Bayesian methods, Data augmentation, Importance sampling, MCMC, Rejection sampling
Source: Global Epidemiology
Bayesian inference has become popular. It offers several pragmatic approaches to account for uncertainty in inference decision-making. Various estimation methods have been introduced to implement Bayesian methods. Although these algorithms are powerful, they are not always easy to grasp for non-statisticians. This paper aims to provide an intuitive framework of four essential Bayesian computational methods for epidemiologists and other health researchers. We do not cover an extensive mathematical discussion of these approaches, but instead offer a non-quantitative description of these algorithms and provide some illuminating examples.
Bayesian computational methods, namely importance sampling, rejection sampling, Markov chain Monte Carlo, and data augmentation are presented.
The substantial amount of research published on Bayesian inference has highlighted its popularity among researchers, while the basic concepts are not always straightforward for interested learners. We show that alternative approaches such as a weighted prior approach, which are intuitively appealing and easy-to-understand, work well in the case of low-dimensional problems and appropriate prior information. Otherwise, MCMC is a trouble-free tool in those cases.
Keywords: Bayesian methods, Data augmentation, Importance sampling, MCMC, Rejection sampling
The growth in the application of Bayesian analysis in the sciences has created a need to present its complex concepts in a language understandable for non-statisticians. To many interested users, its computationally-intensive simulation approaches have the appearance of a “black box.” This article aims to describe in an intuitive manner the pragmatic approaches used in Bayesian methods for posterior simulation. Four Bayesian computational methods are Importance Sampling (IS), Rejection Sampling (RS), Markov Chain Monte Carlo (MCMC) and Data Augmentation (DA). In this study, we aim to provide a comprehensive, yet convenient, explanation of these techniques. Several examples are presented to illustrate the algorithms and the concepts they embody. R software code and other information is available in the Appendix.
In statistical inference, it is not always possible to analytically compute the estimators associated with a given paradigm such as maximum likelihood, Bayes, etc. Numerical approaches are proposed as solutions. Historically, one of the first methods of simulation by a computer was what is known as Monte Carlo integration (MCI). Its computational algorithm relies on repeated random sampling from distributions of interest to obtain numerical results. For example, MCI is used to calculate the expectation of a distribution function when its estimation through integration is impossible. In other words, empirical experimental summation is substituted for analytical (or possibly numerical) integration. This technique is conceptually simple, but not always very efficient if the target probability function is not well-behaved, such as the “fat-tailed” nature of return distributions. Attempts to deal with this problem led to a gradual evolution, first to IS and next to RS, which are regarded as improved versions of MCI. However, in some more challenging cases, these have also proven inadequate to the task. Over the past 30 years or so, MCMC methods have revolutionized statistical computing and permitted ever more complex problems to be handled. Finally, we also discuss DA, which follows a different approach than Bayesian to calculate posterior estimates.
MCI as a computational method was first initiated to solve the integration problem in estimating expectations. Later it was applied to the simulation of Bayesian posterior (BP) distributions. It has a transparent Generate random samples from a distribution function, say “Target,” then numerically calculate the integral by calculating the mean of the values. The expectation obtained by MCI is referred to as the Empirical Average (Appendix Formula-1). For example, imagine we want to estimate the expectation of the function h(x) = exp (x), where the random variable (X) follows a Normal distribution with mean = 0 and SD = 1 (target). Firstly, generate 10,000 random samples from the target distribution, then obtain the values of h(x) from each generated sample, and then calculate numerically the mean and variance of this generated sample (Appendix R Code-1). This random sampling method is not cost-effective when the target distribution is diffuse. In this case, a large sample size is required to obtain acceptable precision.
MCI was improved by IS, a variance reduction technique first presented in statistical physics [1,2]. Contrary to MCI in which samples are treated evenly, IS does not emphasize the procedure of treating all parts of the distribution equally. Instead, it concentrates on those where estimation is critical. In this respect, an alternative function, say “Proposal,” close to the target, is suggested by making an educated guess. A weight which shows the importance of a generated sample is allocated to each generated sample through the “importance” function. For each sample, one compares the likelihood of obtaining that sample from the target distribution proportion to the “likelihood” of sampling it from the proposal distribution. After the sampling process is finished, the obtained relative likelihoods are normalized so that they sum to one. In this way, each point has its likelihood of occurrence as a discrete probability distribution. The expectation obtained by IS method is called Weighted Average (Appendix Formula-2).
Consider a Normal (0, 0.05) and a t-student (df = 1) distribution as the two different proposal distributions for Normal (0, 5). Estimate importance weights through the importance functions target:N05proposal:N0,0.05 and N05t1 for each generated point. After that, normalize the weights as N05N0,0.05∑N05N0,0.05 and N05t1∑N05t1. Then, we have a discrete distribution function for which the properties, such as the mean and the variance, are easily estimable. Estimated means for our generated samples are approximate to what is obtained from MCI, while their variances are considerably lower than the MCI approach (Appendix R Code-2). Using alternative distributions can improve the variance of the samples. A wide proposal distribution leads to inefficient estimates by variance and large sampling numbers (Fig. 1). Choosing an appropriate proposal distribution that looks similar to the target is ideal, though difficult to obtain at times. Unbiased estimates of parameters are obtained for large samples by IS. It also works well when the “importance” function is not highly variable. Indeed, an appropriate proposal distribution leads to lower variances and higher accuracy of the approximation. Robert and Casella provided an example illustrating the use of Normal (0,1) as a proposal to resemble sampling from a Cauchy C(0,1) target distribution caused infinite variance of importance weights [3]. This occurred due to attaching high importance weights to a few points and providing inefficient estimates in terms of variance. A better fit is guaranteed in substituting heavy tails distributions such as t-student.
Fig. 1 Comparison of target and proposals distribution.
In this section, we show how IS plays a role in BP simulation. According to the Bayes rule, the posterior distribution f(θ| y) is proportional to the likelihood of observed data f(y| θ) and the prior distribution of the parameter f(θ) (Appendix Formula-3).
It emphasizes that the likelihood can play the role of “importance” function for priors or for a prior. Using this approach, the expectation of a posterior distribution function can be calculated. In addition, to paint a clearer picture of the posterior density and its properties – such as percentiles – it suffices to generate weighted samples from the prior.
It follows a straightforward
This sampling method, referred to as Weighted Prior (WP), is equivalent to IS. As such, the WP method provides an alternative way to exact simulation from the posterior (Appendix Formula-4).
It is worth mentioning that the bias estimated through this approach is of an order 1/n; hence, large “samples” provide unbiased estimates. It should be considered that the contribution of generated samples to the estimation of the final posterior depends on how much it is supported by data. As such, this means that prior distributions should not be far off from the likelihood, or else many generated points receive small importance weights. We depicted a case when this method fails (Fig. 2).
Fig. 2 As an example, where prior is not supported by data, imagine prior (lower) as Normal (μ=10, SD = 2) and likelihood (upper) Normal (0,1) is far from each other, most of the drawn samples from the prior obtain very small weights. The probability of sampling from −3σ and lower (the parts that obtain higher weights) is nearly 1% so out of 10,000 we expected 100 non-zero weights. Posterior mean estimated equals 1.5 while we expected 0.01 (Posterior ~ Normal (0.01, 0.01). Therefore, this approach proved inefficient in terms of accuracy of estimate and number of sampling.
In contrast to IS, which does not restrict the choice of proposal distribution, RS is more restrictive. Imagine that our objective is to generate samples from a function. The essential condition is that proposal distributions should cover the target. Consider g(θ) as a proposal function from which samples were generated (Appendix Fig. 1. rectangle), f(θ) as target and K as a known constant (Appendix Fig. 1. zigzag triangle).
The algorithm is generate sample θ0 from a proposal, and calculate the likelihood of target and proposal functions at this point. If the likelihood ratio of the target to proposal is greater than a random sample generated from uniform (0, 1), we accept θ0 as a sample from the target; otherwise, we reject it and continue the algorithm. In a Bayesian context, this idea could extend to estimating the posterior distribution by considering K as the maximum likelihood of θ and the prior as the proposal. The probability of each generated sample is compared to the value of the maximum likelihood LMax,so the proposal resembles the posterior distribution (Appendix Formula-5) and (Appendix Fig. 2).
This method fails in high-dimensional space parameters due to a decrease in the acceptance rate. In addition, when Kg(θ) has considerable distance compared to f(θ), the acceptance rate would decrease drastically as well. Likewise, in a Bayesian context when LMax is not practically possible to calculate, this method fails. Armitage presented an example where, by this approach, the posterior estimation of parameters of a linear regression model had been estimated via various priors as the proposals [4], in which using a diffuse prior caused a high rejection rate and rendered the approach futile.
There are several advantages to switching from the frequentist methodology to the Bayesian. This approach can accommodate unobserved variables such as an individual's true disease status in the presence of diagnostic error. Prior distributions provide a powerful tool to apply information from previous literature and to control confounding variables. The posterior distribution, which is equivalent to the combination of prior information and the likelihood of data, is applied as an alternative to the “p-value.” Current advancements in Markov chain Monte Carlo methodology aid Bayesian analyses of complex data sets holding missing observations and multidimensional endpoints. The emergence of the MCMC approach in the 1990s led to a rapid evolution of methods to simulate posterior distributions. It is regarded as MCI through Markov chains. In contrast to approaches having a static mechanism, the MCMC follows a dynamic mechanism whereby samples are generated via a gentle transition through a target distribution function by considering a proposal, eventually converging on a stationary distribution. There are two popular Metropolis-Hastings first introduced by Metropolis in 1953, and then its special case the Gibbs sampler, discussed by Geman and Geman in 1984. Recent developments in this area have generated extensive literature. Armitage provided a neat catalogue of the references and summaries [4]. Here, let us review some technical jargon in MCMC. A chain with the property of being Markov is applied to generate dependent samples. A transition matrix illustrates the probability of moving from one state in the chain to another. It is worth mentioning that the Markov chain should have some properties that are guaranteed to produce samples from a stationary distribution. By stationary, it is meant that subsequent samples are from the same distribution. Firstly, it should be irreducible, meaning that the Markov chain can get to any state from any state within a finite number of iterations; persistent or recurrent means returning to the state at least once, and non-null means the finite mean number of transitions. These three conditions lead to a property referred to as ergodicity by which one can ensure the generation of samples from a stationary distribution. Therefore, the variance of the generated samples is estimable; otherwise, the chain may behave badly and be effectively useless. It also proves the consistency of estimates. Failure of this method occurs when there are convergence issues. For instance, in the case of a non-persistent chain, convergence to a stationary distribution never occurs. Being symmetric is another property that influences the acceptance probability of sampling. It means that g(θ ′ ′| θ′) = g(θ′| θ ′ ′) (Appendix Fig. 3).
The Metropolis-Hastings algorithm to produce a chain of samples by iterative mechanism is defined by the following steps;
When the proposal is symmetric, it reduces the proportion of the values of the target on the candidate and initial points (proposal eliminated from nominator and denominator. Now proceed to decide whether to accept the candidate as the next sample or not. So, generate a random number u from uniform (0, 1) distribution.
Step rβnewβt−1=PosteriorβnewPosteriorβt−1=Beta110.4×Binomial1040.4Beta110.5×Binomial1040.5=1.19
Step Acceptance probability α (βnew, βt−1) = min (1, r(βnew, βt−1)) = min (1, 1.19) = 1.
Step Draw a random number, u, from a Uniform (0, 1), here u = 0.345
Step If u is less than the acceptance probability, the proposed value of βnew will be accepted. Otherwise, we reject βnew and keep,βt−1. Here we accept it.
Step rβnewβt−1=PosteriorβnewPosteriorβt−1=Beta110.2×Binomial1040.2Beta110.3×Binomial1040.3=0.43
Step Acceptance probability α (βnew, βt−1) = min (1, r(βnew, βt−1)) = min (1, 0.43) = 0.43
Step3: Draw a random number, u, from a Uniform (0, 1), here u = 0.675
Step4: Since u > r, we reject it with the probability 43%.
How samples are generated by MCMC is a main question of interest. Imagine we have a uniform distribution (0, 5) as the proposal (dashed line) and a target (zigzag pattern) (Appendix Fig. 4). To show how to estimate this transition matrix, we considered a discrete target distribution as well as two numbers of states involved. Therefore, we need a 2 × 2 transition matrix whose elements are the probability of movement. To estimate the probabilities of movement, if the Markov chain is at state x, element (x, y) in the transition matrix is the probability of moving to y. For instance, consider a Markov chain that has only two possible states, {1,2}. Then the transition matrix be P(x, y)=p11p12p21p22, So there are four elements including, P(x=1, y=1)= Pr(y=1|x=1), P(x=1, y=2)= Pr(y=2|x=1), P(x=2, y=1)= Pr(y=1|x=2) and P(x=2, y=2)= Pr(y=2|x=2). According to the proposal, the probabilities of being in state 1 and 2 are 15and45, respectively. Therefore, estimated transition matrix is (x, y)=13151521545 Table 1. If the process of sampling is repeated n times, based on stochastic process the n^th^-step transition matrix will bePn=1α+pααpp⇒α=15&p=21535352525. It is proven this sampling mechanism converges to a stationary form which shows the target distribution (Appendix Fig. 4). Based on the stochastic process theorem, from this convergence time onwards each sample generated is independent of previous states [5]. From the Bayesian point of view, considering the proposal as prior distribution and target as prior × likelihood, this process of sampling is defined as Metropolis-Hasting which eventually converges to a stationary posterior.
Methods introduced in the past followed the Monte Carlo approach for computing the posterior distribution. On the contrary, DA allows approximate Bayesian analysis with a standard maximum likelihood function. Its philosophy is to translate prior information as equivalent data and add this external information to the observed study data, and then conventional methods of frequentist statistics are applied. No specific tools are required to compute the posterior mean and the variance; inverse-variance weighted averaging is a rule of thumb for estimation [6]. This technique provides an effective remedy to treat bias estimation caused by data sparseness [[7], [8], [9], [10], [11], [12]]. It considers the prior information as a penalty for maximum likelihood estimates and approximates the posterior mode and variance.
DA is a remedial tool for sparse data issues to provide an unbiased estimation of parameters. To illustrate the DA mechanism, we considered how to estimate posterior properties via inverse-variance weighting and then show the influence of prior and likelihood components on estimating posterior distributions. Finally, we depict how to construct data from the prior information to observe its role as equivalent data augmented to actual data. A hypothetical study was reported with Risk Ratio (RR) = 6.3, variance ln(RR) = 0.84 and 95% confidence limits RR = (1.02, 37.3) which was subjected to sparse data due to wide 95% limits.
Suppose that, pooled RR (95% CL) is reported in a meta-analysis study we considered as prior information as such RR with 95% limits between 13 and 3. Mean and variance of prior for Ln(RR) are estimated as
Inverse variances equaling 10.1=10 and 10.84=1.2 illustrate prior information dominated the data information by nearly 8 times. Posterior mean and variance for Ln(RR) could be estimated as the following weighted averaging rule of thumb;
Posterior RR and its 95% CL through DA provided an unbiased estimation of RR with more reasonable values of RR = 1.21 and narrower 95% CL: (1.01, 1.44). In addition, the value of the posterior mean which is closer to prior means showed the influence of the prior as well.
Various prior ranges for Ln(RR), estimated posterior 95% CL, as well as the influence of prior and data were illustrated in Table 2. Prior information weight and data information weight were estimated for each prior. This shows that for prior (16, 6) data and prior had the same influence (equal weights) while for (110, 10) it was data-dominated.
Here, compatibility of prior and data is a key issue as well. DA fails in the case of incompatibility causes misleading results [13]. For our example, Ln6.3−00.84+0.1012=1.9⇒Normal tablePvalue=0.057 showed that the compatibility hypothesis is not rejected.
DA is a kind of penalized-likelihood estimation in which the prior makes the program produce a penalty function that forces the prior constraints. DA is a straightforward and tractable approach, while the MCMC needs the convergence of the Markov Chains to the posterior distribution. It is a matter of concern, especially in the case of small samples.
Appendix Table 2 presents the advantages, failures and remedies of the approaches.
A significant body of literature published on Bayesian inference has proved its popularity among researchers. However, Bayesian inference may not be easy for non-statisticians to understand. [14]. The purpose of our paper was to provide a comprehensive framework and accompanying examples to shed light on the concepts of several Bayesian mechanisms of sampling. We showed that alternative approaches to MCMC, which are intuitively appealing and easy-to-understand, work well in the case of low-dimensional problems and appropriate prior information such as a weighted prior; otherwise, MCMC is a trouble-free tool in the case of high-dimensional problems and inappropriate priors. However, the concept of MCMC approach and its mechanisms are not easy to understand, especially for complex situations. Different studies tried to cover the Bayesian statistical approach as a need, specifically across many sciences such as Psychology, Pharmacology and so on [[15], [16], [17], [18], [19], [20], [21], [22], [23], [24], [25], [26]].
The DA method as an alternative approach to Bayesian gives researchers a more tangible sense of the role of the prior information and the data for inference making. The posterior calculation is simple using DA method as well.
In conclusion, we tried to consider the adequate approaches of the Bayesian simulation with clear examples and provide preliminary work of Bayesian foundation for the epidemiologists and other individual learners. R software code is available in the Appendix.
This research received no specific grant from any funding agency in the public, commercial, or not-for-profit sectors.
Not applicable
Not applicable
All data generated or analyzed during the current study are included in this published article.
Dr. RBY, Prof. KM, Prof. MAM and Prof. HZ had significant contributions to the conception, design, acquisition, analysis and interpretation of the information. Methodological concepts were considered Dr. RBY, Prof. KM, Prof. MAM, Prof. HZ. All authors worked on the drafting and agreed on final approval of the version to be published. Also, the agreement to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.
The authors declare that they have no conflict of interest. Mohammad Mansournia is a Senior Associate Editor on Global Epidemiology and has no access to the peer review of this manuscript.
Hojjat Zeraati, Email: zeraatih@tums.ac.ir.
Kazem Mohammad, Email: mohamadk@sina.tums.ac.ir.
We presented an example which depicted the application of Bayesian method. Three approaches of WP as representative of IS, MCMC and DA were applied. Different Prior distributions were also proposed to compare the results.
Greenland and Mansournia [11] generated observations from a binomial distribution when the number of events is y = 8 and the number of observations n = 10, and then fitted an intercept logistic model for this sparse data defined as logit(π) = β. They considered different priors for intercept β and applied DA. We used WP and MCMC with n = 1000 to estimate β as well. Results showed WP failed when Prior had very wide range (non-informative with large variance), (Table 1-R Code-3)
Fig. 1 Illustration of Rejection Sampling; The top Rectangle is proposal Kg(θ) function which samples generated from, and the Zigzag triangle is target distributionf(θ)confirming the condition f(θ)< K g(θ)for all θ. Right figure showed acceptance rate 1/4 (triangle area divided by rectangle area) and left 1/2. It depends on the choice of g(θ).
Fig. 2 Mechanism of posterior Simulation via Acceptance-Rejection sampling.
Fig. 3 Illustration of Symmetric distribution g(θ) as g (θ'′∣θ') = g (θ'∣θ'′).
Fig. 4 Proposal (dash line) versus target distribution (Zigzag) in Metropolis-Hastings algorithm of simulation with two states and symmetric proposal.
Rcodes-1: MCI
set.seed(1986)
x < −rnorm(10,000,0,1)
y < − 1/x
c(mean(y),var.(y))
R codes-2:
#Figures:
x < − seq(−10, 10, length = 100)
hx < − dnorm(x,0,0.05)
labels <− c(“Target = Norm(0,5)”,”Proposal = Norm(0,0.05)”,”Proposal = t(1)”)
sd < −c(5)
plot(x, hx, type = “l”, lty = 2, xlab = “x value”,
ylab = “Density”)
lines(x, dnorm(x,0,5), lwd = 2)
lines(x, dt(x,1),lwd = 4,lty = 3)
legend(“topright”, inset = 0.05, title = “Distributions”,
labels, lwd = 2, lty = c(1,2,3))
#Importance Sampling:
#Normal
set.seed(1986)
w < − function(x) dnorm(x, 0, 5)/dnorm(x,0,0.05)
f < − function(x) 1/x
x = rnorm(10,000,0,0.05)
W < -w(x)/sum(dnorm(x, 0, 5)/dnorm(x,0,0.05))
Y=W*f(x)
c(mean(Y), var.(Y))
w < − function(x) dnorm(x, 0, 5)/dt(x,1)
f < − function(x) 1/x
x = rt.(10,000,df = 1)
W < -w(x)/sum(dnorm(x, 0, 5)/dt(x,1))
Y = W*f(x)
c(mean(Y), var.(Y))
R codes-3
#Single parameter An Intercept model with different priors in case of data sparseness
#Weighted Prior:
X < -rep(1,10)
Y < -c(1,1,1,1,1,1,1,1,0,0)
data<−cbind(Y,X)
numSimulations<−1000
n < −10,000
output.matrix<− matrix(NA,numSimulations,7)
colnames(output.matrix) < −c(“Mean”,”median”, paste(“Quantile”,c(0.025,0.25,0.5,0.75,0.975)))
set.seed(1986)
for(i in numSimulations){
Beta<−as.matrix(rnorm(n,0,0.71))#normal(mu,sd)
logitfunc<− function(Beta) {
A < -exp(Beta%*%t(X))
w < − ifelse(Y==1, A/(1 + A), 1/(1 + A))
likelihood<−prod(w)
}
WI < -apply(Beta,1,logitfunc)
wi < −round(100n(WI/sum(WI)))
d2 < − rep(Beta,wi)
multi.fun<− function(x) {
c(mean = mean(x), median = median(x),quantile = quantile(x, c(0.025,0.25,0.5,0.75,0.975)))
}
multi.fun(d2)
output.matrix[i,] < − multi.fun(d2)
}
mean(output.matrix[,1]) #mean
mean(output.matrix[,2]) #median
mean(output.matrix[,3]) #Q2.5
mean(output.matrix[,7]) #Q97.5
sd(output.matrix[,1])/numSimulations #SE mean
plot(Beta,WI,main = “Posterior distribution via proposed method for prior Normal(0,var. = 0.5)”)
All data generated or analyzed during the current study are included in this published article.