Multivariate longitudinal data arise frequently in different application fields, where multiple outcomes are measured repeatedly from the same subject. In this paper, we first propose a two-stage weighted least squares estimation procedure for the regression coefficients when the random error follows an irregular autoregressive process, and establish asymptotic normality properties of the resulting estimators. We then apply the SCAD variable selection approach to determine the order of the AR error process. We further propose a test statistic to check whether multiple responses are correlated at the same observation time, and derive the asymptotic distribution of the proposed test statistic. Several simulated examples and a real data analysis are presented to illustrate the finite sample performance of the proposed method.