Commit c698aae5 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

progress on optimiz stan and jags

parent dbd12966
......@@ -18,12 +18,29 @@ model.files.jags.Tf.3 <-
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)
}
return(a)
}
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.call.jags.and.save <- function(jags.model,
list.jags,
path.out,
var.sample,
iter = 50000,
warmup = 5000,
iter = 5000,
warmup = 500,
chains = 3,
thin = 50){
start <- Sys.time()
......@@ -31,8 +48,11 @@ fun.call.jags.and.save <- function(jags.model,
cat(jags.model$bug, file = file.path(path.out,"model.file.bug")
, sep=" ", fill = FALSE, labels = NULL, append = FALSE)
browser()
inits <- fun.generate.init(jags.model, list.jags$N_biomes, chains)
### SEND to jags
jags.output <- jags(data=list.jags,
jags.output <- jags(data=list.jags, inits = inits,
model.file = file.path(path.out,"model.file.bug"),
parameters.to.save = jags.model$pars,
n.chains = chains ,
......@@ -45,6 +65,7 @@ fun.call.jags.and.save <- function(jags.model,
print(jags.output)
saveRDS(jags.output, file = file.path(path.out, paste(var.sample,
"results.jags.rds", sep = '.')))
return(jags.output)
}
run.jags <- function (model.file, trait, init.TF,
......
......@@ -254,7 +254,7 @@ biomes.id <- factor(data.tree[['biomes']])
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTfBn <- fun.center.and.standardized.var(data.tree[[tf]])*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(logG = logG,
......@@ -269,7 +269,7 @@ return(data.frame(logG = logG,
plot.id = plot.id,
sumTnTfBn.diff = fun.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs = fun.standardized.var(sumTnTfBn.abs),
Tf = fun.standardized.var(data.tree[[tf]]),
Tf = fun.center.and.standardized.var(data.tree[[tf]]),
sumTnBn = fun.standardized.var(sumTnBn),
sumTfBn = fun.standardized.var(sumTfBn),
sumBn = fun.standardized.var(sumBn),
......@@ -301,14 +301,14 @@ sumTnBn.A_EV<- data.tree[[vec.CWM.tn[1]]]
sumTnBn.A_D <- data.tree[[vec.CWM.tn[2]]]
sumTnBn.C <- data.tree[[vec.CWM.tn[3]]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTfBn <- fun.center.and.standardized.var(data.tree[[tf]])*data.tree[[BATOT]]
sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *sumTfBn
sumBn <- data.tree[[BATOT]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *fun.center.and.standardized.var(data.tree[[tf]])
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *fun.center.and.standardized.var(data.tree[[tf]])
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *fun.center.and.standardized.var(data.tree[[tf]])
return(
data.frame(logG = logG,
logD = logD,
......
......@@ -5,15 +5,14 @@ load.model <- function(){
list(
name = "LOGLIN.ER.AD.Tf.fixed.biomes",
## parameters to save
pars = c('intercept' , 'mean_logD', 'mean_Tf',
'mean_sumBn', 'mean_sumTfBn', 'mean_sumTnBn',
'mean_sumTnTfBn_abs', 'sigma_inter_species',
pars = c('intercept' ,
'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_inter_tree', 'sigma_logD_species',
'sigma_sumBn_species',
'sigma', 'param_Tf_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumTnBn_biomes', 'param_sumTnTfBn_abs_biomes'),
'sigma', 'vec_Tf_biomes',
'vec_sumBn_biomes', 'vec_sumTfBn_biomes',
'vec_sumTnBn_biomes', 'vec_sumTnTfBn_abs_biomes'),
bug =
"#######################################################
######################## growth model with JAGS #######
......@@ -24,12 +23,12 @@ logG[i] ~ dnorm(theo_g[i],tau)
theo_g[i] <- inter[i] + DD[i] +TTf[i] + SUMBBN[i] +SUMTFBBN[i] + SUMTNBBN[i] + SUMABSBBN[i]
inter[i] <- intercept + intercept_species[species_id[i]] + intercept_tree[tree_id[i]]*tree_01[i] + intercept_plot[plot_id[i]] + intercept_set[set_id[i]]
DD[i] <- (mean_logD + param_logD_species[species_id[i]])*logD[i]
TTf[i] <- (param_Tf_biomes[biomes_id[i]])*Tf[i]
SUMBBN[i] <-(param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
SUMTNBBN[i] <- (param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
SUMABSBBN[i] <- (param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
DD[i] <- (mean_logD + vec_logD_species[species_id[i]])*logD[i]
TTf[i] <- (vec_Tf_biomes[biomes_id[i]])*Tf[i]
SUMBBN[i] <-(vec_sumBn_biomes[biomes_id[i]] + vec_sumBn_species[species_id[i]])*sumBn[i]
SUMTFBBN[i] <- (vec_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
SUMTNBBN[i] <- (vec_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
SUMABSBBN[i] <- (vec_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
}
################################################
......@@ -38,36 +37,36 @@ SUMABSBBN[i] <- (param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
### species random param
for (n in 1:N_species)
{
param_logD_species[n] ~ dnorm(0,tau_logD_species)T(-10, 10)
param_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)T(-10, 10)
intercept_species[n] ~ dnorm(0,tau_inter_species)T(-10, 10)
vec_logD_species[n] ~ dnorm(0,tau_logD_species)T(-7, 7)
vec_sumBn_species[n] ~ dnorm(0,tau_sumBn_species)T(-7, 7)
intercept_species[n] ~ dnorm(0,tau_inter_species)T(-7, 7)
}
## plot effect
for (j in 1:N_plot)
{
intercept_plot[j] ~ dnorm(0,tau_inter_plot)T(-10, 10)
intercept_plot[j] ~ dnorm(0,tau_inter_plot)T(-7, 7)
}
## tree effect
for (j in 1:N_tree)
{
intercept_tree[j] ~ dnorm(0,tau_inter_tree)T(-10, 10)
intercept_tree[j] ~ dnorm(0,tau_inter_tree)T(-7, 7)
}
## set effect
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)T(-10, 10)
intercept_set[j] ~ dnorm(0,tau_inter_set)T(-7, 7)
}
## biomes
for (j in 1:N_biomes)
{
param_Tf_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
param_sumBn_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
param_sumTfBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
param_sumTnBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
param_sumTnTfBn_abs_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
vec_Tf_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
vec_sumBn_biomes[j] ~ dnorm(0, tau0)T(-5, 5)
vec_sumTfBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
vec_sumTnBn_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
vec_sumTnTfBn_abs_biomes[j] ~ dnorm(0 ,tau0)T(-5, 5)
}
###############################################
########### Non-hierarchical parameters ########
......@@ -75,11 +74,9 @@ for (j in 1:N_biomes)
# constants for prior
tau0 <- 1.0E-4
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
#
intercept ~ dnorm(0,tau0)T(-5, 5)
mean_logD ~ dnorm(0,tau0)T(-5, 5)
# variance error
tau <- pow(sigma,-2)
......@@ -99,11 +96,11 @@ sigma_sumBn_species ~ dunif(0.0001,5)
#transformed variables
mean_Tf <- mean(param_Tf_biomes)
mean_sumBn <- mean(param_sumBn_biomes)
mean_sumTfBn <- mean(param_sumTfBn_biomes)
mean_sumTnBn <- mean(param_sumTnBn_biomes)
mean_sumTnTfBn_abs <- mean(param_sumTnTfBn_abs_biomes)
mean_Tf <- mean(vec_Tf_biomes)
mean_sumBn <- mean(vec_sumBn_biomes)
mean_sumTfBn <- mean(vec_sumTfBn_biomes)
mean_sumTnBn <- mean(vec_sumTnBn_biomes)
mean_sumTnTfBn_abs <- mean(vec_sumTnTfBn_abs_biomes)
} # End of the jags model
......
......@@ -43,17 +43,17 @@ parameters{
real<lower=0.0001,upper=3> sigma_logD_species;
real<lower=0.0001,upper=3> sigma_sumBn_species;
real<lower=0.0001,upper=3> sigma;
real<lower=-5,upper=5> param_Tf_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_species[N_species];
real<lower=-5,upper=5> param_sumTfBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnTfBn_abs_biomes[N_biomes];
real<lower=-5,upper=5> param_logD_species[N_species];
real<lower=-5,upper=5> intercept_species[N_species];
real<lower=-5,upper=5> intercept_plot[N_plot];
real<lower=-5,upper=5> intercept_set[N_set];
real<lower=-5,upper=5> intercept_tree[N_tree];
real<lower=-15,upper=15> param_Tf_biomes[N_biomes];
real<lower=-15,upper=15> param_sumBn_biomes[N_biomes];
real<lower=-15,upper=15> param_sumBn_species[N_species];
real<lower=-15,upper=15> param_sumTfBn_biomes[N_biomes];
real<lower=-15,upper=15> param_sumTnBn_biomes[N_biomes];
real<lower=-15,upper=15> param_sumTnTfBn_abs_biomes[N_biomes];
real<lower=-15,upper=15> param_logD_species[N_species];
real<lower=-15,upper=15> intercept_species[N_species];
real<lower=-15,upper=15> intercept_plot[N_plot];
real<lower=-15,upper=15> intercept_set[N_set];
real<lower=-15,upper=15> intercept_tree[N_tree];
}
transformed parameters {
# vector for prediction
......@@ -106,10 +106,6 @@ model {
# slope and intercept
intercept ~ normal(0,sigma0);
mean_logD ~ normal(0,sigma0);
mean_Tf ~ normal(0,sigma0);
mean_sumBn ~ normal(0,sigma0);
mean_sumTfBn ~ normal(0,sigma0);
mean_sumTnBn ~ normal(0,sigma0);
# sigma
sigma_inter_species~ uniform(0,6);
sigma_inter_set~ uniform(0,6);
......
load.model <- function(){
list(name = "LOGLIN.ER.AD.Tf.r.biomes",
pars = c('intercept' , 'mean_logD', 'mean_Tf',
'mean_sumBn', 'mean_sumTfBn', 'mean_sumTnBn',
'mean_sumTnTfBn_abs', 'sigma_inter_species',
'sigma_inter_set', 'sigma_inter_plot',
'sigma_inter_tree', 'sigma_logD_species',
'sigma_Tf_biomes', 'sigma_sumBn_biomes',
'sigma_sumBn_species', 'sigma_sumTfBn_biomes',
'sigma_sumTnBn_biomes', 'sigma_sumTnTfBn_abs_biomes',
'sigma', 'param_Tf_biomes',
'param_sumBn_biomes', 'param_sumTfBn_biomes',
'param_sumTnBn_biomes', 'param_sumTnTfBn_abs_biomes'),## parameters to save
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_biomes;
int<lower=0> N_set;
int<lower=0> N_tree;
int species_id[N_indiv];
int plot_id[N_indiv];
int biomes_id[N_indiv];
int set_id[N_indiv];
int tree_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
real Tf[N_indiv];
real sumTfBn[N_indiv];
real sumTnBn[N_indiv];
real sumTnTfBn_abs[N_indiv];
real tree_01[N_indiv];
}
parameters{
real<lower=-5,upper=5> intercept;
real<lower=-5,upper=5> mean_logD;
real<lower=-5,upper=5> mean_Tf;
real<lower=-5,upper=5> mean_sumBn;
real<lower=-5,upper=5> mean_sumTfBn;
real<lower=-5,upper=5> mean_sumTnBn;
real<lower=-5,upper=5> mean_sumTnTfBn_abs;
real<lower=0.0001,upper=3> sigma_inter_species;
real<lower=0.0001,upper=3> sigma_inter_set;
real<lower=0.0001,upper=3> sigma_inter_plot;
real<lower=0.0001,upper=3> sigma_inter_tree;
real<lower=0.0001,upper=3> sigma_logD_species;
real<lower=0.0001,upper=3> sigma_Tf_biomes;
real<lower=0.0001,upper=3> sigma_sumBn_biomes;
real<lower=0.0001,upper=3> sigma_sumBn_species;
real<lower=0.0001,upper=3> sigma_sumTfBn_biomes;
real<lower=0.0001,upper=3> sigma_sumTnBn_biomes;
real<lower=0.0001,upper=3> sigma_sumTnTfBn_abs_biomes;
real<lower=0.0001,upper=3> sigma;
real<lower=-5,upper=5> param_Tf_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumBn_species[N_species];
real<lower=-5,upper=5> param_sumTfBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnBn_biomes[N_biomes];
real<lower=-5,upper=5> param_sumTnTfBn_abs_biomes[N_biomes];
real<lower=-5,upper=5> param_logD_species[N_species];
real<lower=-5,upper=5> intercept_species[N_species];
real<lower=-5,upper=5> intercept_plot[N_plot];
real<lower=-5,upper=5> intercept_set[N_set];
real<lower=-5,upper=5> intercept_tree[N_tree];
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_species[species_id[i]]
+ intercept_tree[tree_id[i]]*tree_01[i]
+ intercept_plot[plot_id[i]]
+ intercept_set[set_id[i]]
+ (mean_logD + param_logD_species[species_id[i]])*logD[i]
+ (mean_Tf +param_Tf_biomes[biomes_id[i]])*Tf[i]
+ (mean_sumBn + param_sumBn_biomes[biomes_id[i]] + param_sumBn_species[species_id[i]])*sumBn[i]
+ (mean_sumTfBn + param_sumTfBn_biomes[biomes_id[i]])*sumTfBn[i]
+ (mean_sumTnBn + param_sumTnBn_biomes[biomes_id[i]])*sumTnBn[i]
+ (mean_sumTnTfBn_abs + param_sumTnTfBn_abs_biomes[biomes_id[i]])*sumTnTfBn_abs[i]
;
}
}
model {
# constants for prior
real sigma0;
################################################################################
######################## growth model with STAN #############################
###############################################
########### Hierarchical parameters ########
### species random param
param_logD_species ~ normal(0,sigma_logD_species);
param_sumBn_species ~ normal(0,sigma_sumBn_species);
intercept_species ~ normal(0, sigma_inter_species);
### plot random param
intercept_plot ~ normal(0,sigma_inter_plot);
### tree random effect
intercept_tree ~ normal(0, sigma_inter_tree);
intercept_set ~ normal(0, sigma_inter_set);
### biomes random param
param_Tf_biomes ~ normal(0, sigma_Tf_biomes);
param_sumBn_biomes ~ normal(0,sigma_sumBn_biomes);
param_sumTfBn_biomes ~ normal(0 ,sigma_sumTfBn_biomes);
param_sumTnBn_biomes ~ normal(0 ,sigma_sumTnBn_biomes);
param_sumTnTfBn_abs_biomes ~ normal(0 , sigma_sumTnTfBn_abs_biomes);
###############################################
########### Non-Hierarchical parameters ########
# constant for prior
sigma0 <- 10;
# slope and intercept
intercept ~ normal(0,sigma0);
mean_logD ~ normal(0,sigma0);
mean_Tf ~ normal(0,sigma0);
mean_sumBn ~ normal(0,sigma0);
mean_sumTfBn ~ normal(0,sigma0);
mean_sumTnBn ~ normal(0,sigma0);
# sigma
sigma_inter_species~ uniform(0,6);
sigma_inter_set~ uniform(0,6);
sigma_inter_plot~ uniform(0,6);
sigma_inter_tree~ uniform(0,6);
sigma_logD_species~ uniform(0,6);
sigma_sumBn_species~ uniform(0,6);
sigma_Tf_biomes~ uniform(0,6);
sigma_sumBn_biomes~ uniform(0,6);
sigma_sumTfBn_biomes~ uniform(0,6);
sigma_sumTnBn_biomes~ uniform(0,6);
sigma_sumTnTfBn_abs_biomes~ uniform(0,6);
sigma~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
}
......@@ -38,6 +38,7 @@ fun.call.stan.and.save <- function(stan.model, list.stan, path.out,
print(stan.output)
saveRDS(stan.output, file = file.path(path.out, paste(var.sample,
"results.stan", chain_id, "rds", sep = '.')))
return(stan.output)
}
run.stan <- function (model.file, trait, init.TF,
......
......@@ -19,8 +19,8 @@ res <- run.multiple.model.for.set.one.trait(model.files.jags.Tf.1,
type.filling='species',
sample.size = 10000,
var.sample = 'ecocode.id',
iter = 5000,
warmup = 500, chains = 3,
iter = 50,
warmup = 10, chains = 2,
thin = 5,
init.TF = FALSE))
......
......@@ -16,13 +16,12 @@ run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[4], run.stan,
iter = 200, warmup = 20, chains = 2, chain_id = '1',
init.TF = TRUE)
t1000 <- system.time(
res <- run.multiple.model.for.set.one.trait(model.files.stan.Tf.1[4], run.stan,
trait = 'SLA',type.filling='species',
sample.size = 1000,
var.sample = 'ecocode.id',
iter = 200, warmup = 20, chains = 2, chain_id = '1',
init.TF = FALSE))
init.TF = FALSE)
......
#######################################
#######################################
## % Plot trends in communities structure in global data set
##############################
##### load scripts
source("R/process.data/process-fun.R")
source("R/process.data/test.tree.CWM-fun.R")
source("R/utils/plot.R")
source("R/analysis/lmer.run.R")
##########
## load summarise data
data.summarise <- readRDS('output/data.summarise.rds')
# reorder factor
data.summarise$biomes <- factor(data.summarise$biomes, c("Tundra",
"Boreal forest",
"Temperate forest",
"Woodland shrubland",
"Temperate rain forest",
"Temperate grassland desert",
"Subtropical desert",
"Tropical forest savanna",
"Tropical rain forest"))
## shorter name for biomes
names.biomes <- levels(data.summarise$biomes)
names.biomes[2] <- 'Boreal'
names.biomes[3] <- 'Temp forest'
names.biomes[4] <- 'Medit'
names.biomes[5] <- 'Temp rain forest'
names.biomes[6] <- 'Temp grassland'
names.biomes[7] <- 'Subtrop'
names.biomes[8] <- 'Trop forest'
names.biomes[8] <- 'Trop rain forest'
## DO PLOT OF PATTERN
## TODO IN PLOT REMOVE BIOMES Tundra and Subtropical desert
data.temp <- data.summarise[!data.summarise$biomes %in% c('Tundra',
'Temperate grassland desert',
'Subtropical desert'), ]
data.temp$biomes <- factor(data.temp$biomes)
## shorter name for biomes
names.biomes.t <- levels(factor(data.temp$biomes))
names.biomes.t[1] <- 'Boreal'
names.biomes.t[2] <- 'Temp forest'
names.biomes.t[3] <- 'Medit'
names.biomes.t[4] <- 'Temp rain forest'
names.biomes.t[5] <- 'Trop forest'
names.biomes.t[6] <- 'Trop rain forest'
## need to look at what is Temperate grassland desert
## Total basal area for each plots
m <- matrix(c(1:2), 1, 2)
layout(m, widths=c(4,1.),
heights=c(1))
par.old <- par()
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(BATOT ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('Basal area ' ,(cm^2/m^2))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names.biomes.t,
fill = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
bty = 'n', cex = 1.1)
par(par.old)
## species diversity each plots
m <- matrix(c(1:2,3,3), 2, 2)
layout(m, widths=c(4,1),
heights=c(1,1))
par.old <- par()
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(n_sp~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('Number of species' )),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(exp(shannon)~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('Number of species' )),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names.biomes.t,
fill = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
bty = 'n', cex = 1)
par(par.old)
## traits means
m <- matrix(c(1:4), 2, 2)
layout(m, widths=c(4,4),
heights=c(1,1))
par.old <- par()
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(mean.SLA ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('SLA ' ,(mm^2/mg))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(mean.Wood.density ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('Wood density ' ,( mg/mm^3))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(mean.Max.height ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('Max height ' ,(m))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names.biomes.t,
fill = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
bty = 'n', cex = 1.2)
par(par.old)
## traits sd
m <- matrix(c(1:4), 2, 2)
layout(m, widths=c(4,4),
heights=c(1,1))
par.old <- par()
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(sd.SLA ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('SD SLA ' ,(mm^2/mg))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(sd.Wood.density ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('SD Wood density ' ,( mg/mm^3))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mgp=c(1.5,0.5,0), mar= c(.5, 3.1, 0.5, 0.))
boxplot(sd.Max.height ~ biomes, data.temp,
col = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
ylab = expression(paste('SD Max height ' ,(m))),
outline = FALSE, xaxt = 'n', cex.y = 1.5,cex.lab=1.2)
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names.biomes.t,
fill = fun.col.pch.biomes()$col.vec[levels(data.temp$biomes)],
bty = 'n', cex = 1.2)
par(par.old)