Applications of Hierarchical Bayesian Methods to Answer Multilayer Questions with Limited Data

Abstract

There are many types of problems that include variables that are not well defined. Seeking answers to complex problems that involve many variables becomes mathematically challenging. Instead, many investigators use methods like principal component analysis to reduce the number of variables, or linear or logistic regression to rank the impact of the variables and eliminating those with the limited impact. However, eliminating variables can create a loss of integrity, especially for variables that might be associated with low likelihood but have high impact events. The use of hierarchical Bayesian methods resolves this issue by utilizing the benefits of information theory to help answer questions by incorporating a series of prior distributions for a number of variables used to solve an equation. The concept is to create distributions for the range and likelihood for each variable, and then create additional distributions to define the mean and shape values. At least three levels of analysis are required, but the hierarchical solution can include added levels beyond the initial variables (i.e., distributions related to the priors for the shape parameters). The results incorporate uncertainty, variability, and the ability to update the confidence in the values of the variables based on the receipt of new data.

1. Introduction

Suppose you have a complex question with numerous variables that are not well understood. The challenge that confronts us all is that such situations are not unusual—there are many examples of situations where there are numerous variables that can contribute to occurrences in our world, whether these occurrences involve medical issues, infrastructure issues, or science questions. Statistical methods have been developed to address limited information, but most do not permit the incorporation of new information except Bayesian methods.

The development of Bayesian methods that include priors that can be updated with new data or can respond as a result of added data overcomes initial limitations of most models. Bayesian methods developed as a result of information theory, assuming that the absolute or unconditional probability density function p(x) on X is the underlying distribution found through curve-fitting. Priors can be determined based on any combination of subjective or numeric information in the absence of real data, or as data are collected, including parameters such as the mean, variance, and range. Utilization of the observations from the prior data leads to the posterior probability function, which incorporates observations from x in the sample space S, although revealing additional information about the true content of the sample space S is subject to the influence of the proper prior distribution assumptions for x [1, 2, 3].

Many Bayesian practitioners stop with the posterior function, but the ability to develop true statistical inference requires further effort to create a predictive Bayesian solution. The Bayesian posterior methods were used in prior studies [2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12]. Predictive Bayesian methods are an extension of traditional Bayesian approaches, in which unconditional distributions for the quantity of interest are found by integrating over probabilities of parameters of the distribution for the quantity of interest, incorporating both uncertainty and variability in the quantity of interest. They have been termed “believed probabilities.”

Press [13] noted that there are advantages to the predictive Bayesian approach. Practical experience and subjectivity can be accounted for explicitly by fitting known or subjective data to a probability function that can be updated as added information becomes available [13]. Predictive Bayesian methods continually improve the statistical inference based on increased amounts of data (hence more data should advance the understanding of statistical relationships and provide greater confidence in the prior and therefore the solution). The predicted distributions are also important for checking goodness of fit of the resulting predictive model to actual data. However, the analysis can become problematic when the information is so scarce that the analysis yields nothing useful [11].

The use of Monte Carlo methods makes the solutions easier. Through randomized sampling, the resulting predictions are simulated from the posterior predictive distribution, which is the distribution of the unobserved future results based on prior observed data. The more confidence that exists with the priors, the more likely results of likely outcomes can be derived.

However, many times Bayesian methods have been limited to situations where there are one or two variables that contribute to an outcome. While high-quality answers can be derived, the design of the algorithm often oversimplifies the real world where many variables may contribute to the outcome. The ability to incorporate many variables that are unknown or uncertain makes the calculations intractable in the traditional Bayesian processes.

Equally important is the ability to study the events, where multiple variables may impact the probability and impact of a given consequence. Regression models are often used, along with principal component analysis to address such situations. However, both rely on complete data sets, and for many situations, the lack of complete data may be extensive. Examples include much of the public or municipal infrastructure that we rely on so heavily for a functioning society, health risks or impacts, natural disaster risk, and extreme event prediction. This is where hierarchical application to Bayesian methods has value.

2. Methods

Multiple variables create greater challenges. The concept is to create a distribution with shape and location parameters for each variable, then create distributions to define the mean and shape values (priors). Each of those distributions can be solved with MCMC methods to create a predictive distribution. That predictive distribution can be sampled and multiplied (or added) to the results of all other variables solved and sampled similarly. The assignment of distributions need not stop with the priors for shape location and location. Those parameters also can be assigned a distribution, and likewise those (and so on). Model parameters create a structure or a series of levels that look like a hierarchy, whereby the priors of a given solution are dependent on priors to those variable distributions and integrate that information across levels simultaneously [12], thereby separating the observed variability into parts attributable to both random and true differences [13, 14, 15]. Each investigation will have a different series of variables, each with different associated variables, priors, and priors of the priors.

Allenby et al. [16] stated that hierarchical Bayesian models are really “the combination of two things: i) a model written in hierarchical form that is ii) estimated using Bayesian methods.” Shaddick et al. [17] consider there are to be at least three levels: (1) the observation or measurement level, (2) the underlying process level, and (3) the parameter level. Kruschke and Vanpaemel [18] noted that hierarchical Bayesian data analysis involves “describing data by meaningful mathematical models and allocating credibility to parameter values that are consistent with the data and with prior knowledge.” Using Bayesian methods, hierarchical Bayesian models can yield estimates of the true effects at each level of the hierarchy [14, 19]. By considering the results across all levels, hierarchical Bayesian models can be used to rigorously integrate information with a complex underlying structure [14], resulting in a tendency to shrink differences when multiple variables are incorporated [20]. An important aspect of the hierarchical approach is that the model is usually a flexible version of a base model [21], and if needed, the models allow for adding extra levels depending on the hyperparameters [22].

Applying Bayesian prediction and weighting in a unified approach to Bayesian regression models can account for complex design features under the framework of multilevel regression and poststratification [23, 24, 25]. Weighting in a hierarchical model can be used as an extension of linear or logistic regression models. Methods for hierarchical functional data typically require that all curves are observed over or standardized to fall in the same region [26, 27, 28]. While classical weighting usually relies on many user-defined choices for regression that is difficult to codify [29], the Bayesian approach allows prior information to be incorporated and the distributions automatically adjusted [30, 31].

MCMC allows the user to approximate aspects of posterior distributions that cannot be directly calculated (e.g., random samples from the posterior, posterior means, etc.). There are examples of current applications of this approach. Draper [32] considered Bayesian hierarchical Poisson regression models, Wang et al. [33] created hierarchical Bayesian model developed for predicting monthly residential per capita electricity consumption at the state level across the United States, and Maddala et al. [34] studied the relationship of income elasticity on energy demand in the United States by applying a dynamic linear regression model under Bayesian framework. Roman et al. [35] and Neil and Fenton [36] used hierarchical Bayesian model for evaluation of treatments for Covid-19.

There are limits to hierarchical Bayesian model. The first is the underlying resulting model assumption may be wrong [14, 31]. Thais et al. [37] noted that “ill behaved likelihoods” at the lower levels in the hierarchy may create either excessive concentration about a mean or noninformative results. Rouder et al. [38] note that although hierarchical linear models are suitable in several domains, they rarely make good models of psychological process. Heller and Gharamani [39] note the algorithm provides no guide to choosing the “correct” number of clusters.

As a result, Shaddick et al. [18] note that Bayesian hierarchical models are an extremely useful and flexible framework in which to model complex relationships and dependencies in data, while Kruschke and Vanpaemel [17] suggest they provide flexibility in designing models that are appropriate for describing the data at hand that can provide a complete representation of parameter uncertainty (i.e., the posterior distribution) that can be directly interpreted. Some examples on how to create a hierarchical Bayesian model are helpful to demonstrate the process.

3. Applications

3.1 Example 1: Drake equation

Since 1959, SETI has yet to find an alien signal. Two questions arise as a result—what is the probability of there being life in the galaxy, and why have not we received a response to our transmissions? In 1959, a US astronomer, Frank D. Drake, a NASA employee who carried out the first SETI radio telescope experiments, outlined an equation for finding communicable civilizations [40, 41]:

$$N=RL$$E1

$\text{Which was later expanded to:}\mathit{R}{f}_p{n}_e{f}_i{f}_t{f}_c\mathit{L}$.

where N is the number of communicable civilizations, R is the rate at which stars are born in the galaxy, f_p is the fraction of stars with planetary systems, n_e is the number of planets that might hold life, f_t is the fraction of planets with life, f_i is the fraction of planets with life that have evolved, f_c is the number of civilizations of evolved civilizations with the ability to communicate, and L is the length of time over which the communication is possible. An additional factor named C is a recent suggestion for colonization [42]. Bloetscher [3] suggested that the factor n_e actually comprises four factors: planet size (PS), presence of a moon (M), location within the “Goldilocks” or habitable zone (HZ), and the correct star type (ST), creating four unknowns from one. None of the 7–12 factors is fully known, so no specific answer on the likelihood of intelligent life on another planet communicating with Earth is possible. However, hierarchical Bayesian methods can be used to investigate the probability of intelligent life on another planet communicating with Earth. This approach involves the assignment of probability distributions to the underlying factors and using those to develop an MCMC protocol to determine the final predictive solution [3]. Subjective data, shown in Table 1, were included when little or no data were available to specify the parameters of these distributions [3]. Then probability distributions were assigned to the prior parameters within the initial distributions to determine the location and scale parameters of the factor distribution (Table 2). The subjective information serves to create these prior distributions until such time as real data are developed or become available [3].

Variable description	Source of variable	Variable	Drake and Sobel (1991) estimate	Diehl et al. (2006) or Maccone (2010) estimate	Bloetscher mean from literature	Range of potential values	Distribution type	Prior1 distribution	Prior 1 estimate	Prior2 distribution	Prior 2 estimate	Data points
Number of Starts forming each year	Drake and Soble, 1991	R	10	7	20.0	0 to infinity	Gamma	Gamma	4	Gamma	0.2	Numerous, but actually only 1
Percent of stars with planets	Drake and Soble, 1991	fp	0.5	0.5	0.7	0 to 1	Beta	Uniform	1	Uniform	1.4	1
Number of Planets with life	Drake and Soble, 1991	ne	2	1	n/a	0 to infinity	Gamma	Gamma	1	Gamma	0.125	1
Percent of stars of right type	Bloetscher, 2019	ST	n/a	n/a	0.7	0 to 1	Beta	Uniform	1.4	Uniform	2	Thousands, and we find them with planets
Percent of planets of the right size	Bloetscher, 2019	PS	n/a	n/a	0.25	0 to 1	Beta	Uniform	1	Uniform	4	8
Percent of planets in Habitable zone	Bloetscher, 2019	HZ	n/a	n/a	0.375	0 to 1	Beta	Uniform	1	Uniform	2.33	8
Percent of planets with moons to create motion	Bloetscher, 2019	M	n/a	n/a	0.5	0 to 1	Beta	Uniform	1	Uniform	2	Dozens
Percent of planets where life forms	Drake and Soble, 1991	fl	1	0.5	0.5	0 to 1	Beta	Uniform	0.5	Uniform	0.5	1
Percent of planets with life that becomes intelligence	Drake and Soble, 1991	fi	0.01	0.2	0.5	0 to 1	Uniform	n/a		n/a		1
Percent of planets that develop technology to communicate	Drake and Soble, 1991	fc	0.01	0.2	0.2	0 to 1	Beta	Uniform	0.5	Uniform	0.5	1
Lifetime of civilization that wants to communicate	Drake and Soble, 1991	L	10,000	10,000	4600	0 to infinity	Exponential	Gamma	1	n/a		1
Percent of civilizations with a desire to colonize	Walters et al. 1980	C	n/a	n/a	0.2	0 to 1	Beta	Uniform	1	Uniform	5	1
	Mean of N		10	700	21

Table 1.

Summary of Drake parameters used in Bayesian calculations and comparison to prior estimates.

Parameter	Mean (10,000-samples)
R	31.89
fp	0.4392
ST	0.4433
PS	0.2586
HZ	0.3170
M	0.3768
fl	0.5013
fi	0.5016
fc	0.4956
L	4641
C	0.2219
MCMC | Average for N	24.64

Table 2.

Monte Carlo results for parameters used in MCMC.

For the Drake equations, a distribution for N was developed through using the Hierarchical Monte Carlo distributions for the factors of the equation run 10,000 times (see Figure 1). A series of Hierarchical Monte Carlo algorithms were developed for each parameter, and the means were inserted into a Monte Carlo Markov Chain program that uses a Metropolis-Hastings algorithm with a Gibbs sampler to develop a final probabilistic result [43, 44, 45, 46, 47, 48, 49] that was solved for N. Based on suggestions by Glade et al. [50] and Maccone [51], the target MCMC distribution was proposed to be log-normal. Given the uncertainly involved, the standard deviation used for the target distribution was assumed to be the square root of 6, after Wu et al. [52]. Given a multivariate distribution, like the example above, Gibbs sampling breaks down the problem by drawing samples for each parameter directly from that parameter’s conditional distribution or the probability distribution of a parameter given a specific value of another parameter [53].

Figure 1.

Development of HPB for Drake.

The solution is shown in Figure 2 (red data points). Of importance, there is nearly a 50% probability that we are alone in the galaxy. The graph indicates that there is a 95% probability that there are less than 100 communicating civilizations concurrent with Earth, and a 99% that there are 1000 such civilizations.

Figure 2.

Results for the probability of N using the Monte Carlo calculation of the factors placed into a MCMC predictive Bayesian MATLAB calculation. Note that the number of other planets is constantly small despite three options as noted above (reproduced from Bloetcher [3]).

3.2 Example 2: dose response

The use of predictive Bayesian methods for dose–response relationships has also been investigated by a number of authors [10, 11, 12, 54, 55, 56]. Beaudequin et al. [57] developed QMRA with the use of hierarchical Bayesian networks to address the data paucity, combine quantitative and qualitative information including expert opinion, and the ability to offer a systems approach to characterize complexity. They outlined how the Bayesian networks are the current method of choice for determining the risk to human health from exposure to pathogens because of their ability to separate risk and uncertainly, predict outcomes, and deal with poorer quality data [57]. Hence, as subjective data are incorporated, the prior distributions self-adjust [58]. Bloetscher et al. [9] used six sets of Cryptosporidium data to show how the dose–response function changes with new, additional data. As a result of new data, the dose–response is expected to improve, demonstrating that the process can be applied to other organisms. In addition, the paper creates a Predictive Bayesian MCMC solution for the Pareto II distribution with two uncertain parameters.

Given the unlikelihood of reinfection during a single incident (due to a short period of time), the likelihood of infection can be described by the binomial distribution. As such, a binomial function is used to represent the probability of exposure. Figure 3 shows the conceptual model with three levels of probability distributions. Because of the intrinsic difficulty in solving a predictive Bayesian equation with multiple embedded distributions through double integration, an analytical mathematical solution is not achievable. Instead, a probabilistic solution was developed using a Markov Chain Monte Carlo (MCMC) program developed in MATLAB18® with uncertain values for a and k. Six different models of 10,000 iterations were run, each model including an additional dataset and the prior for α increased to account for the additional data.

Figure 3.

Conceptual risk model for cryptosporidium.

When compared with the beta-Poisson models developed by Haas et al. [59], the predictive Bayesian equation derived in this study is less conservative by a factor of over 10 than the beta-Poisson model used by Haas et al. [59] (see Figure 4). However, the beta-Poisson does not accept new data, and therefore cannot be updated, is the likely explanation for the difference.

Figure 4.

Dose–response function for C. Parvum—All six datasets compared with beta-Poisson. The beta-Poisson is the most conservative of the sets. Once the last three datasets (TAMU, UC, and More) are added; all provide the same graph (reproduced from Bloetscher, et al. 2020).

3.3 Example 3: infrastructure

Public water and sewer utility systems are created to develop safe, reliable, and financially self-supporting potable water and sanitary sewage systems, which will meet the water and sewerage needs of the areas served by the utility, to ensure that existing and future utility facilities are constructed, operated, and managed with high reliability and are compatible with the area’s future growth. To gain efficiencies in operation, these new facilities must be developed in accordance with the latest technical and professional standards to protect the health, safety, and welfare of the citizens served now or in the future.

Public infrastructure has been poorly rated by the American Society of Civil Engineers [60, 61, 62, 63, 64, 65], and most public officials acknowledge the deterioration of the infrastructure we rely on daily. Part of the challenge is that many jurisdictions have limited information about their systems, and little data to use to justify spending of specific project. Hence the infrastructure tends to deteriorate further each year as local officials opt to limit budgets in the absence of good needs data. As a result, state and local governments currently spend about 1.8% of its GNP on infrastructure, as compared with 3.1% in 1970 [66]. Twice as much was spent 40 years ago, and a large portion of today’s costs are for growth as opposed to repair and replacement. Asset management is supposed to help this meet this challenge.

An asset management program consists of determining the selected area of study, type of system, and the quality of data used for evaluation. The question is how to collect data that might be useful to a utility that does not involve a lot of destructive testing on buried infrastructure that is costly and inconvenient. When creating an asset management plan, missing data are perceived to be a huge problem, especially when the event data (breaks in pipe as an example) are not tracked. The lack of tracking makes it difficult to determine which factors are the critical ones. Many utilities lack the resources for examining buried infrastructure, so other methods of data collection are needed. The concept in Bloetscher et al. [67] was to develop a means to acquire data on the assets for a condition assessment (buried pipe is not visible and cannot really be assessed). What was found was that for buried infrastructure, much more information was known than anticipated. For one thing, most utilities have a pretty good idea about the pipe materials. Employee memory can be very useful, even if not completely accurate. In most cases, the depth of pipe is fairly similar—the deviations may be known. Soil conditions may be useful—there is an indication that aggressive soil causes more corrosion in ductile iron pipe, and most soil information is readily available. Groundwater is usually known, and if a saltwater interface of a pollution plume exists, it can be mapped and evaluated for impact on pipe. Tree roots will wrap around water and sewer pipes, so their presence is detrimental. Trees are easily noted from aerial photographs. Roads with heavy truck traffic create more vibrations in the soil, causing rocks to move toward the pipe and joints to flex. So, with a little research, there are at least six variables known.

All variable information can be compiled into tables. There is also a need to track events or consequences—breaks, flooding etc.—that would indicate a failure, which is required for predicting future maintenance needs and the most at-risk assets. Finally, the data along with the consequence can be used to predict where the breaks might occur in the future based on past experience. If the break history for a water system, flood records for a stormwater system, or sewer pipe condition from televising is known, the impact of these factors can be developed via a linear regression algorithm. For logistic of linear regression, XLStat® can be used for the statistical analysis. The linear regression algorithm can then be used as a predictive tool to help identify assets that are mostly likely to become a problem.

Data need to be kept up as things change, but exact data are not needed. An example of this type of effort is shown for a medium-sized city in Florida in Figures 5–7. The City’s GIS system was mined for the purposes of this project. Data were retrieved and reviewed to address missing data and clear errors. Nearly 10,000 pipe segments remained. Categorical information on trees, vibrations, soil type, and pipe type is added. Noncategorical data for pipe size, length, and age were also entered. Note that with 10,000 pipe sections and less than 600 breaks, many pipes have no breaks in their history. The linear regression function for XLStat® was used to create equation to identify the factors associated with each variable and the amount of influence that each exerts (see Figure 5). In this case, the equation was:

Figure 5.

Impact of factors on leaks.

Figure 6.

Comparison of predictive and actual breaks over 10 years (correlation desirable).

Figure 7.

Pipe risk—Red pipe is the highest risk for this community, while blue pipes are the lowest risk.

$$\text{Breaks}=-3.5442\text{7E}-03-6.5187\mathrm{E}-03\ast \mathrm{DIA}+2.60\text{7E}-03\ast \mathrm{Age}$$E2

It should be noted that this utility has three main types of pipe, installed at three completely different eras. Because the correlation between pipe type and age was high, and likewise pipe type and diameters, other factors that might impact leaks in other communities were not obvious, so other communities would need to recreate this analysis for their situation.

Figure 6 outlines how the predictive equation correlated for the City’s potable water distribution system (well within one standard deviation). Figure 7 is a GIS map of pipe vulnerability based on the data. Red pipe is the highest priority to schedule for replacement.

The concept should apply to any utility, although the results and factors of concern will be slightly different for each utility. Also, in smaller communities, many variables (ductile iron pipe, PVC pipe, soil condition…) may be so similar that attempts to differentiate factors may be unproductive.

The analysis indicated two things—that age and AC pipe were correlated.

But what if none of this information is fully known? Many of the indicators of failure can be tracked through the information that is required to be included in the as-built drawings, but what if they are not available? Loss on institutional knowledge through retirements can cost the utility much information on actual pipe diameter, pipe depth, age, and breaks given many utilities do not have extensive work order systems. Other information that might be useful is condition that maintenance crews may have knowledge of. A hierarchical Bayesian model could be developed to address these concerns. Where the pipe is actually known, the categorical variable would be set to 1. Otherwise, a beta distribution could be developed with a “confidence mean”—we think it is ductile iron, but it might be PVC or cast iron. The same with pipe diameters, etc. As new variables are developed, confidence could be added and priors adjusted. Criticality could be a distribution as well. Figure 8 shows what an infrastructure assessment Hierarchical Predictive Bayesian model might look like (realizing it might extend far more widely). Currently, research is underway to develop such models, but the data required to create and utilize the models are often lacking even in the most sophisticated organizations.

Figure 8.

Partial diagram of potential hierarchical predictive Bayesian infrastructure risk model.

4. Conclusions

Going back to the beginning of the chapter, “What to do when you have a complex question with numerous variables that are not well understood?” It would appear that the use of hierarchical predictive Bayesian models is a solution to address the challenge. While there may be circumstances where these methods may not work (psychology), for issues such as infrastructure and completely unknown questions such as the Drake equation, the methods seem ideally situated to shed light on the solution in a probabilistic form. The outcome of these methods provides a probability of a given answer—not a specific answer—at different levels of confidence. Uncertainty and variability are by the nature of being a probabilistic answer already included in the solutions. This is why added data can improve the likelihood of a given solution while reducing the potential for less likely solutions. Like the Drake equation solution, the dose–response example is an example of such a process. Infrastructure would be as well—as more data are created, the solutions become more robust and uncertainty is reduced. The results permit us to make better decisions as the data improve our understanding.