Commit 6c62313a authored by kunstler's avatar kunstler
Browse files

fixe error in resampling with only one census

parent d71616db
......@@ -22,7 +22,7 @@ if(grepl('sigma',names(pars)[i])){
}else{
if(grepl('param',names(pars)[i])){
a.t <- -0.3+(chain-1)*0.2
a <- rnorm(list.jags[[pars[i]]],mean = a.t)
a <- rnorm(list.jags[[pars[i]]],mean = a.t, sd = 0.1)
}else{
a <- -0.3+(chain-1)*0.2
}
......@@ -55,11 +55,12 @@ fun.call.jags.and.save <- function(jags.model,
chains = 3,
thin = 50,
method = 'rjags'){
require(runjags)
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]],
print(getwd())
jags.output <- run.jags(data=list.jags[[1]],
inits = inits,
model = file.path(path.out,"model.file.bug"),
monitor = names(jags.model$pars),
......@@ -80,7 +81,7 @@ inits <- fun.generate.init(jags.model, list.jags[[1]], chains)
}
run.jags.b <- function (model.file, trait, init.TF,
data.type='simple',
data.type='simple',
sample.size = NA,
ecocode.var = 'wwf',
Multi.type = 'a',
......@@ -91,8 +92,8 @@ run.jags.b <- function (model.file, trait, init.TF,
require(runjags)
source(model.file, local = TRUE)
model <- load.model()
path.out <- output.dir('lmer', model$name, trait, data.type)
print('OK')
path.out <- output.dir('jags', model$name, trait, data.type)
print(path.out)
dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
......@@ -151,7 +152,6 @@ fun.load.data.list <- function(trait, data.type ,
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]]))
......
......@@ -158,35 +158,36 @@ load.data.for.lmer <- function(trait, data.type,
sp.name))
if(merge.biomes.TF) df <- merge.biomes(df)
if(data.type != 'all.census') {
df <- select.one.census.per.plot(df)
}
if(!is.na(select.biome.)) {
df <- dplyr::filter(df, biomes == select.biome.)
}
if(!is.na(select.set.)) {
df <- dplyr::filter(df, set == select.set.)
}
if(!is.na(sample.size.)){
if(!sample.vec.TF.){
if(!is.na(sample.size.) ){
if(sample.size. < length(unique(df$plot))){
if(!sample.vec.TF.){
if(sample.size. > length(unique(df$plot))){
sample.size. <- length(unique(df$plot))}
plot.selected <- resample.plot.by.var(df, var.sample., sample.size.)
sample.vec <- (1:nrow(df))[df$plot %in% plot.selected]
df <- df[sample.vec, ]
saveRDS(sample.vec,file = file.path(base.dir,paste('sample.vec',
plot.selected <- resample.plot.by.var(df, var.sample., sample.size.)
sample.vec <- (1:nrow(df))[df$plot %in% plot.selected]
df <- df[sample.vec, ]
saveRDS(sample.vec,file = file.path(base.dir,paste('sample.vec',
trait, data.type,'rds',
sep = '.')))
print(paste('sub-sampled by',var.sample.))
}else{
sample.vec <- readRDS(file = file.path(base.dir,paste('sample.vec', trait,
print(paste('sub-sampled by',var.sample.))
}else{
sample.vec <- readRDS(file = file.path(base.dir,paste('sample.vec', trait,
data.type,'rds',
sep = '.')))
df <- df[sample.vec, ]
print(paste('sub-sampled by',var.sample., 'by loading sample.vec'))
df <- df[sample.vec, ]
print(paste('sub-sampled by',var.sample., 'by loading sample.vec'))
}
}
}else{
}
}
if(data.type != 'all.census') {
df <- select.one.census.per.plot(df)
}
list.sd <- get.sd.lmer(df, trait, min.obs = 10)
res <- format.data.for.lmer(df, trait,
data.type = data.type,
......
......@@ -25,12 +25,10 @@ samplesize=$1
for trait in "'SLA'" "'Wood.density'" "'Max.height'"; do
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/jags.run.R'); run.multiple.model.for.set.one.trait(c(model.files.jags.Tf.1[1]), run.jags.b,trait = $trait,data.type='simple', sample.size = NA, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, method = 'parallel');print('done')\"" > Rscript_temp/speciesjags${trait}.sh
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/jags.run.R'); run.multiple.model.for.set.one.trait(model.files.jags.Tf.1[1], run.jags.b, trait = $trait,data.type='simple', sample.size = NA, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, method = 'parallel');print('done')\"" > Rscript_temp/speciesjags${trait}.sh
qsub Rscript_temp/speciesjags${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=4,mem=16gb -N "jags${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/jags.run.R'); run.multiple.model.for.set.one.trait(c(model.files.jags.Tf.1[2]), run.jags.b,trait = $trait,data.type='simple', sample.size = NA, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, merge.biomes.TF = TRUE, method = 'parallel');print('done')\"" > Rscript_temp/speciesjags2${trait}.sh
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/jags.run.R'); run.multiple.model.for.set.one.trait(model.files.jags.Tf.1[2], run.jags.b, trait = $trait,data.type='simple', sample.size = NA, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, merge.biomes.TF = TRUE, method = 'parallel');print('done')\"" > Rscript_temp/speciesjags2${trait}.sh
qsub Rscript_temp/speciesjags2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=4,mem=16gb -N "jags2${trait}" -q opt32G -j oe
done
......@@ -11,6 +11,14 @@ library(runjags)
sample.size = 5000, iter = 200,
warmup = 20, chains = 4,
thin = 2)
run.multiple.model.for.set.one.trait(model.files.jags.Tf.1[1], run.jags.b, trait = 'Wood.density',data.type='simple', sample.size = 1000, var.sample = 'wwf',iter = 200, warmup = 100, chains = 3, thin = 1, init.TF = FALSE, method = 'parallel')
Rscript -e "source('R/analysis/jags.run.R'); library(runjags); sessionInfo(); run.multiple.model.for.set.one.trait(c(model.files.jags.Tf.1[1]), run.jags.b, trait = 'Wood.density',data.type='simple', sample.size = 1000, var.sample = 'wwf',iter = 200, warmup = 100, chains = 3, thin = 1, init.TF = FALSE, method = 'parallel');print('done')"
echo "source('R/analysis/jags.run.R'); library(runjags); sessionInfo(); run.multiple.model.for.set.one.trait(c(model.files.jags.Tf.1[1]), run.jags.b, trait = 'Wood.density',data.type='simple', sample.size = 1000, var.sample = 'wwf',iter = 200, warmup = 100, chains = 3, thin = 1, init.TF = FALSE, method = 'parallel');print('done')" > Rscript_temp/speciesjagsWood.density.R
R CMD BATCH Rscript_temp/speciesjagsWood.density.R
res1 <- run.jags.b(model.files.jags.Tf.1[1], trait = 'Wood.density',
init.TF = FALSE, var.sample = 'wwf',
......
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