# SODF1 example BMIPR Chapter 23 rm (list = ls (all = TRUE)) # CAUTION: cleans environment library (rstan) # rstan (Version 2.19.3, GitRev: 2e1f913d3ca3) # For execution on a local, multicore CPU with excess RAM we recommend calling options(mc.cores = parallel::detectCores()) # To avoid recompilation of unchanged Stan programs, we recommend calling rstan_options(auto_write = TRUE) # For improved execution time, we recommend calling Sys.setenv(LOCAL_CPPFLAGS = '-march=corei7 -mtune=corei7') # although this causes Stan to throw an error on a few processors. library (tidyverse) # Read data # 7 time points: 1,2,5,10,15,20,30 minute # (6batch)x(12unit/batch)x(7disso/unit) = 504disso SODF1 <- tibble(read.csv("SODF1.csv")) SODF1$batch <- as.factor(SODF1$batch) SODF1$unit <- as.factor(SODF1$unit) # Plot data ggplot(SODF1, aes(y=disso,x=minute,shape=unit,color=unit)) + scale_shape_manual(values=rep(1:12,6)) + scale_color_manual(values=rep(1:6,each=12)) + geom_point() + geom_line() + facet_wrap(~batch) # Organize data for stan dataList = list( K=3, N=length(SODF1$disso), B=max(as.numeric(SODF1$batch)), U=max(as.numeric(SODF1$unit)), b=as.numeric(SODF1$batch), u=as.numeric(SODF1$unit), y=SODF1$disso, x=SODF1$minute ) # compile model comp <- stan_model('SODF1d.stan') StanFit <- sampling(comp, data=dataList, iter=8000, chains=4, pars=c("sigma","sigma_b","sigma_u","mu"))