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

fix Makefile and Swiss

parent f1d3e6cc
......@@ -3,7 +3,7 @@ D2 := output/formatted
D3 := output/processed
D4 := figs/test.format.tree
D5 := figs/test.traits
D5 := figs/test.CWM
D6 := figs/test.CWM
sites:= Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS US
D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) ))
D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) ))
......@@ -22,186 +22,192 @@ TRY: $(D2)/TRY/data.TRY.std.rds
$(D2)/TRY/data.TRY.std.rds:
Rscript R/format.data/TRY.R
#-------------------------------------------------------
GLOBAL: $(D3)/Done.txt
$(D3)/Done.txt: R/process.data/global.R R/process.data/process-fun.R $(D2)/traits.std.global.csv $(D5)/Done.txt
Rscript $<
#-------------------------------------------------------
BCI: $(D3)/BCI/Done.txt
$(D3)/BCI/Done.txt: R/process.data/process-fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local'); process_bigplot_dataset('BCI', Rlim=15,std.traits='global');"
$(D3)/BCI/Done.txt: R/process.data/process-fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local');"
$(D2)/BCI/traits.csv: R/find.trait/BCI.R R/find.trait/trait-fun.R $(D2)/BCI/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/BCI/tree.csv: R/format.data/BCI.R $(shell find $(D1)/BCI -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Japan: $(D3)/Japan/Done.txt
$(D3)/Japan/Done.txt: R/process.data/process-fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local'); process_bigplot_dataset('Japan', Rlim=15,std.traits='global');"
$(D3)/Japan/Done.txt: R/process.data/process-fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local');"
$(D2)/Japan/traits.csv: R/find.trait/Japan.R R/find.trait/trait-fun.R $(D2)/Japan/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Japan/tree.csv: R/format.data/Japan.R $(shell find $(D1)/Japan -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Luquillo: $(D3)/Luquillo/Done.txt
$(D3)/Luquillo/Done.txt: R/process.data/process-fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='no');"
$(D3)/Luquillo/Done.txt: R/process.data/process-fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local');"
$(D2)/Luquillo/traits.csv: R/find.trait/Luquillo.R R/find.trait/trait-fun.R $(D2)/Luquillo/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Mbaiki: $(D3)/Mbaiki/Done.txt
$(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='local'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global');"
$(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='local');"
$(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait-fun.R $(D2)/Mbaiki/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Canada: $(D3)/Canada/Done.txt
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local'); process_inventory_dataset('Canada',std.traits='global');"
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local');"
$(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait-fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
France: $(D3)/France/Done.txt
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local'); process_inventory_dataset('France',std.traits='global');"
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local');"
$(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait-fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Fushan: $(D3)/Fushan/Done.txt
$(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='local'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');"
$(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='local');"
$(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait-fun.R $(D2)/Fushan/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
NSW: $(D3)/NSW/Done.txt
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local'); process_inventory_dataset('NSW',std.traits='global');"
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local');"
$(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait-fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
NVS: $(D3)/NVS/Done.txt
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local'); process_inventory_dataset('NVS',std.traits='global');"
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local');"
$(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait-fun.R $(D2)/NVS/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Paracou: $(D3)/Paracou/Done.txt
$(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='local'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global');"
$(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='local');"
$(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait-fun.R $(D2)/Paracou/tree.csv
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Spain: $(D3)/Spain/Done.txt
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local'); process_inventory_dataset('Spain',std.traits='global');"
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local');"
$(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait-fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Sweden: $(D3)/Sweden/Done.txt
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local'); process_inventory_dataset('Sweden',std.traits='global');"
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local');"
$(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait-fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
Swiss: $(D3)/Swiss/Done.txt
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local'); process_inventory_dataset('Swiss',std.traits='global');"
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local');"
$(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait-fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
US: $(D3)/US/Done.txt
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local'); process_inventory_dataset('US',std.traits='global');"
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local');"
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $<
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt
$(D2)/US/tree.csv: R/format.data/US.R $(shell find $(D1)/US -type f)
Rscript $<
Rscript $< ; rm -f $(D4)/Done.txt
#-------------------------------------------------------
TEST.TREE: $(D4)/Done.txt
......
load.model <- function(){
list(
name = "LOGLIN.RE",
## inits function
initial.values = function(){
list("intercept.trait"=runif(1,min=-1,max=1),
"slope.trait.e"=runif(1,min=-2,max=2),
"slope.trait.r"=runif(1,min=-2,max=2),
"mean.B"=runif(1,min=-0.5,max=1),
"sigma.B"=runif(1,min=0.0001,max=1.5),
"sigma"=runif(1,min=0.0001,max=4),
"sigma.alpha"=runif(1,min=0.0001,max=4),
"sigma.plot"=runif(1,min=0.0001,max=4))},
## parameters to save
parameters.to.save =c("alpha","B1",
"slope.trait.e","slope.trait.r","intercept.trait",
"sigma","sigma.plot","mean.B","sigma.B",
"sigma.alpha"),
bug = "################################################################################
######################## growth model with JAGS #############################
model {
############ Likelihood ###################
for (i in 1:N.indiv) {
logG[i] ~ dnorm(theo_g[i],tau)
theo_g[i] <- intercept + intercept_sp[species_id[i]]
+ intercept_plot[plot_id[i]]
+ intercept_ecocode[ecocode_id[i]]
+ param_logD_sp[species_id[i]]*logD[i]
+ param_Tf_ecocode[ecocode_id[i]]*Tf[i]
+ param_sumBn_ecocode[ecocode_id[i]]*sumBn[i]
+ param_sumTfBn_ecocode[ecocode_id[i]]*sumTfBn[i]
+ param_sumTnBn_ecocode[ecocode_id[i]]*sumTnBn[i]
}
################################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N_species)
{
intercept_sp[n] ~ dnorm(0,tau_intercept_sp)
param_logD_sp[n] ~ dnorm(0,tau_logD_sp)
}
## plot effect
for (j in 1:N_plot.species)
{
intercept_plot[j] ~ dnorm(0,tau_intercept_plot)
}
## ecocode effect
for (n in 1:N_ecocode)
{
intercept_ecocode[n] ~ dnorm(0,tau_intercept_ecocode)
param_Tf_ecocode[n] ~ dnorm(0,tau_Tf_ecocode)
param_sumBn_ecocode[n] ~ dnorm(0,tau_sumBn_ecocode)
param_sumTfBn_ecocode[n] ~ dnorm(0,tau_sumTfBn_ecocode)
param_sumTnBn_ecocode[n] ~ dnorm(0,tau_sumTnBn_ecocode)
}
###############################################
########### Non-hierarchical parameters ########
##### TODO NEED TO CHANGE FROM HERE
# constants for prior
tau0 <- 1.0E-4
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau.plot <- pow(sigma.plot,-2)
sigma.plot ~ dunif(0.00001,5)
## traits parameters
slope.trait.e ~ dnorm(0,tau0)
slope.trait.r ~ dnorm(0,tau0)
intercept.trait ~ dnorm(0,tau0)
## prior for hyper parameter mean and sd
for (l in 1:1)
{
mean.B[l]~ dnorm(0,tau0)
tau.B[l] <- pow(sigma.B[1],-2)
sigma.B[l] ~ dunif(0.000001,5)
}
tau.alpha <- pow(sigma.alpha,-2)
sigma.alpha ~ dunif(0.000001,5)
} # End of the jags model
")
}
load.model <- function(){
list(
name = "LOGLIN.RE",
## parameters to save
parameters.to.save <- c('intercept','mean_logD','mean_Tf','mean_sumBn',
'mean_sumTfBn','mean_sumTnBn','sigma_inter_sp',
'sigma_inter_ecocode','sigma_inter_plot','sigma_logD_sp',
'sigma_Tf_ecocode','sigma_sumBn_ecocode',
'sigma_sumTfBn_ecocode','sigma_sumTnBn_ecocode',
'sigma',
'param_Tf_ecocode',
'param_sumBn_ecocode','param_sumTfBn_ecocode',
'param_sumTnBn_ecocode','intercept_ecocode')
,
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_ecocode;
int species_id[N_indiv];
int plot_id[N_indiv];
int ecocode_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];
}
parameters{
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean_logD;
real<lower=-10,upper=10> mean_Tf;
real<lower=-10,upper=10> mean_sumBn;
real<lower=-10,upper=10> mean_sumTfBn;
real<lower=-10,upper=10> mean_sumTnBn;
real<lower=0.0001,upper=5> sigma_inter_sp;
real<lower=0.0001,upper=5> sigma_inter_ecocode;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_logD_sp;
real<lower=0.0001,upper=5> sigma_Tf_ecocode;
real<lower=0.0001,upper=5> sigma_sumBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTfBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTnBn_ecocode;
real<lower=0.0001,upper=5> sigma;
vector[N_species] param_logD_sp;
vector[N_species] intercept_sp;
vector[N_plot] intercept_plot;
vector[N_ecocode] param_Tf_ecocode;
vector[N_ecocode] param_sumBn_ecocode;
vector[N_ecocode] param_sumTfBn_ecocode;
vector[N_ecocode] param_sumTnBn_ecocode;
vector[N_ecocode] intercept_ecocode;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_sp[species_id[i]]
+ intercept_plot[plot_id[i]]
+ intercept_ecocode[ecocode_id[i]]
+ param_logD_sp[species_id[i]]*logD[i]
+ param_Tf_ecocode[ecocode_id[i]]*Tf[i]
+ param_sumBn_ecocode[ecocode_id[i]]*sumBn[i]
+ param_sumTfBn_ecocode[ecocode_id[i]]*sumTfBn[i]
+ param_sumTnBn_ecocode[ecocode_id[i]]*sumTnBn[i]
;
}
}
model {
# constants for prior
real sigma0;
################################################################################
######################## growth model with STAN #############################
###############################################
########### Hierarchical parameters ########
### species random param
param_logD_sp ~ normal(mean_logD,sigma_logD_sp);
intercept_sp ~ normal(0,sigma_inter_sp);
### plot random param
intercept_plot ~ normal(0,sigma_inter_plot);
### ecocode random param
param_Tf_ecocode ~ normal(mean_Tf,sigma_Tf_ecocode);
param_sumBn_ecocode ~ normal(mean_sumBn,sigma_sumBn_ecocode);
param_sumTfBn_ecocode ~ normal(mean_sumTfBn,sigma_sumTfBn_ecocode);
param_sumTnBn_ecocode ~ normal(mean_sumTnBn,sigma_sumTnBn_ecocode);
intercept_ecocode ~ normal(0,sigma_inter_ecocode);
###############################################
########### 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_sp~ uniform(0,6);
sigma_inter_ecocode~ uniform(0,6);
sigma_inter_plot~ uniform(0,6);
sigma_logD_sp~ uniform(0,6);
sigma_Tf_ecocode~ uniform(0,6);
sigma_sumBn_ecocode~ uniform(0,6);
sigma_sumTfBn_ecocode~ uniform(0,6);
sigma_sumTnBn_ecocode~ uniform(0,6);
sigma~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
}
load.model <- function () {
list(name="lmer.LOGLIN.AD.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTnTfBn.abs+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.E.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTnBn+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)") )
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf",
lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs"),
lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumTfBn+sumBn+sumTnBn+sumTnTfBn.abs"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumTfBn+sumBn+sumTnBn+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)"))
}
#TODO CHANEG THAT !!
load.model <- function () {
list(name="lmer.LOGLIN.HD.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTnTfBn.diff+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTnTfBn.diff-1|ecocode.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.R.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.nocomp.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+(logD-1|species.id)+(Tf-1|ecocode.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.simplecomp.Tf",
lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+(logD-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)"))
}
load.model <- function(){
list(
name = "LOGLIN.R",
## parameters to save
parameters.to.save <- c('intercept','mean_logD','mean_Tf','mean_sumBn',
'mean_sumTfBn','sigma_inter_sp',
'sigma_inter_ecocode','sigma_inter_plot','sigma_logD_sp',
'sigma_Tf_ecocode','sigma_sumBn_ecocode',
'sigma_sumTfBn_ecocode',
'sigma',
'param_Tf_ecocode',
'param_sumBn_ecocode','param_sumTfBn_ecocode',
'intercept_ecocode')
,
stan = "
data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_ecocode;
int species_id[N_indiv];
int plot_id[N_indiv];
int ecocode_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
real Tf[N_indiv];
real sumTfBn[N_indiv];
}
parameters{
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean_logD;
real<lower=-10,upper=10> mean_Tf;
real<lower=-10,upper=10> mean_sumBn;
real<lower=-10,upper=10> mean_sumTfBn;
real<lower=0.0001,upper=5> sigma_inter_sp;
real<lower=0.0001,upper=5> sigma_inter_ecocode;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_logD_sp;
real<lower=0.0001,upper=5> sigma_Tf_ecocode;
real<lower=0.0001,upper=5> sigma_sumBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTfBn_ecocode;
real<lower=0.0001,upper=5> sigma;
vector[N_species] param_logD_sp;
vector[N_species] intercept_sp;
vector[N_plot] intercept_plot;
vector[N_ecocode] param_Tf_ecocode;
vector[N_ecocode] param_sumBn_ecocode;
vector[N_ecocode] param_sumTfBn_ecocode;
vector[N_ecocode] intercept_ecocode;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_sp[species_id[i]]