You could take one of the splines packages to create basis functions, add them to a design matrix, and fit with MixedModels, with random effects acting as a smoothing penalty, or even using MCMC. Look at Tutorial in biostatistics: spline smoothing with linear mixed models..