Complex structured data settings are studied where outcomes are multivariate and multilevel and are collected longitudinally. Multivariate outcomes include both continuous and discrete responses. In addition, the data contain a large number of covariates but only some of them are important in explaining the dynamic features of the responses. To delineate the complex association structures of the responses, a model with correlated random effects is proposed. To handle the large dimensionality of covariates, a simultaneous variable selection and parameter estimation method is developed. To implement the method, a computationally feasible algorithm is described. The proposed method is evaluated empirically by simulation studies and illustrated by analyzing the data arising from the Waterloo Smoking Prevention Project.