Power simulation

sample_people = c(25, 50, 70, 100, 150, 250, 500)
nr_of_days = c(2, 10, 30, 60, 90)
dayspans = list("1:38" = 1:38, "17-19,4-6" = c(4:6, 17:19))
fertility_effects = c(0, 0.05, 0.1, 0.2, 0.3, 0.5, 0.8)
trait_effects = c(0, 0.3)
miss_windows = c(0, 1, 2, 3, 4, 6, 8)
nrep <- 1:1000
predictors = c("prc_stirn_b_m", "fertile_broad_m", "fertile_narrow_m")

# res0 <- fit_model(simulate_fertility_effect())

# time = system.time({ for (i in nrep) { fit_model(simulate_fertility_effect()) }})
# length(sample_people) *
# length(nr_of_days) *
# length(dayspans)
# 
# length(fertility_effects) * # nodes
# length(trait_effects) *
# length(miss_windows) *
# length(predictors) *
#   time[3]  / 60 / 60 # hours per node

## prefer to use i,j,k,... for indices
resdf = foreach(i = seq_along(sample_people), .combine = "rbind") %:%
  foreach(j = seq_along(nr_of_days), .combine = "rbind") %:%
    foreach(k = seq_along(dayspans), .combine = "rbind") %:%
    foreach(l = seq_along(fertility_effects), .combine = "rbind") %dopar% {
      source("0_simulation_functions.R")
      
    foreach(m = seq_along(trait_effects), .combine = "rbind") %do% {
        foreach(n = seq_along(miss_windows), .combine = "rbind", .errorhandling = "remove") %do% {
         foreach(repetition = nrep, .combine = "rbind", .maxcombine = 1000, .errorhandling = "remove") %do% {
           cat(i,j,k,l,m,n, "\n") 
           
           simulated_data = simulate_fertility_effect(
             nr_of_people = sample_people[i],
             nr_days = nr_of_days[j],
             dayspan = dayspans[[k]],
             miss_window = miss_windows[n],
             effects = list(
               trait = trait_effects[m],
               fertility_effect = fertility_effects[l],
               noise = 0.5
             )
           )
           foreach(o = seq_along(predictors), .combine = "rbind") %do% {
            
            ## if you want to see progress in the output file 
            result = data.frame(
                nr_of_people = sample_people[i],
                nr_days = nr_of_days[j],
                dayspan = names(dayspans)[k],
                miss_window = miss_windows[n],
                trait_effect = trait_effects[m],
                fertility_effect = fertility_effects[l],
                predictor = predictors[o],
                estimate = NA_real_, 
                std.error = NA_real_, 
                statistic = NA_real_, 
                conf.low = NA_real_, 
                conf.high = NA_real_, 
                p.value = NA_real_, 
                p.value_KR = NA_real_, 
                usable_days = NA_real_
            )
            try({
                result[, c("estimate", "std.error", "statistic", "conf.low", "conf.high", 
                     "p.value", "p.value_KR", "usable_days")] = fit_model(
                        simulated_data
                                  , predictor = predictors[o])
            }, silent = T)
            result
          }
        }
        }
  }
}
saveRDS(resdf,file = "cycle_sims.rds")


closeCluster(cl)
mpi.quit()