Commit d478e536 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

start adde lmer all

parent 40271ad6
......@@ -11,13 +11,13 @@ run.models.for.set.all.traits <- function(model.file,fun.model, traits =
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, type.filling=type.filling,...)
}
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, type.filling, ecoregions = get.ecoregions.for.set(set), ...){
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, type.filling, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait, type.filling=type.filling, ...))
(run.model.for.set.one.trait (m, fun.model,trait, type.filling=type.filling, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), ...){
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling,...))
}
......@@ -73,7 +73,7 @@ run.lmer <- function (model.file, trait,
cat("Ok data with Nobs",nrow(df.lmer),
"\n")
#= Run model
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out,std=std)
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out)
}
......@@ -91,7 +91,9 @@ load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
data.tree.tot <- read.csv(file.path(base.dir,"data.all.csv"),
library(data.table)
data.tree.tot <- fread(file.path(base.dir,"data.all.csv"),
stringsAsFactors = FALSE)
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
......@@ -105,9 +107,10 @@ data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["BA.G"]]>min.BA.G & data.tree[["BA.G"]]<max.BA.G &
data.tree[["D"]]>9.9)
## select species with no missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]]))
data.tree <- data.tree[select.temp,]
# select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor(data.tree[["sp"]]))>min.obs])
......@@ -123,75 +126,70 @@ return((x-mean(x))/sd(x))
### TODO NEED TO CONTIMNUE TO CHANGE FROM HERE
### get variables for lmer
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
fun.get.the.variables.for.lmer.tree.id <- function(data.tree,BATOT, CWM.tn,
abs.CWM.tntf, tf,
min.BA.G=50){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- log(data.tree[["D"]])
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(logG=logG,
logD=logD,
species.id=species.id,
tree.id=tree.id,
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
sumBn=sumBn))
}
species.id <- unclass(factor(paste(data.tree[['set']],
data.tree[["sp"]],sep='.')))
tree.id <- unclass(factor(paste(data.tree[['set']],data.tree[["tree.id"]])))
plot.species.id <- unclass(factor(paste(data.tree[["set"]],
data.tree[["plot"]],
data.tree[["sp"]],sep="")))
fun.get.the.variables.for.lmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- log(data.tree[["D"]])
species.id <- unclass(factor(data.tree[["sp"]]))
tree.id <- unclass(factor(data.tree[["tree.id"]]))
plot.species.id <- unclass(factor(paste(data.tree[["plot"]],data.tree[["sp"]],sep="")))
#= multiply CWMs by BATOT
set.id <- unclass(factor(data.tree[["set"]]))
ecocode.id <- unclass(factor(paste(data.tree[["set"]],data.tree[['ecocode']])))
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
sumTfBn <- data.tree[[tf]]*data.tree[[BATOT]]
sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]]
return(data.frame(logG=logG,
logD=logD,
species.id=species.id,
plot.species.id=plot.species.id,
sumTnTfBn.diff=sumTnTfBn.diff,
sumTnTfBn.abs=sumTnTfBn.abs,
Tf=data.tree[[tf]],
sumTnBn=sumTnBn,
sumTfBn=sumTfBn,
sumBn=sumBn))
return(data.frame(logG= logG,
logD= logD,
species.id= species.id,
tree.id= tree.id,
set.id= set.id,
ecocode.id= ecocode.id,
plot.species.id= plot.species.id,
sumTnTfBn.diff= sumTnTfBn.diff,
sumTnTfBn.abs= sumTnTfBn.abs,
Tf= data.tree[[tf]],
sumTnBn= sumTnBn,
sumTfBn= sumTfBn,
sumBn= sumBn))
}
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
fun.data.for.lmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
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 SLA Wood.density or Max.height ")
fun.data.for.lmer <- function(data.tree, trait, min.obs=10,
type.filling='species') {
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
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=".")
tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".")
data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs)
#= DATA LIST FOR JAGS
if (length(table(table(data.tree[["tree.id"]])))>1){
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
if (length(table(table(data.tree[["tree.id"]])))<2){
lmer.data <- fun.get.the.variables.for.lmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
}
data.tree <- fun.select.data.for.analysis(data.tree,
abs.CWM.tntf,
perc.neigh,
BATOT,
min.obs)
#= DATA LIST FOR LMER
lmer.data <- fun.get.the.variables.for.lmer.tree.id(data.tree,
BATOT,
CWM.tn,
abs.CWM.tntf,
tf)
return(lmer.data)
}
......
##### SCRIPT TO TEST BEFORE TO SEND ON CLUSTER
source("R/analysis/lmer.run.nolog.all.R")
run.lmer(model.files.lmer.Tf.1[3],'SLA',type.filling='species')
run.models.for.set.all.traits(model.files.lmer.Tf.1[3], run.lmer,type.filling='species')
##### SCRIPT TO TEST BEFORE TO SEND ON CLUSTER
source("R/analysis/lmer.run.nolog.R")
sets <- c('France','Spain','Sweden','Swiss','BCI','Fushan','Luquillo','NVS','Japan','Canada','Paracou','Mbaiki','US')
library(doParallel)
cl <- makeCluster(7)
......
......@@ -22,30 +22,25 @@ if(dim(data.all)[1] != sum(mat.perc[['num.obs']]))
########
### TODO
#- compare loading time
#- look at BATOT ALONG MAP AND MAT (log scale)
#- how to compute FD on plot with diferent size
#- pattern of CWM
#- pattern ED angio/conif
## data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
## trim.negative.growth(dbh1=data.all[['D']]*10,
## dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
## data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
## trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
## remove outlier following Condit approach
data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
(data.all[['G']] > -50 & !is.na(data.all[['G']])))] <- NA
trim.negative.growth(dbh1=data.all[['D']]*10,
dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
(data.all[['G']] > -50 & !is.na(data.all[['G']])))] <- NA
trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
### compute mean BATOT, number of species, traits and VAR OF TRAITS per cluster
system.time(data.summarise <- fun.compute.all.var.cluster(data.all))
### NEED TO CHECK WHY JAPAN REACH 300 of BATOT
plot(data.summarise$MAP,data.summarise$BATOT,
,col=col.vec[data.summarise$set])
......@@ -82,33 +77,14 @@ ggplot(data.summarise[,], aes( x=n_sp, y=sd.Wood.density) ) +facet_wrap(~set)+
geom_point(size=1 ) + geom_density2d()
ggplot(mydata100, aes(pretest, posttest)) +
> geom_point()
## plot
geom_point()
fun.plot.hist.trait.per.set(data.all)
to.pdf(fun.hist.var.set(data.all,var='BATOT',cex=0.6),
filename="figs/test.processed/fig.BATOT.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='G',cex=0.6),
filename="figs/test.processed/fig.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='BA.G',cex=0.6),
filename="figs/test.processed/fig.BA.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='D',cex=0.6),
filename="figs/test.processed/fig.D.set.pdf")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.D.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='G',cex=0.6,col=col.vec[data.all$set]),dev=png,
filename="figs/test.processed/fig.xy.D.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.G.set.png")
fun.plot.xy.trait.per.set(data.all,var.x='BATOT',var.y='BA.G')
......
......@@ -24,6 +24,34 @@ data.all.B <-data.all.B[,names(data.all.I)]
data.all <- rbind(data.all.I,data.all.B)
rm(list.all.I,list.all.B,data.all.I,data.all.B)
## remove outlier following Condit approach
data.all[['G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(dbh1=data.all[['D']]*10,
dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
data.all[['BA.G']][!(trim.positive.growth(data.all[['G']]) &
trim.negative.growth(data.all[['D']]*10,dbh2=data.all[['D']]*10 +data.all[['year']]*data.all[['G']]))] <- NA
write.csv(data.all,file=file.path(filedir, "data.all.csv"), row.names=FALSE)
rm(data.all)
gc()
## plots
fun.plot.hist.trait.per.set(data.all)
to.pdf(fun.hist.var.set(data.all,var='BATOT',cex=0.6),
filename="figs/test.processed/fig.BATOT.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='G',cex=0.6),
filename="figs/test.processed/fig.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='BA.G',cex=0.6),
filename="figs/test.processed/fig.BA.G.set.pdf")
to.pdf(fun.hist.var.set(data.all,var='D',cex=0.6),
filename="figs/test.processed/fig.D.set.pdf")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='BA.G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.D.BA.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='D',var.y='G',cex=0.6,col=col.vec[data.all$set]),dev=png,
filename="figs/test.processed/fig.xy.D.G.set.png")
to.dev(fun.plot.xy.set(data.all,var.x='BATOT',var.y='G',cex=0.6),dev=png,
filename="figs/test.processed/fig.xy.BATOT.G.set.png")
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