When analyzing multivariate longitudinal binary data, we estimate the effects on the responses of the covariates while accounting for three types of complex correlations present in the data. These include the correlations within separate responses over time, cross‐correlations between different responses at different times, and correlations between different responses at each time point. The number of parameters thus increases quadratically with the dimension of the correlation matrix, making parameter estimation difficult; the estimated correlation matrix must also meet the positive definiteness constraint. The correlation matrix may additionally be heteroscedastic; however, the matrix structure is commonly considered to be homoscedastic and constrained, such as exchangeable or autoregressive with order one. These assumptions are overly strong, resulting in skewed estimates of the covariate effects on the responses. Hence, we propose probit linear mixed models for multivariate longitudinal binary data, where the correlation matrix is estimated using hypersphere decomposition instead of the strong assumptions noted above. Simulations and real examples are used to demonstrate the proposed methods. An open source R package, BayesMGLM, is made available on GitHub at https://github.com/kuojunglee/BayesMGLM/ with full documentation to produce the results.