We propose a new method for constructing a spatiotemporal statistical model of the distribution of anatomical landmarks (LMs) of a human embryo. This method exhibits potential for the quantitative assessment of the extent of anomalies and is important in the research of congenital malformations. However, a few of the LMs might not be observed at a specific developmental stage because large morphological deformations exist during the early stages of development. It is difficult for conventional statistical shape analysis methods to handle missing LMs in the training dataset. The basic concept of the proposed method is to conduct statistical analyses by predicting and completing the coordinates of the missing LMs. We demonstrated the proposed method in the context of spatiotemporal statistical modeling of 10 LMs on the brain surface using 37 embryonic subjects with Carnegie stages of 19–22. We conducted a comparative study of the spatiotemporal statistical models between four different prediction methods, and we found that deformable surface mapping was the best prediction method in terms of model generalization and specificity.