This study considers analysis of bivariate longitudinal binary data. We propose a model based on marginalized multilevel model framework. The proposed model consists of two levels such that the first level associates the marginal mean of responses with covariates through a logistic regression model and the second level includes subject/time specific random intercepts within a probit regression model. The covariance matrix of multiple correlated time-specific random intercepts for each subject is assumed to represent the within-subject association. The subject-specific random effects covariance matrix is further decomposed into its dependence and variance components through modified Cholesky decomposition method and then the unconstrained version of resulting parameters are modelled in terms of covariates with low-dimensional regression parameters. This provides better explanations related to dependence and variance parameters and a reduction in the number of parameters to be estimated in random effects covariance matrix to avoid possible identifiability problems. Marginal correlations between responses of subjects and within the responses of a subject are derived through a Taylor series-based approximation. Data cloning computational algorithm is used to compute the maximum likelihood estimates and standard errors of the parameters in the proposed model. The validity of the proposed model is assessed through a Monte Carlo simulation study, and results are observed to be at acceptable level. Lastly, the proposed model is illustrated through Mother's Stress and Children's Morbidity study data, where both population-averaged and subject-specific interpretations are drawn through Emprical Bayes estimation of random effects.