######################### ## Functions to Run JAGS source('R/analysis/lmer.run.R') #source('R/analysis/stan.run.R') ### TODO CHANGE TO FIT WITH NEW LMER DATA AND ### FUNCTION TO GENERATE THE GOOD FORMAT OF DATA model.files.jags.Tf.1 <- c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.species.R", "R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.species.fixed.biomes.R") model.files.jags.Tf.2 <- c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.sepcies.tree.R") fun.generate.init.one.param <- function(i, pars, list.jags,chain){ if(grepl('sigma',names(pars)[i])){ a <- 0.1+(chain-1)*0.1 }else{ if(grepl('param',names(pars)[i])){ a.t <- -0.3+(chain-1)*0.2 a <- rnorm(list.jags[[pars[i]]],mean = a.t) }else{ a <- -0.3+(chain-1)*0.2 } } return(a) } fun.generate.init.all.param <- function(chain, list.jags, jags.model){ list.init <- lapply(seq_len(length(jags.model$pars)), fun.generate.init.one.param, jags.model$pars, list.jags, chain) names(list.init) <- names(jags.model$pars) return(list.init) } fun.generate.init <- function(jags.model, list.jags, chains){ chains.v <- seq_len(chains) inits <- lapply(chains.v,fun.generate.init.all.param, list.jags, jags.model) return(inits) } fun.call.jags.and.save <- function(jags.model, list.jags, path.out, var.sample, iter = 5000, warmup = 500, chains = 3, thin = 50, method = 'rjags'){ require(runjags) start <- Sys.time() inits <- fun.generate.init(jags.model, list.jags[[1]], chains) ### SEND to jags jags.output <- run.jags(data=list.jags[[1]], inits = inits, model = file.path(path.out,"model.file.bug"), monitor = names(jags.model$pars), n.chains = chains , sample = iter, burnin = warmup, adapt = 2000, thin = thin, modules= 'glm on', method = method) end <- Sys.time() print(end -start) print(jags.output) saveRDS(list(output = jags.output, list.sd = list.jags[[2]]), file = file.path(path.out, paste(var.sample, "results.jags.rds", sep = '.'))) return(jags.output) } run.jags.b <- function (model.file, trait, init.TF, data.type='simple', sample.size = NA, ecocode.var = 'wwf', Multi.type = 'a', var.sample = NA, select.set = NA, merge.biomes.TF = FALSE, ...) { require(runjags) source(model.file, local = TRUE) model <- load.model() path.out <- output.dir('lmer', model$name, trait, data.type) print(path.out) dir.create(path.out, recursive = TRUE, showWarnings = FALSE) cat(model$bug, file = file.path(path.out,"model.file.bug") , sep=" ", fill = FALSE, labels = NULL, append = FALSE) cat("run jags for model", model.file, " for trait", trait, "\n") if (init.TF) { jags.list <- fun.load.data.list(trait, data.type, sample.size. = sample.size, ecocode.var. = ecocode.var, Multi.type. = Multi.type, var.sample. = var.sample, select.set. = select.set, sample.vec.TF = FALSE, merge.biomes.TF = merge.biomes.TF) cat("Ok data with Nobs", jags.list[[1]]$N_indiv, "\n") } if (!init.TF) { jags.list <- fun.load.data.list(trait, data.type, sample.size. = sample.size, ecocode.var. = ecocode.var, Multi.type. = Multi.type, var.sample. = var.sample, select.set. = select.set, sample.vec.TF = TRUE, merge.biomes.TF = merge.biomes.TF) cat("Ok data with Nobs", jags.list[[1]]$N_indiv, 'run real', "\n") res <- fun.call.jags.and.save(jags.model = model, list.jags = jags.list, path.out = path.out, var.sample = var.sample, ...) return(res) } gc() } #======================================================================== fun.load.data.list <- function(trait, data.type , base.dir = 'output/processed', sample.size., ecocode.var., Multi.type., var.sample. = 'ecocode', select.set. = NA, select.biome. = NA, sample.vec.TF., merge.biomes.TF = FALSE){ list.df.lmer <- load.data.for.lmer(trait, data.type , base.dir, sample.size. = sample.size., ecocode.var. = ecocode.var., Multi.type. = Multi.type., var.sample. = var.sample., select.biome. = select.biome., sample.vec.TF. = sample.vec.TF., merge.biomes.TF = merge.biomes.TF) stan.list <- fun.turn.in.list.for.jags.stan(list.df.lmer[[1]], data.type) return(list(stan.list, list.df.lmer[[2]])) } fun.turn.in.list.for.jags.stan <- function(df, data.type){ if (data.type == 'Multi'){ stan.list <- fun.turn.in.list.for.jags.stan.Multi(df) }else{ stan.list <- fun.turn.in.list.for.jags.stan.no.Multi(df) } return(stan.list) } fun.turn.in.list.for.jags.stan.no.Multi <- function(df){ stan.list <- list(N_indiv = nrow(df), N_tree = nlevels(unclass(factor(df$tree.id))), N_species = nlevels(unclass(factor(df$species.id))), N_plot = nlevels(unclass(factor(df$plot.id))), N_set = nlevels(unclass(factor(df$set.id))), N_biomes = nlevels(unclass(factor(df$biomes.id))), N_ecocode = nlevels(unclass(factor(df$ecocode.id))), species_id = unclass(factor(df$species.id)), plot_id = unclass(factor(df$plot.id)), set_id = unclass(factor(df$set.id)), tree_id = unclass(factor(df$tree.id)), biomes_id = unclass(factor(df$biomes.id)), ecocode_id = unclass(factor(df$ecocode.id)), logG = df$logG, logD = df$logD, Tf = df$Tf, sumBn = df$sumBn, sumTfBn = df$sumTfBn, sumTnBn = df$sumTnBn, sumTnTfBn_abs = df$sumTnTfBn.abs) return(stan.list) } fun.turn.in.list.for.jags.stan.Multi <- function(df){ stan.list <- as.list(df) stan.list$N_indiv <- nrow(df) stan.list$N_tree <- nlevels(unclass(factor(df$tree.id))) stan.list$N_species <- nlevels(unclass(factor(df$species.id))) stan.list$N_plot <- nlevels(unclass(factor(df$plot.id))) stam.list$N_set <- nlevels(unclass(factor(df$set.id))) stan.list$N_biomes <- nlevels(unclass(factor(df$biomes.id))) stan.list$N_ecocode <- nlevels(unclass(factor(df$ecocode.id))) stan.list$species_id <- unclass(factor(df$species.id)) stan.list$plot_id <- unclass(factor(df$plot.id)) stan.list$set_id <- unclass(factor(df$set.id)) stan.list$tree_id <- unclass(factor(df$tree.id)) stan.list$biomes_id <- unclass(factor(df$biomes.id)) stan.list$ecocode_id <- unclass(factor(df$ecocode.id)) return(stan.list) }