Chen and Wehrly [1] have recently proposed intra-correlation, inter-correlation, and total correlation coefficients to assess the correlations under different circumstances for multivariate clustered data with mixed outcomes under multivariate generalized linear mixed models (MGLMMs). However, they restricted their theoretical derivations and numerical computations to joint modeling of clustered data with Gaussian and non-Gaussian continuous outcomes (e.g., exponential and gamma variables) and count outcomes (e.g., Poisson variables) because they were able to derive closed-form solution for marginal variance of outcomes and that for marginal covariance and correlation between outcomes. This was due to the good collaboration of identity and log link functions with normally distributed random effects within MGLMMs framework.We would like to bring to the authors' attention the binary responses which are also commonly observed in clustered data studies. As pointed out by McCulloch [2] and Inan [3], when the link function associating binary responses with covariates is the logit link function and the distribution for the random effects is the normal distribution within generalized linear mixed models, closed-form expression for marginal variance, covariance, and correlation functions cannot be obtained. , we can show that these functions can be approximated through a first-order Taylor series expansion.To illustrate, we follow the conventional notation used in Chen and Wehrly [1] and assume that Y ijt , which is the binary outcome of subject i at measurement t for method j, is associated with some set of covariates x ijt = (x 1,ijt , … , x p j ,ijt ) T through the logit link function. Then, a MGLMM for bivariate clustered binary data with subject and response-specific random-intercepts can be specified as follows [7]:where method j = 1, 2, t = 1, 2, … , T i , ijt is the conditional mean of Y ijt , j = ( j1 , … , jp j ) T is the vector of response-specific regression parameters. The term i = (b i1 , b i2 ) is the vector of subjectspecific random-intercepts coming from a standard multivariate normal distribution, where the covariance matrix is = ) .