Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Biometrics. Author manuscript; available in PMC 2010 June 23.
Published in final edited form as:
PMCID: PMC2890259

Association Models for Clustered Data with Binary and Continuous Responses


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.

Keywords: Bivariate binary and continuous responses, Bridge distribution, Logit link, Probability integral transformation

1. Introduction

In many biomedical studies, the basic sampling unit is a cluster of subjects, such as siblings in a family or rats in a litter. Further, each member of the cluster may have multiple outcomes, including both discrete and continuous endpoints. When the bivariate responses are both continuous and approximately Gaussian, the traditional random effects linear model framework can be used for inference. However, the presence of both categorical/ordinal as well as continuous outcomes, called mixed responses, complicates the situation due to the lack of a discrete multivariate analog of the multivariate Gaussian distribution (Fitzmaurice and Laird, 1995). For example, the developmental toxicity study of ethylene glycol (EG), a high-volume industrial chemical (Price et al., 1985), involves measurements of both binary and continuous responses, namely, (a) fetal malformation and (b) fetal weight of 1028 fetuses from 94 litters. In this study, one of the four levels of EG was administered to 94 pregnant mice over the period of major organogenesis, beginning just after implantation. In particular, 25 pregnant mice received 0 g/kg/day, 24 received 0.75 g/kg/day, 22 received 1.5 g/kg/day, and 23 received 3 g/kg/day. We expect a within cluster (“litter”) association, along with the direct association between the fetal malformation (binary) and fetal weight (continuous) of same fetus. We are primarily interested in modeling the effects of dose on the binary and continuous outcomes, as interpretable regression functions, as well as modeling the within subject and within litter associations between the binary and continuous outcomes.

A number of joint models have been proposed for developmental toxicity studies involving both discrete and continuous responses from each member of a cluster (Geys, Molenberghs, and Ryan, 1999; Regan and Catalano, 1999; Dunson, Chen, and Harry, 2003; Gueorguieva and Agresti, 2001). Catalano and Ryan (1992) suggested a model assuming that the discrete outcome is a discretization of an unobserved continuous latent variable and this latent variable and the continuous response within same subject have a bivariate normal distribution. Joint distribution of the two responses is then formed by a product of marginal distribution of the continuous response and conditional distribution of the binary response conditioned on the continuous response. The other type of models consider conditioning on the binary response. Fitzmaurice and Laird (1995) assumed a logit link function for marginal probability of binary response and a normal distribution for the continuous response given the binary response. Cox and Wermuth (1992) compared different joint distribution models for analyzing data with quantitative and qualitative responses. For more discussions related to conditional models, see Little and Schluchter (1985), Cox and Wermuth (1994), Krzanowski (1988), and Schafer (1997).

The model proposed by Fitzmaurice and Laird (1995), henceforth FL, focuses on estimating and interpreting the marginal regression functions of two responses, while treating the association between the discrete and continuous responses as a nuisance. FL models have the attractive feature that the marginal regression parameters of both responses have good physical interpretations. However, this model only considers the cluster effect on the continuous response given the binary responses from all cluster members. In this model, the specification of the regression model conditional on the cluster-specific random effects is completely unknown. The magnitude of the marginal regression parameters, particularly for the binary response, gets attenuated from that of the conditional regression parameters due to the heterogeneity of clusters. This degree of attenuation is important for future study design as well as understanding the expected effect of covariates on other populations and on individual subjects.

In this article, we avoid the above shortcomings of FL modeling. We propose a model with practical physical interpretations of the marginal as well as conditional (given the cluster-specific random effects) regression functions associated with both discrete and continuous responses. We specify a bivariate distribution of two cluster-specific random effects—presented as separate random intercepts for two types of responses. We use a separate parameter to account for the direct within-subject association between the two responses. Hence, association between the discrete and continuous responses from different subjects within the same cluster is induced by the two correlated random effects. Association between the binary and continuous outcomes on the same subject is induced by correlations between the two within-cluster random effects, as well as the direct dependence of the continuous response on the binary response from the same subject. We use a logit link for the conditional regression function of the binary response given cluster-specific random effect, and an identity link for the regression function of the continuous response conditioned on the binary response and the cluster random effects. The joint distribution of the discrete and continuous responses given the random effects is expressed by the product of the marginal distribution of the binary outcome and the conditional distribution of the continuous response given the binary response. In the FL model, they considered maximum likelihood (henceforth ML) estimation as quite complicated and intractable, and hence they proposed the generalized estimating equations (GEE) method for computational convenience. Contrary to that, our likelihood is much simpler to evaluate and amenable to maximization through routine nonadaptive Gaussian quadrature techniques, because conditional on the random effects, all outcomes are independent.

One drawback of the random effects models for clustered mixed responses is that the marginal regression model of the binary response (after integrating out the random effects) usually is not of logistic form when the subject-specific model conditional on the random effect is of logistic form. In our model, the bridge distribution (Wang and Louis, 2003) of the random effect for the binary response allows the marginal probability of the binary response, integrated over the random intercept, to have a logistic structure with an odds ratio interpretation of the marginal regression effect. The regression parameters in the marginal logistic regression model are proportional to the corresponding regression parameters in the subject-specific conditional logistic model. This is a major advantage over other existing models including the FL model.

The rest of the article is organized as follows. In Section 2, we begin with our random effects model of clustered data with bivariate mixed responses within each subject. We illustrate how different levels of the association between two random effects induce different Kendall's τ values for association between binary and continuous responses from same cluster. In Section 3, we extend this conditional model of the clustered mixed responses to accommodate a direct within subject association between two responses (even after adjusting for the within cluster association). We also give the likelihood representation of our model. In Section 4, we demonstrate the ML estimation of the regression parameters of our models via analyzing the data from a developmental toxicity experiment of EG on mice. Finally, Section 5 discusses the need and advantages of our model. We also study the robustness of our method to model and random effect distribution misspecification through simulation studies.

2. Random Effects Model for Mixed Responses

For subject j = 1,, Mi from cluster i = 1,, N, we have a binary response Y1ij, where Y1ij = 1 when subject j has response 1 and 0 otherwise, and a continuous response Y2ij. Each subject has a K × 1 covariate vector Xij. In this section, we use the usual linear model with normal random (cluster) effect Wi for the continuous outcomes Y2ij from cluster i. Furthermore, for the binary outcomes Y1ij from cluster i, we let the cluster random effects Bi have the bridge density of Wang and Louis (2003). The key to the new approach proposed here is a bivariate density for the joint distribution of the normal random effect Wi and the bridge random effect Bi. This bivariate distribution will induce an association between the vector of continuous and the vector of binary outcomes within the same cluster.

For the time being, we assume that given the cluster-specific random effects (Bi, Wi), the responses Ykij for k = 1, 2; j = 1,, Mi within cluster i are mutually independent. We assume the conditional distribution of binary Y1ij given the random effects (Bi, Wi) to be Bernoulli with success probability pij, with


where the marginal distribution of the random intercept Bi of binary response is a bridge distribution (Wang and Louis, 2003) with unknown parameter ϕ (0 < ϕ < 1) and density


The bridge density is symmetric (mean of Bi is 0), and the variance of Bi is σb2=π2(φ21)/3. It has a slightly heavier tail and is more peaked than the normal density. The marginal probability P[Y1ij = 1 |Xij] (when integrated over the random intercepts (Bi, Wi)) still has logistic form. The parameter ϕ of the bridge distribution being the same across all clusters assures the exchangeability of cluster random effects on the binary responses. The marginal regression parameter of the binary response is proportional to the conditional regression parameter β multiplied by the attenuation parameter ϕ, i.e.,


Small ϕ indicates the important facts: high correlation of the binary responses within a cluster; high heterogeneity among clusters related to the binary response; and a corresponding high degree of attenuation of β in the population due to cluster heterogeneity.

The continuous response Y2ij given the cluster random effect (Bi, Wi) is modeled by a normally distributed linear mixed effects model with mean


and variance σ02. The marginal distribution of the random intercept Wi of the continuous response is normal with mean 0 and variance σ12. The marginal distribution of Y2ij integrated over (Bi, Wi) is again normal with mean


In this model for the clustered mixed outcomes, both the conditional (on the cluster-specific random effects) and marginal regression functions in equations (1)–(4) have useful physical interpretations. As we describe below, the association between a binary and a continuous outcome within same cluster is induced by formulating a joint bivariate distribution for the cluster random effects (Bi, Wi) while preserving the desired forms of the marginal densities of Bi and Wi. This bivariate distribution is constructed from transformations of a pair of latent random variables (Z1i, Z2i) with a bivariate normal distribution with mean zero, unit standard deviations, and correlation ρ.

To formulate the bivariate distribution for (Bi, Wi), we use the probability integral transformation


where Φ(·) is the cumulative distribution function (c.d.f.) of standard normal, and Fb1 (·) is the inverse cumulative distribution


of the bridge density for 0 < u < 1. This assures that the marginal density of Bi will be the bridge distribution with c.d.f.


Because both the bridge cumulative and inverse cumulative distributions have closed-form expressions, they can easily be implemented in standard statistical software and simulations are easy to perform. In order to obtain WiN(0,σ12), we let Wi = σ1Z2i.

In an extensive simulation study, we find the correlation between the two random intercepts Bi and Wi are approximately the same as the correlation between Z1i and Z2i. To understand the associations within (Bi, Wi) and the association within (Y1ij, Y2ij) induced via the correlation ρ of (Z1i, Z2i), we compute Kendall's τ (Hougaard, 2000) for (Bi, Wi) and Kendall's τ for (Y1ij, Y2ij′) for different values of ρ. The transformations Bi=Fb1 (Φ(Z1i)) and Wi = σ1Z2i are both monotone increasing. Kendall's τ for bivariate normal random variables (Z1i, Z2i) with correlation ρ is τ = 2 arcsin(ρ) with −1 ≤ τ ≤ 1. This is also Kendall's τ for (Bi, Wi) because Kendall's τ is invaraint to increasing transformations (Hougaard, 2000). Kendall's τ for (Y1ij, Y2ij′), defined as E[sign{(Y1ijY1i* j*)(Y2ijY2i* j**)}], was simulated via Monte Carlo methods using equations (1), (3), (5), (6), and Wi = σ1Z2i. Figure 1 shows values of Kendall's τ of (Bi, Wi) for different values of ρ. Figure 2 shows Kendall's τ of the binary and continuous responses (Y1ij, Y2ij′) from the same cluster for different values of ρ and a fixed value of ϕ.

Figure 1
Kendall's τ of Bi and Wi for different ρ.
Figure 2
Kendall's τ of Y1ij and Y2ij with ϕ = 0.735 for different ρ.

3. Models for Bivariate Mixed Responses with Clustering

In Section 2, models (1) and (3) induce an exchangeable association structure among binary and continuous responses in the same cluster because given cluster-specific (Bi, Wi), responses within a subject are considered independent. However, we might expect the association between a binary Y1ij and a continuous Y2ij from the same subject to be stronger than the association between Y1ij and Y2ij′ from different subjects j and j′ within cluster i. In this section, we extend our model to accommodate the possibility of the direct association between binary and continuous outcomes from the same subject, even after adjusting for the within cluster association. To this end, we modify the model (3) by incorporating a direct effect of subject-specific binary Y1ij on the continuous Y2ij, such that


where π1ij=E[Y1ij|Xij]=exp{Xij'(φβ)}/[1+exp{Xij(φβ)}], while keeping the regression model for Y1ij identical to that of equation (1). Hence, the association between Y1ij and Y2ij from the same subject is induced by the direct dependence of Y2ij on Y1ij measured by γ as well as the association parameter ρ of the random effects (Bi, Wi).

The conditional expectation of Y2ij integrated over Wi is


because E[Wi | Y1ij, Xij] = EE[Wi | Bi, Y1ij, Xij] = EE [Wi | Bi] = 0. Hence, marginal expectation E[Y2ij|Xij]=EE[Y2ij|Y1ij,Xij]=Xijα is the linear function of Xij, same as in equation (4).

From equation (7), the covariance between responses Y1ij and Y2ij from the same subject is


The covariance of two responses from different subjects j and j′ within the same cluster i is


where Var(Y1ij) = π1ij(1 − π1ij). Our model ensures the anticipated result | Cov[Y1ij, Y2ij] | >| Cov[Y1ij, Y2ij] because Var(Y1ijY1ij) = 2Var(Y1ij) − 2Cov[Y1ij, Y1ij] > 0. For covariate vectors Xij and Xij, we have


where Eb and Ebw, respectively, are the expectation with respect to the marginal density of Bi and that with respect to the joint density of (Bi, Wi).

Now, the likelihood contribution from cluster i is proportional to


where fY1ij | Bi, Wi, Xij (y1ij | bi,wi, xij) is equal to exp{(xij'β+bi)y1ij}/{1+exp(xijβ+bi)}, and fY2ij | Bi, Wi, Y1ij (y2ij | bi, wi, y 1ij, xij) is the normal density with variance σ02 and mean xijα+γ(y1ijπ1ij)+wi.

Unlike the FL model, our likelihood is tractable and maximization is straightforward. In order to use the ML method for estimation, the likelihood is obtained via taking the expectation with respect to the joint density of (Bi, Wi). These integrals are implemented using nonadaptive Gaussian quadrature techniques in PROC NLMIXED of SAS (v9.1) along with optimization tools such as Newton–Raphson with ridging.

4. Example: EG Data

We apply the method described above to analyze data from the development toxicity study using 94 pregnant mice randomized to receive one of the four different doses of EG, 0, 0.75, 1.5, or 3 g/kg/day. Fetal malformation and fetal weight of 1028 fetuses from the 94 litters/clusters, with cluster size ranging from 1 to 16, were recorded. A summary of the fetal malformations and fetal weights is given in Table 1. We clearly observe that the proportion of malformation increases from 0.34% to 57.08% and sample mean of weight decreases from 0.972 to 0.704 as dose level increases from 0 to 3. For more details on the specifics of this data set, see Catalano and Ryan (1992), Fitzmaurice and Laird (1995), and also the monograph by Aerts et al. (2002).

Table 1
Descriptive statistics of data from developmental toxicity study

Treating dose Xi(= Xij) assigned to the ith mice/litter as a covariate, the binary Y1ij (fetal malformation) and the continuous Y2ij (fetal weight) in model 1 are modeled as logit (E [Y1ij | Bi, Wi, Xi]) = β0 + β1Xi + Bi and E[Y2ij | Bi, Wi, Y1ij, Xi] = α0 + α1Xi + γ (Y1ijπ1i) + Wi, where π1i = exp{ϕ(β0 + β1Xi)}/[1 + exp{ϕ(β0 + β1Xi)}], and Y2ij given (Bi, Wi, Y1ij) is normally distributed with variance σ02. The random effects (Bi, Wi) are assumed to have a bivariate density with the marginal density of Bi being the bridge with parameter ϕ, and the marginal density of Wi being normal with mean 0 and variance σ12.

Table 2 shows the regression parameter estimates for the two responses from model 1. The positive estimate for the dose effect on fetal malformation and the negative estimate on fetal weight are consistent with what we observe from Table 1: the malformation rate increases and weight decreases as dose level increases. All of the marginal parameter estimates have slightly different point estimates than the FL model (Table 2 in Fitzmaurice and Laird, 1995), but the estimated standard errors of the parameter estimates are smaller than that of the FL model, especially for the binary response (approximately 20%), leading to tighter 95% confidence intervals.

Table 2
ML estimates of parameters from model 1 and model 2

The estimate of the attenuation parameter ϕ = 0.6618 confirms that heterogeneity of the clusters causes moderate attenuation of the estimated marginal dose effects. The intracluster correlation ρY1 between the binary responses is estimated by 1 − ϕ = 0.3382 (Wang and Louis, 2003). The intracluster correlation ρY2 between the continuous responses are estimated by the between-cluster variability divided by the sum of the within- and between-cluster variabilities σ12/(σ02+σ12)=0.6262. These estimates are somewhat greater than the estimates of the intracluster correlations obtained from the FL model. The negative value −0.6433 for estimate of association parameter ρ between the two cluster-specific random effects Bi and Wi, as well as −0.03369 for estimate of the within subject dependence coefficient γ of Y2ij on Y1ij, suggest reverse association among the two responses fetal malformation and fetal weight. For comparison with the results from model 1, we consider a more complex correlation structure. We change the parameter γ of model 1 to (γ1 + γ2Xi) to verify whether the dose effect on fetal weight is the same for both the normal and malformed fetus. We obtain the p-value for γ2 as 0.7154, which indicates that the direct effect of malformation on weight does not depend on dose.

In smoothed residual plots of model 1 obtained via the LOESS procedure in SAS (Cleveland, 1979), we observe some quadratic trends of dose for both the malformation and weight responses. We reanalyze the data by adding quadratic terms of dose effect for both responses (model 2). Estimates of the parameters of model 2 are in Table 2, where ϕβ2 and α2 are the marginal parameters of squared dose effect for the two responses. No further trends are detected from corresponding residual plots. The estimates of model 2 also suggest that the quadratic dose effect terms are necessary.

Our method allows fractional polynomial regression equations (Royston and Altman, 1994; Geys et al., 1999) for E(Y1ij| Xi) and E(Y2ij| Xi). We have fitted several such models including models using (Xi,Xi),(1/Xi,1/Xi) to our data; however, the corresponding Akaike's information criterion (AIC) and Bayesian information criterion (BIC) values are bigger than those of model 2. For experiments with more dose levels than we have, a fractional polynomial model may become more useful for the dose-response effect.

To see the impact of association between fetal malformation Y1ij and fetal weight Y2ij on the estimates of the regression parameters, we fit separate and independent regression models for Y1ij and Y2ij, i.e., ρ = 0 and γ = 0 (model 0a). We modify model 1 with γ = 0 (model 0b) to accommodate the special case no direct within-subject association between Y1ij and Y2ij. We also consider separate and independent regression models for Y1ij and Y2ij (ρ = 0 and γ = 0) with covariate vector (Xi,Xi2) (model 0c). No noticeable changes in point estimates are observed comparing results of model 0a and 0b with model 1, and model 0c with model 2. However, for model 0a, we observe approximately 35% and 25% larger estimated standard errors for parameters of the binary response Y1ij and continuous response Y2ij, respectively, compared to corresponding estimates from model 1. This indicates that substantially larger sample sizes are required while using model 0a to obtain similar standard errors as model 1. Similarly, 14% and 3% larger standard errors are observed when comparing results for the regression parameters of Y1ij and Y2ij from model 0b to model 1. This is also the case for model 0c compared to model 2. Particularly for this data set, the p-values for the association ρ between the two random effects and direct association γ of Y1ij on Y2ij are both less than 0.0001 implying statistically significant association between the two responses. AIC and BIC of the five fitted models are shown in Table 3. Model 2 has by far the smallest AIC and BIC values (with differences larger than 10 compared to model 1), which indicates it is the best fitted model for our data set.

Table 3
Fit statistics for five models

Using estimates from model 2 and expressions (8)–(11), we obtain from Monte Carlo simulation that, for dose level 0, the estimated correlation between two responses from the same fetus is Corr^[Y1ij,Y2ij]=0.1481, and the estimated correlation between two responses from different fetuses within the same cluster is Corr^[Y1ij,Y2ij]=0.1245. This shows that the former correlation is approximately 19% higher than the latter correlation.

5. Discussion

In this article, we have developed a regression model for bivariate discrete and continuous outcomes motivated by an application in developmental toxicology where the experiment measures the fetal malformation and fetal weight for different dose levels administered to different clusters. Our model is attractive in that it provides both conditional and marginal interpretation of the dose-response modeling, which happens to be the cornerstone of quantitative risk assessment. Similar to a generalized linear mixed model setup, the binary responses are related to the covariates through a logit link whereas the conditional distribution of the continuous responses given the binary responses are related to the covariates by an identity link. Our model can efficiently handle the effects of clustering through introduction of cluster-specific random effects in a regression setup, as well as the association between mixed outcomes on the same (different) fetus.

Our method of estimation is ML as opposed to the GEE methods adapted in the FL approach. With a correctly specified model, provided estimation techniques are computationally feasible, the ML method reigns supreme over any other estimation methods for its fidelity to the likelihood principle as well as for its asymptotic properties. There are clearly several advantages of our ML method over the GEE method used in the FL approach. Firstly, although ML estimation appears to be complicated in a clustered mixed responses setting, our model is easily amenable to ML estimation through routine optimization techniques readily available in standard software like SAS (v9.1). Secondly, ML estimates are asymptotically efficient (Lehmann and Casella, 1998). Thus, using our method we achieve smaller standard errors of the parameter estimates and consequently tighter confidence intervals over the estimates obtained by the GEE method. Thirdly, GEE methods provide only marginal models, thus we are unable to measure the attenuation of the dose effect to the population due to the heterogeneity induced by clustering (litter effect). Our model has both conditional and marginal interpretation in this context.

In theory, the only advantage GEE has over our full-likelihood approach is when the marginal regression model is correctly specified, but the full likelihood is misspecified; in this case, the GEE estimate of the marginal regression parameters will be asymptotically unbiased, but those from our full likelihood could be biased. To investigate the robustness of our ML approach, we performed small-scale simulation studies under misspecification of regression function as well as misspecification of the random effect distribution for binary response. Our current model (model 1) assumes that the parameter (γ) measuring the direct association between the two responses within the same subject does not depend on dose (or other covariates). We generate 50 data sets from a modification of this model with γ replaced by γ1 + γ2Xi assuming γ1 = γ2 for simplicity, and using estimated values from results of model 1. Each data set contains 100 clusters with fixed cluster size 10. Both our likelihood approach and FL's approach using a GEE method (one-step) are applied to the simulated data. The results for estimates of the marginal dose effect are shown in Table 4. We also fit our model to data simulated from an extension of model 1, in which the distribution of Bi of the binary response is standard normal instead of bridge. The average standard error (with corresponding relative bias) of the estimated dose effects on fetal malformation (conditional) and fetal weight (marginal) from model 1 are 0.1543 (0.0748) and 0.0079 (0.0997), respectively. In this small-scale simulation, we observe that our full likelihood approach is robust to regression model misspecifications as well as random effects distributional misspecification compared to the FL approach as far as the bias of regression estimates are concerned. This makes our approach an extremely attractive alternative to the GEE based FL approach.

Table 4
Simulation results of misspecification of regression function
Table thumbnail

Supplementary Material



The authors are grateful for the support provided by grants ROI-CA69222 and ROI-AI60373 from the National Institutes of Health.


Supplementary Materials: Web Appendix derivations of equations (10) and (11) for Xij = 0, as well as the EG data analyzed in Section 4 and an SAS macro implementing the ML method are available under the Paper Information link at the Biometrics website


  • Aerts M, Geys H, Molenberghs G, Ryan LM. Topics in Modelling of Clustered Data. Boca Raton, Florida: Chapman and Hall; 2002.
  • Catalano P, Ryan LM. Bivariate latent variable models for clustered discrete and continuous outcomes. Journal of the American Statistical Association. 1992;87:651–658.
  • Cleveland WS. Robust locally weighted regression and smoothing scatterplots. Journal of the American Statistical Association. 1979;74:829–836.
  • Cox DR, Wermuth N. Response models for mixed binary and quantitative variables. Biometrika. 1992;79:441–461.
  • Cox DR, Wermuth N. Multivariate Dependencies: Models, Analysis and Interpretation. London: Chapman and Hall; 1994.
  • Dunson DB, Chen Z, Harry J. A Bayesian approach for joint modeling of cluster size and subunit-specific outcomes. Biometrics. 2003;59:521–530. [PubMed]
  • Fitzmaurice GM, Laird NM. Regression models for a bivariate discrete and continuous outcome with clustering. Journal of the American Statistical Association. 1995;90:845–852.
  • Geys H, Molenberghs G, Ryan LM. Pseudolikelihood modeling of multivariate outcomes in developmental toxicology. Journal of the American Statistical Association. 1999;94:734–745.
  • Gueorguieva R, Agresti A. A correlated probit model for joint modeling of clustered binary and continuous responses. Journal of the American Statistical Association. 2001;96:1102–1112.
  • Hougaard P. Analysis of Multivariate Survival Data. New York: Springer; 2000.
  • Krzanowski WJ. Principles of Multivariate Analysis. Oxford: Clarendon Press; 1988.
  • Lehmann EL, Casella G. Theory of Point Estimation. 2nd edition. New York: Springer; 1998.
  • Little RJA, Schluchter MD. Maximum likelihood estimation for mixed continuous and categorical data with missing values. Biometrika. 1985;72:497–512.
  • Price CJ, Kimmel CA, Tyl RW, Marr MC. The developmental toxicity of ethylene glycol in rats and mice. Toxicological Applications in Pharmacology. 1985;81:113–127. [PubMed]
  • Regan MM, Catalano PJ. Likelihood models for clustered binary and continuous outcomes: Applications to developmental toxicology. Biometrics. 1999;55:760–768. [PubMed]
  • Royston P, Altman DG. Regression using fractional polynomials of continuous covariates: Parsimonious parametric modelling. Applied Statistics. 1994;43:429–467.
  • Schafer JL. Analysis of Incomplete Multivariate Data. London: Chapman and Hall; 1997.
  • Wang Z, Louis TA. Matching conditional and marginal shapes in binary mixed-effects models using a bridge distribution function. Biometrika. 2003;90:765–775.