We model sparse functional data from multiple subjects with a mixed-effects regression spline. In this model, the expected values for any subject (conditioned on the random effects) can be written as the sum of a population curve and a subject-specific deviate from this population curve. The population curve and the subject-specific deviates are both modeled as free-knot b-splines with k and k′ knots located at tk and tk′, respectively. To identify the number and location of the “free” knots, we sample from the posterior p (k, tk, k′, tk′|y) using reversible jump MCMC methods. Sampling from this posterior distribution is complicated, however, by the flexibility we allow for the model’s covariance structure. No restrictions (other than positive definiteness) are placed on the covariance parameters ψ and σ2 and, as a result, no analytical form for the likelihood p (y|k, tk, k′, tk′) exists. In this paper, we consider two approximations to p(y|k, tk, k′, tk′) and then sample from the corresponding approximations to p(k, tk, k′, tk′|y). We also sample from p(k, tk, k′, tk′, ψ, σ2|y) which has a likelihood that is available in closed form. While sampling from this larger posterior is less efficient, the resulting marginal distribution of knots is exact and allows us to evaluate the accuracy of each approximation. We then consider a real data set and explore the difference between p(k, tk, k′, tk′, ψ, σ2|y) and the more accurate approximation to p(k, tk, k′, tk′|y).