The analysis of anatomy that undergoes rapid changes, such as neuroimaging of the early developing brain, greatly benefits from spatio-temporal statistical analysis methods to represent population variations but also subject-wise characteristics over time. Methods for spatio-temporal modeling and for analysis of longitudinal shape and image data have been presented before, but, to our knowledge, not for diffusion weighted MR images (DW-MRI) fitted with higher-order diffusion models. To bridge the gap between rapidly evolving DW-MRI methods in longitudinal studies and the existing frameworks, which are often limited to the analysis of derived measures like fractional anisotropy (FA), we propose a new framework to estimate a population trajectory of longitudinal diffusion orientation distribution functions (dODFs) along with subject-specific changes by using hierarchical geodesic modeling. The dODF is an angular profile of the diffusion probability density function derived from high angular resolution diffusion imaging (HARDI) and we consider the dODF with the square-root representation to lie on the unit sphere in a Hilbert space, which is a well-known Riemannian manifold, to respect the nonlinear characteristics of dODFs. The proposed method is validated on synthetic longitudinal dODF data and tested on a longitudinal set of 60 HARDI images from 25 healthy infants to characterize dODF changes associated with early brain development.