I had a similar problem a while ago and I too lamented the lack of something as easy as a predict() type solution. If you’re not yet familiar with it, generating posterior predicted values using McElreath’s rethinking R package is just as easy. The rethinking package is a front-end for stan and I can’t say enough positive things about the package or McElreath’s book. If you pursue this route, see the link() function in chapter 4.
Until there’s a similar solution available directly in rstan, I think you’ve got to do this manually. But rather than using expose_stan_functions(), I think it’s easier and more transparent to do this entirely outside of stan and in R. Here’s an example using the rats example that comes with rstan. In that example, assume you’ve fit the very same model I just linked to in stan and that you have a resulting stanfit object called “rats1”. Assume also that your interest is in predicting average rat weight for new data–days 21, 22, and 23.
# generate posterior samples of the stan model coefficients
param.sample <- as.data.frame(rats1)
# assume data2 is your new df. This df has just 3 time periods one day apart in which the rats were weighed
data2 <- data.frame(x=c(21, 22, 23))
data2$xbar <- mean(data2$x)
# posterior predicted values for rat weight using the new data frame (data2):
result.sample<- t(apply(param.sample, 1, function(x) x["mu_alpha"] + x["mu_beta"]*(data2$xbar - data2$x) ))
result.summary <- apply(result.sample, 2, mean)
result.summary.std <- apply(result.sample, 2, sd)
I’m guessing your interest would be in result.summary and in the standard errors around it.