Mixed effects models with bivariate and univariate association parameters for longitudinal bivariate binary response data

When two binary responses are measured for each study subject across time, it may be of interest to model how the bivariate associations and marginal univariate risks involving the two responses change across time. To achieve such a goal, marginal models with bivariate log odds ratio and univariate logit components are extended to include random effects for all components. Specifically, separate normal random effects are specified on the log odds ratio scale for bivariate responses and on the logit scale for univariate responses. Assuming conditional independence given the random effects facilitates the modeling of bivariate associations across time with missing at random incomplete data. We fit the model to a dataset for which such structures are feasible: a longitudinal randomized trial of a cardiovascular educational program where the responses of interest are change in hypertension and hypercholestemia status. The proposed model is compared to a naive bivariate model that assumes independence between time points and univariate mixed effects logit models.     

Association Models for Clustered Data with Binary and Continuous Responses

Summary .  We consider analysis of clustered data with mixed bivariate responses, i.e., where each member of the cluster has a binary and a continuous outcome. We propose a new bivariate random effects model that induces associations among the binary outcomes within a cluster, among the continuous outcomes within a cluster, between a binary outcome and a continuous outcome from different subjects within a cluster, as well as the direct association between the binary and continuous outcomes within the same subject. For the ease of interpretations of the regression effects, the marginal model of the binary response probability integrated over the random effects preserves the logistic form and the marginal expectation of the continuous response preserves the linear form. We implement maximum likelihood estimation of our model parameters using standard software such as PROC NLMIXED of SAS . Our simulation study demonstrates the robustness of our method with respect to the misspecification of the regression model as well as the random effects model. We illustrate our methodology by analyzing a developmental toxicity study of ethylene glycol in mice.     

The Analysis of Correlated Binary Outcomes Using Multivariate Logistic Regression

A method for analyzing correlated binary outcomes when the responses are distinct measurements made simultaneously on a single individual is presented. This extension of univariate logistic regression allows us to model the dependence of the responses on a set of covariates while estimating the degree of association among them. For the case of two dichotomous outcomes, a form of the cumulative bivariate logistic distribution proposed by Gumbel is used to characterize their joint probabilities in terms of logistic marginal probabilities and the correlation coefficient of the responses. The model is then extended to accommodate three or more dichotomous outcomes. A two-step approximation to fitting the multivariate logistic model is also described.     

Modeling continuous covariates with a “spike” at zero: Bivariate approaches

In epidemiology and clinical research, predictors often take value zero for a large amount of observations while the distribution of the remaining observations is continuous. These predictors are called variables with a spike at zero. Examples include smoking or alcohol consumption. Recently, an extension of the fractional polynomial (FP) procedure, a technique for modeling nonlinear relationships, was proposed to deal with such situations. To indicate whether or not a value is zero, a binary variable is added to the model. In a two stage procedure, called FP‐spike, the necessity of the binary variable and/or the continuous FP function for the positive part are assessed for a suitable fit. In univariate analyses, the FP‐spike procedure usually leads to functional relationships that are easy to interpret. This paper introduces four approaches for dealing with two variables with a spike at zero (SAZ). The methods depend on the bivariate distribution of zero and nonzero values. Bi‐Sep is the simplest of the four bivariate approaches. It uses the univariate FP‐spike procedure separately for the two SAZ variables. In Bi‐D3, Bi‐D1, and Bi‐Sub, proportions of zeros in both variables are considered simultaneously in the binary indicators. Therefore, these strategies can account for correlated variables. The methods can be used for arbitrary distributions of the covariates. For illustration and comparison of results, data from a case‐control study on laryngeal cancer, with smoking and alcohol intake as two SAZ variables, is considered. In addition, a possible extension to three or more SAZ variables is outlined. A combination of log‐linear models for the analysis of the correlation in combination with the bivariate approaches is proposed.     

Testing for ordered group effects in binary and continuous outcomes

A likelihood ratio test is proposed for the detection of an ordered group effect on bivariate responses where one response is binary and the other is continuous. The procedure is based on a conditional logistic model for the binary response given the continuous outcome. We also develop a likelihood ratio test for simultaneously determining the goodness of fit of the ordering assumption on both responses. Our approach is motivated by a particular toxicity study application involving laboratory animals that focused on the effect of a food color additive on the development of reticuloendothelial (RE) tumors. A brief discussion on extensions to the methodology introduced here is also given, along with a comparison of the approach with a marginal strategy where the presence of an ordered group effect is assessed independently for each of the two responses.     

Designs for Discrimination Between Bivariate Binary Response Models

We propose a test statistic for discrimination between alternative bivariate binary response models and the optimal design procedure which is an extension of T-optimality. Under certain conditions we prove that the maximum value of the power can be-obtained when the degrees of freedom of the test statistic is one. The conclusion is the same as that in discrimination between alternative univariate separate models. However the test statistics are different.     

Bayesian dose-finding in phase I/II clinical trials using toxicity and efficacy odds ratios

A Bayesian adaptive design is proposed for dose-finding in phase I/II clinical trials to incorporate the bivariate outcomes, toxicity and efficacy, of a new treatment. Without specifying any parametric functional form for the drug dose-response curve, we jointly model the bivariate binary data to account for the correlation between toxicity and efficacy. After observing all the responses of each cohort of patients, the dosage for the next cohort is escalated, deescalated, or unchanged according to the proposed odds ratio criteria constructed from the posterior toxicity and efficacy probabilities. A novel class of prior distributions is proposed through logit transformations which implicitly imposes a monotonic constraint on dose toxicity probabilities and correlates the probabilities of the bivariate outcomes. We conduct simulation studies to evaluate the operating characteristics of the proposed method. Under various scenarios, the new Bayesian design based on the toxicity-efficacy odds ratio trade-offs exhibits good properties and treats most patients at the desirable dose levels. The method is illustrated with a real trial design for a breast medical oncology study.     

Modeling animals' behavioral response by Markov chain models for capture-recapture experiments

A bivariate Markov chain approach that includes both enduring (long-term) and ephemeral (short-term) behavioral effects in models for capture-recapture experiments is proposed. The capture history of each animal is modeled as a Markov chain with a bivariate state space with states determined by the capture status (capture/noncapture) and marking status (marked/unmarked). In this framework, a conditional-likelihood method is used to estimate the population size and the transition probabilities. The classical behavioral model that assumes only an enduring behavioral effect is included as a special case of the bivariate Markovian model. Another special case that assumes only an ephemeral behavioral effect reduces to a univariate Markov chain based on capture/noncapture status. The model with the ephemeral behavioral effect is extended to incorporate time effects; in this model, in contrast to extensions of the classical behavioral model, all parameters are identifiable. A data set is analyzed to illustrate the use of the Markovian models in interpreting animals' behavioral response. Simulation results are reported to examine the performance of the estimators.     

Multivariate character process models for the analysis of two or more correlated function-valued traits

Various methods, including random regression, structured antedependence models, and character process models, have been proposed for the genetic analysis of longitudinal data and other function-valued traits. For univariate problems, the character process models have been shown to perform well in comparison to alternative methods. The aim of this article is to present an extension of these models to the simultaneous analysis of two or more correlated function-valued traits. Analytical forms for stationary and nonstationary cross-covariance functions are studied. Comparisons with the other approaches are presented in a simulation study and in an example of a bivariate analysis of genetic covariance in age-specific fecundity and mortality in Drosophila. As in the univariate case, bivariate character process models with an exponential correlation were found to be quite close to first-order structured antedependence models. The simulation study showed that the choice of the most appropriate methodology is highly dependent on the covariance structure of the data. The bivariate character process approach proved to be able to deal with quite complex nonstationary and nonsymmetric cross-correlation structures and was found to be the most appropriate for the real data example of the fruit fly Drosophila melanogaster.     

An assessment of non-randomized medical treatment of long-term schizophrenia relapse using bivariate binary-response transition models

The analyses of observational longitudinal studies involving concurrent changes in treatment and medical conditions present difficulties because of the multitude of directions of potential relationships: past medication influences current symptoms; past symptoms influence current medication; and current medication is associated with current symptoms. In the context of a long-term study of non-randomized pharmacological treatment of schizophrenic relapse, we present an analysis of bivariate discrete-time transitional data with binary responses in an attempt to understand the transitional and concurrent relationships between schizophrenia relapse and medication use. A naive analysis does not show any association between previous medication and current relapse. However, we provide evidence suggesting that current treatment may impact current relapse for those who have previously taken medication, but not for those who haven't taken medication in the past. When univariate models are specified to assess these associations, the bivariate nature of the problem requires a choice of which response, relapse or medication, should be the dependent variable. In this case, the choice of relapse or medication as a dependent variable does matter. Hence, our results derive from models where both relapse and medication are treated as dependent variables. Specifically, we specify a bivariate log odds ratio for current relapse and current medication use and a separate univariate logit component for each of these outcomes. Each of these components contains transitional associations with previous relapse and medication. Such models represent extensions of univariate transitional association models (e.g. Diggle et al. (1994)) and correspond to bivariate transitional models (e.g. Zeger and Liang (1991)). We incorporate changes in transitional associations into the full-data parametric model for final inference, and investigate if these temporal changes are due to learning effects or the impact of drop-out. We also perform residual analyses and sensitivity analyses in the context of missing data patterns.     

An adaptive bivariate kernel smoothing method for determining instars of Austrosimulium tillyardianum (Diptera: Simuliidae) larvae

1. In insects, instar determination is generally based on the frequency distribution of sclerotised body part measurements. Commonly used univariate methods, such as histograms and univariate kernel smoothing, are not sufficient to reflect the distribution of the measurements, because development of sclerotised body parts is multidimensional. 2. This study used an adaptive bivariate kernel smoothing method, based on 10 pairs of separating variables, to differentiate instars of Austrosimulium tillyardianum (Diptera: Simuliidae) larvae in two‐dimensional space. A variable bandwidth matrix was used and separation lines between instars were defined. Using the Crosby growth ratio, Brooks' rule and the new standard recently proposed, larvae were separated into nine instars. It was found that, using the bivariate kernel smoothing method, the clustering accuracy and determination of separation lines as instar class limits were higher than those associated with the univariate kernel smoothing method. With the exceptions of the paired separating variables, head capsule length and antennal segment 3 length (AS3L), the mean probabilities of correct classifications was > 85%. The pair of separating variables that yielded the greatest classification accuracy comprised mandible length (ML) and AS3L, which had mean probabilities of 0.8984. The clustering accuracy was higher for early‐ and late‐instar larvae, but lower for instars 6 and 7. The adaptive bivariate kernel smoothing method was better than univariate methods for instar determination, especially in the detection of divisions between instars and identification of a larval instar.     

A semiparametric mixture model for analyzing clustered competing risks data

A very general class of multivariate life distributions is considered for analyzing failure time clustered data that are subject to censoring and multiple modes of failure. Conditional on cluster-specific quantities, the joint distribution of the failure time and event indicator can be expressed as a mixture of the distribution of time to failure due to a certain type (or specific cause), and the failure type distribution. We assume here the marginal probabilities of various failure types are logistic functions of some covariates. The cluster-specific quantities are subject to some unknown distribution that causes frailty. The unknown frailty distribution is modeled nonparametrically using a Dirichlet process. In such a semiparametric setup, a hybrid method of estimation is proposed based on the i.i.d. Weighted Chinese Restaurant algorithm that helps us generate observations from the predictive distribution of the frailty. The Monte Carlo ECM algorithm plays a vital role for obtaining the estimates of the parameters that assess the extent of the effects of the causal factors for failures of a certain type. A simulation study is conducted to study the consistency of our methodology. The proposed methodology is used to analyze a real data set on HIV infection of a cohort of female prostitutes in Senegal.     

A comparison of bivariate and univariate QTL mapping in livestock populations

This study presents a multivariate, variance component-based QTL mapping model implemented via restricted maximum likelihood (REML). The method was applied to investigate bivariate and univariate QTL mapping analyses, using simulated data. Specifically, we report results on the statistical power to detect a QTL and on the precision of parameter estimates using univariate and bivariate approaches. The model and methodology were also applied to study the effectiveness of partitioning the overall genetic correlation between two traits into a component due to many genes of small effect, and one due to the QTL. It is shown that when the QTL has a pleiotropic effect on two traits, a bivariate analysis leads to a higher statistical power of detecting the QTL and to a more precise estimate of the QTL''s map position, in particular in the case when the QTL has a small effect on the trait. The increase in power is most marked in cases where the contributions of the QTL and of the polygenic components to the genetic correlation have opposite signs. The bivariate REML analysis can successfully partition the two components contributing to the genetic correlation between traits.     

Quantitative risk assessment for multivariate continuous outcomes with application to neurotoxicology: the bivariate case

The neurotoxic effects of chemical agents are often investigated in controlled studies on rodents, with multiple binary and continuous endpoints routinely collected. One goal is to conduct quantitative risk assessment to determine safe dose levels. Such studies face two major challenges for continuous outcomes. First, characterizing risk and defining a benchmark dose are difficult. Usually associated with an adverse binary event, risk is clearly definable in quantal settings as presence or absence of an event; finding a similar probability scale for continuous outcomes is less clear. Often, an adverse event is defined for continuous outcomes as any value below a specified cutoff level in a distribution assumed normal or log normal. Second, while continuous outcomes are traditionally analyzed separately for such studies, recent literature advocates also using multiple outcomes to assess risk. We propose a method for modeling and quantitative risk assessment for bivariate continuous outcomes that address both difficulties by extending existing percentile regression methods. The model is likelihood based; it allows separate dose-response models for each outcome while accounting for the bivariate correlation and overall characterization of risk. The approach to estimation of a benchmark dose is analogous to that for quantal data without the need to specify arbitrary cutoff values. We illustrate our methods with data from a neurotoxicity study of triethyl tin exposure in rats.     

Latent model for correlated binary data with diagnostic error

We propose a methodology for modeling correlated binary data measured with diagnostic error. A shared random effect is used to induce correlations in repeated true latent binary outcomes and in observed responses and to link the probability of a true positive outcome with the probability of having a diagnosis error. We evaluate the performance of our proposed approach through simulations and compare it with an ad hoc approach. The methodology is illustrated with data from a study that assessed the probability of corneal arcus in patients with familial hypercholesterolemia.     

Nonparametric estimation of the bivariate recurrence time distribution

This article considers statistical models in which two different types of events, such as the diagnosis of a disease and the remission of the disease, occur alternately over time and are observed subject to right censoring. We propose nonparametric estimators for the joint distribution of bivariate recurrence times and the marginal distribution of the first recurrence time. In general, the marginal distribution of the second recurrence time cannot be estimated due to an identifiability problem, but a conditional distribution of the second recurrence time can be estimated non-parametrically. In the literature, statistical methods have been developed to estimate the joint distribution of bivariate recurrence times based on data on the first pair of censored bivariate recurrence times. These methods are inefficient in the model considered here because recurrence times of higher orders are not used. Asymptotic properties of the proposed estimators are established. Numerical studies demonstrate the estimators perform well with practical sample sizes. We apply the proposed method to the South Verona, Italy, psychiatric case register (PCR) data set for illustration of the methods and theory.     

Predicting events in clinical trials using two time‐to‐event outcomes

In clinical trials with time‐to‐event outcomes, it is of interest to predict when a prespecified number of events can be reached. Interim analysis is conducted to estimate the underlying survival function. When another correlated time‐to‐event endpoint is available, both outcome variables can be used to improve estimation efficiency. In this paper, we propose to use the convolution of two time‐to‐event variables to estimate the survival function of interest. Propositions and examples are provided based on exponential models that accommodate possible change points. We further propose a new estimation equation about the expected time that exploits the relationship of two endpoints. Simulations and the analysis of real data show that the proposed methods with bivariate information yield significant improvement in prediction over that of the univariate method.     

A unification of multivariate methods for meta-analysis of genetic association studies

Methods for multivariate meta-analysis of genetic association studies are reviewed, summarized and presented in a unified framework. Modifications of standard models are described in detail in order to be applied in genetic association studies. The model based on summary data is uniformly defined for both discrete and continuous outcomes and analytical expressions for the covariance of the two jointly modeled outcomes are derived for both cases. The models based on the binary nature of the data are fitted using both prospective and retrospective likelihood. Furthermore, formal tests for assessing the genetic model of inheritance are developed based on standard normal theory. The general model is compared to the recently proposed genetic model-free bivariate approach (either using summary or binary data), and it is clearly shown that the estimates provided by this approach are nearly identical to the estimates derived by the general bivariate model using the aforementioned tests for the genetic model. The methods developed here as well as the tests, are easily implemented in all major statistical packages, escaping the need of self written software. The methods are applied in several already published meta-analyses of genetic association studies (with both discrete and continuous outcomes) and the results are compared against the widely used univariate approach as well as against the genetic model free approaches. Illustrative examples of code in Stata are given in the appendix. It is anticipated that the methods developed in this work will be widely applied in the meta-analysis of genetic association studies.     

Analysis of bivariate flow karyotypes

Bivariate flow karyotype analysis is performed using data from chromosomes stained with two fluorescent dyes, typically chromomycin A3 and Hoechst-33258, and measured in a flow cytometer or cell sorter (Carrano et al.: Proceedings of the National Academy of Sciences of the United States of America 76:1382-1384, 1979; Gray et al.: Proceedings of the National Academy of Sciences of the United States of America 72:1231-1234, 1975; Langlois et al.: Proceedings of the National Academy of Sciences of the United States of America 79:7876-7880, 1982). In the resulting bivariate histogram, most chromosome types appear as individual peaks. In sorting of chromosomes to purify a specific chromosomal type, its corresponding peak in the bivariate histogram is delineated by a rectangular region which surrounds it. All events (objects) that fall within this region trigger the sorting process. In most cases, peaks for different chromosomal types overlap to some extent, and in addition there is always an underlying background due to chromosome fragments and clumps. Thus the sorted population will not be pure; it may include more than one chromosome type and will include debris. To determine the purity of a sort, i.e., the percentage of the sorted material that is of the actual chromosomal type desired, two methods of mathematical analysis have been developed. In the more general method, the bivariate data within an analysis region that includes the sort region, are fit with a series of bivariate Gaussian functions, one for each peak. In a simplified method, the data within the analysis region are transformed into a univariate distribution of either chromomycin A3 or Hoechst-33258 fluorescence. The peaks in these univariate distributions are fit with univariate Gaussian functions. In both methods the purity is determined mathematically. The results of both methods agree well with independent methods of analysis.     

Analysis of multivariate recurrent event data with time‐dependent covariates and informative censoring

Multivariate recurrent event data are usually encountered in many clinical and longitudinal studies in which each study subject may experience multiple recurrent events. For the analysis of such data, most existing approaches have been proposed under the assumption that the censoring times are noninformative, which may not be true especially when the observation of recurrent events is terminated by a failure event. In this article, we consider regression analysis of multivariate recurrent event data with both time‐dependent and time‐independent covariates where the censoring times and the recurrent event process are allowed to be correlated via a frailty. The proposed joint model is flexible where both the distributions of censoring and frailty variables are left unspecified. We propose a pairwise pseudolikelihood approach and an estimating equation‐based approach for estimating coefficients of time‐dependent and time‐independent covariates, respectively. The large sample properties of the proposed estimates are established, while the finite‐sample properties are demonstrated by simulation studies. The proposed methods are applied to the analysis of a set of bivariate recurrent event data from a study of platelet transfusion reactions.