Commit 8e54c6f4 authored by kunstler's avatar kunstler
Browse files

new launch jags

parent 37afb3ae
......@@ -36,6 +36,7 @@ list.init <- lapply(seq_len(length(jags.model$pars)),
list.jags,
chain)
names(list.init) <- names(jags.model$pars)
return(list.init)
}
......@@ -53,24 +54,11 @@ fun.call.jags.and.save <- function(jags.model,
warmup = 500,
chains = 3,
thin = 50,
parallel = FALSE){
method = 'rjags'){
require(runjags)
start <- Sys.time()
inits <- fun.generate.init(jags.model, list.jags[[1]], chains)
### SEND to jags
if(!parallel){
jags.output <- run.jags(data=list.jags[[1]],
inits = inits,
model = file.path(path.out,"model.file.bug"),
monitor = names(jags.model$pars),
n.chains = chains ,
sample = iter,
burnin = warmup,
adapt = 2000,
thin = thin,
modules= 'glm on',
method = 'rjags')
}else{
jags.output <- run.jags(data=list.jags[[1]],
inits = inits,
model = file.path(path.out,"model.file.bug"),
......@@ -81,8 +69,7 @@ inits <- fun.generate.init(jags.model, list.jags[[1]], chains)
adapt = 2000,
thin = thin,
modules= 'glm on',
method = 'parallel')
}
method = method)
end <- Sys.time()
print(end -start)
print(jags.output)
......
......@@ -87,9 +87,12 @@ mat.param <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
mat.param.sd <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.param.sd, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
mat.R2 <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
mat.R2c <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2c, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
mat.R2m <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2m, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
bold.index <- which(((mat.param - 1.96*mat.param.sd) >0 & mat.param > 0) |
((mat.param + 1.96*mat.param.sd) <0 & mat.param <0),
arr.ind = TRUE)
......@@ -98,17 +101,19 @@ mat.param.mean.sd <- matrix(paste0(round(mat.param, 3),
round(mat.param.sd, 3),
')'), ncol = 3)
mat.param <- rbind(mat.param.mean.sd,
round(mat.R2, 4))
round(mat.R2m, 4),
round(mat.R2c, 4))
colnames(mat.param) <- c('Wood density', 'SLA', 'Maximum height')
row.names(mat.param) <- c('$\\gamma$', '$m_1$', '$\\alpha_0$',
'$\\alpha_i$', '$\\alpha_r$',
'$\\alpha_s$', '$R^2$*')
'$\\alpha_s$', '$R^2_m$*', '$R^2_c$*')
##+ Table2_Effectsize, echo = FALSE, results='asis', message=FALSE
pandoc.table(mat.param, caption = "Standaridized parameters estimates and standard error (in bracket) estimated for each traits and $R^2$* of models. See Fig 1. in main text for explanation of parameters",
digits = 3, justify = c('left', rep('right', 3)),
emphasize.strong.cells = bold.index, split.tables = 200)
## \* We report the conditional $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133–142 (2013).
## \* We report the conditional and marginal $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133–142 (2013), modified by Johnson, P. C. D. Extension of Nakagawa and Schielzeth’s R2GLMM to random slopes models. Methods in Ecology and Evolution 5, 944–946 (2014).
.
......@@ -91,9 +91,12 @@ mat.param <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
mat.param.sd <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.param.sd, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
mat.R2 <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
mat.R2c <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2c, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
mat.R2m <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2m, list.res = list.all.results,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
bold.index <- which(((mat.param - 1.96*mat.param.sd) >0 & mat.param > 0) |
((mat.param + 1.96*mat.param.sd) <0 & mat.param <0),
arr.ind = TRUE)
......@@ -102,11 +105,12 @@ mat.param.mean.sd <- matrix(paste0(round(mat.param, 3),
round(mat.param.sd, 3),
')'), ncol = 3)
mat.param <- rbind(mat.param.mean.sd,
round(mat.R2, 4))
round(mat.R2m, 4),
round(mat.R2c, 4))
colnames(mat.param) <- c('Wood density', 'SLA', 'Maximum height')
row.names(mat.param) <- c('$\\gamma$', '$m_1$', '$\\alpha_0$',
'$\\alpha_i$', '$\\alpha_r$',
'$\\alpha_s$', '$R^2$*')
'$\\alpha_s$', '$R^2_m$*', '$R^2_c$*')
```
``` {r Table2_Effectsize, echo = FALSE, results='asis', message=FALSE}
......@@ -115,6 +119,9 @@ pandoc.table(mat.param, caption = "Standaridized parameters estimates and standa
emphasize.strong.cells = bold.index, split.tables = 200)
```
\* We report the conditional $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133142 (2013).
\* We report the conditional and marginal $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133142 (2013), modified by Johnson, P. C. D. Extension of Nakagawa and Schielzeths R2GLMM to random slopes models. Methods in Ecology and Evolution 5, 944946 (2014).
``` {r }
.
```
......@@ -19,11 +19,11 @@
-------------------------------------------------------------------------------------------------------------
set # of trees # of species # of plots/quadrats % of angiosperm % of evergreen
------------------------ ------------ -------------- --------------------- ----------------- ----------------
Sweden 202469 26 22552 27.0 72.9
Sweden 202519 26 22552 27.0 73.0
New Zealand 53775 117 1415 94.0 99.1
US 1369815 493 59840 63.3 37.2
US 1370481 493 59840 63.4 37.2
Canada 495008 75 14983 34.4 64.9
......@@ -31,21 +31,21 @@ Australia 906 101 63
France 184316 127 17611 74.1 28.5
Switzerland 28207 65 2597 36.1 55.5
Switzerland 28301 61 2597 36.4 55.2
Spain 418805 122 36462 34.7 81.6
Panama 26994 238 2033 99.8 77.7
Panama 27028 238 2033 99.8 77.6
Paracou 46382 714 2157 100.0 83.5
Paracou 46513 716 2157 100.0 83.5
Japan 4637 138 318 72.2 70.8
Japan 4619 140 318 72.7 70.4
Taiwan 14701 72 623 92.0 75.3
Puerto Rico 14011 82 399 100.0 99.0
Central African Republic 17602 203 989 99.5 72.5
Central African Republic 17591 204 989 99.5 72.4
-------------------------------------------------------------------------------------------------------------
Table: Data description, with number of individual tree, species and plot in NFI data and quadrat in LPP data, and percentage of angiosperm and evergreen species.
......@@ -97,23 +97,33 @@ Table: Traits coverage in each sites. Percentage of species with species level t
-------------------------------------------------------------------------
&nbsp; Wood density SLA Maximum height
---------------- ------------------ ------------------ ------------------
**$\gamma$** **0.426 (0.014)** **0.409 (0.014)** **0.422 (0.013)**
**$\gamma$** **0.441 (0.014)** **0.409 (0.014)** **0.429 (0.013)**
**$m_1$** **-0.125 (0.042)** **0.108 (0.048)** 0.058 (0.043)
**$m_1$** **-0.121 (0.039)** **0.115 (0.051)** 0.067 (0.044)
**$\alpha_0$** **-0.323 (0.056)** **-0.246 (0.077)** **-0.336 (0.068)**
**$\alpha_0$** **-0.308 (0.058)** **-0.259 (0.069)** **-0.349 (0.066)**
**$\alpha_i$** -0.022 (0.017) **0.082 (0.027)** -0.023 (0.027)
**$\alpha_i$** -0.019 (0.017) **0.076 (0.025)** -0.025 (0.026)
**$\alpha_r$** **0.058 (0.021)** 0.018 (0.039) **-0.08 (0.037)**
**$\alpha_r$** **0.056 (0.022)** -0.007 (0.037) **-0.088 (0.039)**
**$\alpha_s$** **0.041 (0.016)** 0.055 (0.029) **0.062 (0.016)**
**$\alpha_s$** **0.042 (0.015)** **0.065 (0.028)** **0.064 (0.016)**
**$R^2$*** 0.7072 0.7514 0.7156
**$R^2_m$*** 0.1242 0.1355 0.1203
**$R^2_c$*** 0.717 0.7373 0.7168
-------------------------------------------------------------------------
Table: Standaridized parameters estimates and standard error (in bracket) estimated for each traits and $R^2$* of models. See Fig 1. in main text for explanation of parameters
\* We report the conditional $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133–142 (2013).
\* We report the conditional and marginal $R^2$ of the models using the methods of Nakagawa, S. & Schielzeth, H. A general and simple method for obtaining R2 from generalized linear mixed-effects models. Methods in Ecology and Evolution 4, 133–142 (2013), modified by Johnson, P. C. D. Extension of Nakagawa and Schielzeth’s R2GLMM to random slopes models. Methods in Ecology and Evolution 5, 944–946 (2014).
```r
.
```
```
## Error in eval(expr, envir, enclos): object '.' not found
```
......@@ -4,8 +4,7 @@
Competition is very important to understand and predict the dynamics
of plant community composition. In terrestrial
vegetation, where plants strongly modify the environment in their
immediate neighbourhood, competition is conspusious and influences growth and survival of
individuals. Predicting competition via phenotypic traits may be
immediate neighbourhood, competition is conspusious but our ability to predict its consequences on plant performances is extremely limited. Predicting competition via phenotypic traits may be
simpler and more general than via competition coefficients between
each pair of species, an intractable approach at global scale. Here
for the first time we show how three functional traits - wood density,
......
......@@ -25,10 +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 = $samplesize, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, parallel = TRUE);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(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, parallel = TRUE);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 = $samplesize, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = FALSE, merge.biomes.TF = TRUE, parallel = TRUE);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(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, parallel = TRUE);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
......
......@@ -24,7 +24,7 @@ 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), run.jags.b,trait = $trait,data.type='simple', sample.size = $samplesize, var.sample = 'wwf',iter = $iter, warmup = $warmup, chains = 3, thin = $thin, init.TF = TRUE);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(c(model.files.jags.Tf.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 = TRUE);print('done')\"" > Rscript_temp/speciesjags${trait}.sh
qsub Rscript_temp/speciesjags${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "jags${trait}" -q opt32G -j oe
......
......@@ -173,6 +173,14 @@ plot.growth.ba(traits = c('Wood.density', 'SLA', 'Max.height'),
labels.leg = c('High trait value', 'Low trait value'))
dev.off()
mat.R2c <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2c, list.res = list.all.results.set,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
mat.R2m <- do.call('cbind', lapply(c('Wood.density', 'SLA', 'Max.height'),
extract.R2m, list.res = list.all.results.set,
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species'))
##### PLOTS FOR SLIDES
pdf('figs/fig_res_fix_slide1.pdf', height = 7, width = 14)
plot.param(list.all.results.set , data.type = "simple",
......
......@@ -16,13 +16,13 @@ res1 <- run.jags.b(model.files.jags.Tf.1[1], trait = 'Wood.density',
init.TF = FALSE, var.sample = 'wwf',
sample.size = 20000, iter = 2000,
warmup = 1000, chains = 3,
thin = 2, parallel = TRUE)
thin = 2, method= 'parallel')
res1 <- run.jags.b(model.files.jags.Tf.1[2], trait = 'Wood.density',
init.TF = FALSE, var.sample = 'wwf',
sample.size = 10000, iter = 2000,
warmup = 1000, chains = 4,
thin = 2, parallel = TRUE, merge.biomes.TF = TRUE)
thin = 2, method = 'parallel', merge.biomes.TF = TRUE)
### PROCESS OUTPUT !! MESSY
......
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