Correcting for nonrandom ascertainment in generalized linear mixed models (GLMMs), fitted using Gibbs sampling

Gibbs sampling-based generalized linear mixed models (GLMMs) provide a convenient and flexible way to extend variance components models for multivariate normally distributed continuous traits to other classes of phenotype. This includes binary traits and right-censored failure times such as age-at-onset data. The approach has applications in many areas of genetic epidemiology. However, the required GLMMs are sensitive to nonrandom ascertainment. In the absence of an appropriate correction for ascertainment, they can exhibit marked positive bias in the estimated grand mean and serious shrinkage in the estimated magnitude of variance components. To compound practical difficulties, it is currently difficult to implement a conventional adjustment for ascertainment because of the need to undertake repeated integration across the distribution of random effects. This is prohibitively slow when it must be repeated at every iteration of the Markov chain Monte Carlo (MCMC) procedure. This paper motivates a correction for ascertainment that is based on sampling random effects rather than integrating across them and can therefore be implemented in a general-purpose Gibbs sampling environment such as WinBUGS. The approach has the characteristic that it returns ascertainment-adjusted parameter estimates that pertain to the true distribution of determinants in the ascertained sample rather than in the general population. The implications of this characteristic are investigated and discussed. This paper extends the utility of Gibbs sampling-based GLMMs to a variety of settings in which family data are ascertained nonrandomly.     

A Bayesian semiparametric Markov regression model for juvenile dermatomyositis

Juvenile dermatomyositis (JDM) is a rare autoimmune disease that may lead to serious complications, even to death. We develop a 2‐state Markov regression model in a Bayesian framework to characterise disease progression in JDM over time and gain a better understanding of the factors influencing disease risk. The transition probabilities between disease and remission state (and vice versa) are a function of time‐homogeneous and time‐varying covariates. These latter types of covariates are introduced in the model through a latent health state function, which describes patient‐specific health over time and accounts for variability among patients. We assume a nonparametric prior based on the Dirichlet process to model the health state function and the baseline transition intensities between disease and remission state and vice versa. The Dirichlet process induces a clustering of the patients in homogeneous risk groups. To highlight clinical variables that most affect the transition probabilities, we perform variable selection using spike and slab prior distributions. Posterior inference is performed through Markov chain Monte Carlo methods. Data were made available from the UK JDM Cohort and Biomarker Study and Repository, hosted at the UCL Institute of Child Health.     

Bayesian adaptive group lasso with semiparametric hidden Markov models

This paper presents a Bayesian adaptive group least absolute shrinkage and selection operator method to conduct simultaneous model selection and estimation under semiparametric hidden Markov models. We specify the conditional regression model and the transition probability model in the hidden Markov model into additive nonparametric functions of covariates. A basis expansion is adopted to approximate the nonparametric functions. We introduce multivariate conditional Laplace priors to impose adaptive penalties on regression coefficients and different groups of basis expansions under the Bayesian framework. An efficient Markov chain Monte Carlo algorithm is then proposed to identify the nonexistent, constant, linear, and nonlinear forms of covariate effects in both conditional and transition models. The empirical performance of the proposed methodology is evaluated via simulation studies. We apply the proposed model to analyze a real data set that was collected from the Alzheimer's Disease Neuroimaging Initiative study. The analysis identifies important risk factors on cognitive decline and the transition from cognitive normal to Alzheimer's disease.     

Bayesian inference for generalized linear mixed models with predictors subject to detection limits: an approach that leverages information from auxiliary variables

This paper is motivated from a retrospective study of the impact of vitamin D deficiency on the clinical outcomes for critically ill patients in multi‐center critical care units. The primary predictors of interest, vitamin D2 and D3 levels, are censored at a known detection limit. Within the context of generalized linear mixed models, we investigate statistical methods to handle multiple censored predictors in the presence of auxiliary variables. A Bayesian joint modeling approach is proposed to fit the complex heterogeneous multi‐center data, in which the data information is fully used to estimate parameters of interest. Efficient Monte Carlo Markov chain algorithms are specifically developed depending on the nature of the response. Simulation studies demonstrate the outperformance of the proposed Bayesian approach over other existing methods. An application to the data set from the vitamin D deficiency study is presented. Possible extensions of the method regarding the absence of auxiliary variables, semiparametric models, as well as the type of censoring are also discussed. Copyright © 2015 John Wiley & Sons, Ltd.     

A chain of evidence with mixed comparisons: models for multi-parameter synthesis and consistency of evidence

Multi-parameter evidence synthesis is a generalization of meta-analysis in which several parameters are estimated jointly. Here this approach is applied to a common form of data structure in which mixed pairwise comparisons are made between treatments in a 'chain of evidence' structure. The data set investigated here, which first appeared in the confidence profile method literature, features three types of data relating to thrombolytic treatment following acute myocardial infarction. There is information on reperfusion (coronary patency) following treatment, information on survival in reperfused and non-reperfused patients (conditional survival), and independent information on the effect of treatment on survival without specifying reperfusion state (overall survival). The objective of this study is to explore models for combining these three types of evidence within a single model, and some of the evidence consistency issues that arise. Bayesian Markov chain Monte Carlo methods are used to fit two models. The first, proposed in the confidence profile literature, assumes there are no differences between studies in baseline effects (equal study effects). The second model assumes fixed treatment effects but allows for study differences in baseline reperfusion and conditional survival rates by assuming random study effects. The equal study effects model fits the data poorly, and cross-validation shows that the overall survival evidence is not consistent with the reperfusion and conditional survival evidence under this model. The evidence appears to be consistent within the more loosely structured random study effects model, but results are sensitive to prior assumptions about the between-study variance. The inconsistency between the three evidence types relates to overall survival, rather than to the effect of treatment on survival. Multi-parameter evidence synthesis, with suitable checks for evidence consistency, can be used to help determine whether chains of evidence 'add up' epidemiologically within an overall model, and to construct parameter distributions for an entire model simultaneously, consistent with all the available data.     

Bayesian inference for multivariate meta‐analysis Box–Cox transformation models for individual patient data with applications to evaluation of cholesterol‐lowering drugs

In this paper, we propose a class of Box–Cox transformation regression models with multidimensional random effects for analyzing multivariate responses for individual patient data in meta‐analysis. Our modeling formulation uses a multivariate normal response meta‐analysis model with multivariate random effects, in which each response is allowed to have its own Box–Cox transformation. Prior distributions are specified for the Box–Cox transformation parameters as well as the regression coefficients in this complex model, and the deviance information criterion is used to select the best transformation model. Because the model is quite complex, we develop a novel Monte Carlo Markov chain sampling scheme to sample from the joint posterior of the parameters. This model is motivated by a very rich dataset comprising 26 clinical trials involving cholesterol‐lowering drugs where the goal is to jointly model the three‐dimensional response consisting of low density lipoprotein cholesterol (LDL‐C), high density lipoprotein cholesterol (HDL‐C), and triglycerides (TG) (LDL‐C, HDL‐C, TG). Because the joint distribution of (LDL‐C, HDL‐C, TG) is not multivariate normal and in fact quite skewed, a Box–Cox transformation is needed to achieve normality. In the clinical literature, these three variables are usually analyzed univariately; however, a multivariate approach would be more appropriate because these variables are correlated with each other. We carry out a detailed analysis of these data by using the proposed methodology. Copyright © 2013 John Wiley & Sons, Ltd.     

Level‐adjusted funnel plots based on predicted marginal expectations: an application to prophylactic antibiotics in gallstone surgery

Funnel plots are widely used to visualize grouped data, for example, in institutional comparison. This paper extends the concept to a multi‐level setting, displaying one level at a time, adjusted for the other levels, as well as for covariates at all levels. These level‐adjusted funnel plots are based on a Markov chain Monte Carlo fit of a random effects model, translating the estimated model parameters to predicted marginal expectations. Working within the estimation framework, we accommodate outlying institutions using heavy‐tailed random effects distributions. We also develop computer‐efficient methods to compute predicted probabilities in the case of dichotomous outcome data and various random effect distributions. We apply the method to a data set on prophylactic antibiotics in gallstone surgery. Copyright © 2014 John Wiley & Sons, Ltd.     

The effect of provider-level ascertainment bias on profiling nursing homes

Profiling health care providers for the purpose of public reporting and quality improvement has become commonplace. Recently, the Centers for Medicare and Medicaid Services (CMS) began publishing measures of quality for every Medicare/Medicaid-certified nursing home in the country. The facility-specific quality indicators (QIs) reported by CMS are based on quarterly measures from the minimum data set (MDS). However, some QIs from the MDS are potentially subject to ascertainment bias. Ascertainment bias would occur if there was variation in the way items that make up QIs are measured by nurses from each facility. This is potentially a problem for difficult-to-measure items such as pain and pressure ulcers. To assess the impact of ascertainment bias on profiling, we utilize data from a reliability study of nursing homes from six states. We develop methods for profiling providers in situations where the data consist of a response variable for each subject based on assessments from an internal rater, and, for a subset of subjects in each facility, a response variable based on assessments from an independent (external) rater. The internal assessments are potentially subject to provider-level ascertainment bias, whereas the independent assessments are considered the 'gold standard'. Our methods extend popular Bayesian approaches for profiling by using the paired observations from the subset of subjects with error-prone and error-free assessments to adjust for ascertainment bias. We apply the methods to MDS merged with the reliability data, and compare the bias-corrected profiles with those of standard approaches.     

Hierarchical Bayes Analysis of Multilevel Health Services Data: A Canadian Neonatal Mortality Study

This article presents a comprehensive analysis of multilevel health services data in order to provide detailed information about important sources of outcome variation. We present a study of neonatal mortality variation among 17 neonatal intensive care units across Canada and briefly introduce the Neonatal Health Services in Canada Project, under which this work was carried out. We illustrate the potential of hierarchical Bayes modeling as a general strategy for systematically evaluating response-covariate associations at each level of the hierarchy, examining cross-level interaction, and attributing unexplained variability. We engage in an extensive discussion on relevant issues that are important for institutional comparison, provider profiling, and quality improvement.     

Genetic variance components analysis for binary phenotypes using generalized linear mixed models (GLMMs) and Gibbs sampling

The common complex diseases such as asthma are an important focus of genetic research, and studies based on large numbers of simple pedigrees ascertained from population‐based sampling frames are becoming commonplace. Many of the genetic and environmental factors causing these diseases are unknown and there is often a strong residual covariance between relatives even after all known determinants are taken into account. This must be modelled correctly whether scientific interest is focused on fixed effects, as in an association analysis, or on the covariances themselves. Analysis is straightforward for multivariate Normal phenotypes, but difficulties arise with other types of trait. Generalized linear mixed models (GLMMs) offer a potentially unifying approach to analysis for many classes of phenotype including multivariate Normal traits, binary traits, and censored survival times. Markov Chain Monte Carlo methods, including Gibbs sampling, provide a convenient framework within which such models may be fitted. In this paper, Bayesian inference Using Gibbs Sampling (a generic Gibbs sampler; BUGS) is used to fit GLMMs for multivariate Normal and binary phenotypes in nuclear families. BUGS is easy to use and readily available. We motivate a suitable model structure for Normal phenotypes and show how the model extends to binary traits. We discuss parameter interpretation and statistical inference and show how to circumvent a number of important theoretical and practical problems that we encountered. Using simulated data we show that model parameters seem consistent and appear unbiased in smaller data sets. We illustrate our methods using data from an ongoing cohort study. Genet. Epidemiol. 17:118–140, 1999. © 1999 Wiley‐Liss, Inc.     

Reducing Monte Carlo error in the Bayesian estimation of risk ratios using log‐binomial regression models

In cohort studies, binary outcomes are very often analyzed by logistic regression. However, it is well known that when the goal is to estimate a risk ratio, the logistic regression is inappropriate if the outcome is common. In these cases, a log‐binomial regression model is preferable. On the other hand, the estimation of the regression coefficients of the log‐binomial model is difficult owing to the constraints that must be imposed on these coefficients. Bayesian methods allow a straightforward approach for log‐binomial regression models and produce smaller mean squared errors in the estimation of risk ratios than the frequentist methods, and the posterior inferences can be obtained using the software WinBUGS. However, Markov chain Monte Carlo methods implemented in WinBUGS can lead to large Monte Carlo errors in the approximations to the posterior inferences because they produce correlated simulations, and the accuracy of the approximations are inversely related to this correlation. To reduce correlation and to improve accuracy, we propose a reparameterization based on a Poisson model and a sampling algorithm coded in R. Copyright © 2015 John Wiley & Sons, Ltd.     

Bayesian modeling and inference for clinical trials with partial retrieved data following dropout

In randomized clinical trials, it is common that patients may stop taking their assigned treatments and then switch to a standard treatment (standard of care available to the patient) but not the treatments under investigation. Although the availability of limited retrieved data on patients who switch to standard treatment, called off‐protocol data, could be highly valuable in assessing the associated treatment effect with the experimental therapy, it leads to a complex data structure requiring the development of models that link the information of per‐protocol data with the off‐protocol data. In this paper, we develop a novel Bayesian method to jointly model longitudinal treatment measurements under various dropout scenarios. Specifically, we propose a multivariate normal mixed‐effects model for repeated measurements from the assigned treatments and the standard treatment, a multivariate logistic regression model for those stopping the assigned treatments, logistic regression models for those starting a standard treatment off protocol, and a conditional multivariate logistic regression model for completely withdrawing from the study. We assume that withdrawing from the study is non‐ignorable, but intermittent missingness is assumed to be at random. We examine various properties of the proposed model. We develop an efficient Markov chain Monte Carlo sampling algorithm. We analyze in detail via the proposed method a real dataset from a clinical trial. Copyright © 2013 John Wiley & Sons, Ltd.     

Hidden Markov models for zero-inflated Poisson counts with an application to substance use

Paradigms for substance abuse cue-reactivity research involve pharmacological or stressful stimulation designed to elicit stress and craving responses in cocaine-dependent subjects. It is unclear as to whether stress induced from participation in such studies increases drug-seeking behavior. We propose a 2-state Hidden Markov model to model the number of cocaine abuses per week before and after participation in a stress-and cue-reactivity study. The hypothesized latent state corresponds to 'high' or 'low' use. To account for a preponderance of zeros, we assume a zero-inflated Poisson model for the count data. Transition probabilities depend on the prior week's state, fixed demographic variables, and time-varying covariates. We adopt a Bayesian approach to model fitting, and use the conditional predictive ordinate statistic to demonstrate that the zero-inflated Poisson hidden Markov model outperforms other models for longitudinal count data.     

Tutorial in biostatistics: spline smoothing with linear mixed models

The semi-parametric regression achieved via penalized spline smoothing can be expressed in a linear mixed models framework. This allows such models to be fitted using standard mixed models software routines with which many biostatisticians are familiar. Moreover, the analysis of complex correlated data structures that are a hallmark of biostatistics, and which are typically analysed using mixed models, can now incorporate directly smoothing of the relationship between an outcome and covariates. In this paper we provide an introduction to both linear mixed models and penalized spline smoothing, and describe the connection between the two. This is illustrated with three examples, the first using birth data from the U.K., the second relating mammographic density to age in a study of female twin-pairs and the third modelling the relationship between age and bronchial hyperresponsiveness in families. The models are fitted in R (a clone of S-plus) and using Markov chain Monte Carlo (MCMC) implemented in the package WinBUGS.     

A monotone data augmentation algorithm for multivariate nonnormal data: With applications to controlled imputations for longitudinal trials

An efficient monotone data augmentation (MDA) algorithm is proposed for missing data imputation for incomplete multivariate nonnormal data that may contain variables of different types and are modeled by a sequence of regression models including the linear, binary logistic, multinomial logistic, proportional odds, Poisson, negative binomial, skew-normal, skew-t regressions, or a mixture of these models. The MDA algorithm is applied to the sensitivity analyses of longitudinal trials with nonignorable dropout using the controlled pattern imputations that assume the treatment effect reduces or disappears after subjects in the experimental arm discontinue the treatment. We also describe a heuristic approach to implement the controlled imputation, in which the fully conditional specification method is used to impute the intermediate missing data to create a monotone missing pattern, and the missing data after dropout are then imputed according to the assumed nonignorable mechanisms. The proposed methods are illustrated by simulation and real data analyses. Sample SAS code for the analyses is provided in the supporting information     

Bayesian approach for flexible modeling of semicompeting risks data

Semicompeting risks data arise when two types of events, non‐terminal and terminal, are observed. When the terminal event occurs first, it censors the non‐terminal event, but not vice versa. To account for possible dependent censoring of the non‐terminal event by the terminal event and to improve prediction of the terminal event using the non‐terminal event information, it is crucial to model their association properly. Motivated by a breast cancer clinical trial data analysis, we extend the well‐known illness–death models to allow flexible random effects to capture heterogeneous association structures in the data. Our extension also represents a generalization of the popular shared frailty models that usually assume that the non‐terminal event does not affect the hazards of the terminal event beyond a frailty term. We propose a unified Bayesian modeling approach that can utilize existing software packages for both model fitting and individual‐specific event prediction. The approach is demonstrated via both simulation studies and a breast cancer data set analysis. Copyright © 2014 John Wiley & Sons, Ltd.     

Constrained S-estimators for linear mixed effects models with covariance components

Linear mixed effects (LME) models are increasingly used for analyses of biological and biomedical data. When the multivariate normal assumption is not adequate for an LME model, then a robust estimation approach is preferable to the maximum likelihood one. M-estimators were considered before for robust estimation of the LME models, and recently a constrained S-estimator was proposed. This S-estimator cannot be applied directly to LME models with correlated error terms and vector random effects with correlated dimensions. Therefore, a modification is proposed, which extends application of the constrained S-estimator to the LME models for multivariate responses with correlated dimensions and to longitudinal data. Also, a new computational algorithm is developed for computing constrained S-estimators. Performance of the S-estimators based on the original Tukey's biweight and translated biweight is evaluated in a small simulation study with repeated multivariate responses with correlated dimensions. The proposed methodology is applied to jointly analyze repeated measures on three cholesterol components, HDL, LDL, and triglycerides.     

A discrete random effects probit model with application to the demand for preventive care

I have developed a random effects probit model in which the distribution of the random intercept is approximated by a discrete density. Monte Carlo results show that only three to four points of support are required for the discrete density to closely mimic normal and chi-squared densities and provide unbiased estimates of the structural parameters and the variance of the random intercept. The empirical application shows that both observed family characteristics and unobserved family-level heterogeneity are important determinants of the demand for preventive care.     

Predicting time to prostate cancer recurrence based on joint models for non-linear longitudinal biomarkers and event time outcomes

Biological markers that are both sensitive and specific for tumour regrowth or metastasis are increasingly becoming available and routinely monitored during the regular follow-up of patients treated for cancer. Obtained by a simple blood test, these markers provide an inexpensive non-invasive means for the early detection of recurrence (or progression). Currently, the longitudinal behaviour of the marker is viewed as an indicator of early disease progression, and is applied by a physician in making clinical decisions. One marker that has been studied for use in both population screening for early disease and for detection of recurrence in prostate cancer patients is PSA. The elevation of PSA levels is known to precede clinically detectable recurrence by 2 to 5 years, and current clinical practice often relies partially on multiple recent rises in PSA to trigger a change in treatment. However, the longitudinal trajectory for individual markers is often non-linear; in many cases there is a decline immediately following radiation therapy or surgery, a plateau during remission, followed by an exponential rise following the recurrence of the cancer. The aim of this article is to determine the multiple aspects of the longitudinal PSA biomarker trajectory that can be most sensitive for predicting time to clinical recurrence. Joint Bayesian models for the longitudinal measures and event times are utilized based on non-linear hierarchical models, implied by unknown change-points, for the longitudinal trajectories, and a Cox proportional hazard model for progression times, with functionals of the longitudinal parameters as covariates in the Cox model. Using Markov chain Monte Carlo sampling schemes, the joint model is fit to longitudinal PSA measures from 676 patients treated at Massachusetts General Hospital between the years 1988 and 1995 with follow-up to 1999. Based on these data, predictive schemes for detecting cancer recurrence in new patients based on their longitudinal trajectory are derived.     

A Bayesian toolkit for genetic association studies

We present a range of modelling components designed to facilitate Bayesian analysis of genetic-association-study data. A key feature of our approach is the ability to combine different submodels together, almost arbitrarily, for dealing with the complexities of real data. In particular, we propose various techniques for selecting the "best" subset of genetic predictors for a specific phenotype (or set of phenotypes). At the same time, we may control for complex, non-linear relationships between phenotypes and additional (non-genetic) covariates as well as accounting for any residual correlation that exists among multiple phenotypes. Both of these additional modelling components are shown to potentially aid in detecting the underlying genetic signal. We may also account for uncertainty regarding missing genotype data. Indeed, at the heart of our approach is a novel method for reconstructing unobserved haplotypes and/or inferring the values of missing genotypes. This can be deployed independently or, alternatively, it can be fully integrated into arbitrary genotype- or haplotype-based association models such that the missing data and the association model are "estimated" simultaneously. The impact of such simultaneous analysis on inferences drawn from the association model is shown to be potentially significant. Our modelling components are packaged as an "add-on" interface to the widely used WinBUGS software, which allows Markov chain Monte Carlo analysis of a wide range of statistical models. We illustrate their use with a series of increasingly complex analyses conducted on simulated data based on a real pharmacogenetic example.