Run model on cluster

library(dplyr); library(brms)
setwd("/usr/users/rarslan/updated_data/")

args = commandArgs()
dataset = args[6]
uptobyear = args[7]

load(paste0(dataset, ".rdata"))

if (dataset == "swed_subset_children") {
    swed_subset_children.1$age = swed_subset_children.1$age/10
}

model_data = get(paste0(dataset, ".1")) %>% tbl_df %>% 
    filter(byear < uptobyear) %>%
    mutate(
        age = age*10,
        age = ifelse(is.na(age), "unknown", 
                                cut(round(age,1), c(0, 25, 50, max(age,na.rm = TRUE)), include.lowest = TRUE) %>% as.character()
        ) %>% factor() %>% relevel("(25,50]")
        ) %>%
    select(children, birth_cohort, male, maternalage.factor, paternalage.mean, paternalage, paternal_loss, maternal_loss, older_siblings, nr.siblings, last_born, idParents, age) %>% 
    na.omit()


model_formula = children ~ paternalage + age + birth_cohort + male + maternalage.factor + paternalage.mean + paternal_loss + maternal_loss + older_siblings + nr.siblings + last_born + (1 | idParents)


if (dataset == "swed_subset_children") {
    model_prior = c(set_prior("normal(0,5)", class = "b"),
                                    set_prior("student_t(3, 0, 5)", class = "sd"))
    model_family = poisson()
} else {
    model_prior = c(set_prior("normal(0,5)", class = "b"),
                                    set_prior("normal(0,5)", class = "b", nlpar = "hu"),
                                    set_prior("student_t(3, 0, 5)", class = "sd"))
    model_formula_hu = update(model_formula,  hu ~ . )
    model_formula = bf(model_formula, model_formula_hu)
    model_family = hurdle_poisson()
}

model = brm(    model_formula,
                         prior = model_prior, 
                         family = model_family, data = model_data, 
                         chains = 6, iter = 1000, warmup = 500, cores = 6, ranef = FALSE)

summary(model)

saveRDS(model,file = paste0("coefs/", dataset, "/s1_control_age.rds"))