Application of Gibbs sampling for inference in a mixed major gene-polygenic inheritance model in animal populations

The application of Gibbs sampling is considered for inference in a mixed inheritance model in animal populations. Implementation of the Gibbs sampler on scalar components, as used for human populations, appeared not to be efficient, and an approach with blockwise sampling of genotypes was proposed for use in animal populations. The blockwise sampling of genotypes was proposed for use in animal populations. The blockwise sampling by which genotypes of a sire and its final progeny were sampled jointly was effective in improving mixing, although further improvements could be looked for. Posterior densities of parameters were visualised from Gibbs samples; from the former highly marginalised Bayesian point and interval estimates can be obtained.     

A bayesian approach to the multistate Jolly-Seber capture-recapture model

This article considers a Bayesian approach to the multistate extension of the Jolly-Seber model commonly used to estimate population abundance in capture-recapture studies. It extends the work of George and Robert (1992, Biometrika79, 677-683), which dealt with the Bayesian estimation of a closed population with only a single state for all animals. A super-population is introduced to model new entrants in the population. Bayesian estimates of abundance are obtained by implementing a Gibbs sampling algorithm based on data augmentation of the missing data in the capture histories when the state of the animal is unknown. Moreover, a partitioning of the missing data is adopted to ensure the convergence of the Gibbs sampling algorithm even in the presence of impossible transitions between some states. Lastly, we apply our methodology to a population of fish to estimate abundance and movement.     

Comparison of REML and Gibbs sampling estimates of multi-trait genetic parameters in Scots pine

Multi-trait (co)variance estimation is an important topic in plant and animal breeding. In this study we compare estimates obtained with restricted maximum likelihood (REML) and Bayesian Gibbs sampling of simulated data and of three traits (diameter, height and branch angle) from a 26-year-old partial diallel progeny test of Scots pine (Pinus sylvestris L.). Based on the results from the simulated data we can conclude that the REML estimates are accurate but the mode of posterior distributions from the Gibbs sampling can be overestimated depending on the level of the heritability. The mean and median of the posteriors were considerably higher than the expected values of the heritabilities. The confidence intervals calculated with the delta method were biased downwardly. The highest probablity density (HPD) interval provides a better interval estimate, but could be slightly biased at the lower level. Similar differences between REML and Gibbs sampling estimates were found for the Scots pine data. We conclude that further simulation studies are needed in order to evaluate the effect of different priors on (co)variance components in the genetic individual model.     

Conjugate Gibbs sampling for Bayesian phylogenetic models.

We propose a new Markov Chain Monte Carlo (MCMC) sampling mechanism for Bayesian phylogenetic inference. This method, which we call conjugate Gibbs, relies on analytical conjugacy properties, and is based on an alternation between data augmentation and Gibbs sampling. The data augmentation step consists in sampling a detailed substitution history for each site, and across the whole tree, given the current value of the model parameters. Provided convenient priors are used, the parameters of the model can then be directly updated by a Gibbs sampling procedure, conditional on the current substitution history. Alternating between these two sampling steps yields a MCMC device whose equilibrium distribution is the posterior probability density of interest. We show, on real examples, that this conjugate Gibbs method leads to a significant improvement of the mixing behavior of the MCMC. In all cases, the decorrelation times of the resulting chains are smaller than those obtained by standard Metropolis Hastings procedures by at least one order of magnitude. The method is particularly well suited to heterogeneous models, i.e. assuming site-specific random variables. In particular, the conjugate Gibbs formalism allows one to propose efficient implementations of complex models, for instance assuming site-specific substitution processes, that would not be accessible to standard MCMC methods.     

Predicting and correcting bias caused by measurement error in line transect sampling using multiplicative error models

Line transect sampling is one of the most widely used methods for animal abundance assessment. Standard estimation methods assume certain detection on the transect, no animal movement, and no measurement errors. Failure of the assumptions can cause substantial bias. In this work, the effect of error measurement on line transect estimators is investigated. Based on considerations of the process generating the errors, a multiplicative error model is presented and a simple way of correcting estimates based on knowledge of the error distribution is proposed. Using beta models for the error distribution, the effect of errors and of the proposed correction is assessed by simulation. Adequate confidence intervals for the corrected estimates are obtained using a bootstrap variance estimate for the correction and the delta method. As noted by Chen (1998, Biometrics 54, 899-908), even unbiased estimators of the distances might lead to biased density estimators, depending on the actual error distribution. In contrast with the findings of Chen, who used an additive model, unbiased estimation of distances, given a multiplicative model, lead to overestimation of density. Some error distributions result in observed distance distributions that make efficient estimation impossible, by removing the shoulder present in the original detection function. This indicates the need to improve field methods to reduce measurement error. An application of the new methods to a real data set is presented.     

Random effects in censored ordinal regression: latent structure and Bayesian approach

This paper discusses random effects in censored ordinal regression and presents a Gibbs sampling approach to fit the regression model. A latent structure and its corresponding Bayesian formulation are introduced to effectively deal with heterogeneous and censored ordinal observations. This work is motivated by the need to analyze interval-censored ordinal data from multiple studies in toxicological risk assessment. Application of our methodology to the data offers further support to the conclusions developed earlier using GEE methods yet provides additional insight into the uncertainty levels of the risk estimates.     

Bayesian estimation in animal breeding using the Dirichlet process prior for correlated random effects

In the case of the mixed linear model the random effects are usually assumed to be normally distributed in both the Bayesian and classical frameworks. In this paper, the Dirichlet process prior was used to provide nonparametric Bayesian estimates for correlated random effects. This goal was achieved by providing a Gibbs sampler algorithm that allows these correlated random effects to have a nonparametric prior distribution. A sampling based method is illustrated. This method which is employed by transforming the genetic covariance matrix to an identity matrix so that the random effects are uncorrelated, is an extension of the theory and the results of previous researchers. Also by using Gibbs sampling and data augmentation a simulation procedure was derived for estimating the precision parameter M associated with the Dirichlet process prior. All needed conditional posterior distributions are given. To illustrate the application, data from the Elsenburg Dormer sheep stud were analysed. A total of 3325 weaning weight records from the progeny of 101 sires were used.     

A hierarchical Bayesian model for combining multiple 2 x 2 tables using conditional likelihoods

This paper introduces a hierarchical Bayesian model for combining multiple 2 x 2 tables that allows the flexibility of different odds ratio estimates for different tables and at the same time allows the tables to borrow information from each other. The proposed model, however, is different from a full Bayesian model in that the nuisance parameters are eliminated by conditioning instead of integration. The motivation is a more robust model and a faster and more stable Gibbs algorithm. We work out a Gibbs scheme using the adaptive rejection sampling for log concave density and an algorithm for the mean and variance of the noncentral hypergeometric distribution. The model is applied to a multicenter ulcer clinical trial.     

Hierarchical Bayes estimation of hunting success rates with spatial correlations

A Bayesian hierarchical generalized linear model is used to estimate hunting success rates at the subarea level for postseason harvest surveys. The model includes fixed week effects, random geographic effects, and spatial correlations between neighboring subareas. The computation is done by Gibbs sampling and adaptive rejection sampling techniques. The method is illustrated using data from the Missouri Turkey Hunting Survey in the spring of 1996. Bayesian model selection methods are used to demonstrate that there are significant week differences and spatial correlations of hunting success rates among counties. The Bayesian estimates are also shown to be quite robust in terms of changes of hyperparameters.     

A comparison of alternative methods to compute conditional genotype probabilities for genetic evaluation with finite locus models

An increased availability of genotypes at marker loci has prompted the development of models that include the effect of individual genes. Selection based on these models is known as marker-assisted selection (MAS). MAS is known to be efficient especially for traits that have low heritability and non-additive gene action. BLUP methodology under non-additive gene action is not feasible for large inbred or crossbred pedigrees. It is easy to incorporate non-additive gene action in a finite locus model. Under such a model, the unobservable genotypic values can be predicted using the conditional mean of the genotypic values given the data. To compute this conditional mean, conditional genotype probabilities must be computed. In this study these probabilities were computed using iterative peeling, and three Markov chain Monte Carlo (MCMC) methods – scalar Gibbs, blocking Gibbs, and a sampler that combines the Elston Stewart algorithm with iterative peeling (ESIP). The performance of these four methods was assessed using simulated data. For pedigrees with loops, iterative peeling fails to provide accurate genotype probability estimates for some pedigree members. Also, computing time is exponentially related to the number of loci in the model. For MCMC methods, a linear relationship can be maintained by sampling genotypes one locus at a time. Out of the three MCMC methods considered, ESIP, performed the best while scalar Gibbs performed the worst.     

Indirect Methods Produce Higher Estimates of Fine Root Production and Turnover Rates than Direct Methods

The production and turnover of fine roots play substantial roles in the biogeochemical cycles of terrestrial ecosystems. However, the disparity among the estimates of both production and turnover, particularly due to technical limitations, has been debated for several decades. Here, we conducted a meta-analysis to compare published estimates of fine root production and turnover rates derived from different methods at the same sites and at the same sampling time. On average, the estimates of fine root production and turnover rates were 87% and 124% higher, respectively, by indirect methods than by direct methods. The substantially higher fine root production and turnover estimated by indirect methods, on which most global carbon models are based, indicate the necessity of re-assessing the global carbon model predictions for atmospheric carbon sequestration in soils as a result of the production and turnover of fine roots.     

HMM sampling and applications to gene finding and alternative splicing

The standard method of applying hidden Markov models to biological problems is to find a Viterbi (maximal weight) path through the HMM graph. The Viterbi algorithm reduces the problem of finding the most likely hidden state sequence that explains given observations, to a dynamic programming problem for corresponding directed acyclic graphs. For example, in the gene finding application, the HMM is used to find the most likely underlying gene structure given a DNA sequence. In this note we discuss the applications of sampling methods for HMMs. The standard sampling algorithm for HMMs is a variant of the common forward-backward and backtrack algorithms, and has already been applied in the context of Gibbs sampling methods. Nevetheless, the practice of sampling state paths from HMMs does not seem to have been widely adopted, and important applications have been overlooked. We show how sampling can be used for finding alternative splicings for genes, including alternative splicings that are conserved between genes from related organisms. We also show how sampling from the posterior distribution is a natural way to compute probabilities for predicted exons and gene structures being correct under the assumed model. Finally, we describe a new memory efficient sampling algorithm for certain classes of HMMs which provides a practical sampling alternative to the Hirschberg algorithm for optimal alignment. The ideas presented have applications not only to gene finding and HMMs but more generally to stochastic context free grammars and RNA structure prediction.     

Study on mapping Quantitative Trait Loci for animal complex binary traits using Bayesian-Markov chain Monte Carlo approach

It is a challenging issue to map Quantitative Trait Loci (QTL) underlying complex discrete traits, which usually show discontinuous distribution and less information, using conventional statistical methods. Bayesian-Markov chain Monte Carlo (Bayesian-MCMC) approach is the key procedure in mapping QTL for complex binary traits, which provides a complete posterior distribution for QTL parameters using all prior information. As a consequence, Bayesian estimates of all interested variables can be obtained straightforwardly basing on their posterior samples simulated by the MCMC algorithm. In our study, utilities of Bayesian-MCMC are demonstrated using simulated several animal outbred full-sib families with different family structures for a complex binary trait underlied by both a QTL and polygene. Under the Identity-by-Descent-Based variance component random model, three samplers basing on MCMC, including Gibbs sampling, Metropolis algorithm and reversible jump MCMC, were implemented to generate the joint posterior distribution of all unknowns so that the QTL parameters were obtained by Bayesian statistical inferring. The results showed that Bayesian-MCMC approach could work well and robust under different family structures and QTL effects. As family size increases and the number of family decreases, the accuracy of the parameter estimates will be improved. When the true QTL has a small effect, using outbred population experiment design with large family size is the optimal mapping strategy.     

Estimating population density from indirect sign: track counts and the Formozov–Malyshev–Pereleshin formula

For many purposes it is often desirable to estimate animal population densities over large areas. Where total counts are not possible and sightings are relatively rare, a range of methods exists to estimate densities from indirect sign. Such methods are frequently unreliable and usually require independent calibration or confirmation. We present an analytical method for estimating population density from track counts. The method, widely known in the Russian Federation but not in the English language scientific literature, requires counts of tracks of known age, together with estimates of animal daily travel distances. We use simulations to verify the theoretical basis of the approach and to indicate potential precision that may be achieved. We illustrate application of the approach using a large data set on ungulate track counts in the Russian Far East. We suggest that under most circumstances, nonparametric bootstrapping will be the most appropriate method for deriving estimates of confidence intervals about density estimates. As with other approaches to estimating density from indirect sign, the method that we discuss is vulnerable to violations of an array of underlying assumptions. However, it is easily applied and could represent an important method by which the relationship between indices of abundance and absolute density can be understood.     

A bivariate quantitative genetic model for a linear Gaussian trait and a survival trait

With the increasing use of survival models in animal breeding to address the genetic aspects of mainly longevity of livestock but also disease traits, the need for methods to infer genetic correlations and to do multivariate evaluations of survival traits and other types of traits has become increasingly important. In this study we derived and implemented a bivariate quantitative genetic model for a linear Gaussian and a survival trait that are genetically and environmentally correlated. For the survival trait, we considered the Weibull log-normal animal frailty model. A Bayesian approach using Gibbs sampling was adopted. Model parameters were inferred from their marginal posterior distributions. The required fully conditional posterior distributions were derived and issues on implementation are discussed. The two Weibull baseline parameters were updated jointly using a Metropolis-Hasting step. The remaining model parameters with non-normalized fully conditional distributions were updated univariately using adaptive rejection sampling. Simulation results showed that the estimated marginal posterior distributions covered well and placed high density to the true parameter values used in the simulation of data. In conclusion, the proposed method allows inferring additive genetic and environmental correlations, and doing multivariate genetic evaluation of a linear Gaussian trait and a survival trait.     

Animal Density and Track Counts: Understanding the Nature of Observations Based on Animal Movements

Counting animals to estimate their population sizes is often essential for their management and conservation. Since practitioners frequently rely on indirect observations of animals, it is important to better understand the relationship between such indirect indices and animal abundance. The Formozov-Malyshev-Pereleshin (FMP) formula provides a theoretical foundation for understanding the relationship between animal track counts and the true density of species. Although this analytical method potentially has universal applicability wherever animals are readily detectable by their tracks, it has long been unique to Russia and remains widely underappreciated. In this paper, we provide a test of the FMP formula by isolating the influence of animal travel path tortuosity (i.e., convolutedness) on track counts. We employed simulations using virtual and empirical data, in addition to a field test comparing FMP estimates with independent estimates from line transect distance sampling. We verify that track counts (total intersections between animals and transects) are determined entirely by density and daily movement distances. Hence, the FMP estimator is theoretically robust against potential biases from specific shapes or patterns of animal movement paths if transects are randomly situated with respect to those movements (i.e., the transects do not influence animals’ movements). However, detectability (the detection probability of individual animals) is not determined simply by daily travel distance but also by tortuosity, so ensuring that all intersections with transects are counted regardless of the number of individual animals that made them becomes critical for an accurate density estimate. Additionally, although tortuosity has no bearing on mean track encounter rates, it does affect encounter rate variance and therefore estimate precision. We discuss how these fundamental principles made explicit by the FMP formula have widespread implications for methods of assessing animal abundance that rely on indirect observations.     

Estimation of additive, dominance and epistatic variance components using finite locus models implemented with a single-site Gibbs and a descent graph sampler

In a previous contribution, we implemented a finite locus model (FLM) for estimating additive and dominance genetic variances via a Bayesian method and a single-site Gibbs sampler. We observed a dependency of dominance variance estimates on locus number in the analysis FLM. Here, we extended the FLM to include two-locus epistasis, and implemented the analysis with two genotype samplers (Gibbs and descent graph) and three different priors for genetic effects (uniform and variable across loci, uniform and constant across loci, and normal). Phenotypic data were simulated for two pedigrees with 6300 and 12,300 individuals in closed populations, using several different, non-additive genetic models. Replications of these data were analysed with FLMs differing in the number of loci. Simulation results indicate that the dependency of non-additive genetic variance estimates on locus number persisted in all implementation strategies we investigated. However, this dependency was considerably diminished with normal priors for genetic effects as compared with uniform priors (constant or variable across loci). Descent graph sampling of genotypes modestly improved variance components estimation compared with Gibbs sampling. Moreover, a larger pedigree produced considerably better variance components estimation, suggesting this dependency might originate from data insufficiency. As the FLM represents an appealing alternative to the infinitesimal model for genetic parameter estimation and for inclusion of polygenic background variation in QTL mapping analyses, further improvements are warranted and might be achieved via improvement of the sampler or treatment of the number of loci as an unknown.     

Modelling and optimal control of immune response of renal transplant recipients

We consider the increasingly important and highly complex immunological control problem: control of the dynamics of immunosuppression for organ transplant recipients. The goal in this problem is to maintain the delicate balance between over-suppression (where opportunistic latent viruses threaten the patient) and under-suppression (where rejection of the transplanted organ is probable). First, a mathematical model is formulated to describe the immune response to both viral infection and introduction of a donor kidney in a renal transplant recipient. Some numerical results are given to qualitatively validate and demonstrate that this initial model exhibits appropriate characteristics of primary infection and reactivation for immunosuppressed transplant recipients. In addition, we develop a computational framework for designing adaptive optimal treatment regimes with partial observations and low-frequency sampling, where the state estimates are obtained by solving a second deterministic optimal tracking problem. Numerical results are given to illustrate the feasibility of this method in obtaining optimal treatment regimes with a balance between under-suppression and over-suppression of the immune system.     

Analysis of rabbit doe longevity using a semiparametric log-Normal animal frailty model with time-dependent covariates

Data on doe longevity in a rabbit population were analysed using a semiparametric log-Normal animal frailty model. Longevity was defined as the time from the first positive pregnancy test to death or culling due to pathological problems. Does culled for other reasons had right censored records of longevity. The model included time dependent covariates associated with year by season, the interaction between physiological state and the number of young born alive, and between order of positive pregnancy test and physiological state. The model also included an additive genetic effect and a residual in log frailty. Properties of marginal posterior distributions of specific parameters were inferred from a full Bayesian analysis using Gibbs sampling. All of the fully conditional posterior distributions defining a Gibbs sampler were easy to sample from, either directly or using adaptive rejection sampling. The marginal posterior mean estimates of the additive genetic variance and of the residual variance in log frailty were 0.247 and 0.690.     

Modelling and optimal control of immune response of renal transplant recipients

We consider the increasingly important and highly complex immunological control problem: control of the dynamics of immunosuppression for organ transplant recipients. The goal in this problem is to maintain the delicate balance between over-suppression (where opportunistic latent viruses threaten the patient) and under-suppression (where rejection of the transplanted organ is probable). First, a mathematical model is formulated to describe the immune response to both viral infection and introduction of a donor kidney in a renal transplant recipient. Some numerical results are given to qualitatively validate and demonstrate that this initial model exhibits appropriate characteristics of primary infection and reactivation for immunosuppressed transplant recipients. In addition, we develop a computational framework for designing adaptive optimal treatment regimes with partial observations and low-frequency sampling, where the state estimates are obtained by solving a second deterministic optimal tracking problem. Numerical results are given to illustrate the feasibility of this method in obtaining optimal treatment regimes with a balance between under-suppression and over-suppression of the immune system.