Growth_curve_2.R

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).