One major challenge of using the phylogenetic comparative method (PCM) is the analysis of the evolution of interrelated continuous and discrete traits in a single multivariate statistical framework. In addition, more intricate parameters such as branch-specific directional selection have rarely been integrated into such multivariate PCM frameworks. Here, originally motivated to analyze the complex evolutionary trajectories of group size (continuous variable) and social systems (discrete variable) in African subterranean rodents, we develop a flexible approach using approximate Bayesian computation (ABC). Specifically, our multivariate ABC-PCM method allows the user to flexibly model an underlying latent evolutionary function between continuous and discrete traits. The ABC-PCM also simultaneously incorporates complex evolutionary parameters such as branch-specific selection. This study highlights the flexibility of ABC-PCMs in analyzing the evolution of phenotypic traits interrelated in a complex manner.