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

vs.data <- list(I=I, T=T, Tmax=Tmax,
                Obs=Obs, 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)

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

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

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

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

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

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

vs.data5 <- list(I=I, T=T, Tmax=Tmax,
                 Obs=Obs5, Time=Time, Dose=Dose5,
                 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)

vs.fit5 <- stan(file="varSlopes2.stan", data=vs.data5, iter=2000, chains=5)

print(vs.fit5,pars=c("var_innov","treat_eff","slope_mu",
                    "slope_std","slope_r2","lp__"))