In radiotherapy for prostate cancer irradiation of neighboring organs at risk may lead to undesirable side-effects. Given this setting, the bladder presents the largest inter-fraction shape variations hampering the computation of the actual delivered dose vs. planned dose. This paper proposes a population model, based on longitudinal data, able to estimate the probability of bladder presence during treatment, using only the planning computed tomography (CT) scan as input information. As in previously-proposed principal component analysis (PCA) population-based models, we have used the data to obtain the dominant eigenmodes that describe bladder geometric variations between fractions. However, we have used a longitudinal analysis along each mode in order to properly characterize patient's variance from the total population variance. We have proposed is a mixed-effects (ME) model in order to separate intra- and inter-patient variability, in an effort to control confounding cohort effects. Other than using PCA, bladder shapes are represented by using spherical harmonics (SPHARM) that additionally enables data compression without information lost. Based on training data from repeated CT scans, the ME model was thus implemented following dimensionality reduction by means of SPHARM and PCA. We have evaluated the model in a leave-one-out cross validation framework on the training data but also using independent data. Probability maps (PMs) were thus generated with several draws from the learnt model as predicted regions where the bladder will likely move and deform. These PMs were compared with the actual regions using metrics based on mutual information distance and misestimated voxels. The prediction was also compared with two previous population PCA-based models. The proposed model was able to reduce the uncertainties in the estimation of the probable region of bladder motion and deformation. This model can thus be used for tailoring radiotherapy treatments.