new koppen ecocode

R/analysis/jags.output-fun.R

@@ -16,16 +16,44 @@ summarise.jags.output <- function(x,
 	     deviance = quantile(x$BUGSoutput$sims.matrix[, "deviance"], probs = c(0.025, 0.5, 0.975)))	
 }
 
-summarise.jags.output.table <- function(f, 	probs = c(0.05, 0.1, 0.5, 0.9, 0.95) ){
+summarise.jags.output.list <- function(f, 	probs = c(0.05, 0.1, 0.5, 0.9, 0.95) ){
 	list(file.names= files.details(f),
 	     list.summary= summarise.jags.output(read.jags.output(f)))
 }
 
  
+
 files.details <- function(x){
-	s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
-	names(s)  <- c("d1", "d2", "model", "trait", "set", "ecocode", "file" )
-	s
+	s <- data.frame(t(strsplit(x, "/", fixed = TRUE)[[1]]),
+                        stringsAsFactors =  FALSE)
+	names(s)  <- c("d1", "d2", "data.type", "trait",
+                        "model", "file" )
+	s[-(1:2)]
+}
+
+
+format.all.output.jags <- function(file.name = "wwf.results.jags.rds",
+                               list.file.name,
+                               models,
+                               traits = c("SLA", "Wood.density", "Max.height"),
+                               data.type = "simple"){
+files <-  c()
+for (trait in traits){
+    for(model in models){
+       source(model, local = TRUE)
+       model.obj <- load.model()
+       pathout <- output.dir('jags', model.obj$name, trait, data.type)
+       files <- c(files,file.path(pathout,file.name))
+      }
+    }
+
+out <- lapply(files, summarise.jags.output.list)
+names(out) <- lapply(lapply(files,files.details),
+                     function(x) paste(as.vector(x[names(x) != 'file']),
+                                       collapse="_"))
+saveRDS(out,file=file.path('output' , list.file.name))
+rm(out)
+gc()
 }
 
 

R/analysis/lmer.run.R

@@ -77,7 +77,7 @@ call.lmer.and.save <- function(formula, df.lmer, path.out,
 
 run.lmer <- function (model.file, trait,
                       min.obs = 10, sample.size = NA,
-                      ecocode.var = 'wwf',
+                      ecocode.var = 'koppen',
                       data.type ='simple',
                       var.sample = NA,
                       select.biome = NA,
@@ -153,7 +153,7 @@ load.data.for.lmer <- function(trait, data.type,
                                base.dir = "output/processed",
                                sample.size. ,
                                ecocode.var.,
-                               var.sample. = 'wwf',
+                               var.sample. = 'koppen',
                                select.biome. = NA,
                                select.set. = NA,
                                sample.vec.TF. = FALSE,
@@ -330,7 +330,7 @@ return(as.vector(scale(x,
 
 
 get.variables <- function(data.tree, BATOT, CWM.tn,
-                          abs.CWM.tntf, tf, ecocode.var = 'wwf',
+                          abs.CWM.tntf, tf, ecocode.var = 'koppen',
                          min.BA.G = 40, min.MAT = 10){
 logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
 logD <- scale.d(log(data.tree[["D"]]))
@@ -370,7 +370,7 @@ return(data.frame(logG =  logG,
 
 
 get.variables.intra <- function(data.tree, BATOT, CWM.tn,
-                          abs.CWM.tntf, tf, ecocode.var = 'wwf',
+                          abs.CWM.tntf, tf, ecocode.var = 'koppen',
                          min.BA.G = 40, min.MAT = 10){
 logG <- scale.d(log(data.tree[["BA.G"]] + min.BA.G+1))
 logD <- scale.d(log(data.tree[["D"]]))
@@ -416,7 +416,7 @@ return(data.frame(logG =  logG,
 # Function to prepare data for lmer with CWM data
 #============================================================
 format.data.for.lmer <-  function(data.tree, trait, min.obs = 10,
-                                  data.type = 'simple', ecocode.var = 'wwf') {
+                                  data.type = 'simple', ecocode.var = 'koppen') {
 if(! trait %in%  c("SLA", "Leaf.N", "Seed.mass",
                    "SLA", "Wood.density", "Max.height"))
     stop("need trait to be in SLA Leaf.N Seed.mass

docs/paper/SupplMat.md

@@ -20,11 +20,6 @@ and $N_i$ is the number of species in the local neighbourhood of the tree $i$.
 
 Two main data type were used: national forest inventories data -- NFI, large permanent plots data -- LPP.
 
-
-```
-## Warning: package 'plyr' was built under R version 3.1.3
-```
-
 ### Panama
 
 - Data set name: Panama
@@ -259,6 +254,14 @@ Two main data type were used: national forest inventories data -- NFI, large per
 	- Kooyman, R.M. and Westoby, M. (2009) Costs of height gain in rainforest saplings: main stem scaling, functional traits and strategy variation across 75 species. Annals of Botany 104: 987-993.
 	- Kooyman, R.M., Rossetto, M., Allen, C. and Cornwell, W. (2012) Australian tropical and sub-tropical rainforest: phylogeny, functional biogeography and environmental gradients. Biotropica 44: 668-679.
 
+## Age distribution of some forest data analysed.
+
+
+![Age distribution of forest area in 20-year age class for France, Switzerland and Sweden, estimated by Vilen et al.[@Vilen-2012]. The last class plotted at 150 years is for age > 140 years (except for Sweden where the last class 110 is age > 100 years).](../../figs/age_europe.pdf)
+
+![Age distribution of forest area in 20-year age class for North America (USA and Canada), estimated by Pan et al.[@Pan-2011]. The last class plotted at 150 years is for age > 140 years.](../../figs/age_na.pdf)
+
+
 # Supplementary discussion
 
 ## Trait effects and potential mechanisms

docs/paper/Suppl_Mat.Rmd

@@ -53,6 +53,14 @@ list.t <- dlply(dat, 1, paste_name_data)
 writeLines(unlist(list.t[dat[["Country"]]]))
 ```
 
+## Age distribution of some forest data analysed.
+
+
+![Age distribution of forest area in 20-year age class for France, Switzerland and Sweden, estimated by Vilen et al.[@Vilen-2012]. The last class plotted at 150 years is for age > 140 years (except for Sweden where the last class 110 is age > 100 years).](../../figs/age_europe.pdf)
+
+![Age distribution of forest area in 20-year age class for North America (USA and Canada), estimated by Pan et al.[@Pan-2011]. The last class plotted at 150 years is for age > 140 years.](../../figs/age_na.pdf)
+
+
 # Supplementary discussion
 
 ## Trait effects and potential mechanisms

launch.cluster/launch_all_lmer.bash

@@ -64,11 +64,11 @@ samplesize=$1
 
 
 # # # # ecocode 3
-#   	echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO${trait}.sh
-# 	qsub Rscript_temp/allECO${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO${trait}" -q opt32G -j oe
+  	echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO${trait}.sh
+	qsub Rscript_temp/allECO${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO${trait}" -q opt32G -j oe
 
-#   	echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[2], run.lmer,merge.biomes.TF = TRUE,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO2${trait}.sh
-# 	qsub Rscript_temp/allECO2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO2${trait}" -q opt32G -j oe
+  	echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[2], run.lmer,merge.biomes.TF = TRUE,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO2${trait}.sh
+	qsub Rscript_temp/allECO2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO2${trait}" -q opt32G -j oe
 
 # # # ecocode TP
   	echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.4[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECOTP${trait}.sh

scripts/analysis/jags.all.output.R

+# save all output in a list
+source("R/analysis/jags.output-fun.R")
+source("R/analysis/jags.run.R")
+
+
+
+## ## save list of all output for  NA
+
+format.all.output.jags(file.name = "wwf.results.jags.rds",
+                   list.file.name = 'list.jags.out.all.NA.simple.set.sample.rds',
+                   models = c(model.files.jags.Tf.1[1]),
+                   traits = c("SLA", "Wood.density", "Max.height"))
+
+

scripts/analysis/results.fig.R

@@ -155,9 +155,12 @@ plot.param.mean.and.biomes.fixed(list.all.results.ecocode ,
                xlim = c(-0.28, 0.27))
 dev.off()
 
-pdf('figs/figres1.ecocode.pdf', height = 14, width = 16)
-plot.param(list.all.results.set , data.type = 'simple',
-               model = c('lmer.LOGLIN.ER.AD.Tf.r.ecocode.species'),
+
+pdf('figs/figres12.ecocode.TP.pdf', height = 14, width = 16)
+plot.param.mean.and.biomes.fixed(list.all.results.ecocode ,
+                                 data.type = "simple",
+               models = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species',
+                   'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.fixed.biomes.species'),
                traits = c('Wood.density' , 'SLA', 'Max.height'),
                param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
                    "sumBn", "Tf"),
@@ -170,11 +173,13 @@ plot.param(list.all.results.set , data.type = 'simple',
                               expression("Size "(gamma %*% log('D'))) ) ,
                col.vec =  fun.col.pch.biomes()$col.vec,
                pch.vec =  fun.col.pch.biomes()$pch.vec,
-               xlim = c(-0.25, 0.25))
+               names.bio = names.biomes ,
+               xlim = c(-0.28, 0.27))
 dev.off()
 
 
 
+
 ### PLOT TRADE OFF param
 
 pdf('figs/figres4.pdf', width = 8, height = 7)

scripts/format.data/age.europe.R

@@ -27,11 +27,10 @@ theme_simple <-  function(){
   theme(axis.line = element_line(color = 'black'))
 }
 
-pdf('figs/age.europe.pdf')
+pdf('figs/age_europe.pdf')
 ggplot(df, aes(x = age.c, y = density)) +
     geom_bar(stat = "identity", position = "stack",
              colour = "black",  fill="grey") +
-#    scale_x_log10(limits=c(1, 1000)) +    
-    theme_simple() + facet_wrap(~ country)+ xlab("age (yr.)")
+    theme_simple()  + facet_wrap(~ country)+ xlab("age (yr.)")
 dev.off()
 

scripts/format.data/age.na.R

@@ -9,6 +9,16 @@ data.age <- as.data.frame(age_raster, na.rm = TRUE)
 vec.age <- data.age[!is.na(data.age[,1]),]
 df.age <- data.frame(age= vec.age)
 
+library(dplyr)
+
+df <-  mutate(df.age,
+              age.class = cut(age, breaks = c(seq(0,140, by = 20), 860),
+                              include.lowest = TRUE))
+levels(df$age.class) <- seq(10,150, by = 20)
+N.age <- table(df$age.class)
+P.age <- as.vector(N.age)/sum(N.age)
+
+df.age.class <- data.frame(age.m = seq(10,150, by = 20), density = P.age)
 
 # plot age dist
 library(ggplot2)
@@ -30,8 +40,9 @@ theme_simple <-  function(){
   theme(axis.line = element_line(color = 'black'))
 }
 
-pdf('figs/age.na.pdf')
-ggplot(df.age, aes(x = age)) + geom_histogram(aes(y = ..density..), colour = "black",  fill="grey") +
-    geom_density(, col = 'red')+  scale_x_log10(limits=c(1, 1000))  +
-        theme_simple() + ggtitle("North America (Pan et al . 2011)") + xlab("age (yr.)")
+pdf('figs/age_na.pdf')
+ggplot(df.age.class, aes(x = age.m, y = density)) +
+    geom_bar(stat = "identity", position = "stack",
+             colour = "black",  fill="grey") +
+    theme_simple() + xlab("age (yr.)") +  ggtitle("North America")
 dev.off()