  Functional mixed models are widely useful for regression analysis with
dependent functional data, including longitudinal functional data with scalar
predictors. However, existing algorithms for Bayesian inference with these
models only provide either scalable computing or accurate approximations to the
posterior distribution, but not both. We introduce a new MCMC sampling strategy
for highly efficient and fully Bayesian regression with longitudinal functional
data. Using a novel blocking structure paired with an orthogonalized basis
reparametrization, our algorithm jointly samples the fixed effects regression
functions together with all subject- and replicate-specific random effects
functions. Crucially, the joint sampler optimizes sampling efficiency for these
key parameters while preserving computational scalability. Perhaps
surprisingly, our new MCMC sampling algorithm even surpasses state-of-the-art
algorithms for frequentist estimation and variational Bayes approximations for
functional mixed models -- while also providing accurate posterior uncertainty
quantification -- and is orders of magnitude faster than existing Gibbs
samplers. Simulation studies show improved point estimation and interval
coverage in nearly all simulation settings over competing approaches. We apply
our method to a large physical activity dataset to study how various
demographic and health factors associate with intraday activity.
