Commit 6fd5b3c4 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

run with cat

parent dedc93b2
...@@ -4,7 +4,7 @@ D3 := output/processed ...@@ -4,7 +4,7 @@ D3 := output/processed
D4 := figs/test.format.tree D4 := figs/test.format.tree
D5 := figs/test.traits D5 := figs/test.traits
D6 := figs/test.CWM D6 := figs/test.CWM
sites:= Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS US sites:= Fushan Paracou BCI Mbaiki Luquillo Japan Spain Sweden Canada France Swiss NSW NVS
D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) )) D3Done := $(addsuffix /Done.txt,$(addprefix $(D3)/, $(sites) ))
D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) )) D2traits := $(addsuffix /traits.csv,$(addprefix $(D2)/, $(sites) ))
D2tree := $(addsuffix /tree.csv,$(addprefix $(D2)/, $(sites) )) D2tree := $(addsuffix /tree.csv,$(addprefix $(D2)/, $(sites) ))
...@@ -71,7 +71,7 @@ $(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -t ...@@ -71,7 +71,7 @@ $(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -t
Mbaiki: $(D3)/Mbaiki/Done.txt Mbaiki: $(D3)/Mbaiki/Done.txt
$(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv $(D3)/Mbaiki/Done.txt: R/process.data/process-fun.R $(D2)/Mbaiki/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='no'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global');" Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='global');"
$(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait-fun.R $(D2)/Mbaiki/tree.csv $(D2)/Mbaiki/traits.csv: R/find.trait/Mbaiki.R R/find.trait/trait-fun.R $(D2)/Mbaiki/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -84,7 +84,7 @@ $(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f) ...@@ -84,7 +84,7 @@ $(D2)/Mbaiki/tree.csv: R/format.data/Mbaiki.R $(shell find $(D1)/Mbaiki -type f)
Canada: $(D3)/Canada/Done.txt Canada: $(D3)/Canada/Done.txt
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv $(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='no'); process_inventory_dataset('Canada',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='global');"
$(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait-fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/Canada/traits.csv: R/find.trait/Canada.R R/find.trait/trait-fun.R $(D2)/Canada/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -96,7 +96,7 @@ $(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f) ...@@ -96,7 +96,7 @@ $(D2)/Canada/tree.csv: R/format.data/Canada.R $(shell find $(D1)/Canada -type f)
France: $(D3)/France/Done.txt France: $(D3)/France/Done.txt
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv $(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='no'); process_inventory_dataset('France',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='global');"
$(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait-fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/France/traits.csv: R/find.trait/France.R R/find.trait/trait-fun.R $(D2)/France/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -109,7 +109,7 @@ $(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f) ...@@ -109,7 +109,7 @@ $(D2)/France/tree.csv: R/format.data/France.R $(shell find $(D1)/France -type f)
Fushan: $(D3)/Fushan/Done.txt Fushan: $(D3)/Fushan/Done.txt
$(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv $(D3)/Fushan/Done.txt: R/process.data/process-fun.R $(D2)/Fushan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='no'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');" Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='global');"
$(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait-fun.R $(D2)/Fushan/tree.csv $(D2)/Fushan/traits.csv: R/find.trait/Fushan.R R/find.trait/trait-fun.R $(D2)/Fushan/tree.csv
...@@ -123,7 +123,7 @@ $(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f) ...@@ -123,7 +123,7 @@ $(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
NSW: $(D3)/NSW/Done.txt NSW: $(D3)/NSW/Done.txt
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv $(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='no'); process_inventory_dataset('NSW',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='global');"
$(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait-fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/NSW/traits.csv: R/find.trait/NSW.R R/find.trait/trait-fun.R $(D2)/NSW/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -136,7 +136,7 @@ $(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f) ...@@ -136,7 +136,7 @@ $(D2)/NSW/tree.csv: R/format.data/NSW.R $(shell find $(D1)/NSW -type f)
NVS: $(D3)/NVS/Done.txt NVS: $(D3)/NVS/Done.txt
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv $(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='no'); process_inventory_dataset('NVS',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='global');"
$(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait-fun.R $(D2)/NVS/tree.csv $(D2)/NVS/traits.csv: R/find.trait/NVS.R R/find.trait/trait-fun.R $(D2)/NVS/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -149,7 +149,7 @@ $(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f) ...@@ -149,7 +149,7 @@ $(D2)/NVS/tree.csv: R/format.data/NVS.R $(shell find $(D1)/NVS -type f)
Paracou: $(D3)/Paracou/Done.txt Paracou: $(D3)/Paracou/Done.txt
$(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv $(D3)/Paracou/Done.txt: R/process.data/process-fun.R $(D2)/Paracou/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='no'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global');" Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='global');"
$(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait-fun.R $(D2)/Paracou/tree.csv $(D2)/Paracou/traits.csv: R/find.trait/Paracou.R R/find.trait/trait-fun.R $(D2)/Paracou/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -162,7 +162,7 @@ $(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type ...@@ -162,7 +162,7 @@ $(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type
Spain: $(D3)/Spain/Done.txt Spain: $(D3)/Spain/Done.txt
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv $(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='no'); process_inventory_dataset('Spain',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='global');"
$(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait-fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/Spain/traits.csv: R/find.trait/Spain.R R/find.trait/trait-fun.R $(D2)/Spain/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -175,7 +175,7 @@ $(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f) ...@@ -175,7 +175,7 @@ $(D2)/Spain/tree.csv: R/format.data/Spain.R $(shell find $(D1)/Spain -type f)
Sweden: $(D3)/Sweden/Done.txt Sweden: $(D3)/Sweden/Done.txt
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv $(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='no'); process_inventory_dataset('Sweden',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='global');"
$(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait-fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/Sweden/traits.csv: R/find.trait/Sweden.R R/find.trait/trait-fun.R $(D2)/Sweden/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -188,7 +188,7 @@ $(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f) ...@@ -188,7 +188,7 @@ $(D2)/Sweden/tree.csv: R/format.data/Sweden.R $(shell find $(D1)/Sweden -type f)
Swiss: $(D3)/Swiss/Done.txt Swiss: $(D3)/Swiss/Done.txt
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv $(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='no'); process_inventory_dataset('Swiss',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='global');"
$(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait-fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/Swiss/traits.csv: R/find.trait/Swiss.R R/find.trait/trait-fun.R $(D2)/Swiss/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -201,7 +201,7 @@ $(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f) ...@@ -201,7 +201,7 @@ $(D2)/Swiss/tree.csv: R/format.data/Swiss.R $(shell find $(D1)/Swiss -type f)
US: $(D3)/US/Done.txt US: $(D3)/US/Done.txt
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv $(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='no'); process_inventory_dataset('US',std.traits='global');" Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='global');"
$(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds $(D2)/US/traits.csv: R/find.trait/US.R R/find.trait/trait-fun.R $(D2)/US/tree.csv $(D2)/TRY/data.TRY.std.rds
Rscript $< ; rm -f $(D5)/Done.txt Rscript $< ; rm -f $(D5)/Done.txt
...@@ -218,13 +218,13 @@ $(D4)/Done.txt: R/format.data/test.tree.R $(D3tree) ...@@ -218,13 +218,13 @@ $(D4)/Done.txt: R/format.data/test.tree.R $(D3tree)
#------------------------------------------------------- #-------------------------------------------------------
TEST.TRAITS: $(D5)/Done.txt TEST.TRAITS: $(D5)/Done.txt
$(D5)/Done.txt: R/find.trait/test.traits.R $(D4)/Done.txt $(D3traits) $(D5)/Done.txt: R/find.trait/test.traits.R $(D3traits)
Rscript $< Rscript $<
#------------------------------------------------------- #-------------------------------------------------------
TEST.CWM: $(D6)/Done.txt TEST.CWM: $(D6)/Done.txt
$(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D5)/Done.txt $(D3Done) $(D6)/Done.txt: R/process.data/test.tree.CWM.R $(D3Done)
Rscript $< Rscript $<
#------------------------------------------------------- #-------------------------------------------------------
......
...@@ -58,7 +58,8 @@ model.files.glmer.Tf.4 <- ...@@ -58,7 +58,8 @@ model.files.glmer.Tf.4 <-
fun.call.glmer.and.save <- function(formula,df.lmer,path.out){ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time() Start <- Sys.time()
glmer.output <- glmer(formula=formula,data=df.lmer,family=binomial) glmer.output <- glmer(formula=formula,data=df.lmer,
family=binomial(link = 'cloglog'))
end <- Sys.time() end <- Sys.time()
print(end - Start) print(end - Start)
print(summary(glmer.output)) print(summary(glmer.output))
...@@ -67,12 +68,14 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){ ...@@ -67,12 +68,14 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
run.glmer <- function (model.file, trait, run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA, min.obs=5, sample.size=NA,
type.filling) { type.filling, fname = 'data.all.no.std.csv') {
require(lme4) require(lme4)
source(model.file, local = TRUE) source(model.file, local = TRUE)
model <- load.model() model <- load.model()
#= Path for output #= Path for output
path.out <- output.dir.glmer(model$name, trait, if(fname == 'data.all.no.std.csv') dir <- 'all.no.log'
if(fname == 'data.all.csv') dir <- 'all.log'
path.out <- output.dir.glmer(dir, model$name, trait,
type.filling=type.filling) type.filling=type.filling)
print(path.out) print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE) dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
...@@ -81,7 +84,8 @@ run.glmer <- function (model.file, trait, ...@@ -81,7 +84,8 @@ run.glmer <- function (model.file, trait,
trait,"\n") trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait, df.glmer <- load.and.prepare.data.for.glmer(trait,
min.obs, sample.size, min.obs, sample.size,
type.filling=type.filling) # return a DF type.filling=type.filling,
fname = fname) # return a DF
cat("Ok data with Nobs",nrow(df.glmer),"\n") cat("Ok data with Nobs",nrow(df.glmer),"\n")
#= Run model #= Run model
fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out) fun.call.glmer.and.save(formula=model$glmer.formula,df.glmer,path.out)
...@@ -89,8 +93,25 @@ run.glmer <- function (model.file, trait, ...@@ -89,8 +93,25 @@ run.glmer <- function (model.file, trait,
#======================================================================== #========================================================================
output.dir.glmer <- function (model, trait, type.filling) { output.dir.glmer <- function (dir, model, trait, type.filling) {
file.path("output/glmer", 'all',trait,type.filling,model) file.path("output/glmer", dir,trait,type.filling,model)
}
fun.load.data.all <- function(base.dir,fname){
data.all.sample <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
intern = TRUE))
data.tree.tot <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE,
nrows = nrows,
colClasses = classes)
return(data.tree.tot)
} }
...@@ -99,37 +120,37 @@ output.dir.glmer <- function (model, trait, type.filling) { ...@@ -99,37 +120,37 @@ output.dir.glmer <- function (model, trait, type.filling) {
#============================================================ #============================================================
load.and.prepare.data.for.glmer <- function(trait, load.and.prepare.data.for.glmer <- function(trait,
min.obs, sample.size, type.filling, min.obs, sample.size, type.filling,
fname = 'data.all.no.std.csv',
base.dir = "output/processed/"){ base.dir = "output/processed/"){
### load data ### load data
data.all.sample <- read.csv(file=file.path(base.dir, "data.all.csv"), data.tree.tot <- fun.load.data.all(base.dir,fname)
stringsAsFactors=FALSE,nrows=1000)
classes <- sapply(data.all.sample,class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',
intern=TRUE))
data.tree.tot <- read.csv(file=file.path(base.dir, "data.all.csv"),
stringsAsFactors=FALSE,
nrows=nrows,
colClasses=classes)
fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling) fun.data.for.glmer(data.tree.tot,trait,type.filling=type.filling)
} }
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh, fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs, BATOT,min.obs,
min.perc.neigh=0.80){ min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95){
## remove tree with NA ## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) & data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) ) (!is.na(data.tree[["D"]])) &
(!is.na(data.tree[['year']])))
## remove tree with zero ## remove tree with zero
data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9) data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits ## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) & data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),] !is.na(data.tree[[BATOT]])),]
# select obs abov min perc.neigh # select obs abov min perc.neigh species
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
!is.na(data.tree[[perc.neigh]])) min.perc.neigh.sp &
!is.na(data.tree[[perc.neigh.sp]]))
# select obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs & !is.na(data.tree[[perc.neigh.gs]]))
# select species with minimum obs # select species with minimum obs
data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in% data.tree <- subset(data.tree,subset=data.tree[["sp"]] %in%
names(table(factor(data.tree[["sp"]])))[table(factor( names(table(factor(data.tree[["sp"]])))[table(factor(
...@@ -146,6 +167,7 @@ fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn, ...@@ -146,6 +167,7 @@ fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,
abs.CWM.tntf,tf,min.BA.G=50){ abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]] dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(log(data.tree[["D"]])) logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
logyear <- log(data.tree[["year"]])
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]]) MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]]) MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
species.id <- factor(data.tree[["sp"]]) species.id <- factor(data.tree[["sp"]])
...@@ -161,6 +183,7 @@ sumTnTfBn.diff <- sumTnBn-sumTfBn ...@@ -161,6 +183,7 @@ sumTnTfBn.diff <- sumTnBn-sumTfBn
sumBn <- data.tree[[BATOT]] sumBn <- data.tree[[BATOT]]
return(data.frame(dead=dead, return(data.frame(dead=dead,
logD=logD, logD=logD,
logyear=logyear,
MAT = MAT, MAT = MAT,
MAP = MAP, MAP = MAP,
species.id=species.id, species.id=species.id,
...@@ -192,9 +215,13 @@ CWM.tn <- paste(trait,"CWM",'fill',sep=".") ...@@ -192,9 +215,13 @@ CWM.tn <- paste(trait,"CWM",'fill',sep=".")
abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".") abs.CWM.tntf <- paste(trait,"abs.CWM",'fill',sep=".")
tf <- paste(trait,"focal",sep=".") tf <- paste(trait,"focal",sep=".")
BATOT <- "BATOT" BATOT <- "BATOT"
perc.neigh <- paste(trait,"perc",type.filling,sep=".") perc.neigh.sp <- paste(trait, "perc", 'species', sep = ".")
data.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf, perc.neigh.gs <- paste(trait, "perc", 'genus', sep = ".")
perc.neigh,BATOT,min.obs) data.tree <- fun.select.data.for.analysis(data.tree,
abs.CWM.tntf,
perc.neigh.sp,
perc.neigh.gs,
BATOT,min.obs)
#= DATA LIST FOR JAGS #= DATA LIST FOR JAGS
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf) glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
return(glmer.data) return(glmer.data)
......
This diff is collapsed.
...@@ -12,17 +12,17 @@ run.models.for.set.all.traits <- function(model.file, fun.model, traits = ...@@ -12,17 +12,17 @@ run.models.for.set.all.traits <- function(model.file, fun.model, traits =
} }
run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait, run.multiple.model.for.set.one.trait <- function(model.files, fun.model, trait,
type.filling, ...){ type.filling, cat.TF, ...){
for (m in model.files) for (m in model.files)
(run.model.for.set.one.trait (m, fun.model, trait, (run.model.for.set.one.trait (m, fun.model, trait,
type.filling = type.filling, ...)) type.filling = type.filling, cat.TF = cat.TF, ...))
} }
run.model.for.set.one.trait <- function(model.file, fun.model, trait, run.model.for.set.one.trait <- function(model.file, fun.model, trait,
type.filling, ...){ type.filling, cat.TF, ...){
fun.model <- match.fun(fun.model) fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling = type.filling, ...)) try(fun.model(model.file, trait, type.filling = type.filling,cat.TF = cat.TF, ...))
} }
...@@ -62,26 +62,31 @@ fun.call.lmer.and.save <- function(formula, df.lmer, path.out){ ...@@ -62,26 +62,31 @@ fun.call.lmer.and.save <- function(formula, df.lmer, path.out){
run.lmer <- function (model.file, trait, run.lmer <- function (model.file, trait,
min.obs = 10, sample.size = NA, min.obs = 10, sample.size = NA,
type.filling) { type.filling, cat.TF, fname = 'data.all.no.std.csv') {
require(lme4) require(lme4)
source(model.file, local = TRUE) source(model.file, local = TRUE)
model <- load.model() model <- load.model()
#= Path for output #= Path for output
path.out <- output.dir.lmer(model$name, trait, 'all', if(fname == 'data.all.no.std.csv') dir <- 'all.no.log'
if(fname == 'data.all.csv') dir <- 'all.log'
path.out <- output.dir.lmer(model$name, trait, dir,
type.filling = type.filling) type.filling = type.filling)
print(path.out) print(path.out)
dir.create(path.out, recursive = TRUE, showWarnings = FALSE) dir.create(path.out, recursive = TRUE, showWarnings = FALSE)
cat("run lmer for model", model.file, " for trait", cat("run lmer for model", model.file, " for trait",
trait, "\n") trait, "\n")
df.lmer <- load.and.prepare.data.for.lmer(trait, df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size, min.obs, sample.size,
type.filling = type.filling) type.filling = type.filling,
fname = fname, cat.TF = cat.TF)
# return a DF # return a DF
cat("Ok data with Nobs", nrow(df.lmer), cat("Ok data with Nobs", nrow(df.lmer),
"\n") "\n")
#= Run model #= Run model
fun.call.lmer.and.save(formula = model$lmer.formula.tree.id, fun.call.lmer.and.save(formula = model$lmer.formula.tree.id,
df.lmer, path.out) df.lmer, path.out)
rm(df.lmer)
gc()
} }
...@@ -91,36 +96,38 @@ output.dir.lmer <- function (model, trait, set, type.filling) { ...@@ -91,36 +96,38 @@ output.dir.lmer <- function (model, trait, set, type.filling) {
file.path("output/lmer", set, trait, type.filling, model) file.path("output/lmer", set, trait, type.filling, model)
} }
fun.load.data.all <- function(base.dir,fname){
#============================================================ data.all.sample <- read.csv(file = file.path(base.dir, fname),
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
## library(data.table)
## data.tree.tot <- fread(file.path(base.dir, "data.all.csv"),
## stringsAsFactors = FALSE)
data.all.sample <- read.csv(file = file.path(base.dir, "data.all.csv"),
stringsAsFactors = FALSE, nrows = 1000) stringsAsFactors = FALSE, nrows = 1000)
classes <- sapply(data.all.sample, class) classes <- sapply(data.all.sample, class)
classes[classes=='integer'] <- "numeric" classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv', nrows <- as.numeric(system(paste('wc -l < ',file.path(base.dir, fname)),
intern = TRUE)) intern = TRUE))
data.tree.tot <- read.csv(file = file.path(base.dir, "data.all.csv"), data.tree.tot <- read.csv(file = file.path(base.dir, fname),
stringsAsFactors = FALSE, stringsAsFactors = FALSE,
nrows = nrows, nrows = nrows,
colClasses = classes) colClasses = classes)
return(data.tree.tot)
}
fun.data.for.lmer(data.tree.tot, trait, type.filling = type.filling) #============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
cat.TF,
fname = 'data.all.no.std.csv',
base.dir = "output/processed"){
data.tree.tot <- fun.load.data.all(base.dir,fname)
fun.data.for.lmer(data.tree.tot, trait, type.filling = type.filling, cat.TF = cat.TF)
} }
fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp, fun.select.data.for.analysis <- function(data.tree, abs.CWM.tntf, perc.neigh.sp,
perc.neigh.gs, BATOT, min.obs, perc.neigh.gs, BATOT, min.obs,
min.perc.neigh.sp = 0.90, min.perc.neigh.gs = 0.98, min.BA.G = -60, min.perc.neigh.sp = 0.90,
min.perc.neigh.gs = 0.95,
min.BA.G = -60,
max.BA.G = 500){ max.BA.G = 500){
## remove tree with NA ## remove tree with NA
data.tree <- subset(data.tree, subset = (!is.na(data.tree[["BA.G"]])) & data.tree <- subset(data.tree, subset = (!is.na(data.tree[["BA.G"]])) &
...@@ -134,10 +141,12 @@ select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) & ...@@ -134,10 +141,12 @@ select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])) !is.na(data.tree[[BATOT]]))
data.tree <- data.tree[select.temp, ] data.tree <- data.tree[select.temp, ]
# select obs abov min perc.neigh species # select obs abov min perc.neigh species
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] > min.perc.neigh.sp & data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.sp]] >
min.perc.neigh.sp &
!is.na(data.tree[[perc.neigh.sp]])) !is.na(data.tree[[perc.neigh.sp]]))
# select obs abov min perc.neigh genus # select obs abov min perc.neigh genus
data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] > min.perc.neigh.gs & data.tree <- subset(data.tree, subset = data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs &
!is.na(data.tree[[perc.neigh.gs]])) !is.na(data.tree[[perc.neigh.gs]]))
# select species with minimum obs # select species with minimum obs
data.tree <- subset(data.tree, subset = data.tree[["sp"]] %in% data.tree <- subset(data.tree, subset = data.tree[["sp"]] %in%
...@@ -154,7 +163,7 @@ return((x-mean(x))/sd(x)) ...@@ -154,7 +163,7 @@ return((x-mean(x))/sd(x))
fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn, fun.get.the.variables.for.lmer.tree.id <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf, abs.CWM.tntf, tf,
min.BA.G = 60){ min.BA.G = 60){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G)) logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]] + min.BA.G))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]])) logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]]) MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]]) MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
...@@ -189,13 +198,75 @@ return(data.frame(logG = logG, ...@@ -189,13 +198,75 @@ return(data.frame(logG = logG,
sumBn = fun.center.and.standardized.var(sumBn))) sumBn = fun.center.and.standardized.var(sumBn)))
} }
fun.get.the.variables.for.lmer.tree.id.cat <- function(data.tree, BATOT, CWM.tn,
abs.CWM.tntf, tf,
min.BA.G = 60){
logG <- fun.center.and.standardized.var(log(data.tree[["BA.G"]]+min.BA.G))
logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
MAT <- fun.center.and.standardized.var(data.tree[["MAT"]])
MAP <- fun.center.and.standardized.var(data.tree[["MAP"]])
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])
#get the three cwm per cat
vec.CWM.tn <- paste(CWM.tn, c('A_EV', 'A_D', 'C'), sep = '.'