More Accurate Dose-Response Estimation Using Monte-Carlo Uncertainty Analysis: The Data Cube Approach

The traditional q₁ ^* methodology for constructing upper confidence limits (UCLs) for the low-dose slopes of quantal dose-response functions has two limitations: (i) it is based on an asymptotic statistical result that has been shown via Monte Carlo simulation not to hold in practice for small, real bioassay experiments (Portier and Hoel, 1983); and (ii) it assumes that the multistage model (which represents cumulative hazard as a polynomial function of dose) is correct. This paper presents an uncertainty analysis approach for fitting dose-response functions to data that does not require specific parametric assumptions or depend on asymptotic results. It has the advantage that the resulting estimates of the dose-response function (and uncertainties about it) no longer depend on the validity of an assumed parametric family nor on the accuracy of the asymptotic approximation. The method derives posterior densities for the true response rates in the dose groups, rather than deriving posterior densities for model parameters, as in other Bayesian approaches (Sielken, 1991), or resampling the observed data points, as in the bootstrap and other resampling methods. It does so by conditioning constrained maximum-entropy priors on the observed data. Monte Carlo sampling of the posterior (constrained, conditioned) probability distributions generate values of response probabilities that might be observed if the experiment were repeated with very large sample sizes. A dose-response curve is fit to each such simulated dataset. If no parametric model has been specified, then a generalized representation (e.g., a power-series or orthonormal polynomial expansion) of the unknown dose-response function is fit to each simulated dataset using “model-free” methods. The simulation-based frequency distribution of all the dose-response curves fit to the simulated datasets yields a posterior distribution function for the low-dose slope of the dose-response curve. An upper confidence limit on the low-dose slope is obtained directly from this posterior distribution. This “Data Cube” procedure is illustrated with a real dataset for benzene, and is seen to produce more policy-relevant insights than does the traditional q₁ ^* methodology. For example, it shows how far apart are the 90%, 95%, and 99% limits and reveals how uncertainty about total and incremental risk vary with dose level (typically being dominated at low doses by uncertainty about the response of the control group, and being dominated at high doses by sampling variability). Strengths and limitations of the Data Cube approach are summarized, and potential decision-analytic applications to making better informed risk management decisions are briefly discussed.     

Estimation and testing with overdispersed proportions using the beta-logistic regression model of Heckman and Willis

Methods are presented for modeling dose-related effects in proportion data when extra-binomial variability is a concern. Motivation is taken from experiments in developmental toxicology, where similarity among conceptuses within a litter leads to intralitter correlations and to overdispersion in the observed proportions. Appeal is made to the well-known beta-binomial distribution to represent the overdispersion. From this, an exponential function of the linear predictor is used to model the dose-response relationship. The specification was introduced previously for econometric applications by Heckman and Willis; it induces a form of logistic regression for the mean response, together with a reciprocal biexponential model for the intralitter correlation. Large-sample, likelihood-based methods for estimating and testing the joint proportion-correlation response are studied. A developmental toxicity data set illustrates the methods.     

Bayesian predictive inference for time series count data

Correlated count data arise often in practice, especially in repeated measures situations or instances in which observations are collected over time. In this paper, we consider a parametric model for a time series of counts by constructing a likelihood-based version of a model similar to that of Zeger (1988, Biometrika 75, 621-629). The model has the advantage of incorporating both overdispersion and autocorrelation. We consider a Bayesian approach and propose a class of informative prior distributions for the model parameters that are useful for prediction. The prior specification is motivated from the notion of the existence of data from similar previous studies, called historical data, which is then quantified into a prior distribution for the current study. We derive the Bayesian predictive distribution and use a Bayesian criterion, called the predictive L measure, for assessing the predictions for a given time series model. The distribution of the predictive L measure is also derived, which will enable us to compare the predictive ability for each model under consideration. Our methodology is motivated by a real data set involving yearly pollen counts, which is examined in some detail.     

Bayesian regularization for flexible baseline hazard functions in Cox survival models

Fully Bayesian methods for Cox models specify a model for the baseline hazard function. Parametric approaches generally provide monotone estimations. Semi‐parametric choices allow for more flexible patterns but they can suffer from overfitting and instability. Regularization methods through prior distributions with correlated structures usually give reasonable answers to these types of situations. We discuss Bayesian regularization for Cox survival models defined via flexible baseline hazards specified by a mixture of piecewise constant functions and by a cubic B‐spline function. For those “semi‐parametric” proposals, different prior scenarios ranging from prior independence to particular correlated structures are discussed in a real study with microvirulence data and in an extensive simulation scenario that includes different data sample and time axis partition sizes in order to capture risk variations. The posterior distribution of the parameters was approximated using Markov chain Monte Carlo methods. Model selection was performed in accordance with the deviance information criteria and the log pseudo‐marginal likelihood. The results obtained reveal that, in general, Cox models present great robustness in covariate effects and survival estimates independent of the baseline hazard specification. In relation to the “semi‐parametric” baseline hazard specification, the B‐splines hazard function is less dependent on the regularization process than the piecewise specification because it demands a smaller time axis partition to estimate a similar behavior of the risk.     

On symmetric semiparametric two-sample problem

We consider a two-sample problem where data come from symmetric distributions. Usual two-sample data with only magnitudes recorded, arising from case-control studies or logistic discriminant analyses, may constitute a symmetric two-sample problem. We propose a semiparametric model such that, in addition to symmetry, the log ratio of two unknown density functions is modeled in a known parametric form. The new semiparametric model, tailor-made for symmetric two-sample data, can also be viewed as a biased sampling model subject to symmetric constraint. A maximum empirical likelihood estimation approach is adopted to estimate the unknown model parameters, and the corresponding profile empirical likelihood ratio test is utilized to perform hypothesis testing regarding the two population distributions. Symmetry, however, comes with irregularity. It is shown that, under the null hypothesis of equal symmetric distributions, the maximum empirical likelihood estimator has degenerate Fisher information, and the test statistic has a mixture of χ²-type asymptotic distribution. Extensive simulation studies have been conducted to demonstrate promising statistical powers under correct and misspecified models. We apply the proposed methods to two real examples.     

A Bayesian nonparametric testing procedure for paired samples

We propose a Bayesian hypothesis testing procedure for comparing the distributions of paired samples. The procedure is based on a flexible model for the joint distribution of both samples. The flexibility is given by a mixture of Dirichlet processes. Our proposal uses a spike-slab prior specification for the base measure of the Dirichlet process and a particular parametrization for the kernel of the mixture in order to facilitate comparisons and posterior inference. The joint model allows us to derive the marginal distributions and test whether they differ or not. The procedure exploits the correlation between samples, relaxes the parametric assumptions, and detects possible differences throughout the entire distributions. A Monte Carlo simulation study comparing the performance of this strategy to other traditional alternatives is provided. Finally, we apply the proposed approach to spirometry data collected in the United States to investigate changes in pulmonary function in children and adolescents in response to air polluting factors.     

A Variational Bayes Approach to the Analysis of Occupancy Models

Detection-nondetection data are often used to investigate species range dynamics using Bayesian occupancy models which rely on the use of Markov chain Monte Carlo (MCMC) methods to sample from the posterior distribution of the parameters of the model. In this article we develop two Variational Bayes (VB) approximations to the posterior distribution of the parameters of a single-season site occupancy model which uses logistic link functions to model the probability of species occurrence at sites and of species detection probabilities. This task is accomplished through the development of iterative algorithms that do not use MCMC methods. Simulations and small practical examples demonstrate the effectiveness of the proposed technique. We specifically show that (under certain circumstances) the variational distributions can provide accurate approximations to the true posterior distributions of the parameters of the model when the number of visits per site (K) are as low as three and that the accuracy of the approximations improves as K increases. We also show that the methodology can be used to obtain the posterior distribution of the predictive distribution of the proportion of sites occupied (PAO).     

Synthesis of evidence from epidemiological studies with interval-censored exposure due to grouping

We describe a method for assessing dose-response effects from a series of case-control and cohort studies in which the exposure information is interval censored. The interval censoring of the exposure variable is dealt with through the use of retrospective models in which the exposure is treated as a multinomial response and disease status as a binary covariate. Polychotomous logistic regression models are adopted in which the dose-response relationship between exposure and disease may be modeled in a discrete or continuous fashion. Partial conditioning is possible to eliminate some of the nuisance parameters. The methods are applied to the motivating study of the relationship between chorionic villus sampling and the occurrence of terminal transverse limb reduction.     

Semiparametric estimation in copula models for bivariate sequential survival times

Sequentially observed survival times are of interest in many studies but there are difficulties in analyzing such data using nonparametric or semiparametric methods. First, when the duration of followup is limited and the times for a given individual are not independent, induced dependent censoring arises for the second and subsequent survival times. Non-identifiability of the marginal survival distributions for second and later times is another issue, since they are observable only if preceding survival times for an individual are uncensored. In addition, in some studies a significant proportion of individuals may never have the first event. Fully parametric models can deal with these features, but robustness is a concern. We introduce a new approach to address these issues. We model the joint distribution of the successive survival times by using copula functions, and provide semiparametric estimation procedures in which copula parameters are estimated without parametric assumptions on the marginal distributions. This provides more robust estimates and checks on the fit of parametric models. The methodology is applied to a motivating example involving relapse and survival following colon cancer treatment.     

A Generalized Logistic Model for Quantal Response Bioassay

In bioassay, where different levels of the stimulus may represent different doses of a drug, the binary response is the death or survival of an individual receiving a specified dose. In such applications, it is common to model the probability of a positive response P at the stimulus level x by P = F(x′β), where F is a cumulative distribution function and β is a vector of unknown parameters which characterize the response function. The two most popular models used for modelling binary response bioassay involve the probit model [BLISS (1935), FINNEY (1978)], and the logistic model [BERKSON (1944), BROWN (1982)]. However, these models have some limitations. The use of the probit model involves the inverse of the standard normal distribution function, making it rather intractable. The logistic model has a simple form and a closed expression for the inverse distribution function, however, neither the logistic nor the probit can provide a good fit to response functions which are not symmetric or are symmetric but have a steeper or gentler incline in the central probability region. In this paper we introduce a more realistic model for the analysis of quantal response bioassay. The proposed model, which we refer to it as the generalized logistic model, is a family of response curves indexed by shape parameters m₁ and m₂. This family is rich enough to include the probit and logistic models as well as many others as special cases or limiting distributions. In particular, we consider the generalized logistic three parameter model where we assume that m₁ = m, m is a positive real number, and m₂ = 1. We apply this model to various sets of data, comparing the fit results to those obtained previously by other dose-response curves such as the logistic and probit, and showing that the fit can be improved by using the generalized logistic.     

Combining multiple comparisons and modeling techniques in dose-response studies

The analysis of data from dose-response studies has long been divided according to two major strategies: multiple comparison procedures and model-based approaches. Model-based approaches assume a functional relationship between the response and the dose, taken as a quantitative factor, according to a prespecified parametric model. The fitted model is then used to estimate an adequate dose to achieve a desired response but the validity of its conclusions will highly depend on the correct choice of the a priori unknown dose-response model. Multiple comparison procedures regard the dose as a qualitative factor and make very few, if any, assumptions about the underlying dose-response model. The primary goal is often to identify the minimum effective dose that is statistically significant and produces a relevant biological effect. One approach is to evaluate the significance of contrasts between different dose levels, while preserving the family-wise error rate. Such procedures are relatively robust but inference is confined to the selection of the target dose among the dose levels under investigation. We describe a unified strategy to the analysis of data from dose-response studies which combines multiple comparison and modeling techniques. We assume the existence of several candidate parametric models and use multiple comparison techniques to choose the one most likely to represent the true underlying dose-response curve, while preserving the family-wise error rate. The selected model is then used to provide inference on adequate doses.     

Testing threshold and hormesis in a random effects dose-response model applied to developmental toxicity data

Here we describe a random effects threshold dose-response model for clustered binary-response data from developmental toxicity studies. For our model we assume that a hormetic effect occurs in addition to a threshold effect. Therefore, the dose-response curve is based on two components: relationships below the threshold (hormetic u-shaped model) and those above the threshold (logistic model). In the absence of hormesis and threshold effects, the estimation procedure is straightforward. We introduce score tests that are derived from a random effects hormetic-threshold dose-response model. The model and tests are applied to clustered binary data from developmental toxicity studies of animals to test for hormesis and threshold effects. We also compare the score test and likelihood ratio test to test for hormesis and threshold effects in a simulated study.     

Statistical optimization of parametric accelerated failure time model for mapping survival trait loci

Most existing statistical methods for mapping quantitative trait loci (QTL) are not suitable for analyzing survival traits with a skewed distribution and censoring mechanism. As a result, researchers incorporate parametric and semi-parametric models of survival analysis into the framework of the interval mapping for QTL controlling survival traits. In survival analysis, accelerated failure time (AFT) model is considered as a de facto standard and fundamental model for data analysis. Based on AFT model, we propose a parametric approach for mapping survival traits using the EM algorithm to obtain the maximum likelihood estimates of the parameters. Also, with Bayesian information criterion (BIC) as a model selection criterion, an optimal mapping model is constructed by choosing specific error distributions with maximum likelihood and parsimonious parameters. Two real datasets were analyzed by our proposed method for illustration. The results show that among the five commonly used survival distributions, Weibull distribution is the optimal survival function for mapping of heading time in rice, while Log-logistic distribution is the optimal one for hyperoxic acute lung injury.     

Quantitative Exposure-Response Assessment Approaches to Evaluate Acute Inhalation Toxicity of Phosgene

Phosgene has been a long-term subject of toxicological research due to its widespread use, high toxicity, and status as a model of chemically induced lung injury. To take advantage of the abundant data set for the acute inhalation toxicity of phosgene, methods for exposure-response analysis that use more data than the traditional no-observed-adverse-effect level approach were used to perform an exposure-response assessment for phosgene. Categorical regression is particularly useful for acute exposures due to the ability to combine studies of various exposure durations, and thus provide estimates of effect severity for a range of both exposure concentrations and durations. Results from the categorical regression approach were compared to those from parametric curve fitting models (i.e., benchmark concentration models) that make use of information from an entire dose-response, but only for one exposure duration. While categorical regression analysis provided results that were comparable to benchmark concentration results, categorical regression provides an improvement over that technique by accounting for the effects of both exposure concentration and duration on response. The other major advantage afforded by categorical regression is the ability to combine studies, allowing the quantitative use of a larger data set, which increases confidence in the final result.     

A Gibbs sampling approach to linkage analysis.

We present a Monte Carlo approach to estimation of the recombination fraction theta and the profile likelihood for a dichotomous trait and a single marker gene with 2 alleles. The method is an application of a technique known as 'Gibbs sampling', in which random samples of each of the unknowns (here genotypes, theta and nuisance parameters, including the allele frequencies and the penetrances) are drawn from their posterior distributions, given the data and the current values of all the other unknowns. Upon convergence, the resulting samples derive from the marginal distribution of all the unknowns, given only the data, so that the uncertainty in the specification of the nuisance parameters is reflected in the variance of the posterior distribution of theta. Prior knowledge about the distribution of theta and the nuisance parameters can be incorporated using a Bayesian approach, but adoption of a flat prior for theta and point priors for the nuisance parameters would correspond to the standard likelihood approach. The method is easy to program, runs quickly on a microcomputer, and could be generalized to multiple alleles, multipoint linkage, continuous phenotypes and more complex models of disease etiology. The basic approach is illustrated by application to data on cholesterol levels and an a low-density lipoprotein receptor gene in a single large pedigree.     

Confronting different models of community structure to species-abundance data: a Bayesian model comparison

Species abundances are undoubtedly the most widely available macroecological data, but can we use them to distinguish among several models of community structure? Here we present a Bayesian analysis of species‐abundance data that yields a full joint probability distribution of each model's parameters plus a relatively parameter‐independent criterion, the posterior Bayes factor, to compare these models. We illustrate our approach by comparing three classical distributions: the zero‐sum multinomial (ZSM) distribution, based on Hubbell's neutral model, the multivariate Poisson lognormal distribution (MPLN), based on niche arguments, and the discrete broken stick (DBS) distribution, based on MacArthur's broken stick model. We give explicit formulas for the probability of observing a particular species‐abundance data set in each model, and argue that conditioning on both sample size and species count is needed to allow comparisons between the two distributions. We apply our approach to two neotropical communities (trees, fish). We find that DBS is largely inferior to ZSM and MPLN for both communities. The tree data do not allow discrimination between ZSM and MPLN, but for the fish data ZSM (neutral model) overwhelmingly outperforms MPLN (niche model), suggesting that dispersal plays a previously underestimated role in structuring tropical freshwater fish communities. We advocate this approach for identifying the relative importance of dispersal and niche‐partitioning in determining diversity of different ecological groups of species under different environmental conditions.     

A Bayesian semiparametric joint hierarchical model for longitudinal and survival data

This article proposes a new semiparametric Bayesian hierarchical model for the joint modeling of longitudinal and survival data. We relax the distributional assumptions for the longitudinal model using Dirichlet process priors on the parameters defining the longitudinal model. The resulting posterior distribution of the longitudinal parameters is free of parametric constraints, resulting in more robust estimates. This type of approach is becoming increasingly essential in many applications, such as HIV and cancer vaccine trials, where patients' responses are highly diverse and may not be easily modeled with known distributions. An example will be presented from a clinical trial of a cancer vaccine where the survival outcome is time to recurrence of a tumor. Immunologic measures believed to be predictive of tumor recurrence were taken repeatedly during follow-up. We will present an analysis of this data using our new semiparametric Bayesian hierarchical joint modeling methodology to determine the association of these longitudinal immunologic measures with time to tumor recurrence.     

Bayesian semiparametric modeling for matched case-control studies with multiple disease states

We present a Bayesian approach to analyze matched "case-control" data with multiple disease states. The probability of disease development is described by a multinomial logistic regression model. The exposure distribution depends on the disease state and could vary across strata. In such a model, the number of stratum effect parameters grows in direct proportion to the sample size leading to inconsistent MLEs for the parameters of interest even when one uses a retrospective conditional likelihood. We adopt a semiparametric Bayesian framework instead, assuming a Dirichlet process prior with a mixing normal distribution on the distribution of the stratum effects. We also account for possible missingness in the exposure variable in our model. The actual estimation is carried out through a Markov chain Monte Carlo numerical integration scheme. The proposed methodology is illustrated through simulation and an example of a matched study on low birth weight of newborns (Hosmer, D. A. and Lemeshow, S., 2000, Applied Logistic Regression) with two possible disease groups matched with a control group.     

A Caution for Monte Carlo Risk Assessment of Long Term Exposures based on Short-Term Exposure Study Data

Monte Carlo risk assessments commonly take as input empirical or parametric exposure distributions from specially designed exposure studies. The exposure studies typically have limited duration, since their design is based on statistical and practical factors (such as cost and respondent burden). For these reasons, the exposure period studied rarely corresponds to the biologic exposure period, which we define as the time at risk that is relevant for quantifying exposure that may result in health effects. Both the exposure period studied and the biologic exposure period will often differ from the exposure interval used in a Monte Carlo analysis. Such time period differences, which are often not accounted for, can have dramatic effects on the ultimate risk assessment. When exposure distributions are right skewed and/or follow a lognormal distribution, exposure will usually be overestimated for percentiles above the median by direct use of exposure study empirical data, since biologic exposure periods are generally longer than the exposure periods in exposure assessment studies. We illustrate the effect that biologic exposure time period and response error can have on exposure distributions, using soil ingestion as an example. Beginning with variance components from lognormally distributed soil ingestion estimates, we illustrate the effect of different modeling assumptions, and the sensitivity of the resulting analyses to these assumptions. We develop a strategy for determining appropriate exposure input distributions for soil ingestion, and illustrate this using data on soil ingestion in children.     

Bayesian analysis of response to selection: a case study using litter size in Danish Yorkshire pigs

Implementation of a Bayesian analysis of a selection experiment is illustrated using litter size [total number of piglets born (TNB)] in Danish Yorkshire pigs. Other traits studied include average litter weight at birth (WTAB) and proportion of piglets born dead (PRBD). Response to selection for TNB was analyzed with a number of models, which differed in their level of hierarchy, in their prior distributions, and in the parametric form of the likelihoods. A model assessment study favored a particular form of an additive genetic model. With this model, the Monte Carlo estimate of the 95% probability interval of response to selection was (0.23; 0.60), with a posterior mean of 0.43 piglets. WTAB showed a correlated response of -7.2 g, with a 95% probability interval equal to (-33.1; 18.9). The posterior mean of the genetic correlation between TNB and WTAB was -0.23 with a 95% probability interval equal to (-0.46; -0.01). PRBD was studied informally; it increases with larger litters, when litter size is >7 piglets born. A number of methodological issues related to the Bayesian model assessment study are discussed, as well as the genetic consequences of inferring response to selection using additive genetic models.