library(rstan) library(doParallel) library(foreach) rstan::rstan_options(auto_write = TRUE) Sys.setenv(LOCAL_CPPFLAGS = '-march=native') options(mc.cores = parallel::detectCores()) library(dplyr) library(tibble) # GENERATE DATA # generate data set.seed(123) ldat1 <- list() ldat2 <- list() ldat3 <- list() nsim <- 2000 for(i in 1:nsim){ N <- 250 x <- rep(0:1, each = N/2) g0 <- 0.5 g1 <- 0.8 p <- (1-x)*g0 + x*g1 y <- rbinom(N, 1, p) idx_ymis <- as.logical(rbinom(N, 1, 0.2)) N_mis <- sum(idx_ymis) d <- tibble(y = y, idx_ymis = idx_ymis, x = x) # setup for no missing data analysis ld1 <- list(N_obs = N - 0, N_mis = 0, y_obs = d$y, X_obs = d$x, X_mis = numeric(0)) # setup for complete case analysis ld2 <- list(N_obs = N - N_mis, N_mis = 0, y_obs = d$y[!d$idx_ymis], X_obs = d$x[!d$idx_ymis], X_mis = numeric(0)) # set up for marginalized missing ld3 <- list(N_obs = N - N_mis, N_mis = N_mis, y_obs = d$y[!d$idx_ymis], X_obs = d$x[!d$idx_ymis], X_mis = d$x[d$idx_ymis]) ldat1[[i]] <- ld1 ldat2[[i]] <- ld2 ldat3[[i]] <- ld3 } m1 <- rstan::stan_model("logistic_003.stan") starttime <- Sys.time() message("Start simulation at ", starttime) # initiate clusters debug = F cl <- NA if(!debug){ cl <- makeCluster(4, outfile="") # cl <- makeCluster(3, outfile="") registerDoParallel(cl) } else { registerDoSEQ() } # predictive probability set.seed(66) res1 <- foreach(i = 1:nsim, .errorhandling = 'pass', .packages=c("rstan") #.options.snow=opts, ) %dopar%{ options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) message(paste0("Iteration ", i)) s1 <- rstan::sampling(object = m1, data = ldat1[[i]], chains = 4, thin = 1, iter = 2000, refresh = 2000, control = list(max_treedepth = 15) ) smry <- summary(s1, probs = c(0.025, 0.975), pars = c("b0", "b1", "p0", "p1"))$summary m <- as.data.frame(cbind(desc = 1, idx = i, smry)) m$par <- rownames(m) return(m) } res2 <- foreach(i = 1:nsim, .errorhandling = 'pass', .packages=c("rstan") #.options.snow=opts, ) %dopar%{ options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) message(paste0("Iteration ", i)) s1 <- rstan::sampling(object = m1, data = ldat2[[i]], chains = 4, thin = 1, iter = 2000, refresh = 2000, control = list(max_treedepth = 15) ) smry <- summary(s1, probs = c(0.025, 0.975), pars = c("b0", "b1", "p0", "p1"))$summary m <- as.data.frame(cbind(desc = 2, idx = i, smry)) m$par <- rownames(m) return(m) } res3 <- foreach(i = 1:nsim, .errorhandling = 'pass', .packages=c("rstan") #.options.snow=opts, ) %dopar%{ options(mc.cores = parallel::detectCores()) rstan_options(auto_write = TRUE) message(paste0("Iteration ", i)) s1 <- rstan::sampling(object = m1, data = ldat3[[i]], chains = 4, thin = 1, iter = 2000, refresh = 2000, control = list(max_treedepth = 15) ) smry <- summary(s1, probs = c(0.025, 0.975), pars = c("b0", "b1", "p0", "p1"))$summary m <- as.data.frame(cbind(desc = 3, idx = i, smry)) m$par <- rownames(m) return(m) } endtime <- Sys.time() difftime(endtime, starttime, units = "hours") # turn of parallel cluster stopCluster(cl) mres1 <- do.call(rbind, res1) mres2 <- do.call(rbind, res2) mres3 <- do.call(rbind, res3) mres <- rbind(mres1, mres2, mres3) lsav <- list(mres = mres, starttime = starttime, endtime = endtime) saveRDS(lsav, "mres.RDS")