Skip to content

GitLab

Commit b5f8f7c4 authored by kunstler's avatar kunstler
Browse files

add TP with intra

parent 87f8641f
Hide whitespace changes
Inline Side-by-side
R/analysis/lmer.output-fun.R
View file @ b5f8f7c4
......@@ -656,10 +656,7 @@ plot.param <- function(list.res,
param.print = 1:5,
traits.print = 1:3,
col.names = fun.col.param(),
data.type = "all.no.log",
col.vec,
pch.vec,
names.bio,
data.type = "simple",
add.param.descrip.TF = TRUE,
...){
m <- matrix(c(1:3), 1, 3)
......@@ -673,6 +670,8 @@ for (i in traits[traits.print]){
"_", model,
sep = '')]]$lmer.summary
param.mean <- list.temp$fixed.coeff.E[param.vec]
param.mean[!names(param.mean) %in% c("Tf", "sumTfBn")] <-
-param.mean[!names(param.mean) %in% c("Tf", "sumTfBn")]
param.std <- list.temp$fixed.coeff.Std.Error
names(param.std) <- names(list.temp$fixed.coeff.E)
param.std <- param.std[param.vec]
......@@ -688,10 +687,10 @@ for (i in traits[traits.print]){
plot(param.mean[param.print], (1:length(param.vec))[param.print],
yaxt = 'n', xlab = NA, ylab = NA,
pch = 16, cex = 2, cex.lab = 1.7, cex.axis = 1.5,
ylim = range(1:length(param.vec)), ...)
ylim = range(1-0.21, length(param.vec)+0.21), ...)
mtext(traits.names[i], side=3, cex =1.7, line = 1)
box(lwd= 2)
lines(c(0, 0), c(1-0.2, length(param.vec)+0.2))
lines(c(0, 0), c(par()$usr[3], par()$usr[4]))
if(i == traits[1]) {lapply(1:length(param.vec),
fun.axis.one.by.one,
side = 2,
......@@ -699,12 +698,12 @@ for (i in traits[traits.print]){
cols.vec = col.names[param.vec])
if(add.param.descrip.TF){
mtext("Max growth", side=2, at = 4.95, cex =1.6,
line = 17.9, col = '#e41a1c')
lines(c(-0.71, -0.71), c(4.5, 5.28), col = '#e41a1c',
line = 16.9, col = '#e41a1c')
lines(c(-0.5, -0.5), c(4.5, 5.28), col = '#e41a1c',
lwd = 2.5)
mtext("Competition", side=2, at = 2.5, cex =1.6,
line = 17.9, col = '#377eb8')
lines(c(-0.71, -0.71), c(0.82, 4.28), col = '#377eb8',
line = 16.9, col = '#377eb8')
lines(c(-0.5, -0.5), c(0.82, 4.28), col = '#377eb8',
lwd = 2.5)
}
}
......@@ -954,27 +953,18 @@ fun.layout <- function(b = border.size()){
layout(m2, widths=c(wid, b$legend.m), heights = 1)
}
fun.param.descrip <- function(add.param.descrip.TF, x.line = -0.63){
if(add.param.descrip.TF == 2){
mtext("Max growth", side=2, at = 4.95, cex =1.6,
line = 17.9, col = '#e41a1c')
lines(c(x.line, x.line), c(4.55, 5.33), col = '#e41a1c',
lwd = 2.5)
mtext("Competition", side=2, at = 2.5, cex =1.6,
line = 17.9, col = '#377eb8')
lines(c(x.line, x.line), c(0.82, 4.25), col = '#377eb8',
lwd = 2.5)
}
if(add.param.descrip.TF == 1){
mtext("Max growth", side=2, at = 4, cex =1.6,
fun.param.descrip <- function(seq.jitter, n.param, x.line = -0.63){
mtext("Max growth", side=2, at = n.param, cex =1.6,
line = 17.9, col = '#e41a1c')
lines(c(x.line, x.line), c(3.7, 4.3), col = '#e41a1c',
lines(c(x.line, x.line),
c(n.param - abs(min(seq.jitter))-0.15, n.param + abs(min(seq.jitter))+0.15),
col = '#e41a1c',
lwd = 2.5)
mtext("Competition", side=2, at = 2, cex =1.6,
mtext("Competition", side=2, at = (n.param - 1)/2 + max(seq.jitter)+abs(min(seq.jitter)), cex =1.6,
line = 17.9, col = '#377eb8')
lines(c(x.line, x.line), c(0.72, 3.28), col = '#377eb8',
lines(c(x.line, x.line),
c(1 - abs(min(seq.jitter)-0.15), n.param-1 + abs(min(seq.jitter))+0.15), col = '#377eb8',
lwd = 2.5)
}
}
fun.legend <- function(biomes.names, biomes.c, col.vec, pch.vec){
......@@ -1027,16 +1017,15 @@ plot.param.mean.and.biomes.fixed <- function(list.res,
'Compet x trait dissimilarity'),
param.print = 1:5,
col.names = fun.col.param() ,
data.type = "all.no.log",
data.type = "simple",
col.vec,
pch.vec,
names.bio,
add.param.descrip.TF = 1,
...){
if(! add.param.descrip.TF %in% 0:2) stop("add.param.descrip.TF need to be in 0 1 2")
...){
col.vec[2] <- col.vec[1]
biomes.c <- as.character(biomes)
Var <- "Trait indep"
intra <- "intra"
fun.layout()
b <- border.size()
##################################
......@@ -1048,8 +1037,8 @@ b <- border.size()
"_", models[1],
sep = '')]]$lmer.summary
param.mean.m <- list.temp.m$fixed.coeff.E[param.vec]
param.mean.m[c("sumTnTfBn.abs", "sumTnBn", "sumBn")] <-
-param.mean.m[c("sumTnTfBn.abs", "sumTnBn", "sumBn")]
param.mean.m[!names(param.mean.m) %in% c("Tf", "sumTfBn")] <-
-param.mean.m[!names(param.mean.m) %in% c("Tf", "sumTfBn")]
param.std.m <- list.temp.m$fixed.coeff.Std.Error
names(param.std.m) <- names(list.temp.m$fixed.coeff.E)
param.std.m <- param.std.m[param.vec]
......@@ -1060,7 +1049,7 @@ b <- border.size()
list.fixed <- fun.get.fixed.biomes(param.vec[n.vars], list.temp,
biomes.vec = biomes)
param.mean <- list.fixed$fixed.biomes
if (param.vec[n.vars] %in% c("sumTnTfBn.abs", "sumTnBn", "sumBn")) param.mean <- -param.mean
if (!param.vec[n.vars] %in% c("Tf", "sumTfBn")) param.mean <- -param.mean
param.std <- list.fixed$fixed.biomes.std
seq.jitter <- seq(25, -25, length.out = length(biomes)+1)/120
if(n.vars == 1){
......@@ -1072,14 +1061,14 @@ b <- border.size()
col = c('black', col.vec[biomes.c]), ...)
mtext('', side=3, cex =1.7, line = 1)
box(lwd= 2)
lines(c(0, 0), c(1-0.38, length(param.vec)+0.38), lwd= 1)
lines(c(0, 0), c(par()$usr[3], par()$usr[4]), lwd= 1)
if(i == traits[1]) {
lapply(1:length(param.vec),
fun.axis.one.by.one,
side = 2,
labels = param.names,
cols.vec = col.names[param.vec])
fun.param.descrip(add.param.descrip.TF)
fun.param.descrip(seq.jitter,length(param.vec))
}
if(i == traits[2]){ mtext('Standardized coefficients', side=1, cex =1.5,
line = 4)}
......
R/analysis/lmer.run.R
View file @ b5f8f7c4
......@@ -34,13 +34,17 @@ model.files.lmer.Tf.1 <- file.path(path.model,
c("model.lmer.LOGLIN.ER.AD.Tf.r.set.species.R",
"model.lmer.LOGLIN.ER.AD.Tf.r.set.fixed.biomes.species.R"))
model.files.lmer.Tf.2 <- 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"))
c("model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R",
"model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.fixed.biomes.species.R"))
model.files.lmer.Tf.3 <- file.path(path.model,
c("model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.species.R"))
model.files.lmer.Tf.4 <- 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.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"))
model.files.lmer.Tf.intra.2 <- file.path(path.model,
c("model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species.R",
"model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species.R"))
......
R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species.R 0 → 100644
View file @ b5f8f7c4
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+biomes.id+MAT+MAP+Tf+Tf:biomes.id+logD+sumBn.inter+sumBn.inter:biomes.id+sumBn.intra+sumBn.intra:biomes.id+sumTfBn+sumTfBn:biomes.id+sumTnBn+sumTnBn:biomes.id+sumTnTfBn.abs+sumTnTfBn.abs:biomes.id +(1+logD+sumBn.inter+sumBn.intra||species.id)+(1|plot.id)+(1+Tf+sumBn.inter+sumBn.intra+sumTfBn+sumTnBn+sumTnTfBn.abs||set.id)"))
}
R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species.R 0 → 100644
View file @ b5f8f7c4
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+Tf+logD+sumBn.inter+MAT+MAP+sumBn.intra+sumTfBn+sumTnBn+sumTnTfBn.abs +(1+logD+sumBn.inter+sumBn.intra||species.id)+(1|plot.id)+(1+Tf+sumBn.inter+sumBn.intra+sumTfBn+sumTnBn+sumTnTfBn.abs||set.id)"))
}
R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.r.ecocode.fixed.biomes.species.R 0 → 100644
View file @ b5f8f7c4
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.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+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)+(sumBn-1|ecocode.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
View file @ b5f8f7c4
......@@ -38,27 +38,38 @@ samplesize=$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[2], merge.biomes.TF = TRUE, run.lmer,$trait);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
# # 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.4[1], run.lmer,$trait);print('done')\"" > Rscript_temp/allTP${trait}.sh
qsub Rscript_temp/allTP${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.TP${trait}" -q opt32G -j oe
# # # 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.3[1], run.lmer,$trait);print('done')\"" > Rscript_temp/allTP${trait}.sh
# qsub Rscript_temp/allTP${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.TP${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], merge.biomes.TF = TRUE, run.lmer,$trait);print('done')\"" > Rscript_temp/allTP2${trait}.sh
qsub Rscript_temp/allTP2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.TP2${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.3[2], merge.biomes.TF = TRUE, run.lmer,$trait);print('done')\"" > Rscript_temp/allTP2${trait}.sh
# qsub Rscript_temp/allTP2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.TP2${trait}" -q opt32G -j oe
# # # INTRA 2
# 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 = 'intra');print('done')\"" > Rscript_temp/allINTRA${trait}.sh
# # # INTRA 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.intra.1[1], run.lmer,$trait, data.type = 'intra');print('done')\"" > Rscript_temp/allINTRA${trait}.sh
# qsub Rscript_temp/allINTRA${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRA${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], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRA2${trait}.sh
# 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.intra.1[2], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRA2${trait}.sh
# qsub Rscript_temp/allINTRA2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRA2${trait}" -q opt32G -j oe
# # INTRA 2
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.intra.2[1], run.lmer,$trait, data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB${trait}.sh
qsub Rscript_temp/allINTRAB${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB${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.intra.2[2], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB2${trait}.sh
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.3[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO${trait}.sh
# 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,$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
done
scripts/analysis/lmer.all.no.log.output.R
View file @ b5f8f7c4
......@@ -9,7 +9,7 @@ source("R/analysis/lmer.run.R")
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.simple.set.rds',
models = c(model.files.lmer.Tf.0, model.files.lmer.Tf.1, model.files.lmer.Tf.2,
models = c(model.files.lmer.Tf.0, model.files.lmer.Tf.1,
model.files.lmer.Tf.3, model.files.lmer.Tf.4),
traits = c("SLA", "Wood.density", "Max.height"))
......
scripts/analysis/results.fig.R
View file @ b5f8f7c4
......@@ -16,12 +16,21 @@ 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.intra <-
readRDS('output/list.lmer.out.all.NA.intra.set.rds')
# check convergence
fun.get.conv(list.all.results.set)
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.intra)
vec.rel.grad.set <- sapply(list.all.results.intra,
function(list.t) { max(abs(list.t[['relgrad']]))})
vec.rel.grad.set < 0.001
## plots
pdf('figs/figres12.pdf', height = 14, width = 16)
......@@ -32,25 +41,80 @@ plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
param.print = 1:5,
param.names = c(expression(paste('Trait sim ',
(alpha[s]))),
expression(paste('Tolerance ',
(alpha[t] ))),
expression(paste('Effect ',
(alpha[e] ))),
expression(paste('Trait indep ',
(alpha[0]))),
expression(paste("Direct trait ",
(m[1]))),
expression(paste("Size ",
(gamma %*% log(D))) )) ,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha[0])),
expression("Direct trait "(m[1])),
expression("Size "(gamma %*% log('D'))) ) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.28, 0.27))
dev.off()
pdf('figs/figres12.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.intra.r.set.species',
'lmer.LOGLIN.ER.AD.Tf.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"),
param.print = 1:6,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha['0 intra'])),
expression('Trait indep'(alpha['0 inter'])),
expression("Direct trait "(m[1])),
expression("Size "(gamma %*% log('D'))) ) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.28, 0.27),
xlim = c(-0.28, 0.31),
add.param.descrip.TF =2)
dev.off()
pdf('figs/figres12.TP.pdf', height = 14, width = 16)
plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
models = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species',
'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.fixed.biomes.species'),
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
param.print = 1:5,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha[0])),
expression("Direct trait "(m[1])),
expression("Size "(gamma %*% log('D'))) ) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
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.MAT.MAP.r.set.species'),
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
param.print = 1:5,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha[0])),
expression("Direct trait "(m[1])),
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))
dev.off()
### PLOT TRADE OFF param
......
scripts/format.data/age.europe.R 0 → 100644
View file @ b5f8f7c4
################################ READ DATA
data.europe <- read.csv("data/raw/Europe_Age.csv",
stringsAsFactors = FALSE)
library(tidyr)
df <- gather(data.europe, country, density, -age.c)
library(dplyr)
# 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.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.)")
dev.off()
x11()
ggplot(data.age, aes(x = age.max)) + geom_histogram(aes(y = ..density..), colour = "black", fill="grey") + geom_density(, col = 'red')+ scale_x_continuous(limits=c(1, 300))+ theme_simple() + ggtitle("France") + xlab("age (yr.)")
scripts/format.data/age.france.R 0 → 100644
View file @ b5f8f7c4
################################ READ DATA
data.france <- read.csv("data/raw/France/dataIFN.FRANCE.csv",
stringsAsFactors = FALSE)
library(dplyr)
data.age <- data.france %>% group_by(idp) %>% summarise(age.max = max(age, 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.france.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("France") + xlab("age (yr.)")
dev.off()
scripts/format.data/age.na.R 0 → 100644
View file @ b5f8f7c4
## READ AND PLOT HIST OF AGE DATA FROM PAN et al. 2011
library(raster)
age_raster=raster("data/raw/US/ca04_usak06_1km.tif")
plot(age_raster)
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)
# 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.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.)")
dev.off()
scripts/format.data/age.quebec.R 0 → 100644
View file @ b5f8f7c4
################################ READ DATA
data.quebec <- read.csv("data/raw/Canada/data.quebec.csv",
stringsAsFactors = FALSE)
library(dplyr)
data.age <- data.quebec %>% group_by(ID_PEP) %>% summarise(age.max = max(AGE.1, 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.quebec.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("Quebec")+ xlab("age (yr.)")
dev.off()
scripts/format.data/age.spain.R 0 → 100644
View file @ b5f8f7c4
################################ READ DATA
data.spain <- read.delim("data/raw/Spain/Export_ageplots.txt",
stringsAsFactors = FALSE)
library(dplyr)
data.age <- data.spain %>% group_by(PlotCode) %>%
summarise(age.max = max(Age1,Age2,Age3, na.rm = TRUE)) %>%
ungroup() %>% filter(age.max>0)
# 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.spain.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("France") + xlab("age (yr.)")
dev.off()
scripts/format.data/age.sweden.R 0 → 100644
View file @ b5f8f7c4
### sweden DATA
data.swe <- read.table("data/raw/Sweden/Swe_NFI_all.txt", header = T,
stringsAsFactors = F, sep = "\t")
data.swe$tree.id <- apply(cbind(data.swe[["Year"]], data.swe[["TractNr"]],
data.swe[["PlotNr"]],
data.swe[["TreeNr"]]),