© 2022 Elsevier Inc.In the analysis of repeated or clustered measurements, it is crucial to determine the dynamics that affect the mean, variance, and correlations of the data, which will be possible using appropriate models. One of these models is the joint mean–covariance model, which is a multivariate heteroscedastic regression model with autoregressive covariance structures. In these models, parameter estimation is usually carried on under normality assumption, but the resulting estimators will be very sensitive to the outliers or non-normality of data. In this study, we propose a robust alternative method and an EM algorithm for estimating the parameters of joint mean–covariance models. Robustification is achieved using a multivariate heavy-tailed distribution with the same number of parameters as the multivariate normal distribution. To simplify the estimation procedure, a modified Cholesky decomposition is adopted to factorize the dependence structure in terms of unconstrained autoregressive and scale innovation parameters. Also, the technique for the prediction of future responses is given. An intensive simulation study and a real data example are provided to demonstrate the performance of the proposed method.