Multivariate Regression with Two Different Types of Response Problem Setting:
I have an interesting question related with longitudinal study and multivariate regression. I found that in lots of biomedical studies, multiple discrete and continuous endpoints are very common. For example, continuous fetal weight $(Y_{1i})$ and binary fetal death $(Y_{2i})$ are observed in toxicological study. We are interested in studying the effects of dose levels on fetal weight and fetal death and will jointly modeling these two outcomes together since they may be correlated. If we assume we have a vector of variables $X_i$ and there is correlation between $Y_{1i}$ and $Y_{2i}$, how can we propose a marginal (population-average) model to jointly modeling $(Y_{1i}, Y_{2i})$ by their means as function of $X_i$?
Solving Attempt:
Since there is correlation between two types response $(Y_{1i}, Y_{2i})$ and we want to obtain marginal model, I think it may be good to use generalized estimating equations for obtaining estimated coefficients and then also provide estimators for the covariance of these estimated coefficients. However, I don't know how to mathematically layout the model since there are two correlated response.

To handle such data via generalized estimating equations (GEE), you would need to specify a variance structure that incorporates both continuous and time-to-event outcomes. This might represent my lack of imagination, but I don't see an obvious way to do that.
Joint modeling of this type is based on random-effect or frailty models, as reviewed by Ibrahim et al., for example. So far as I can tell, all the packages under the R Survival Task View heading "Joint modelling of time-to-event and longitudinal data" use such approaches, as do the frailtypack and JM packages listed elsewhere on that page. Marginal estimates are then obtained from the resulting models.
