Commit 143a6ccf authored by Georges Kunstler's avatar Georges Kunstler
Browse files

added model with no log BA

parent 67ee0df5
#### function to analyse lmer output
library(lme4)
read.lmer.output <- function(file.name){
tryCatch(readRDS(file.name), error=function(cond)return(NULL)) # Choose a return value in case of error
}
summarise.lmer.output <- function(x){
list( nobs = nobs(x),
R2m =Rsquared.glmm.lmer(x)$Marginal,
R2c =Rsquared.glmm.lmer(x)$Conditional,
AIC = AIC(x),
deviance = deviance(x),
conv=x@optinfo$conv,
effect.response.var=variance.fixed.glmm.lmer.effect.and.response(x),
fixed.coeff.E=fixef(x),
fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
fixed.var=variance.fixed.glmm.lmer(x))
}
summarise.lmer.output.list <- function(f ){
output.lmer <- read.lmer.output(f)
if(!is.null(output.lmer)){
res <- list(files.details=files.details(f),
lmer.summary=summarise.lmer.output( output.lmer))
}else{
res <- NULL
}
return(res)
}
files.details <- function(x){
s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
names(s) <- c("d1", "d2", "set", "ecocode", "trait", "filling", "model", "file" )
s[-(1:2)]
}
#### R squred functions
# Function rsquared.glmm requires models to be input as a list (can include fixed-
# effects only models,but not a good idea to mix models of class "mer" with models
# of class "lme") FROM http://jslefche.wordpress.com/2013/03/13/r2-for-linear-mixed-effects-models/
Rsquared.glmm.lmer <- function(i){
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
# Calculate marginal R-squared (fixed effects/total variance)
Rm <- VarF/(VarF+VarRand+VarResid)
# Calculate conditional R-squared (fixed effects+random effects/total variance)
Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid)
# Bind R^2s into a matrix and return with AIC values
Rsquared.mat <- data.frame(Class=class(i),Family="Gaussian",Marginal=Rm,
Conditional=Rc,AIC=AIC(i))
return(Rsquared.mat)
}
variance.fixed.glmm.lmer <- function(i){
# Get variance of for each fixed effects by multiplying coefficients by design matrix
var.vec <- apply(fixef(i) * t(i@pp$X),MARGIN=1,var)
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
names(var.vec) <- paste(names(var.vec),"VAR",sep=".")
return(var.vec)
}
variance.fixed.glmm.lmer.effect.and.response <- function(i){
if(sum(c("sumTfBn","sumTnBn") %in% names(fixef(i)))==2){
# Get variance of for each fixed effects by multiplying coefficients by design matrix
var.vec <- var(as.vector(fixef(i)[c("sumTfBn","sumTnBn")] %*% t(i@pp$X[,c("sumTfBn","sumTnBn")])))
# Get variance of fixed effects by multiplying coefficients by design matrix
VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
# Get variance of random effects by extracting variance components
VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
}else{
var.vec <- NA
}
names(var.vec) <- paste("effect.response","VAR",sep=".")
return(var.vec)
}
## function to turn lmer output from list to DF
fun.format.in.data.frame <- function(list.res,names.param){
dat.t <- data.frame(t(rep(NA,3*length(names.param))))
names(dat.t) <- c(names.param,paste(names.param,"Std.Error",sep=".")
,paste(names.param,"VAR",sep="."))
dat.t[,match(names(list.res$lmer.summary$fixed.coeff.E),names(dat.t))] <-
list.res$lmer.summary$fixed.coeff.E
dat.t[,length(names.param)+match(names(list.res$lmer.summary$fixed.coeff.E),names(dat.t))] <-
list.res$lmer.summary$fixed.coeff.Std.Error
dat.t[,match(names(list.res$lmer.summary$fixed.var),names(dat.t))] <-
list.res$lmer.summary$fixed.var
res <- data.frame(list.res$files.details,list.res$lmer.summary[1:7],dat.t)
return(res)
}
#!/usr/bin/env Rscript
source("R/analysis/lmer.nolog.output-fun.R")
source("R/analysis/lmer.nolog.run.R")
## TODO NEED TO CHANGE THAT TO LOOP OVER FILES AND NOT SCAN ALL FILES
sets <- c('BCI','Canada','France','Fushan','NVS','Paracou',
'Spain','US','Swiss','Sweden','NSW','Mbaiki','Luquillo','Japan')
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (set in sets){
ecoregions <- get.ecoregions.for.set(set)
for (ecoregion in ecoregions){
for (trait in traits){
for (model in c(model.files.lmer.Tf.1,model.files.lmer.Tf.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model=model.obj$name, trait, set, ecoregion,type.filling)
files <- c(files,file.path(pathout,"results.no.std.nolog.rds"))
}
}
}
}
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/lmer.list.out.rds')
names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
saveRDS(DF.results,file='output/lmer.DF.rds')
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
files <- list.files("output/lmer", recursive=TRUE, full.names =TRUE, pattern = "rds")
source("R/analysis/lmer.run.R")
## TODO NEED TO CHANGE THAT TO LOOP OVER FILES AND NOT SCAN ALL FILES
sets <- c('BCI','Canada','France','Fushan','NVS','Paracou',
'Spain','US','Swiss','Sweden','NSW','Mbaiki','Luquillo','Japan')
traits <- c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass")
type.filling <- 'species'
files <- c()
for (set in sets){
ecoregions <- get.ecoregions.for.set(set)
for (ecoregion in ecoregions){
for (trait in traits){
for (model in c(model.files.lmer.Tf.1,model.files.lmer.Tf.2)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model=model.obj$name, trait, set, ecoregion,type.filling)
files <- c(files,file.path(pathout,"results.no.std.rds"))
}
}
}
}
out <- lapply(files, summarise.lmer.output.list)
names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
saveRDS(out,file='output/lmer.list.out.rds')
......
......@@ -53,6 +53,12 @@ table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$trait,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set)
t(table(DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$best.model,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set))/(apply(table(DF.best.and.all.AIC[
DF.best.and.all.AIC$filling=='species',]$best.model,
DF.best.and.all.AIC[DF.best.and.all.AIC$filling=='species',]$set),
MARGIN=2,sum))
## AIC weights
AIC.weights <- do.call('rbind',
......@@ -63,6 +69,20 @@ DF.AIC.weights <- data.frame(DF.best.and.all.AICc[,1],AIC.weights)
names(DF.AIC.weights) <- c('id2',paste('AIC.weight',names(DF.AIC.weights)[-1],sep='.'))
DF.best.and.all.AIC <- merge(DF.best.and.all.AIC,DF.AIC.weights,by='id2')
## assume ER best model if either E or R best model
f <- function(i,DF){
d.AIC <- as.vector(DF$delta.AIC[DF$id2==i] )
best <- DF$model[DF$id2==i][DF$delta.AIC[DF$id2==i]==0]
if(best %in% c('lmer.LOGLIN.E','lmer.LOGLIN.R')) {
d.AIC[ DF$model[DF$id2==i]=='lmer.LOGLIN.ER'] <- 0}
return(d.AIC)
}
d.AIC <- do.call('c',lapply(unique(DF.results$id2),FUN=f,DF.results))
DF.results$delta.AIC <- d.AIC
## compute a global AIC
fun.global.aic <- function(DF.results){
DF.results <- DF.results[DF.results$nobs>1000,] # select set ecocode with more than 1000 obs
......@@ -136,14 +156,14 @@ fun.plot.panel.lmer.res.x.y(models=models,
var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,0.17))
dev.off()
models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
names(models) <- c('Effect/response','Absolute distance')
var.y.l <- list('ER.perc.var','abs.perc.var')
pdf('figs/perc.var.relative.BATOT.MAP.pdf',width=12,height=8)
fun.plot.panel.lmer.res.x.y(models=models,
traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,100))
dev.off()
## models <- c('lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
## names(models) <- c('Effect/response','Absolute distance')
## var.y.l <- list('ER.perc.var','abs.perc.var')
## pdf('figs/perc.var.relative.BATOT.MAP.pdf',width=12,height=8)
## fun.plot.panel.lmer.res.x.y(models=models,
## traits=c('Wood.density','SLA','Leaf.N','Max.height'),DF.results,
## var.x='MAP',var.y.l=var.y.l,ylim=c(-0.015,100))
## dev.off()
......
......@@ -50,13 +50,16 @@ fun.test.if.multi.census <- function(data){
return("tree.id" %in% names(data))
}
fun.call.lmer <- function(formula,df.lmer){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
return(lmer.output)
}
fun.call.lmer.and.save <- function(formula,df.lmer,path.out){
Start <- Sys.time()
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
end <- Sys.time()
print(end - Start)
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.no.std.rds"))
saveRDS(lmer.output,file=file.path(path.out, "results.rds"))
}
run.lmer <- function (model.file, trait, set, ecoregion,
......@@ -82,10 +85,37 @@ run.lmer <- function (model.file, trait, set, ecoregion,
#= Run model
if(test.if.multi.census){
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out)
fun.ci.boot(df.lmer,formula=model$lmer.formula.tree.id,path.out,level=0.95,nsim=500)
}else{
fun.call.lmer.and.save(formula=model$lmer.formula,df.lmer,path.out)
fun.ci.boot(df.lmer,formula=model$lmer.formula,path.out,level=0.95,nsim=500)
}
}
## new function to compute boot ci
fun.ci.boot <- function(df.lmer,formula,path.out,level=0.95,nsim=500){
require(boot)
require(multicore)
bb <- boot(data=df.lmer, statistic=boot.fun,R= nsim,formula=formula)
bci <- lapply(seq_along(bb$t0), boot.out = bb, boot::boot.ci,
type = "perc", conf = level)
citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
a <- (1 - level)/2
a <- c(a, 1 - a)
pct <- paste('CI',round(a, 3),sep='.')
dimnames(citab) <- list(names(bb[["t0"]]), pct)
saveRDS(citab,file=file.path(path.out, "results.ci.rds"))
}
boot.fun <- function(data, indices, formula){
df.lmer <- data[indices,] # select obs. in bootstrap sample
res <- fun.call.lmer(formula=formula,df.lmer)
fixef(res)
}
#========================================================================
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
......@@ -98,7 +128,7 @@ output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
#============================================================
load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
min.obs, sample.size, type.filling,
base.dir = "output/processed/",std=TRUE){
base.dir = "output/processed/",std=FALSE){
### load data
if(std) { data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"),
stringsAsFactors = FALSE)}else{
......
###########################
###########################
### FUNCTION TO RUN LMER ESTIMATION
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, std,...){
for(trait in traits)
run.multiple.model.for.set.one.trait(model.file,fun.model, trait, set, type.filling=type.filling, std=std,...)
}
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), std, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait, set,type.filling=type.filling,std=std, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.filling, ecoregions = get.ecoregions.for.set(set), std, ...){
fun.model <- match.fun(fun.model)
for (e in ecoregions)
try(fun.model(model.file, trait, set, e, type.filling=type.filling,std=std,...))
}
#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.lmer.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.R.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.R")
model.files.lmer.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R")
model.files.lmer.Tf.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.R.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R")
model.files.lmer.Tf.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.HD.Tf.R",
"R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R")
fun.test.if.multi.census <- function(data){
return("tree.id" %in% names(data))
}
fun.call.lmer <- function(formula,df.lmer){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
return(lmer.output)
}
fun.call.lmer.and.save <- function(formula,df.lmer,path.out,std){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE)
print(summary(lmer.output))
if(std) { saveRDS(lmer.output,file=file.path(path.out, "results.no.std.nolog.rds"))
}else{saveRDS(lmer.output,file=file.path(path.out, "results.nolog.rds"))
}
}
run.lmer <- function (model.file, trait, set, ecoregion,
min.obs=10, sample.size=NA,
type.filling,std) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
#= Path for output
path.out <- output.dir.lmer(model$name, trait, set,
ecoregion,type.filling=type.filling)
print(path.out)
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run lmer for model",model.file," for set",
set,"ecoregion",ecoregion,"trait",
trait,"\n")
df.lmer <- load.and.prepare.data.for.lmer(trait, set, ecoregion,
min.obs, sample.size,
type.filling=type.filling,std=std) # return a DF
test.if.multi.census <- fun.test.if.multi.census(df.lmer)
cat("Ok data with Nobs",nrow(df.lmer),
"multiple census", test.if.multi.census ,"\n")
#= Run model
if(test.if.multi.census){
fun.call.lmer.and.save(formula=model$lmer.formula.tree.id,df.lmer,path.out,std=std)
}else{
fun.call.lmer.and.save(formula=model$lmer.formula,df.lmer,path.out,std=std)
}
}
## ## new function to compute boot ci
## fun.ci.boot <- function(df.lmer,formula,path.out,level=0.95,nsim=500){
## require(boot)
## require(multicore)
## bb <- boot(data=df.lmer, statistic=boot.fun,R= nsim,formula=formula)
## bci <- lapply(seq_along(bb$t0), boot.out = bb, boot::boot.ci,
## type = "perc", conf = level)
## citab <- t(sapply(bci, function(x) x[["percent"]][4:5]))
## a <- (1 - level)/2
## a <- c(a, 1 - a)
## pct <- paste('CI',round(a, 3),sep='.')
## dimnames(citab) <- list(names(bb[["t0"]]), pct)
## saveRDS(citab,file=file.path(path.out, "results.ci.no.std.rds"))
## }
## boot.fun <- function(data, indices, formula){
## df.lmer <- data[indices,] # select obs. in bootstrap sample
## res <- fun.call.lmer(formula=formula,df.lmer)
## fixef(res)
## }
#========================================================================
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
file.path("output/lmer", set,ecoregion,trait,type.filling,model)
}
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait, set, ecoregion,
min.obs, sample.size, type.filling,
base.dir = "output/processed/",std){
### load data
if(std) { data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.no.std.csv"),
stringsAsFactors = FALSE)}else{
data.tree.tot <- read.csv(file.path(base.dir, set,ecoregion,"data.tree.tot.csv"),
stringsAsFactors = FALSE)}
fun.data.for.lmer(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){
## remove tree with NA
data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
(!is.na(data.tree[["D"]])) )
## 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 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]]))
return(data.tree)
}
fun.center.and.standardized.var <- function(x){
return((x-mean(x))/sd(x))
}
### 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){
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))
}
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
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))
}
#============================================================
# 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 ")
# 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)
}
return(lmer.data)
}
##### SCRIPT TO TEST plots resuls with regression lines
source("R/analysis/lmer.run.R")
output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
file.path("output/lmer", set,ecoregion,trait,type.filling,model)
}
fun.compute.resid <- function(i){
return(fitted(i) +resid(i) -i@pp$X %*%fixef(i))
}
fun.plot.colors.density <- function(x,y,...){
mat <- cbind(x,y)
data <- as.data.frame(mat)
......@@ -28,12 +19,15 @@ colors.dens <- densCols(mat)
points(x,y, col=colors.dens, pch=20,...)
}
fun.boxplot.breaks <- function(x,y,Nclass=15,...){
fun.boxplot.breaks <- function(x,y,Nclass=15,add=TRUE,...){
x1 <- x[!is.na(x) & !is.na(y)]
y1 <- y[!is.na(x) & !is.na(y)]
x <- x1