Commit 8550a229 authored by kunstler's avatar kunstler
Browse files

issue in process data

parent 5910f7da
......@@ -16,25 +16,54 @@ model.files.jags.Tf.2 <-
fun.generate.init.one.param <- function(i, pars, list.jags,chain){
fun.generate.init.one.param <- function(i, pars, list.jags,chain, chains.v){
if(grepl('sigma',names(pars)[i])){
a <- 0.1+(chain-1)*0.1
a <- 0.05+(chain-min(chains.v))/max(chains.v))*0.5
}else{
if(grepl('param',names(pars)[i])){
a.t <- -0.3+(chain-1)*0.2
if(grepl('logD', names(pars)[i])){
a.t <- 2/3+(chain-2.5)*0.15
a <- rnorm(list.jags[[pars[i]]],mean = a.t, sd = 0.1)
}else{
if(grepl('intercept', names(pars)[i])){
a.t <- 0.2+(chain-2.5)*0.1
a <- rnorm(list.jags[[pars[i]]],mean = a.t, sd = 0.1)
}else{
if(grepl('sumBn', names(pars)[i])){
a.t <- -0.25+(chain-2.5)*0.1
a <- rnorm(list.jags[[pars[i]]],mean = a.t, sd = 0.1)
}else{
a.t <- 0+(chain-2.5)*0.15
a <- rnorm(list.jags[[pars[i]]],mean = a.t, sd = 0.1)
}
}
}
}else{
a <- -0.3+(chain-1)*0.2
if(grepl('logD', names(pars)[i])){
a <- 2/3+(chain-2.5)*0.15
}else{
if(grepl('intercept', names(pars)[i])){
a <- 0.2+(chain-2.5)*0.1
}else{
if(grepl('sumBn', names(pars)[i])){
a <- -0.25+(chain-2.5)*0.1
}else{
a <- 0+(chain-2.5)*0.15
}
}
}
}
}
return(a)
}
fun.generate.init.all.param <- function(chain, list.jags, jags.model){
fun.generate.init.all.param <- function(chain, list.jags, jags.model, chains.v){
list.init <- lapply(seq_len(length(jags.model$pars)),
fun.generate.init.one.param,
jags.model$pars,
list.jags,
chain)
chain,
chains.v)
names(list.init) <- names(jags.model$pars)
return(list.init)
......@@ -42,7 +71,7 @@ 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)
inits <- lapply(chains.v,fun.generate.init.all.param, list.jags, jags.model, chains.v)
return(inits)
}
......@@ -122,6 +151,7 @@ print('OK')
select.set. = select.set,
sample.vec.TF = TRUE,
merge.biomes.TF = merge.biomes.TF)
browser()
cat("Ok data with Nobs", jags.list[[1]]$N_indiv, 'run real',
"\n")
res <- fun.call.jags.and.save(jags.model = model,
......
......@@ -337,6 +337,14 @@ scale.d <- function(x, ...){
return(as.vector(scale(x, ...)))
}
scale.nc <- function(x, center = FALSE){
return(as.vector(scale(x,
center = center,
scale = apply(x,2, sd, na.rm = TRUE))))
}
get.variables <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, ecocode.var = 'wwf',
......@@ -369,11 +377,11 @@ return(data.frame(logG = logG,
biomes.id = biomes.id,
plot.id = plot.id,
tree.id = tree.id,
sumTnTfBn.abs = scale.d(sumTnTfBn.abs, center = FALSE),
sumTnTfBn.abs = scale.nc(sumTnTfBn.abs, center = FALSE),
Tf = scale.d(data.tree[[tf]]),
sumTnBn = scale.d(sumTnBn, center = FALSE),
sumTfBn = scale(sumTfBn, center = FALSE),
sumBn = scale(sumBn, center = FALSE),
sumTnBn = scale.nc(sumTnBn, center = FALSE),
sumTfBn = scale.nc(sumTfBn, center = FALSE),
sumBn = scale.nc(sumBn, center = FALSE),
stringsAsFactors = FALSE))
}
......@@ -396,13 +404,13 @@ biomes.id <- factor(data.tree[['biomes']])
if(Multi.type == 'a') {
traits <- c('SLA', 'Wood.density', 'Max.height')
sumTnTfBn.abs.Multi <-
scale.d(data.tree[['Multi1.abs.CWM.fill']]*
scale.nc(data.tree[['Multi1.abs.CWM.fill']]*
data.tree[[BATOT]], center = FALSE)
}
if(Multi.type == 'b') {
traits <- c('SLA', 'Wood.density', 'Max.height', 'Seed.mass')
sumTnTfBn.abs.Multi <-
scale.d(data.tree[['Multi2.abs.CWM.fill']]*
scale.nc(data.tree[['Multi2.abs.CWM.fill']]*
data.tree[[BATOT]], center = FALSE)
}
......@@ -414,11 +422,11 @@ for (t in traits){
abs.CWM.tntf <- paste0(t, '.abs.CWM.fill')
CWM.tn <- paste0(t, '.CWM.fill')
tf <- paste0(t, '.focal')
list.sumTnTfBn.abs[[t]] <- scale.d(data.tree[[abs.CWM.tntf]]*
list.sumTnTfBn.abs[[t]] <- scale.nc(data.tree[[abs.CWM.tntf]]*
data.tree[[BATOT]], center = FALSE)
list.sumTnBn[[t]] <- scale.d(data.tree[[CWM.tn]]*data.tree[[BATOT]],
list.sumTnBn[[t]] <- scale.nc(data.tree[[CWM.tn]]*data.tree[[BATOT]],
center = FALSE)
list.sumTfBn[[t]] <- scale.d(data.tree[[tf]]*data.tree[[BATOT]],
list.sumTfBn[[t]] <- scale.nc(data.tree[[tf]]*data.tree[[BATOT]],
center = FALSE)
list.Tf[[t]] <- scale.d(data.tree[[tf]])
}
......@@ -493,17 +501,17 @@ data.frame(logG = logG,
plot.id = plot.id,
tree.id = tree.id,
sumTnTfBn.abs =
scale.d(sumTnTfBn.abs, center = FALSE),
scale.nc(sumTnTfBn.abs, center = FALSE),
Tf.A_EV= scale.d(Tf.A_EV),
Tf.A_D= scale.d(Tf.A_D),
Tf.C= scale.d(Tf.C),
sumTnBn.A_EV= scale.d(sumTnBn.A_EV, center = FALSE),
sumTnBn.A_D= scale.d(sumTnBn.A_D, center = FALSE),
sumTnBn.C = scale.d(sumTnBn.C, center = FALSE),
sumTfBn.A_EV= scale.d(sumTfBn.A_EV, center = FALSE),
sumTfBn.A_D= scale.d(sumTfBn.A_D, center = FALSE),
sumTfBn.C = scale.d(sumTfBn.C, center = FALSE),
sumBn = scale.d(sumBn, center = FALSE)))
sumTnBn.A_EV= scale.nc(sumTnBn.A_EV, center = FALSE),
sumTnBn.A_D= scale.nc(sumTnBn.A_D, center = FALSE),
sumTnBn.C = scale.nc(sumTnBn.C, center = FALSE),
sumTfBn.A_EV= scale.nc(sumTfBn.A_EV, center = FALSE),
sumTfBn.A_D= scale.nc(sumTfBn.A_D, center = FALSE),
sumTfBn.C = scale.nc(sumTfBn.C, center = FALSE),
sumBn = scale.nc(sumBn, center = FALSE)))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+Tf+Tf:biomes.id+logD+sumBn+sumBn:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)))
}
......@@ -8,7 +8,7 @@ source("R/analysis/lmer.run.R")
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.simple.set.0.rds',
models = c(model.files.lmer.Tf.1),
models = c(model.files.lmer.Tf.0),
traits = c("SLA", "Wood.density", "Max.height")
)
......@@ -19,7 +19,7 @@ format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
)
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.simple.set.TF.rds',
list.file.name = 'list.lmer.out.all.NA.simple.set.TP.rds',
models = c(model.files.lmer.Tf.3),
traits = c("SLA", "Wood.density", "Max.height")
)
......
......@@ -108,21 +108,32 @@ dev.off()
list.all.results.set <-
readRDS('output/list.lmer.out.all.NA.simple.set.rds')
list.all.results.set.TP <-
readRDS('output/list.lmer.out.all.NA.simple.set.TF.rds')
list.all.results.set.0 <-
readRDS('output/list.lmer.out.all.NA.simple.set.0.rds')
fun.get.conv(list.all.results.set)
vec.rel.grad.set <- sapply(list.all.results.set,
function(list.t) { max(abs(list.t[['relgrad']]))})
vec.rel.grad.set < 0.001
fun.get.conv(list.all.results.set.TP)
vec.rel.grad.set.TP <- sapply(list.all.results.set.TP,
function(list.t) { max(abs(list.t[['relgrad']]))})
vec.rel.grad.set.TP < 0.001
# Figmean
pdf('figs/figres1.pdf', height = 7, width = 14)
plot.param(list.all.results.set , data.type = "simple",
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species',
plot.param(list.all.results.set.TP , data.type = "simple",
model = 'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species',
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
param.print = 1:6,
param.print = 1:5,
param.names = c(expression(paste('Trait sim ',
(alpha[s] %*% abs(t[c] - t[f])))),
expression(paste('Response ',
......
......@@ -43,7 +43,7 @@ for(trait in c('SLA', 'Wood.density', 'Max.height')){
## ### PROCESS OUTPUT !! MESSY
## res <- as.mcmc(res1)
res <- as.mcmc(output.jags[[1]])
## crosscorr.plot(res)
## rejectionRate(res)
## gelman.diag(res)
......@@ -59,46 +59,46 @@ for(trait in c('SLA', 'Wood.density', 'Max.height')){
## ggs_crosscorrelation(res.gg.sel)
## res.df <- data.frame(do.call('rbind', res))
## res.df.sel <- res.df[,select.vars]
## library(IDPmisc)
res.df <- data.frame(do.call('rbind', res))
res.df.sel <- res.df[,select.vars]
library(IDPmisc)
## panel.hist <- function(x, ...)
## {
## usr <- par("usr"); on.exit(par(usr))
## par(usr = c(usr[1:2], 0, 1.5) )
## h <- hist(x, plot = FALSE)
## breaks <- h$breaks; nB <- length(breaks)
## y <- h$counts; y <- y/max(y)
## rect(breaks[-nB], 0, breaks[-1], y, col="blue4", ...)
## }
panel.hist <- function(x, ...)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(usr[1:2], 0, 1.5) )
h <- hist(x, plot = FALSE)
breaks <- h$breaks; nB <- length(breaks)
y <- h$counts; y <- y/max(y)
rect(breaks[-nB], 0, breaks[-1], y, col="blue4", ...)
}
## panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
## {
## usr <- par("usr"); on.exit(par(usr))
## par(usr = c(0, 1, 0, 1))
## r <- abs(cor(x, y, method = "spearman"))
## txt <- format(c(r, 0.123456789), digits=digits)[1]
## txt <- paste(prefix, txt, sep="")
## if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
## text(0.5, 0.5, txt, cex = cex * r)
## }
panel.cor <- function(x, y, digits=2, prefix="", cex.cor)
{
usr <- par("usr"); on.exit(par(usr))
par(usr = c(0, 1, 0, 1))
r <- abs(cor(x, y, method = "spearman"))
txt <- format(c(r, 0.123456789), digits=digits)[1]
txt <- paste(prefix, txt, sep="")
if(missing(cex.cor)) cex <- 0.8/strwidth(txt)
text(0.5, 0.5, txt, cex = cex * r)
}
## betterPairs <- function(YourData){
## return(pairs(YourData, lower.panel=function(...) {par(new=TRUE);ipanel.smooth(...)}, diag.panel=panel.hist, upper.panel=panel.cor))
## }
betterPairs <- function(YourData){
return(pairs(YourData, lower.panel=function(...) {par(new=TRUE);ipanel.smooth(...)}, diag.panel=panel.hist, upper.panel=panel.cor))
}
## # Example
# Example
## betterPairs(res.df.sel[, 1:7])
betterPairs(as.data.frame(res[, 1:18]))
## 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)
res2 <- run.jags.b(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 = FALSE)
......
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