library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores=5)

vsA.data <- list(I=I, T=T, Tmax=Tmax,
                Obs=ObsA, Time=Time, Dose=Dose,
                slope_mu_mu=slope_mu_mu,slope_mu_std=slope_mu_std,
                obs_rel=obs_rel, obs_mean_t1=obs_mean_t1,
                obs_std_t1=obs_std_t1)

vsA.fit <- stan(file="varSlope3.stan", data=vsA.data, iter=2000, chains=5)

print(vsA.fit,pars=c("var_innov","treat_eff","slope_mu",
                    "slope_std","slope_r2","lp__"))
plot(vsA.fit,pars=c("var_innov","treat_eff","slope_mu",
                    "slope_std","slope_r2"))
pairs(vsA.fit,pars=c("var_innov","treat_eff","slope_mu",
                    "slope_std","slope_r2","lp__"))

summary(monitor(extract(vsA.fit,pars="slope",permuted=FALSE),print=FALSE)[,"Rhat"])
summary(monitor(extract(vsA.fit,pars="slope_res",permuted=FALSE),print=FALSE)[,"Rhat"])
summary(monitor(extract(vsA.fit,pars="T0_val",permuted=FALSE),print=FALSE)[,"Rhat"])
summary(monitor(extract(vsA.fit,pars="innov",permuted=FALSE),print=FALSE)[,"Rhat"])

innovA_est <- extract(vsA.fit,pars="innov")$innov
innovA_mean <- apply(innovA_est,c(2,3),mean)
pairs(t(innovA_mean))
cor(t(innovA_mean))
summary(as.vector(cor(t(innovA_mean))))
hist(as.vector(cor(t(innovA_mean))))

slope_est <- extract(vsA.fit,pars="slope")$slope

summary(as.vector(cor(t(innovA_mean),t(Dose))))