Commit 386ee232 authored by kunstler's avatar kunstler
Browse files

progress on jags

parent 5e61d7ad
......@@ -2,49 +2,47 @@
## Functions to Run JAGS
source('R/analysis/lmer.run.R')
source('R/analysis/stan.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.fixed.biomes.R")
model.files.jags.Tf.1b <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.fixed.biomes.b.R")
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.fixed.biomes.no.plot.R")
model.files.jags.Tf.2b <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.fixed.biomes.no.plot.b.R")
model.files.jags.Tf.0 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.fixed.R")
model.files.jags.Tf.00 <-
c("R/analysis/model.jags/model.LOGLIN.fixed.R")
model.files.jags.Tf.3 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.biomes.R")
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.sepcies.tree.R")
model.files.jags.Tf.4 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.ecocode.R")
fun.generate.init0 <- function(x, N){
if(!grepl('vec',x)){
a <- c(0.5)
}else{
a <- rep(0.5, N)
}
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, N, chains){
init0 <- lapply(jags.model$pars,fun.generate.init0, N)
names(init0) <- jags.model$pars
inits <- lapply(1:chains,function(x,init0) lapply(init0,function(x,a)x+a, a = (x-1)*0.1), init0)
return(inits)
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,
......@@ -55,21 +53,15 @@ fun.call.jags.and.save <- function(jags.model,
warmup = 500,
chains = 3,
thin = 50,
biomes.TF = TRUE){
parallel = FALSE){
start <- Sys.time()
model.file <-
cat(jags.model$bug, file = file.path(path.out,"model.file.bug")
, sep=" ", fill = FALSE, labels = NULL, append = FALSE)
if(biomes.TF){
inits <- fun.generate.init(jags.model, list.jags$N_biomes, chains)
}
if(!biomes.TF){
inits <- fun.generate.init(jags.model, 1, chains)
}
inits <- fun.generate.init(jags.model, list.jags[[1]], chains)
### SEND to jags
jags.output <- jags(data=list.jags, inits = inits,
if(!parallel){
jags.output <- jags(data=list.jags[[1]], inits = inits,
model.file = file.path(path.out,"model.file.bug"),
parameters.to.save = jags.model$pars,
parameters.to.save = names(jags.model$pars),
n.chains = chains ,
DIC = FALSE,
n.iter = iter,
......@@ -78,45 +70,66 @@ inits <- fun.generate.init(jags.model, 1, chains)
end <- Sys.time()
print(end -start)
print(jags.output)
saveRDS(jags.output, file = file.path(path.out, paste(var.sample,
"results.jags.rds", sep = '.')))
}else{
## jags.output <- jags.parallel(data=list.jags[[1]], inits = inits,
## model.file = file.path(path.out,"model.file.bug"),
## parameters.to.save = names(jags.model$pars),
## n.chains = chains ,
## DIC = FALSE,
## n.iter = iter,
## n.burnin = warmup,
## n.thin = thin)
## TODO CHANGE TO run.jags
}
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 <- function (model.file, trait, init.TF,
min.obs = 10, sample.size = NA,
type.filling='species', cat.TF,
fname = 'data.all.no.log.rds',
var.sample = 'ecocode.id',
data.type='simple',
sample.size = NA,
ecocode.var = 'wwf',
Multi.type = 'a',
var.sample = NA,
select.set = NA,
merge.biomes.TF = FALSE,
...) {
require(R2jags)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
if(fname == 'data.all.no.log.rds') dir <- 'all.no.log'
if(fname == 'data.all.log.rds') dir <- 'all.log'
path.out <- output.dir(type = 'jags', model = model$name,
trait = trait, set = dir,
type.filling = type.filling)
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.stan(trait,fname = fname, cat.TF = cat.TF,
sample.size = sample.size,
var.sample = var.sample,
sample.vec.TF = FALSE)
cat("Ok data with Nobs", jags.list$N_indiv,
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.stan(trait,fname = fname, cat.TF = cat.TF,
sample.size = sample.size,
var.sample = var.sample,
sample.vec.TF = TRUE)
cat("Ok data with Nobs", jags.list$N_indiv, 'run real',
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,
......@@ -129,3 +142,81 @@ run.jags <- function (model.file, trait, init.TF,
}
#========================================================================
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)
}
......@@ -514,7 +514,7 @@ tf <- paste(trait, "focal", sep = ".")
BATOT <- "BATOT"
#= DATA LIST FOR LMER
if(data.type == 'simple' | data.type == 'all.census') {
print('no cat')
print('no Multi')
lmer.data <- get.variables(data.tree,
BATOT,
CWM.tn,
......
#!/bin/bash
export LD_LIBRARY_PATH=/usr/lib64/R/library
mkdir -p Rscript_temp
pwd
for trait in "Wood.density" "SLA" "Leaf.N" "Max.height" "Seed.mass"; do
Rscript -e "source('R/analysis/lmer.run.R'); load.and.save.data.for.lmer('$trait')"
for trait in "'SLA'" "'Wood.density'" "'Max.height'"; do
Rscript -e "source('R/analysis/lmer.run.R'); load.and.save.data.for.lmer($trait);print('done')"&
Rscript -e "source('R/analysis/lmer.run.R'); load.and.save.data.for.lmer($trait, data.type = 'all.census');print('done')"&
done
Rscript -e "source('R/analysis/lmer.run.R'); load.and.save.data.for.lmer('SLA', data.type = 'Multi')"
Rscript -e "source('R/analysis/lmer.run.R'); load.and.save.data.for.lmer('SLA', data.type = 'Multi');print('done')"&
......@@ -13,17 +13,17 @@ format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
traits = c("SLA", "Wood.density", "Max.height")
)
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.all.census.set.rds',
models = c(model.files.lmer.Tf.2),
traits = c("SLA", "Wood.density", "Max.height"),
data.type = 'all.census'
)
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.simple.set.Multi.rds',
models = c(model.files.lmer.Tf.Multi),
traits = c("Multi"),
data.type = 'Multi'
)
## format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
## list.file.name = 'list.lmer.out.all.NA.all.census.set.rds',
## models = c(model.files.lmer.Tf.2),
## traits = c("SLA", "Wood.density", "Max.height"),
## data.type = 'all.census'
## )
## format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
## list.file.name = 'list.lmer.out.all.NA.simple.set.Multi.rds',
## models = c(model.files.lmer.Tf.Multi),
## traits = c("Multi"),
## data.type = 'Multi'
## )
### test jags
source('R/analysis/jags.run.R')
library(R2jags)
# look at data
data <- readRDS('output/processed/data.Wood.density.no.log.rds')
# number of obs per plot
count.plot <- table(data$plot)
length(count.plot)
hist(count.plot, breaks = 50, prob = TRUE)
sort( table(count.plot)/sum(count.plot))
# if use only one obs per indiv
biomes.count <- table(data$biomes)
biomes.count.unique <- table(data$biomes[!duplicated(data$tree.id)])
cbind(biomes.count, biomes.count.unique)
## biomes.count biomes.count.unique
## Boreal forest 231218 218615
## Subtropical desert 1998 1998
## Temperate forest 1368965 859256
## Temperate grassland desert 43895 43812
## Temperate rain forest 36066 29158
## Tropical forest savanna 98485 42374
## Tropical rain forest 8482 7298
## Tundra 5251 5195
## Woodland shrubland 616223 524105
(biomes.count - biomes.count.unique)/biomes.count
## test different models
## Boreal forest Subtropical desert
## 0.054507002 0.000000000
## Temperate forest Temperate grassland desert
## 0.372331652 0.001890876
## Temperate rain forest Tropical forest savanna
## 0.191537736 0.569741585
## Tropical rain forest Tundra
## 0.139589719 0.010664635
## Woodland shrubland
## 0.149488091
run.jags(model.files.jags.Tf.1[1], trait = 'Wood.density',
init.TF = TRUE, var.sample = 'wwf',
sample.size = 5000, iter = 2000,
warmup = 200, chains = 4,
thin = 2)
# what I need to test:
# - model without tree id
# - model with no biomes effect
# - model without plot effect
# - model without T effect
# - there was a problem with T() need to check if this changes anything. with no at all or less in stan as well.
res1 <- run.jags(model.files.jags.Tf.1[1], trait = 'Wood.density',
init.TF = FALSE, var.sample = 'wwf',
sample.size = 1000, iter = 2000,
warmup = 200, chains = 3,
thin = 2, parallel = FALSE)
# - memory issue how to know what went wrong
# check some results
res <- readRDS("/home/kunstler/Documents/currentwork/Dropbox/westobyLAB/workshop/data.processing2/trait.competition.workshop/output/jags/all.no.log/Wood.density/species/LOGLIN.ER.AD.Tf.fixed.biomes/ecocode.id.results.jags.rds")
traceplot(res)
cor(res$BUGSoutput$sims.matrix)
## test different models
run.multiple.model.for.set.one.trait(model.files.jags.Tf.1, run.jags,
trait = 'Wood.density',
type.filling='species',
sample.size = 100000,
var.sample = 'ecocode.id',
iter = 200, warmup = 20, chains = 3,
init.TF = TRUE)
res2 <- run.jags(model.files.jags.Tf.1[2], trait = 'Wood.density',
init.TF = FALSE, var.sample = 'wwf',
sample.size = 10000, iter = 2000,
warmup = 200, chains = 2,
thin = 2, merge.biomes.TF = TRUE)
res0 <- run.jags(model.files.jags.Tf.0, trait = 'Wood.density',
type.filling = 'species',
init.TF = FALSE, cat.TF = FALSE,
sample.size = 100000, iter = 2000,
warmup = 200, chains = 2,
thin = 2,
biomes.TF = FALSE)
res0
run.multiple.model.for.set.one.trait(c(model.files.jags.Tf.1, model.files.jags.Tf.1b), run.jags,
run.multiple.model.f
or.set.one.trait(c(model.files.jags.Tf.1, model.files.jags.Tf.1b), run.jags,
trait = 'Wood.density',
type.filling='species',
sample.size = 100000,
var.sample = 'ecocode.id',
iter = 200, warmup = 20, chains = 3, thin = 2,
init.TF = FALSE)
init.TF = TRUE)
res1 <- run.jags(model.files.jags.Tf.1, trait = 'Wood.density',
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment