Commit 4f06fce4 authored by kunstler's avatar kunstler
Browse files

add model intra ecocode

parent 2a0742ee
......@@ -23,7 +23,7 @@ run.model.for.set.one.trait <- function(model.file, fun.model, trait,
#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
# Run lmer() (in package lme4) for one trait, one model
#=====================================================================
path.model <- "R/analysis/model.lmer"
......@@ -33,18 +33,12 @@ model.files.lmer.Tf.0 <- file.path(path.model,
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.r.tree.id.set.species.R",
"model.lmer.LOGLIN.ER.AD.Tf.r.tree.id.set.fixed.biomes.species.R"))
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.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"))
c("model.lmer.LOGLIN.ER.AD.Tf.intra.r.ecocode.species.R"))
model.files.lmer.Tf.Multi <- file.path(path.model,
c("model.lmer.LOGLIN.ER.AD.Tf.Multi.r.set.species.R",
"model.lmer.LOGLIN.ER.AD.Tf.Multi.r.set.fixed.biomes.species.R"))
call.lmer.and.save <- function(formula, df.lmer, path.out,
......@@ -56,7 +50,6 @@ call.lmer.and.save <- function(formula, df.lmer, path.out,
relgrad <- with(lmer.output@optinfo$derivs,solve(Hessian,gradient))
print('test convergence of relative gradient of hessian')
print(max(abs(relgrad)) )
## conf.param <- confint(lmer.output)
if(is.na(select.biome)){
name.file <- paste(var.sample,ecocode.var,
"results.nolog.all.rds", sep ='.')
......@@ -74,7 +67,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', Multi.type = 'a',
ecocode.var = 'wwf',
data.type ='simple',
var.sample = NA,
select.biome = NA,
......@@ -93,7 +86,6 @@ run.lmer <- function (model.file, trait,
data.type = data.type,
sample.size. = sample.size,
ecocode.var. = ecocode.var,
Multi.type. = Multi.type,
var.sample. = var.sample,
select.biome. = select.biome,
merge.biomes.TF = merge.biomes.TF)
......@@ -147,22 +139,17 @@ load.data.for.lmer <- function(trait, data.type,
base.dir = "output/processed",
sample.size. ,
ecocode.var.,
Multi.type.,
var.sample. = 'ecocode',
var.sample. = 'wwf',
select.biome. = NA,
select.set. = NA,
sample.vec.TF. = FALSE,
merge.biomes.TF = FALSE){
require(dplyr)
if (! data.type %in% c('simple', 'Multi', 'all.census')) stop('not good data.type')
if (data.type == 'Multi'){
df <- readRDS(file.path(base.dir,paste('data', 'Multi', data.type, 'rds',
sep = '.')))
}else{
if (! data.type %in% c('simple', 'intra', 'all.census'))
stop('not good data.type')
df <- readRDS( file.path(base.dir,paste('data', trait, data.type, 'rds',
sep = '.')))
}
df <- mutate(df,sp.name = ifelse(is.na(sp.name), 'missing.sp',
df <- dplyr::mutate(df,sp.name = ifelse(is.na(sp.name), 'missing.sp',
sp.name))
if(merge.biomes.TF) df <- merge.biomes(df)
......@@ -185,9 +172,11 @@ load.data.for.lmer <- function(trait, data.type,
sep = '.')))
print(paste('sub-sampled by',var.sample.))
}else{
sample.vec <- readRDS(file = file.path(base.dir,paste('sample.vec', trait,
data.type,'rds',
sep = '.')))
sample.vec <- readRDS(file =
file.path(base.dir,paste('sample.vec',
trait,
data.type,'rds',
sep = '.')))
df <- df[sample.vec, ]
print(paste('sub-sampled by',var.sample., 'by loading sample.vec'))
}
......@@ -199,8 +188,7 @@ load.data.for.lmer <- function(trait, data.type,
list.sd <- get.sd.lmer(df, trait, min.obs = 10)
res <- format.data.for.lmer(df, trait,
data.type = data.type,
ecocode.var = ecocode.var.,
Multi.type = Multi.type.)
ecocode.var = ecocode.var.)
return( list(res, list.sd))
}
......@@ -295,38 +283,19 @@ data.tree <- select.data.for.analysis(data.tree,
return(data.tree)
}
select.data.Multi <- function(data.tree, Multi.type = 'a', min.obs = 10) {
if(Multi.type == 'a')
traits <- c("SLA", "Wood.density", "Max.height")
if(Multi.type == 'b')
traits <- c("SLA", "Wood.density", "Max.height", 'Seed.mass')
df <- select.data.trait(data.tree, traits[1], min.obs)
for (t in traits[-1]){
df <- select.data.trait(df, t, min.obs)
}
return(df)
}
load.and.save.data.for.lmer <- function(trait,
min.obs= 10,
data.type = 'simple',
Multi.type = 'a',
base.dir = "output/processed"){
fname <- switch(data.type, 'simple' = 'data.all.no.log.all.census.rds',
'all.census' = 'data.all.no.log.all.census.rds',
'Multi' = 'data.all.global.all.census.rds')
'intra' = 'data.all.no.log.all.census.rds')
data.tree.tot <- readRDS(file.path(base.dir, fname))
if (data.type == 'Multi'){
df <- select.data.Multi(data.tree.tot, Multi.type)
saveRDS(df,file = file.path(base.dir,paste('data', 'Multi', data.type, 'rds',
sep = '.')))
}else{
df <- select.data.trait(data.tree.tot, trait)
saveRDS(df,file = file.path(base.dir,paste('data', trait, data.type, 'rds',
sep = '.')))
}
}
###################
......@@ -386,9 +355,9 @@ return(data.frame(logG = logG,
}
get.variables.Multi <- function(data.tree, BATOT, Multi.type = 'a',
ecocode.var = 'wwf',
min.BA.G = 40, min.MAT = 10){
get.variables.intra <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, ecocode.var = 'wwf',
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"]]))
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
......@@ -396,45 +365,17 @@ MAP <- scale.d(log(data.tree[["MAP"]]))
species.id <- as.character(factor(data.tree[["sp.name"]]))
plot.id <- as.character(factor(data.tree[["plot"]]))
tree.id <- as.character(factor(data.tree[["tree.id"]]))
set.id <- as.character(factor(data.tree[["set"]]))
ecocode.id <- as.character(factor(data.tree[[ecocode.var]]))
biomes.id <- factor(data.tree[['biomes']])
#= multiply CWMs by BATOT for all traits selected
if(Multi.type == 'a') {
traits <- c('SLA', 'Wood.density', 'Max.height')
sumTnTfBn.abs.Multi <-
scale.nc(data.tree[['Multi1.abs.CWM.fill']]*
data.tree[[BATOT]], center = FALSE)
}
if(Multi.type == 'b') {
traits <- c('SLA', 'Wood.density', 'Max.height', 'Seed.mass')
sumTnTfBn.abs.Multi <-
scale.nc(data.tree[['Multi2.abs.CWM.fill']]*
data.tree[[BATOT]], center = FALSE)
}
list.Tf <- vector('list')
list.sumTnTfBn.abs <- vector('list')
list.sumTnBn <- vector('list')
list.sumTfBn <- vector('list')
for (t in traits){
abs.CWM.tntf <- paste0(t, '.abs.CWM.fill')
CWM.tn <- paste0(t, '.CWM.fill')
tf <- paste0(t, '.focal')
list.sumTnTfBn.abs[[t]] <- scale.nc(data.tree[[abs.CWM.tntf]]*
data.tree[[BATOT]], center = FALSE)
list.sumTnBn[[t]] <- scale.nc(data.tree[[CWM.tn]]*data.tree[[BATOT]],
center = FALSE)
list.sumTfBn[[t]] <- scale.nc(data.tree[[tf]]*data.tree[[BATOT]],
center = FALSE)
list.Tf[[t]] <- scale.d(data.tree[[tf]])
}
names(list.sumTnTfBn.abs) <- paste0(traits, '.sumTnTf.abs')
names(list.sumTnBn) <- paste0(traits, '.sumTnBn')
names(list.sumTfBn) <- paste0(traits, '.sumTfBn')
names(list.Tf) <- paste0(traits, '.Tf')
sumBn <- scale.d(data.tree[[BATOT]], center = FALSE)
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*
(data.tree[[BATOT]] - data.tree[['BAintra']])
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumBn.inter <- (data.tree[[BATOT]] - data.tree[['BAintra']])
return(data.frame(logG = logG,
logD = logD,
MAT = MAT,
......@@ -445,87 +386,27 @@ return(data.frame(logG = logG,
biomes.id = biomes.id,
plot.id = plot.id,
tree.id = tree.id,
list.sumTnTfBn.abs,
list.Tf,
list.sumTnBn,
list.sumTfBn,
sumTnTfBn.abs.Multi = sumTnTfBn.abs.Multi,
sumBn = sumBn,
sumTnTfBn.abs = scale.nc(sumTnTfBn.abs, center = FALSE),
Tf = scale.d(data.tree[[tf]]),
sumTnBn = scale.nc(sumTnBn, center = FALSE),
sumTfBn = scale.nc(sumTfBn, center = FALSE),
sumBn.inter = scale.nc(sumBn.inter, center = FALSE),
sumBn.intra = scale.nc(data.tree[['BAintra']],
center = FALSE),
stringsAsFactors = FALSE))
}
get.variables.cat <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, ecocode.var = 'wwf',
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"]]))
MAT <- scale.d(log(data.tree[["MAT"]]+min.MAT))
MAP <- scale.d(log(data.tree[["MAP"]]))
species.id <- factor(data.tree[["sp.name"]])
plot.id <- factor(data.tree[["plot"]])
tree.id <- factor(data.tree[["tree.id"]])
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[[ecocode.var]])
biomes.id <- factor(data.tree[['biomes']])
#get the three cwm per cat
vec.CWM.tn <- paste(CWM.tn, c('A_EV', 'A_D', 'C'), sep = '.')
vec.abs.CWM.tntf <- abs.CWM.tntf
sumTnTfBn.abs<- data.tree[[vec.abs.CWM.tntf]]
sumTnBn.A_EV<- data.tree[[vec.CWM.tn[1]]]
sumTnBn.A_D <- data.tree[[vec.CWM.tn[2]]]
sumTnBn.C <- data.tree[[vec.CWM.tn[3]]]
sumTfBn <- scale.d(data.tree[[tf]])*data.tree[[BATOT]]
sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *sumTfBn
sumBn <- data.tree[[BATOT]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *
scale.d(data.tree[[tf]])
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *
scale.d(data.tree[[tf]])
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *
scale.d(data.tree[[tf]])
return(
data.frame(logG = logG,
logD = logD,
MAT = MAT,
MAP = MAP,
species.id = species.id,
set.id = set.id,
ecocode.id = ecocode.id,
biomes.id = biomes.id,
plot.id = plot.id,
tree.id = tree.id,
sumTnTfBn.abs =
scale.nc(sumTnTfBn.abs, center = FALSE),
Tf.A_EV= scale.d(Tf.A_EV),
Tf.A_D= scale.d(Tf.A_D),
Tf.C= scale.d(Tf.C),
sumTnBn.A_EV= scale.nc(sumTnBn.A_EV, center = FALSE),
sumTnBn.A_D= scale.nc(sumTnBn.A_D, center = FALSE),
sumTnBn.C = scale.nc(sumTnBn.C, center = FALSE),
sumTfBn.A_EV= scale.nc(sumTfBn.A_EV, center = FALSE),
sumTfBn.A_D= scale.nc(sumTfBn.A_D, center = FALSE),
sumTfBn.C = scale.nc(sumTfBn.C, center = FALSE),
sumBn = scale.nc(sumBn, center = FALSE)))
}
#============================================================
# 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',
Multi.type = 'a') {
data.type = 'simple', ecocode.var = 'wwf') {
if(! trait %in% c("SLA", "Leaf.N", "Seed.mass",
"SLA", "Wood.density", "Max.height", "Multi"))
"SLA", "Wood.density", "Max.height"))
stop("need trait to be in SLA Leaf.N Seed.mass
SLA Wood.density or Max.height Multi")
SLA Wood.density or Max.height")
# get vars names
CWM.tn <- paste(trait, "CWM", 'fill', sep = ".")
abs.CWM.tntf <- paste(trait, "abs.CWM", 'fill', sep = ".")
......@@ -533,7 +414,7 @@ tf <- paste(trait, "focal", sep = ".")
BATOT <- "BATOT"
#= DATA LIST FOR LMER
if(data.type == 'simple' | data.type == 'all.census') {
print('no Multi')
print('simple')
lmer.data <- get.variables(data.tree,
BATOT,
CWM.tn,
......@@ -541,11 +422,10 @@ if(data.type == 'simple' | data.type == 'all.census') {
tf,
ecocode.var)
}
if(data.type =='Multi') {
print('Multi')
lmer.data <- get.variables.Multi(data.tree,
if(data.type =='intra') {
print('intra')
lmer.data <- get.variables.intra(data.tree,
BATOT,
Multi.type,
ecocode.var)
}
return(lmer.data)
......@@ -554,10 +434,9 @@ if(data.type =='Multi') {
## SD for plot
get.sd.lmer <- function(data.tree, trait, min.obs = 10) {
if(! trait %in% c("SLA", "Leaf.N", "Seed.mass",
"SLA", "Wood.density", "Max.height", "Multi"))
"SLA", "Wood.density", "Max.height"))
stop("need trait to be in SLA Leaf.N Seed.mass
SLA Wood.density or Max.height Multi")
if(trait != 'Multi'){
SLA Wood.density or Max.height")
# get vars names
CWM.tn <- paste(trait, "CWM", 'fill', sep = ".")
abs.CWM.tntf <- paste(trait, "abs.CWM", 'fill', sep = ".")
......@@ -566,9 +445,6 @@ BATOT <- "BATOT"
#= DATA SD LIST FOR LMER
list.sd <- compute.sd.mean.var.lmer(data.tree, tf, CWM.tn,
abs.CWM.tntf, min.BA.G = 40)
}else{
list.sd <- NA
}
return(list.sd)
}
......@@ -582,6 +458,9 @@ list.sd <- list(
'mean.tf' = mean(data.tree[[tf]] , na.rm = TRUE),
'sd.tf' = sd(data.tree[[tf]] , na.rm = TRUE),
'sd.sumBn' = sd( data.tree[['BATOT']], na.rm = TRUE),
'sd.sumBn.inter' = sd( data.tree[['BATOT']] - data.tree[['BAintra']],
na.rm = TRUE),
'sd.sumBn.intra' = sd( data.tree[['BAintra']], na.rm = TRUE),
'sd.sumTfBn' = sd(data.tree[[tf]]* data.tree[['BATOT']], na.rm = TRUE),
'sd.sumTnTfBn.abs' = sd(data.tree[[abs.CWM.tntf]]*
data.tree[['BATOT']], na.rm = TRUE),
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.fixed.biomes.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.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|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.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|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.tree.id.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|tree.id)+(1|species.id)+(1|plot.id)+Tf+logD+MAT+MAP+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.Multi.r.set.fixed.biomes.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+SLA.Tf+Wood.density.Tf+Max.height.Tf+SLA.Tf:biomes.id+Wood.density.Tf:biomes.id+Max.height.Tf:biomes.id+logD+sumBn+sumBn:biomes.id+SLA.sumTfBn+Wood.density.sumTfBn+Max.height.sumTfBn+SLA.sumTfBn:biomes.id+Wood.density.sumTfBn:biomes.id+Max.height.sumTfBn:biomes.id+SLA.sumTnBn+Wood.density.sumTnBn+Max.height.sumTnBn+SLA.sumTnBn:biomes.id+Wood.density.sumTnBn:biomes.id+Max.height.sumTnBn:biomes.id+sumTnTfBn.abs.Multi+sumTnTfBn.abs.Multi:biomes.id +(logD-1|species.id) +(sumBn-1|species.id)+(SLA.Tf-1|set.id)+(Wood.density.Tf-1|set.id)+(Max.height.Tf-1|set.id)+(sumBn-1|set.id)+(SLA.sumTfBn-1|set.id)+(Wood.density.sumTfBn-1|set.id)+(Max.height.sumTfBn-1|set.id)+(SLA.sumTnBn-1|set.id)+(Wood.density.sumTnBn-1|set.id)+(Max.height.sumTnBn-1|set.id)+(sumTnTfBn.abs.Multi-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.Multi.r.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+SLA.Tf+Wood.density.Tf+Max.height.Tf+logD+sumBn+SLA.sumTfBn+Wood.density.sumTfBn+Max.height.sumTfBn+SLA.sumTnBn+Wood.density.sumTnBn+Max.height.sumTnBn+sumTnTfBn.abs.Multi +(logD-1|species.id) +(sumBn-1|species.id)+(SLA.Tf-1|set.id)+(Wood.density.Tf-1|set.id)+(Max.height.Tf-1|set.id)+(sumBn-1|set.id)+(SLA.sumTfBn-1|set.id)+(Wood.density.sumTfBn-1|set.id)+(Max.height.sumTfBn-1|set.id)+(SLA.sumTnBn-1|set.id)+(Wood.density.sumTnBn-1|set.id)+(Max.height.sumTnBn-1|set.id)+(sumTnTfBn.abs.Multi-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.intra.r.set.fixed.biomes.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+biomes.id+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)(logD-1|species.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.intra.r.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+Tf+logD+sumBn.inter+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)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.ecocode.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs +(1+logD+sumBn||species.id)+(1|plot.id)+(1+Tf+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs||ecocode.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.tree.id.set.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|tree.id)+(1|species.id)+(1|plot.id)++(1|biomes.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)+(sumTfBn-1|biomes.id)+(sumTnBn-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.tree.id.set.fixed.biomes.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|tree.id)+(1|set.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)+(Tf-1|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.tree.id.set.species",
var.BLUP = 'set.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|tree.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|set.id)+(sumBn-1|set.id)+(sumTfBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
......@@ -188,7 +188,6 @@ fun.turn.list.in.DF <- function(sp, res.list) {
fun.extract.format.sp.traits.TRY <- function(sp, sp.syno.table, data) {
### test data sp and sp.syno.table match
require(gdata)
browser()
sp.syno.table[["Latin_name_syn"]] <-
trim(as.character(sp.syno.table[["Latin_name_syn"]]) )
data[["Latin_name"]] <- trim(as.character(data[["Latin_name"]]))
......
......@@ -267,7 +267,8 @@ return(data)
fun.order.colums <- function(data, var.order = c("obs.id", "tree.id", "sp",
"sp.name", "cluster", "plot",
"ecocode", "koppen", "wwf", "D", "G", "BA.G", "year",
"ecocode", "koppen", "wwf", "D", "G",
"BA.G", "year",
"dead", "Lon", "Lat", "census", "MAT",
"MAP", "biomes", "BA.w", "Leaf.N.focal",
"Leaf.N.CWM.fill",
......@@ -293,10 +294,12 @@ fun.order.colums <- function(data, var.order = c("obs.id", "tree.id", "sp",
"Max.height.perc.species",
"Multi1.abs.CWM.fill",
"Multi2.abs.CWM.fill", "BATOT",
"BAintra",
"Phylo.group", "Pheno.T",
"LeafType.T", "cat")){
require(dplyr)
data <- data %>% dplyr::select(-Leaf.N.genus, -Seed.mass.genus, - SLA.genus,
data <- data %>% dplyr::select(-Leaf.N.genus, -Seed.mass.genus,
- SLA.genus,
-Wood.density.genus, -Max.height.genus)
data <- data[, var.order]
......@@ -384,13 +387,13 @@ require(dplyr)
data.plot.sp<- mutate(data, plot.sp = paste(plot.c, sp)) %>% group_by(plot.sp) %>%
dplyr::summarise(BAintra = sum(BA.w)) %>% ungroup()
data <- left_join(data, data.plot.sp, by = 'plot.sp') %>% dplyr::select(-plot.sp)
## remove BA obs tree
data <- data %>%
mutate(
BATOT = BATOT - BA.w,
BAintra = BAintra - BA.w,
BAintra = BAintra - BA.w,
Leaf.N.CWM.fill = (Leaf.N.CWM.fill - BA.w*Leaf.N.mean)/BATOT,
Leaf.N.CWM.fill = ifelse(BATOT == 0, 0, Leaf.N.CWM.fill),
SLA.CWM.fill = (SLA.CWM.fill - BA.w*SLA.mean)/BATOT,
......@@ -506,7 +509,6 @@ fun.CWM.traits.all.plot.census.dplyr <- function(data,data.TRAITS){
by = 'sp')
# compute CWM Tn and BATOT
data <- fun.CWM.Tn(data)
# compute CWM abs
data.CWM.abs <- data %>% group_by(plot.c) %>%
do(fun.CWM.abs.all(.)) %>%
......@@ -544,6 +546,7 @@ BA.n <- areas_by_species_within_neighbourhood(idx=i.t - 1L,
type=as.vector(sp.num) - 1L, n_types=sp.length)
BATOT <- sum(BA.n)
names(BA.n) <- sp.names
BAintra <- sum(BA.n[sp.names[sp.num[obs.id == i]]])
## compute community weighted mean for each traits
res.list <- lapply(traits.list, FUN=format.one.trait.CWM,
sp.names[sp.num][obs.id == i], BA.n)
......@@ -559,7 +562,7 @@ Multi2.CWM.abs.fill <- format.multi.traits.CWM(traits.list,
'Max.height',
'Seed.mass'))
res.vec <- c(i, unlist(res.list), Multi1.CWM.abs.fill, Multi2.CWM.abs.fill,
BATOT)
BATOT, BAintra)
return(res.vec)
}
......@@ -589,7 +592,7 @@ colnames(ans) <- c("obs.id", unlist(lapply(names(list.traits),
y=name.var)),
"Multi1.abs.CWM.fill",
"Multi2.abs.CWM.fill",
"BATOT")
"BATOT", "BAintra")
return(as.data.frame(ans))
}
......
......@@ -22,7 +22,7 @@ load_all_processed_data <- function(set, pathdir){
data[[ecocode.select]] <- load.processed.data(file.path(pathdir, set,
ecocode.select))
}
data
data
}
list_all_processed_data <- function(set, pathdir){
......@@ -44,7 +44,7 @@ fun.test.one.var.q(data[[var.name[i]]], t.q.01.vec[i], t.q.9.vec[i])
}
# test range D G BA.G
# test range D G BA.G
fun.test.range.vars.q <- function(data.tree, set){
vars.name <- c( "D", "dead", "G", "BA.G")
t.q.01.vec <- c(9.9, 0, -2, -80)
......@@ -71,7 +71,7 @@ if(any(vec.test))
fun.test.CWM.plot.r <- function(data, data.TRAITS, data.tree){
rep.count <- 0
repeat{
res <- fun.test.CWM.plot(data.t = data, data.TRAITS, data.tree)
rep.count <- rep.count +1
if(!is.na(res)){
......@@ -132,7 +132,7 @@ fun.test.CWM.plot <- function(data.t, data.TRAITS, data.tree) {
!data.TRAITS[data.TRAITS[["sp"]] == samp.sp, "Wood.density.genus"] &
!is.na(data[["Wood.density.CWM.fill"]][data[["obs.id"]] == samp.id])) {
test.focal <- all.equal(data.TRAITS[data.TRAITS[["sp"]] == samp.sp,
"Wood.density.mean"],
"Wood.density.mean"],
data[["Wood.density.focal"]][data[["obs.id"]]
== samp.id])
res.fill <- fun.compute.CWM.test.plot(samp.id, samp.plot, data.tree.t,
......@@ -158,7 +158,7 @@ rep.count <- 0
if(!is.na(res)){
break
}
if(rep.count>1000) stop('error no Wood density data')
if(rep.count>1000) stop('error no Wood density data')
}
return(res)
}
...