Commit a037c86f authored by kunstler's avatar kunstler
Browse files

progress figure and process new results cluster

parent eba0a4c6
......@@ -193,13 +193,7 @@ return(res)
###########################
### FUNCTIOMNS TO PLOTS RESULTS
fun.plot.error.bar <- function(x, y, sd, small.bar = (max(x) - min(x))/40){
segments( x, unlist(y - 1.96*sd), x , unlist(y +1.96*sd))
segments( x - small.bar, unlist(y - 1.96*sd), x + small.bar , unlist(y -1.96*sd))
segments( x - small.bar, unlist(y + 1.96*sd), x + small.bar, unlist(y +1.96*sd))
}
### FUNCTIONS TO PLOTS RESULTS
fun.plot.error.bar.horiz <- function(x, y, sd, small.bar = (max(x) - min(x))/40, ...){
segments( unlist(x - 1.96*sd), y, unlist(x +1.96*sd), y, ...)
......@@ -216,6 +210,13 @@ return(t.col)
}
fun.traits.names <- function(){
c(Wood.density= 'Wood density',
SLA = 'Specific leaf area',
Max.height = 'Maximum height',
Seed.size = 'Seed size',
Leaf.N = 'Leaf N content')
}
......@@ -461,9 +462,15 @@ for (i in traits){
}
dat.res <- data.frame(do.call('rbind', list.df),
traits = rep(traits, each = 2*N.pred))
dat.res <- change.traits.levels(dat.res)
return(dat.res)
}
change.traits.levels <- function(df,
traits.names = fun.traits.names()){
levels(df$traits) <- traits.names[levels(df$traits)]
return(df)
}
theme_simple <- function(){
theme_bw() +
......@@ -473,7 +480,8 @@ theme_simple <- function(){
panel.grid.minor = element_blank(),
legend.title=element_blank(),
strip.background = element_blank(),
legend.key = element_blank()
legend.key = element_blank(),
legend.position=c(.3,.07)
) +
#draws x and y axis line
theme(axis.line = element_line(color = 'black'))
......@@ -517,14 +525,11 @@ fun.axis.one.by.one <- function(x, side, labels, cols.vec , cex.axis = 2){
}
plot.param.BLUP <- function(list.res,
model = 'lmer.LOGLIN.ER.AD.Tf',
data.type = "all.no.log",
traits = c('Wood.density', 'SLA', 'Max.height'),
traits.names = c(Wood.density= 'Wood density',
SLA = 'Specific leaf area',
Max.height = 'Maximum height'),
traits.names = fun.traits.names() ,
param.vec = c("Tf","sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs"),
param.names = c('Direct trait',
......@@ -766,20 +771,20 @@ for (i in traits){
if(i == traits[length(traits)]) {
par(mai=c(1.2,0,0.6,small.m))
}
seq.jitter <- seq(-15,15, length.out = length(biomes))/120
seq.jitter <- seq(25, -25, length.out = length(biomes))/120
if(n.vars == 1){
plot(param.mean, seq.jitter+n.vars,
yaxt = 'n', xlab = NA, ylab = NA,
, pch = pch.vec[biomes.c] , cex = 2, cex.axis = 1.5, cex.lab = 1.5,
ylim = range(1-0.15, length(param.vec)+0.15),
col = col.vec[biomes.c], ...)
if(i == traits[2]) mtext('Standardized coefficients', side=1, cex =1.4,
if(i == traits[2]) mtext('Standardized coefficients', side=1, cex =1.5,
line = 4)
}
if(n.vars != 1){
points(param.mean, seq.jitter+n.vars,
col = col.vec[biomes.c], pch = pch.vec[biomes.c], cex = 2)}
mtext(traits.names[i], side=3, cex =1.4, line = 1)
mtext(traits.names[i], side=3, cex =1.5, line = 1)
box(lwd= 2)
abline(v = 0)
if(i == traits[1]) {lapply(1:length(param.vec),
......@@ -787,7 +792,7 @@ for (i in traits){
side = 2,
labels = param.names,
cols.vec = col.names[param.vec],
cex.axis = 1.6)
cex.axis = 2)
}
fun.plot.error.bar.horiz(param.mean,
seq.jitter+n.vars,
......
......@@ -8,6 +8,13 @@ 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.3[-2]),
traits = c("SLA", "Wood.density", "Max.height")
)
format.all.output.lmer(file.name = "biomes.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.biomes.simple.set.rds',
models = c(model.files.lmer.Tf.3),
traits = c("SLA", "Wood.density", "Max.height")
)
......@@ -99,7 +99,7 @@ plot.param(list.all.results.set , data.type = "simple",
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.45, 0.51))
xlim = c(-0.45, 0.45))
dev.off()
pdf('figs/figres1b.pdf', height = 7, width = 14)
......@@ -216,7 +216,7 @@ plot.param.BLUP(list.all.results.species.id.set , data.type = "simple",
dev.off()
pdf('figs/figres2.pdf', height = 7, width = 14)
pdf('figs/figres2.pdf', height = 7, width = 16)
plot.param.biomes.fixed(list.all.results.set, data.type = 'simple',
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.fixed.biomes.species',
param.vec = c("sumTnTfBn.abs", "sumTnBn", "sumTfBn",
......@@ -232,7 +232,7 @@ plot.param.biomes.fixed(list.all.results.set, data.type = 'simple',
param.print = 1:4,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
xlim = c(-0.4, 0.3))
xlim = c(-0.35, 0.29))
dev.off()
# logD
......@@ -240,12 +240,12 @@ dev.off()
### PLOT TRADE OFF
pdf('figs/figres3.pdf', width = 7, height = 5.5)
pdf('figs/figres3.pdf', width = 3.5, height = 9)
plot.growth.ba(traits = c('Wood.density', 'SLA', 'Max.height'),
model = 'lmer.LOGLIN.ER.AD.Tf.r.set.species',
type = 'Tabs', data.type = 'simple',
dir.root = '.',
list.res = list.all.results.species.id.set,
list.res = list.all.results.set,
labels.leg = c('High trait value', 'Low trait value'))
dev.off()
......
## REMOVE OUTLIERS
source("R/process.data/process-fun.R")
source("R/process.data/test.tree.CWM-fun.R")
source("R/utils/plot.R")
filedir <- "output/processed"
## remove outlier following Condit approach and reformat data
data.all <- readRDS( 'output/processed/data.all.global.t.rds')
## compute mean and sd for each traits
library(dplyr)
data.sd.traits <- summarise(data.all,
sd.SLA = sd(SLA.focal, na.rm = TRUE),
sd.Leaf.N = sd(Leaf.N.focal, na.rm = TRUE),
sd.Wood.density = sd(Wood.density.focal, na.rm = TRUE),
sd.Max.height = sd(Max.height.focal, na.rm = TRUE),
sd.Seed.mass = sd(Seed.mass.focal, na.rm = TRUE),
mean.SLA = mean(SLA.focal, na.rm = TRUE),
mean.Leaf.N = mean(Leaf.N.focal, na.rm = TRUE),
mean.Wood.density = mean(Wood.density.focal, na.rm = TRUE),
mean.Max.height = mean(Max.height.focal, na.rm = TRUE),
mean.Seed.mass = mean(Seed.mass.focal, na.rm = TRUE)
)
saveRDS(data.sd.traits, 'output/processed/data.sd.traits.global.rds')
data.all1 <- fun.reform.data.and.remove.outlier(data.all, err.limit = 4,
select.one.census.rand = TRUE,
std.traits.TF = FALSE)
saveRDS(data.all1, 'output/processed/data.all.global.rds')
rm(data.all1)
data.all1 <- fun.reform.data.and.remove.outlier(data.all, err.limit = 4,
select.one.census.rand = FALSE,
std.traits.TF = FALSE)
saveRDS(data.all1, 'output/processed/data.all.global.all.census.rds')
rm(data.all1)
gc()
####
## remove outlier following Condit approach and reformat data
data.all <- readRDS( 'output/processed/data.all.no.log.t.rds')
data.sd.traits <- summarise(data.all,
sd.SLA = sd(SLA.focal, na.rm = TRUE),
sd.Leaf.N = sd(Leaf.N.focal, na.rm = TRUE),
sd.Wood.density = sd(Wood.density.focal, na.rm = TRUE),
sd.Max.height = sd(Max.height.focal, na.rm = TRUE),
sd.Seed.mass = sd(Seed.mass.focal, na.rm = TRUE),
mean.SLA = mean(SLA.focal, na.rm = TRUE),
mean.Leaf.N = mean(Leaf.N.focal, na.rm = TRUE),
mean.Wood.density = mean(Wood.density.focal, na.rm = TRUE),
mean.Max.height = mean(Max.height.focal, na.rm = TRUE),
mean.Seed.mass = mean(Seed.mass.focal, na.rm = TRUE)
)
saveRDS(data.sd.traits, 'output/processed/data.sd.traits.no.log.rds')
data.all1 <- fun.reform.data.and.remove.outlier(data.all)
saveRDS(data.all1, 'output/processed/data.all.no.log.rds')
rm(data.all1)
data.all1 <- fun.reform.data.and.remove.outlier(data.all, select.one.census.rand = FALSE)
saveRDS(data.all1, 'output/processed/data.all.no.log.rds')
rm(data.all1)
gc()
cat("finished", file = file.path(filedir, "Done.txt"))
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