run ecocode TP

R/analysis/lmer.run.R

@@ -39,6 +39,9 @@ model.files.lmer.Tf.2 <- file.path(path.model,
 model.files.lmer.Tf.3 <- file.path(path.model,
       c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species.R",
         "model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.fixed.biomes.species.R"))
+model.files.lmer.Tf.4 <- file.path(path.model,
+      c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species.R",
+        "model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.fixed.biomes.species.R"))
 model.files.lmer.Tf.intra.1 <- file.path(path.model,
       c("model.lmer.LOGLIN.ER.AD.Tf.intra.r.set.species.R",
         "model.lmer.LOGLIN.ER.AD.Tf.intra.r.set.fixed.biomes.species.R"))

R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.fixed.biomes.species.R

+load.model <- function () {
+    list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.fixed.biomes.species",
+         var.BLUP = 'ecocode.id',
+         lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+biomes.id+Tf+MAT+MAP+Tf:biomes.id+logD+sumBn+sumBn:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"))
+}
+

R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species.R

+load.model <- function () {
+    list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.ecocode.species",
+         var.BLUP = 'ecocode.id',
+         lmer.formula.tree.id=formula("logG~1+(1|ecocode.id)+(1|species.id)+(1|plot.id)+Tf+MAP+MAT+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"))
+}
+
+
+

launch.cluster/launch_all_lmer.bash

@@ -63,12 +63,19 @@ samplesize=$1
 # 	qsub Rscript_temp/allINTRAB2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB2${trait}" -q opt32G -j oe
 
 
-# # # 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
+# # # # 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[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
+	qsub Rscript_temp/allECOTP${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECOTP${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.4[2], run.lmer,merge.biomes.TF = TRUE,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECOTP2${trait}.sh
+	qsub Rscript_temp/allECOTP2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECOTP2${trait}" -q opt32G -j oe
 
 
 done

launch.cluster/launch_all_lmer_sample.bash

+#!/bin/bash
+# Georges Kunstler 14/05/2014
+# read one variable
+export LD_LIBRARY_PATH=/usr/lib64/R/library
+
+mkdir -p Rscript_temp
+
+
+# test parameter
+nbargs=$#
+echo "number of arguments="$nbargs
+if [ $nbargs -ne 1 ]
+then
+  echo "need one and only one argument"
+  echo " usage :"
+  echo " ./launch_all_lmer.sh  sample.size"
+  exit 100
+fi
+
+
+samplesize=$1
+# "'Seed.mass'" "'Leaf.N'"
+ for trait in "'SLA'" "'Wood.density'" "'Max.height'" ; do
+
+
+# # ALL data 1
+  	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.1[1], run.lmer,$trait, sample.vec.TF = TRUE);print('done')\"" > Rscript_temp/allf${trait}.sh
+	qsub Rscript_temp/allf${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.f${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.1[2], merge.biomes.TF = TRUE, run.lmer,$trait, sample.vec.TF = TRUE);print('done')\"" > Rscript_temp/allf2${trait}.sh
+	qsub Rscript_temp/allf2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.f2${trait}" -q opt32G -j oe
+
+
+
+done
+

scripts/analysis/results.fig.R

@@ -16,6 +16,10 @@ names.biomes <- c('Desert', 'Desert', 'Woodland/shrubland', 'Temperate forest',
 list.all.results.set <-
     readRDS('output/list.lmer.out.all.NA.simple.set.rds')
 
+list.all.results.ecocode <-
+    readRDS('output/list.lmer.out.all.NA.simple.ecocode.rds')
+
+
 list.all.results.intra <-
     readRDS('output/list.lmer.out.all.NA.intra.set.rds')
 
@@ -26,6 +30,13 @@ 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.ecocode)
+vec.rel.grad.ecocode <- sapply(list.all.results.ecocode,
+                         function(list.t) { max(abs(list.t[['relgrad']]))})
+vec.rel.grad.ecocode < 0.001
+
+
 fun.get.conv(list.all.results.intra)
 vec.rel.grad.set <- sapply(list.all.results.intra,
                          function(list.t) { max(abs(list.t[['relgrad']]))})
@@ -102,7 +113,7 @@ dev.off()
 pdf('figs/figres12.TP.intra.pdf', height = 14, width = 16)
 plot.param.mean.and.biomes.fixed(list.all.results.intra , data.type = "intra",
                models = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species',
-                   'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species'),
+               'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species'),
                traits = c('Wood.density' , 'SLA', 'Max.height'),
                param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
                    "sumBn.intra", "sumBn.inter", "Tf"),
@@ -124,7 +135,8 @@ dev.off()
 
 
 pdf('figs/figres12.ecocode.pdf', height = 14, width = 16)
-plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
+plot.param.mean.and.biomes.fixed(list.all.results.ecocode ,
+                                 data.type = "simple",
                models = c('lmer.LOGLIN.ER.AD.Tf.r.ecocode.species',
                    'lmer.LOGLIN.ER.AD.Tf.r.ecocode.fixed.biomes.species'),
                traits = c('Wood.density' , 'SLA', 'Max.height'),

scripts/format.data/age.swiss.R

+######################### READ DATA read individuals tree data
+data.swiss1 <- read.csv("data/raw/Swiss/LFI12.csv",
+                        header = TRUE,
+                        stringsAsFactors = FALSE)
+data.swiss2 <- read.csv("data/raw/Swiss/LFI23.csv",
+                        header = TRUE, stringsAsFactors = FALSE)
+data.swiss3 <- read.csv("data/raw/Swiss/LFI34.csv",
+                        header = TRUE,
+                        stringsAsFactors = FALSE)
+data.swiss <- rbind(data.swiss1, data.swiss2, data.swiss3)
+rm(data.swiss1, data.swiss2, data.swiss3)
+data.swiss <- data.swiss[order(data.swiss$BANR), ]
+
+library(dplyr)
+
+data.age <- data.swiss %>% group_by(CLNR) %>% summarise(age.max = max(ALTERD1, na.rm = TRUE))
+
+
+
+# plot age dist
+library(ggplot2)
+theme_simple <-  function(){
+  theme_bw() +
+  theme(
+    plot.background = element_blank(),
+    panel.grid.major = element_blank(),
+    panel.grid.minor = element_blank(),
+    legend.title=element_blank(),
+    strip.background = element_blank(),
+    legend.key = element_blank(),
+    legend.position=c(.1,.15),
+    strip.text = element_text(size=14),
+   axis.title.x = element_text(size=14),
+   axis.title.y = element_text(size=14)
+  ) +
+  #draws x and y axis line
+  theme(axis.line = element_line(color = 'black'))
+}
+
+pdf('figs/age.swiss.pdf')
+ggplot(data.age, aes(x = age.max)) + geom_histogram(aes(y = ..density..), colour = "black",  fill="grey") +
+    geom_density(, col = 'red')+  scale_x_log10(limits=c(1, 1000))  +theme_simple() +
+    ggtitle("Swiss") + xlab("age (yr.)")
+dev.off()
+
+
+