Bayesian phylogenetic model selection using reversible jump Markov chain Monte Carlo

A common problem in molecular phylogenetics is choosing a model of DNA substitution that does a good job of explaining the DNA sequence alignment without introducing superfluous parameters. A number of methods have been used to choose among a small set of candidate substitution models, such as the likelihood ratio test, the Akaike Information Criterion (AIC), the Bayesian Information Criterion (BIC), and Bayes factors. Current implementations of any of these criteria suffer from the limitation that only a small set of models are examined, or that the test does not allow easy comparison of non-nested models. In this article, we expand the pool of candidate substitution models to include all possible time-reversible models. This set includes seven models that have already been described. We show how Bayes factors can be calculated for these models using reversible jump Markov chain Monte Carlo, and apply the method to 16 DNA sequence alignments. For each data set, we compare the model with the best Bayes factor to the best models chosen using AIC and BIC. We find that the best model under any of these criteria is not necessarily the most complicated one; models with an intermediate number of substitution types typically do best. Moreover, almost all of the models that are chosen as best do not constrain a transition rate to be the same as a transversion rate, suggesting that it is the transition/transversion rate bias that plays the largest role in determining which models are selected. Importantly, the reversible jump Markov chain Monte Carlo algorithm described here allows estimation of phylogeny (and other phylogenetic model parameters) to be performed while accounting for uncertainty in the model of DNA substitution.     

QTL mapping in outbred half-sib families using Bayesian model selection

In this article, we propose a model selection method, the Bayesian composite model space approach, to map quantitative trait loci (QTL) in a half-sib population for continuous and binary traits. In our method, the identity-by-descent-based variance component model is used. To demonstrate the performance of this model, the method was applied to map QTL underlying production traits on BTA6 in a Chinese half-sib dairy cattle population. A total of four QTLs were detected, whereas only one QTL was identified using the traditional least square (LS) method. We also conducted two simulation experiments to validate the efficiency of our method. The results suggest that the proposed method based on a multiple-QTL model is efficient in mapping multiple QTL for an outbred half-sib population and is more powerful than the LS method based on a single-QTL model.     

Comparing vector–host and SIR models for dengue transmission

Various simple mathematical models have been used to investigate dengue transmission. Some of these models explicitly model the mosquito population, while others model the mosquitoes implicitly in the transmission term. We study the impact of modeling assumptions on the dynamics of dengue in Thailand by fitting dengue hemorrhagic fever (DHF) data to simple vector–host and SIR models using Bayesian Markov chain Monte Carlo estimation. The parameter estimates obtained for both models were consistent with previous studies. Most importantly, model selection found that the SIR model was substantially better than the vector–host model for the DHF data from Thailand. Therefore, explicitly incorporating the mosquito population may not be necessary in modeling dengue transmission for some populations.     

Bayesian Estimation of Positively Selected Sites

In protein-coding DNA sequences, historical patterns of selection can be inferred from amino acid substitution patterns. High relative rates of nonsynonymous to synonymous changes (=d _N /d _S ) are a clear indicator of positive, or directional, selection, and several recently developed methods attempt to distinguish these sites from those under neutral or purifying selection. One method uses an empirical Bayesian framework that accounts for varying selective pressures across sites while conditioning on the parameters of the model of DNA evolution and on the phylogenetic history. We describe a method that identifies sites under diversifying selection using a fully Bayesian framework. Similar to earlier work, the method presented here allows the rate of nonsynonymous to synonymous changes to vary among sites. The significant difference in using a fully Bayesian approach lies in our ability to account for uncertainty in parameters including the tree topology, branch lengths, and the codon model of DNA substitution. We demonstrate the utility of the fully Bayesian approach by applying our method to a data set of the vertebrate -globin gene. Compared to a previous analysis of this data set, the hierarchical model found most of the same sites to be in the positive selection class, but with a few striking exceptions.     

Bayesian point estimation of quantitative trait loci

In this article, we consider the problem of the estimation of quantitative trait loci (QTL), those chromosomal regions at which genetic information affecting some quantitative trait is encoded. Generally the number of such encoding sites is unknown, and associations between neutral molecular marker genotypes and observed trait phenotypes are sought to locate them. We consider a Bayesian model for simple experimental designs, and discuss the existing approaches to inference for this problem. In particular, we focus on locating positions of the best candidate markers segregating for the trait, a situation which is of primary interest in comparative mapping. We introduce a loss function for estimating both the number of QTL and their location, and we illustrate its application via simulated and real data.     

bgc: Software for Bayesian estimation of genomic clines

Introgression in admixed populations can be used to identify candidate loci that might underlie adaptation or reproductive isolation. The Bayesian genomic cline model provides a framework for quantifying variable introgression in admixed populations and identifying regions of the genome with extreme introgression that are potentially associated with variation in fitness. Here we describe the bgc software, which uses Markov chain Monte Carlo to estimate the joint posterior probability distribution of the parameters in the Bayesian genomic cline model and designate outlier loci. This software can be used with next‐generation sequence data, accounts for uncertainty in genotypic state, and can incorporate information from linked loci on a genetic map. Output from the analysis is written to an HDF5 file for efficient storage and manipulation. This software is written in C++ . The source code, software manual, compilation instructions and example data sets are available under the GNU Public License at http://sites.google.com/site/bgcsoftware/ .     

On the Bayesian analysis of ring-recovery data

Vounatsou and Smith (1995, Biometrics 51, 687-708) describe the modern Bayesian analysis of ring-recovery data. Here we discuss and extend their work. We draw different conclusions from two major data analyses. We emphasize the extreme sensitivity of certain parameter estimates to the choice of prior distribution and conclude that naive use of Bayesian methods in this area can be misleading. Additionally, we explain the discrepancy between the Bayesian and classical analyses when the likelihood surface has a flat ridge. In this case, when there is no unique maximum likelihood estimate, the Bayesian estimators are remarkably precise.     

On the Bayesian estimation of a closed population size in the presence of heterogeneity and model uncertainty

Summary .   We consider the estimation of the size of a closed population, often of interest for wild animal populations, using a capture–recapture study. The estimate of the total population size can be very sensitive to the choice of model used to fit to the data. We consider a Bayesian approach, in which we consider all eight plausible models initially described by Otis et al. (1978, Wildlife Monographs 62, 1–135) within a single framework, including models containing an individual heterogeneity component. We show how we are able to obtain a model-averaged estimate of the total population, incorporating both parameter and model uncertainty. To illustrate the methodology we initially perform a simulation study and analyze two datasets where the population size is known, before considering a real example relating to a population of dolphins off northeast Scotland.     

A stochastic formulation of the gompertzian growth model for in vitro bactericidal kinetics: parameter estimation and extinction probability

Time-kill curves have frequently been employed to study the antimicrobial effects of antibiotics. The relevance of pharmacodynamic modeling to these investigations has been emphasized in many studies of bactericidal kinetics. Stochastic models are needed that take into account the randomness of the mechanisms of both bacterial growth and bacteria-drug interactions. However, most of the models currently used to describe antibiotic activity against microorganisms are deterministic. In this paper we examine a stochastic differential equation representing a stochastic version of a pharmacodynamic model of bacterial growth undergoing random fluctuations, and derive its solution, mean value and covariance structure. An explicit likelihood function is obtained both when the process is observed continuously over a period of time and when data is sampled at time points, as is the custom in these experimental conditions. Some asymptotic properties of the maximum likelihood estimators for the model parameters are discussed. The model is applied to analyze in vitro time-kill data and to estimate model parameters; the probability of the bacterial population size dropping below some critical threshold is also evaluated. Finally, the relationship between bacterial extinction probability and the pharmacodynamic parameters estimated is discussed.     

Bayesian predictive information criterion for the evaluation of hierarchical Bayesian and empirical Bayes models

The problem of evaluating the goodness of the predictive distributionsof hierarchical Bayesian and empirical Bayes models is investigated.A Bayesian predictive information criterion is proposed as anestimator of the posterior mean of the expected loglikelihoodof the predictive distribution when the specified family ofprobability distributions does not contain the true distribution.The proposed criterion is developed by correcting the asymptoticbias of the posterior mean of the loglikelihood as an estimatorof its expected loglikelihood. In the evaluation of hierarchicalBayesian models with random effects, regardless of our parametricfocus, the proposed criterion considers the bias correctionof the posterior mean of the marginal loglikelihood becauseit requires a consistent parameter estimator. The use of thebootstrap in model evaluation is also discussed.     

Influences of error distributions of net ecosystem exchange on parameter estimation of a process-based terrestrial model: A case of broad-leaved Korean pine mixed forest in Changbaishan, China

Model predictions can be improved by parameter estimation from measurements. It was assumed that measurement errors of net ecosystem exchange (NEE) of CO₂ follow a normal distribution. However, recent studies have shown that errors in eddy covariance measurements closely follow a double exponential distribution. In this paper, we compared effects of different distributions of measurement errors of NEE data on parameter estimation. NEE measurements in the Changbaishan forest were assimilated into a process-based terrestrial ecosystem model. We used the Markov chain Monte Carlo method to derive probability density functions of estimated parameters. Our results showed that modeled annual total gross primary production (GPP) and ecosystem respiration (Re) using the normal error distribution were higher than those using the double exponential distribution by 61–86 gC m^?2 a^?1 and 107–116 gC m^?2 a^?1, respectively. As a result, modeled annual sum of NEE using the normal error distribution was lower by 29–47 gC m^?2 a^?1 than that using the double exponential error distribution. Especially, modeled daily NEE based on the normal distribution underestimated the strong carbon sink in the Changbaishan forest in the growing season. We concluded that types of measurement error distributions and corresponding cost functions can substantially influence the estimation of parameters and carbon fluxes.     

Bayesian analysis of a time series of counts with covariates: an application to the control of an infectious disease

This paper presents a Bayesian analysis of a time series of counts to assess its dependence on an explanatory variable. The time series represented is the incidence of the infectious disease ESBL-producing Klebsiella pneumoniae in an Australian hospital and the explanatory variable is the number of grams of antibiotic (third generation) cephalosporin used during that time. We demonstrate that there is a statistically significant relationship between disease occurrence and use of the antibiotic, lagged by three months. The model used is a parameter-driven model in the form of a generalized linear mixed model. Comparison of models is made in terms of mean square error.     

Stochastic search gene suggestion: a Bayesian hierarchical model for gene mapping

Mapping the genes for a complex disease, such as diabetes or rheumatoid arthritis (RA), involves finding multiple genetic loci that may contribute to the onset of the disease. Pairwise testing of the loci leads to the problem of multiple testing. Looking at haplotypes, or linear sets of loci, avoids multiple tests but results in a contingency table with sparse counts, especially when using marker loci with multiple alleles. We propose a hierarchical Bayesian model for case-parent triad data that uses a conditional logistic regression likelihood to model the probability of transmission to a diseased child. We define hierarchical prior distributions on the allele main effects to model the genetic dependencies present in the human leukocyte antigen (HLA) region of chromosome 6. First, we add a hierarchical level for model selection that accounts for both locus and allele selection. This allows us to cast the problem of identifying genetic loci relevant to the disease into a problem of Bayesian variable selection. Second, we attempt to include linkage disequilibrium as a covariance structure in the prior for model coefficients. We evaluate the performance of the procedure with some simulated examples and then apply our procedure to identifying genetic markers in the HLA region that influence risk for RA. Our software is available on the website http://www.epigenetic.org/Linkage/ssgs-public/.     

Exploratory Bayesian model selection for serial genetics data

Characterizing the process by which molecular and cellular level changes occur over time will have broad implications for clinical decision making and help further our knowledge of disease etiology across many complex diseases. However, this presents an analytic challenge due to the large number of potentially relevant biomarkers and the complex, uncharacterized relationships among them. We propose an exploratory Bayesian model selection procedure that searches for model simplicity through independence testing of multiple discrete biomarkers measured over time. Bayes factor calculations are used to identify and compare models that are best supported by the data. For large model spaces, i.e., a large number of multi-leveled biomarkers, we propose a Markov chain Monte Carlo (MCMC) stochastic search algorithm for finding promising models. We apply our procedure to explore the extent to which HIV-1 genetic changes occur independently over time.     

A comparison of three models of the inositol trisphosphate receptor

The inositol (1,4,5)-trisphosphate receptor (IPR) plays a crucial role in calcium dynamics in a wide range of cell types, and is often a central feature in quantitative models of calcium oscillations and waves. We compare three mathematical models of the IPR, fitting each of them to the same data set to determine ranges for the parameter values. Each of the fits indicates that fast activation of the receptor, followed by slow inactivation, is an important feature of the model, and also that the speed of inositol trisphosphate (IP₃) binding cannot necessarily be assumed to be faster than Ca²⁺ activation. In addition, the model which assumed saturating binding rates of Ca²⁺ to the IPR demonstrated the best fit. However, lack of convergence in the fitting procedure indicates that responses to step increases of [Ca²⁺] and [IP₃] provide insufficient data to determine the parameters unambiguously in any of the models.     

Modeling the euglycemic hyperinsulinemic clamp by stochastic differential equations

The Euglycemic Hyperinsulinemic Clamp (EHC) is the most widely used experimental procedure for the determination of insulin sensitivity. In the present study, 16 subjects with BMI between 18.5 and 63.6 kg/m² have been studied with a long-duration (5 hours) EHC. In order to explain the oscillations of glycemia occurring in response to the hyperinsulinization and to the continuous glucose infusion at varying speeds, we first hypothesized a system of ordinary differential equations (ODEs), with limited success. We then extended the model and represented the experiment using a system of stochastic differential equations (SDEs). The latter allow for distinction between (i) random variation imputable to observation error and (ii) system noise (intrinsic variability of the metabolic system), due to a variety of influences which change over time. The stochastic model of the EHC was fitted to data and the system noise was estimated by means of a (simulated) maximum likelihood procedure, for a series of different hypothetical measurement error values. We showed that, for the whole range of reasonable measurement error values: (i) the system noise estimates are non-negligible; and (ii) these estimates are robust to changes in the likely value of the measurement error. Explicit expression of system noise is physiologically relevant in this case, since glucose uptake rate is known to be affected by a host of additive influences, usually neglected when modeling metabolism. While in some of the studied subjects system noise appeared to only marginally affect the dynamics, in others the system appeared to be driven more by the erratic oscillations in tissue glucose transport rather than by the overall glucose-insulin control system. It is possible that the quantitative relevance of the unexpressed effects (system noise) should be considered in other physiological situations, represented so far only with deterministic models.The work was supported by grants from the Danish Medical Research Council and the Lundbeck Foundation to S. Ditlevsen.