Skip to content

GitLab

Commit 93f19645 authored by Georges Kunstler's avatar Georges Kunstler
Browse files

global analysis lmer and glmer

parent f9153d45
Hide whitespace changes
Inline Side-by-side
Makefile
View file @ 93f19645
......@@ -25,15 +25,15 @@ $(D2)/TRY/data.TRY.std.rds:
#-------------------------------------------------------
GLOBAL: $(D3)/Done.txt
$(D3)/Done.txt: R/process.data/global.R R/process.data/process-fun.R $(D2)/traits.std.global.csv $(D5)/Done.txt
Rscript $< ; Rscript R/process.data/merge.all.processed.data.R
$(D3)/Done.txt: R/process.data/merge.all.processed.data.R sites R/process.data/process-fun.R $(D2)/traits.std.global.csv $(D5)/Done.txt
Rscript $<
#-------------------------------------------------------
BCI: $(D3)/BCI/Done.txt
$(D3)/BCI/Done.txt: R/process.data/process-fun.R $(D2)/BCI/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='local');"
Rscript -e "source('$<'); process_bigplot_dataset('BCI', Rlim=15,std.traits='no'); process_bigplot_dataset('BCI', Rlim=15,std.traits='global');"
$(D2)/BCI/traits.csv: R/find.trait/BCI.R R/find.trait/trait-fun.R $(D2)/BCI/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt
......@@ -46,7 +46,7 @@ $(D2)/BCI/tree.csv: R/format.data/BCI.R $(shell find $(D1)/BCI -type f)
Japan: $(D3)/Japan/Done.txt
$(D3)/Japan/Done.txt: R/process.data/process-fun.R $(D2)/Japan/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='local');"
Rscript -e "source('$<'); process_bigplot_dataset('Japan', Rlim=15,std.traits='no'); process_bigplot_dataset('Japan', Rlim=15,std.traits='global');"
$(D2)/Japan/traits.csv: R/find.trait/Japan.R R/find.trait/trait-fun.R $(D2)/Japan/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt
......@@ -58,7 +58,7 @@ $(D2)/Japan/tree.csv: R/format.data/Japan.R $(shell find $(D1)/Japan -type f)
Luquillo: $(D3)/Luquillo/Done.txt
$(D3)/Luquillo/Done.txt: R/process.data/process-fun.R $(D2)/Luquillo/traits.csv
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='local');"
Rscript -e "source('$<'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='no'); process_bigplot_dataset('Luquillo', Rlim=15,std.traits='global');"
$(D2)/Luquillo/traits.csv: R/find.trait/Luquillo.R R/find.trait/trait-fun.R $(D2)/Luquillo/tree.csv
Rscript $< ; rm -f $(D5)/Done.txt
......@@ -71,7 +71,7 @@ $(D2)/Luquillo/tree.csv: R/format.data/Luquillo.R $(shell find $(D1)/Luquillo -t
Mbaiki: $(D3)/Mbaiki/Done.txt
$(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='local');"
Rscript -e "source('$<'); process_bigplot_dataset('Mbaiki', Rlim=15,std.traits='no'); 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
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)
Canada: $(D3)/Canada/Done.txt
$(D3)/Canada/Done.txt: R/process.data/process-fun.R $(D2)/Canada/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('Canada',std.traits='no'); 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
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)
France: $(D3)/France/Done.txt
$(D3)/France/Done.txt: R/process.data/process-fun.R $(D2)/France/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('France',std.traits='no'); 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
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)
Fushan: $(D3)/Fushan/Done.txt
$(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='local');"
Rscript -e "source('$<'); process_bigplot_dataset('Fushan', Rlim=15,std.traits='no'); 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
......@@ -123,7 +123,7 @@ $(D2)/Fushan/tree.csv: R/format.data/Fushan.R $(shell find $(D1)/Fushan -type f)
NSW: $(D3)/NSW/Done.txt
$(D3)/NSW/Done.txt: R/process.data/process-fun.R $(D2)/NSW/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('NSW',std.traits='no'); 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
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)
NVS: $(D3)/NVS/Done.txt
$(D3)/NVS/Done.txt: R/process.data/process-fun.R $(D2)/NVS/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('NVS',std.traits='no'); 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
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)
Paracou: $(D3)/Paracou/Done.txt
$(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='local');"
Rscript -e "source('$<'); process_bigplot_dataset('Paracou', Rlim=15,std.traits='no'); 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
Rscript $< ; rm -f $(D5)/Done.txt
......@@ -162,7 +162,7 @@ $(D2)/Paracou/tree.csv: R/format.data/Paracou.R $(shell find $(D1)/Paracou -type
Spain: $(D3)/Spain/Done.txt
$(D3)/Spain/Done.txt: R/process.data/process-fun.R $(D2)/Spain/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('Spain',std.traits='no'); 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
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)
Sweden: $(D3)/Sweden/Done.txt
$(D3)/Sweden/Done.txt: R/process.data/process-fun.R $(D2)/Sweden/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('Sweden',std.traits='no'); 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
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)
Swiss: $(D3)/Swiss/Done.txt
$(D3)/Swiss/Done.txt: R/process.data/process-fun.R $(D2)/Swiss/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('Swiss',std.traits='no'); 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
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)
US: $(D3)/US/Done.txt
$(D3)/US/Done.txt: R/process.data/process-fun.R $(D2)/US/traits.csv
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='local');"
Rscript -e "source('$<'); process_inventory_dataset('US',std.traits='no'); 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
Rscript $< ; rm -f $(D5)/Done.txt
......
R/analysis/glmer.all.output.R 0 → 100644
View file @ 93f19645
#!/usr/bin/env Rscript
source("R/analysis/glmer.output-fun.R")
source("R/analysis/glmer.run.all.R")
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (trait in traits){
for(model in c(model.files.lmer.Tf.1, model.files.lmer.Tf.2,
model.files.lmer.Tf.3, model.files.lmer.Tf.4)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.glmer(model.obj$name, trait, 'all',
type.filling=type.filling)
files <- c(files,file.path(pathout,"glmer.results.global.rds"))
}
}
out <- lapply(files, summarise.glmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/glmer.global.list.out.rds')
R/analysis/glmer.nolog.output.figs.R 0 → 100644
View file @ 93f19645
#!/usr/bin/env Rscript
rm(list = ls())
source("R/analysis/lmer.output-fun.R")
source("R/utils/plot.R")
## load results
DF.results <- readRDS('output/glmer.global.std.DF.rds')
table(DF.results$conv)
DF.results$id <- paste(DF.results$set, DF.results$ecocode, sep = ".")
## climate
site.clim.all <- read.csv( file.path("output/processed", "all.sites.clim.csv"),
stringsAsFactors = FALSE)
site.clim.all$id <- paste(site.clim.all$set, site.clim.all$ecocode, sep = ".")
mean.MAT <- tapply(site.clim.all$MAT, INDEX = site.clim.all$id,
FUN = mean, na.rm = TRUE)
mean.MAP <- tapply(site.clim.all$MAP, INDEX = site.clim.all$id,
FUN = mean, na.rm = TRUE)
unique.set <- tapply(site.clim.all$set, INDEX = site.clim.all$id,
FUN = unique, na.rm = TRUE)
data.clim.ecocode <- data.frame(id = names(mean.MAT),
MAT = mean.MAT,
MAP = mean.MAP,
set = unique.set, stringsAsFactors = FALSE)
biomes.vec <- fun.overly.plot.on.biomes(MAP = data.clim.ecocode$MAP/10,
MAT = data.clim.ecocode$MAT,
names.vec = 1:nrow(data.clim.ecocode))
data.clim.ecocode$biomes <- as.character(unlist(biomes.vec))
### merge climate and lmer results
DF.results <- merge(DF.results,
data.clim.ecocode[, c('id','MAT','MAP','biomes')] ,
by = "id")
DF.results$id2 <- paste(DF.results$id, DF.results$trait, DF.results$filling)
## add value for NVS
DF.results$biomes[DF.results$ecocode == 'MixedCool'] <- 'temperate_rainforest'
## ADD percentage traits
site.perc.all <- read.csv('output/processed/all.sites.perc.traits.csv',
stringsAsFactors = FALSE)
site.perc.all$id <- paste(site.perc.all$set, site.perc.all$ecocode, sep = ".")
##
DF.results <- merge(DF.results,
subset(site.perc.all,
select = c('id', 'perc.gymno',
'perc.ev', 'BATOT.m')),
by = "id")
### Analysis of the results
DF.R2m.diff <- do.call("rbind", lapply(1:nrow(DF.results),
fun.compute.criteria.diff,
DF.results, "R2m"))
DF.R2c.diff <- do.call("rbind", lapply(1:nrow(DF.results),
fun.compute.criteria.diff,
DF.results, "R2c"))
DF.AIC.diff <- do.call("rbind", lapply(1:nrow(DF.results),
fun.compute.criteria.diff,
DF.results, "AIC"))
DF.delta.AIC <- do.call("rbind", lapply(1:nrow(DF.results),
fun.compute.delta.AIC, DF.results))
DF.var.fixed <- fun.ratio.var.fixed.effect(DF.results)
DF.results <- cbind(DF.results, DF.R2m.diff, DF.R2c.diff,
DF.AIC.diff, DF.delta.AIC, DF.var.fixed)
## best model
DF.best.and.all.AIC <- do.call('rbind', lapply(unique(DF.results$id2),
FUN = fun.AIC, DF.results))
table(DF.best.and.all.AIC$best.model, DF.best.and.all.AIC$trait)
# global AIC
global.AIC.list <- fun.global.aic(DF.results)
### plots
models <- c('glmer.LOGLIN.ER.Tf', 'glmer.LOGLIN.ER.Tf')
names(models) <- c('Effect/response effect', 'Effect/response response')
list.params <- list(c(Response = 'sumTnBn'),
c(Effect = 'sumTfBn'))
## pdf('figs/parameters.glmer.MAP.ER.all.nolog.pdf', width = 9, height = 7)
fun.plot.panel.lmer.parameters.c(models = models,
traits = c('Wood.density', 'SLA',
'Leaf.N', 'Max.height'),
DF.results, var.x = 'MAP',
list.params = list.params,
small.bar = 0.2, threshold.delta.AIC = 10000000,
ylim = c(-0.2, 0.2), ylim.same = TRUE,
log = 'x', lwd = 1.5)
## dev.off()
models <- c('glmer.LOGLIN.AD.Tf')
names(models) <- c('Absolut distance effect')
list.params <- list(c(Response = 'sumTnTfBn.abs'))
## pdf('figs/parameters.glmer.MAP.ER.all.nolog.pdf', width = 9, height = 7)
fun.plot.panel.lmer.parameters.c(models = models,
traits = c('Wood.density', 'SLA',
'Leaf.N', 'Max.height'),
DF.results, var.x = 'MAP',
list.params = list.params,
small.bar = 0.2, threshold.delta.AIC = 10000000,
ylim = c(-0.005, 0.2), ylim.same = FALSE,
log = 'x', lwd = 1.5)
## dev.off()
models <- c('glmer.LOGLIN.ER.Tf')
names(models) <- c('ER Tf')
list.params <- list(c(Response = 'Tf'))
fun.plot.panel.lmer.parameters.c(models = models,
traits = c('Wood.density', 'SLA',
'Leaf.N', 'Max.height'),
DF.results, var.x = 'MAP',
list.params = list.params,
small.bar = 0.2, threshold.delta.AIC = 10000000,
ylim = c(-1, 1), ylim.same = FALSE, log = 'x')
models <- c('glmer.LOGLIN.ER.Tf')
names(models) <- c('ER Tf')
list.params <- list(c(Response = 'logD'))
fun.plot.panel.lmer.parameters.c(models = models,
traits = c('Wood.density', 'SLA',
'Leaf.N', 'Max.height'),
DF.results, var.x = 'MAP',
list.params = list.params,
small.bar = 0.2, threshold.delta.AIC = 10000000,
ylim = c(-1, 1), ylim.same = FALSE, log = 'x')
R/analysis/glmer.output-fun.R
View file @ 93f19645
......@@ -106,7 +106,7 @@ names(dat.t) <- c(names.param,paste(names.param,"Std.Error",sep=".")
,paste(names.param,"VAR",sep="."))
dat.t[,match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
list.res$glmer.summary$fixed.coeff.E
dat.t[,match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
dat.t[,length(names.param)+match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
list.res$glmer.summary$fixed.coeff.Std.Error
dat.t[,match(names(list.res$glmer.summary$fixed.var),names(dat.t))] <-
list.res$glmer.summary$fixed.var
......
R/analysis/glmer.run.R
View file @ 93f19645
......@@ -4,51 +4,44 @@
library(lme4)
get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
}
run.models.for.set.all.traits <- function(set,model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),type.filling, ...){
run.models.for.set.all.traits <- function(model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),
type.filling, ...){
for(trait in traits)
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, set, type.filling=type.filling, ...)
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, set,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, set,type.filling=type.filling, ...))
try(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,
type.filling, ...){
fun.model <- match.fun(fun.model)
for (e in ecoregions)
try(fun.model(model.file, trait, set, e, type.filling=type.filling,...))
try(fun.model(model.file, trait, type.filling=type.filling,...))
}
#=====================================================================
# Run glmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.glmer.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.R.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.R")
model.files.glmer.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.AD.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.R")
model.files.glmer.Tf.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
fun.test.if.multi.census <- function(data){
return("tree.id" %in% names(data))
}
model.files.glmer.Tf.1 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.R")
fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
......@@ -59,21 +52,21 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
saveRDS(glmer.output,file=file.path(path.out, "glmer.results.global.rds"))
}
run.glmer <- function (model.file, trait, set, ecoregion,
min.obs=10, sample.size=NA,
run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA,
type.filling) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
path.out <- output.dir.glmer(model$name, trait, set,
ecoregion,type.filling=type.filling)
path.out <- output.dir.glmer(model$name, trait,
type.filling=type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run glmer for model",model.file," for set",
set,"ecoregion",ecoregion,"trait",
'all',"trait",
trait,"\n")
df.glmer <- load.and.prepare.data.for.glmer(trait, set, ecoregion,
df.glmer <- load.and.prepare.data.for.glmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
cat("Ok data with Nobs",nrow(df.glmer),"\n")
......@@ -83,24 +76,36 @@ run.glmer <- function (model.file, trait, set, ecoregion,
#========================================================================
output.dir.glmer <- function (model, trait, set, ecoregion,type.filling) {
file.path("output/glmer", set,ecoregion,trait,type.filling,model)
output.dir.glmer <- function (model, trait, type.filling) {
file.path("output/glmer", 'all',trait,type.filling,model)
}
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.glmer <- function(trait, set, ecoregion,
load.and.prepare.data.for.glmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.global.csv"), stringsAsFactors = FALSE)
data.all.sample <- read.csv(file=file.path(base.dir, "data.all.csv"),
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.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,BATOT,min.obs,
min.perc.neigh=0.90,min.BA.G=-50,max.BA.G=150){
fun.select.data.for.analysis <- function(data.tree,abs.CWM.tntf,perc.neigh,
BATOT,min.obs,
min.perc.neigh=0.80){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
(!is.na(data.tree[["D"]])) )
......@@ -109,11 +114,13 @@ data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)
## select species with missing traits
data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]])),]
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh &
!is.na(data.tree[[perc.neigh]]))
# 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])
# select obs abov min perc.neigh
data.tree <- subset(data.tree,subset=data.tree[[perc.neigh]] > min.perc.neigh & !is.na(data.tree[[perc.neigh]]))
names(table(factor(data.tree[["sp"]])))[table(factor(
data.tree[["sp"]]))>min.obs])
return(data.tree)
}
......@@ -125,9 +132,11 @@ return((x-mean(x))/sd(x))
fun.get.the.variables.for.glmer.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(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="")))
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']])
#= multiply CWMs by BATOT
sumTnTfBn.abs <- data.tree[[abs.CWM.tntf]]*data.tree[[BATOT]]
sumTnBn <- data.tree[[CWM.tn]]*data.tree[[BATOT]]
......@@ -138,46 +147,26 @@ return(data.frame(dead=dead,
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))
plot.id=plot.id,
set.id = set.id,
ecocode.id = ecocode.id,
sumTnTfBn.diff= fun.center.and.standardized.var(sumTnTfBn.diff),
sumTnTfBn.abs= fun.center.and.standardized.var(sumTnTfBn.abs),
Tf= fun.center.and.standardized.var(data.tree[[tf]]),
sumTnBn=fun.center.and.standardized.var(sumTnBn),
sumTfBn=fun.center.and.standardized.var(sumTfBn),
sumBn=fun.center.and.standardized.var(sumBn)))
}
fun.get.the.variables.for.glmer.no.tree.id <- function(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf,min.BA.G=50){
dead <- data.tree[["dead"]]
logD <- fun.center.and.standardized.var(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(dead=dead,
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))
}
#============================================================
# Function to prepare data for lmer with new CWM data
# that will be used in a simple linear model
#============================================================
fun.data.for.glmer <- function(data.tree,trait,min.obs=10,type.filling='species') {
if(! trait %in% c("SLA", "Leaf.N","Seed.mass","SLA","Wood.density","Max.height"))
fun.data.for.glmer <- 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
......@@ -186,7 +175,8 @@ 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.tree <- fun.select.data.for.analysis(data.tree,abs.CWM.tntf,
perc.neigh,BATOT,min.obs)
#= DATA LIST FOR JAGS
glmer.data <- fun.get.the.variables.for.glmer.no.tree.id(data.tree,BATOT,CWM.tn,abs.CWM.tntf,tf)
return(glmer.data)
......
R/analysis/glmer.run.all.R 0 → 100644
View file @ 93f19645
###########################
###########################
### FUNCTION TO RUN GLMER ESTIMATION
library(lme4)
run.models.for.set.all.traits <- function(model.file,fun.model, traits =
c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),
type.filling, ...){
for(trait in 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, ...){
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 <- function(model.file,fun.model, trait,
type.filling, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling,...))
}
#=====================================================================
# Run glmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.glmer.Tf.1 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.E.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.R.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.HD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.Tf.R")
model.files.glmer.Tf.2 <-
c("R/analysis/model.glmer.all/model.glmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.ER.AD.Tf.R",
"R/analysis/model.glmer.all/model.glmer.LOGLIN.simplecomp.Tf.R")