In this paper, we propose an algorithm to construct a spatiotemporal statistical shape model (SSM) of the brain surface of a human embryo. Our model is based on level sets and the model is constructed independently for each Carnegie stage (CS), in which feature space is constructed by spatially-weighted principal component analysis (PCA). The statistics of the shape variation were calculated under the assumption of q-Gaussian distribution to reduce the risk of overfitting for datasets of small sizes. To define the statistics of the shape distribution of the intermediate CS, we consider 18 interpolation methods, which are the combinations of three average interpolations (linear, B-spline and information geometry) and six covariance interpolation (rotation, affine-invariant, log-Euclidean, information geometry, Wasserstein geometry and tensor B-spline) techniques, which are exhaustively compared in the experiments. We also propose a method to allow interpolation between distributions with a mixed number of dimensions, i.e., the rank of the covariance matrix. The SSM was constructed and evaluated using 60 sets of brain labels acquired from human embryos from the Kyoto collection. We found that the best interpolation method was a combination of the linear average interpolation method with an information geometry–based covariance interpolation method. We also found that use of the q-Gaussian distribution and selection of the number of dimensions are effective methods for improving the performance of the spatiotemporal SSM of the human embryo.