-
Notifications
You must be signed in to change notification settings - Fork 0
/
run_models_local.R
54 lines (38 loc) · 1.8 KB
/
run_models_local.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
## Packages needed
library(rstan)
library(matrixStats)
library(doParallel)
load('Rout/model_run_setup.RData')
for(i in 1:length(model_settings)){
writeLines('Doing the following job:')
print(model_settings[i, ])
options(mc.cores = model_settings$Nchain[i])
stopifnot(model_settings$Nchain[i]>getDoParWorkers()) # check worker number assigned
mod = stan_model(file = as.character(model_settings$mod[i])) # compile
stan_input_job = stan_inputs[[model_settings$dataset[i]]]
analysis_data_stan = stan_input_job$analysis_data_stan
analysis_data_stan$trt_mat = stan_input_job$Trt_matrix
analysis_data_stan$K_trt = ncol(analysis_data_stan$trt_mat)
x_intercept = stan_input_job$cov_matrices$X_int[[model_settings$cov_matrices[i]]]
if(ncol(x_intercept)==0) x_intercept = array(0, dim=c(nrow(x_intercept),1))
analysis_data_stan$x_intercept = x_intercept
analysis_data_stan$K_cov_intercept= ncol(x_intercept)
x_slope = stan_input_job$cov_matrices$X_slope[[model_settings$cov_matrices[i]]]
if(ncol(x_slope)==0) x_slope = array(0, dim=c(nrow(x_slope),1))
analysis_data_stan$x_slope = x_slope
analysis_data_stan$K_cov_slope=ncol(x_slope)
# sample posterior
out = sampling(mod,
data=c(analysis_data_stan,
all_priors[[model_settings$prior[i]]]),
iter=model_settings$Niter[i],
chain=model_settings$Nchain[i],
thin=model_settings$Nthin[i],
warmup=model_settings$Nwarmup[i],
save_warmup = FALSE,
seed=i,
pars=c('L_Omega','theta_rand_id'), # we don't save these as it takes up vast memory!
include=FALSE)
save(out, file = paste0('Rout/model_fits_',i,'.RData'))# save output
writeLines('Finished job')
}