Longitudinal regression analysis for clinical imaging studies is essential to investigate unknown relationships between subject-wise changes over time and subject-specific characteristics, represented by covariates such as disease severity or a level of genetic risk. Image-derived data in medical image analysis, e.g. diffusion tensors or geometric shapes, are often represented on nonlinear Riemannian manifolds. Hierarchical geodesic models were suggested to characterize subject-specific changes of nonlinear data on Riemannian manifolds as extensions of a linear mixed effects model. We propose a new hierarchical multi-geodesic model to enable analysis of the relationship between subject-wise anatomical shape changes on a Riemannian manifold and multiple subject-specific characteristics. Each individual subject-wise shape change is represented by a univariate geodesic model. The effects of subject-specific covariates on the estimated subject-wise trajectories are then modeled by multivariate intercept and slope models which together form a multi-geodesic model. Validation was performed with a synthetic example on a S2 manifold. The proposed method was applied to a longitudinal set of 72 corpus callosum shapes from 24 autism spectrum disorder subjects to study the relationship between anatomical shape changes and the autism severity score, resulting in statistics for the population but also for each subject. To our knowledge, this is the first longitudinal framework to model anatomical developments over time as functions of both continuous and categorical covariates on a nonlinear shape space.