stanmodel <- " data { int<lower=0> N; real x[N]; real Y[N]; } parameters { real alpha; real beta; real<lower=.5,upper= 1> lambda; // original gamma in the JAGS example real<lower=0> tau; } transformed parameters { real sigma; real m[N]; for (i in 1:N) m[i] = alpha - beta * pow(lambda, x[i]); sigma = 1 / sqrt(tau); } model { // priors alpha ~ normal(0.0, 1000); beta ~ normal(0.0, 1000); lambda ~ uniform(.5, 1); tau ~ gamma(.0001, .0001); // likelihood Y ~ normal(m, sigma); } generated quantities{ real Y_mean[N]; real Y_pred[N]; for(i in 1:N){ // Posterior parameter distribution of the mean Y_mean[i] = alpha - beta * pow(lambda, x[i]); // Posterior predictive distribution Y_pred[i] = normal_rng(Y_mean[i], sigma); } } " library(rstan) rstan_options(auto_write = TRUE) options(mc.cores = parallel::detectCores()) fit <- stan(model_code = stanmodel, model_name = "GrowthCurve", data = dat) Y_mean <- extract(fit, "Y_mean") Y_mean_cred <- apply(Y_mean$Y_mean, 2, quantile, c(0.05, 0.95)) Y_mean_mean <- apply(Y_mean$Y_mean, 2, mean) Y_pred <- extract(fit, "Y_pred") Y_pred_cred <- apply(Y_pred$Y_pred, 2, quantile, c(0.05, 0.95)) Y_pred_mean <- apply(Y_pred$Y_pred, 2, mean) plot(dat$Y ~ dat$x, xlab="x", ylab="Y", ylim=c(1.6, 2.8), main="Non-linear Growth Curve") lines(dat$x, Y_mean_mean) points(dat$x, Y_pred_mean, pch=19) lines(dat$x, Y_mean_cred[1,], col=4) lines(dat$x, Y_mean_cred[2,], col=4) lines(dat$x, Y_pred_cred[1,], col=2) lines(dat$x, Y_pred_cred[2,], col=2) legend(x="bottomright", bty="n", lwd=2, lty=c(NA, NA, 1, 1,1), legend=c("observation", "prediction", "mean prediction", "90% mean cred. interval", "90% pred. cred. interval"), col=c(1,1,1,4,2), pch=c(1, 19, NA, NA, NA)) fit # Inference for Stan model: GrowthCurve. # 4 chains, each with iter=2000; warmup=1000; thin=1; # post-warmup draws per chain=1000, total post-warmup draws=4000. # # mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat # alpha 2.65 0.00 0.07 2.53 2.60 2.65 2.70 2.82 1067 1 # beta 0.97 0.00 0.07 0.83 0.93 0.97 1.02 1.12 1425 1 # lambda 0.86 0.00 0.03 0.79 0.85 0.87 0.89 0.92 935 1 # tau 108.99 0.81 31.55 56.97 86.13 106.61 127.93 178.90 1508 1 # sigma 0.10 0.00 0.02 0.07 0.09 0.10 0.11 0.13 1472 1 # . # Rows omitted # . # lp__ 47.00 0.05 1.57 43.16 46.17 47.35 48.17 48.96 1116 1 # # Samples were drawn using NUTS(diag_e) at Tue Oct 27 07:33:40 2015. # For each parameter, n_eff is a crude measure of effective sample size, # and Rhat is the potential scale reduction factor on split chains (at # convergence, Rhat=1).