Multilevel functional principal components analysis (MFPCA) facilitates estimation of hierarchical covariance structures for functional data produced by wearable sensors, including continuous glucose monitors (CGM), all while accounting for covariate effects. There are several existing methods to efficiently fit these types of models, including the eminent fast covariance estimation and the recently proposed procedure of fitting appropriate localized mixed effects models and smoothing (Xiao et al. 2016, Leroux et al. 2023). However, these methods do not inherently account for uncertainty in the eigenfunctions during the estimation procedure. Most rely on bootstrap or asymptotic analytic results to perform inference after estimation. Towards this end, we fit MFPCA within a fully-Bayesian framework using MCMC, treating the orthogonal eigenfunctions as random. A model constructed in this way automatically accounts for variability in eigenfunction estimation and its interplay with both features of the data and the assumed hierarchical structure. The flexibility of this method also makes it well-suited to exploring the imposition of additional constraints on the eigenfunctions, such as mutual orthogonality across levels. We assess the convergence of our model using Grassmannian distances between the spaces spanned by sampled eigenfunctions at each level. After performing validation using a variety of simulated functional data, we compare the results of our model to the prominent existing approaches using 4-hour windows of CGM data for persons with diabetes centered around known mealtimes.