Joint modelling of mean-covariance structure is an important topic in clustered data analysis. Existing methods, such as those based on modified Cholesky decomposition (MCD), alternative Cholesky decomposition (ACD) and hyperspherical coordinates decomposition (HPC), have two main restrictions. First, they often assume that responses in the same cluster are naturally ordered, for example, by time in longitudinal studies. Second, the existing methods model transformed parameters, for instance, the generalized autoregressive parameters and innovation variances in MCD/ACD, and the hyperspherical coordinates in HPC, making the dependence of correlation coefficients or variances on covariates hardly understandable. As an alternative, a data-driven method that models directly the mean, variances and correlation coefficients for clustered data is proposed. Comparing to the existing methods, the proposed approach not only has no need of natural order in responses but also works on original correlation coefficients and variances. The proposed models are flexible and interpretable, and the parameter estimators in joint generalized estimating equations (GEE) are shown to be consistent and asymptotically normally distributed. Consistent model selection criteria in spirit of quasi-likelihood under independence model criterion (QIC) are considered. The use of the proposed approach is demonstrated by intensive simulation studies and real data analysis.