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

first version of departemental talk

parent 88c42fb8
......@@ -8,6 +8,7 @@ figs
*.xlsx
.Rhistory
*.pdf
*.png
*.doc
*.docx
*.html
......
......@@ -11,7 +11,6 @@ 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,
model.files.lmer.Tf.5,
model.files.lmer.Tf.CAT.1)){
source(model, local = TRUE)
......@@ -26,7 +25,8 @@ for (trait in traits){
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
out <- lapply(files, summarise.lmer.output.all.list)
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
......@@ -34,31 +34,31 @@ names(out) <- lapply(lapply(files,files.details.all),
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.ecocode.id.no.log.rds')
## LOOP OVER FILES AND NOT SCAN ALL FILES
traits <- c("SLA", "Wood.density","Max.height")
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)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
type.filling=type.filling)
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
## traits <- c("SLA", "Wood.density","Max.height")
## 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)){
## source(model, local = TRUE)
## model.obj <- load.model()
## pathout <- output.dir('lmer', model.obj$name, trait, 'all.no.log',
## type.filling=type.filling)
files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
## files <- c(files,file.path(pathout,"species.id.results.nolog.all.rds"))
}
}
## }
## }
out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
names(out) <- lapply(lapply(files,files.details.all),
function(x) paste(as.vector(x[names(x) != 'file']),
collapse="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.species.id.no.log.rds')
## out <- lapply(files, summarise.lmer.output.all.list, random.name = 'biomes.id')
## names(out) <- lapply(lapply(files,files.details.all),
## function(x) paste(as.vector(x[names(x) != 'file']),
## collapse="_"))
## ### remove missing
## out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
## saveRDS(out,file='output/list.lmer.out.all.species.id.no.log.rds')
## ## LOOP OVER FILES AND NOT SCAN ALL FILES
......
......@@ -21,7 +21,7 @@ summarise.lmer.output <- function(x){
fixed.var = variance.fixed.glmm.lmer(x))
}
summarise.lmer.all.output <- function(x, random.name = 'set.id'){
summarise.lmer.all.output <- function(x, random.name ){
list( nobs = nobs(x),
R2m = Rsquared.glmm.lmer(x)$Marginal,
R2c = Rsquared.glmm.lmer(x)$Conditional,
......@@ -49,12 +49,19 @@ summarise.lmer.output.list <- function(f ){
}
summarise.lmer.output.all.list <- function(f, random.name = 'set.id' ){
summarise.lmer.output.all.list <- function(f ){
output.lmer <- read.lmer.output(f)
if(!is.null(output.lmer)){
res <- list(files.details = files.details.all(f),
details <- files.details.all(f)
source(file.path('R','analysis', 'model.lmer',
paste('model',details$model,'R', sep = '.')))
model <- load.model()
res <- list(files.details = details,
lmer.summary = summarise.lmer.all.output( output.lmer,
random.name ),
random.name =
model$var.BLUP),
terms = terms(output.lmer),
vcov = vcov(output.lmer))
}else{
......@@ -82,6 +89,16 @@ files.details.all <- function(x){
s[-(1:2)]
}
BLUP.CI <- function(fm, var.BLUP){ ## NOT WORKING WHEN MULTIPLE term per factor
cV <- ranef(fm, condVar = TRUE, whichel = var.BLUP)
ranvar <- attr(cV[[1]], "postVar")
res <- do.call('rbind',
lapply(1:dim(ranvar)[3],
function(i, array) sqrt(diag(array[ , , i])),
array = ranvar))
colnames(res) <- names(cV[[1]])
return(res)
}
#### R squared functions
......@@ -325,6 +342,12 @@ segments( x - small.bar, unlist(y - 1.96*sd), x + small.bar , unlist(y -1.96*
segments( x - small.bar, unlist(y + 1.96*sd), x + small.bar, unlist(y +1.96*sd))
}
fun.plot.error.bar.horiz <- function(x, y, sd, small.bar = (max(x) - min(x))/40){
segments( unlist(x - 1.96*sd), y, unlist(x +1.96*sd), y)
segments( unlist(x - 1.96*sd), y-small.bar, unlist(x -1.96*sd), y+small.bar)
segments( unlist(x + 1.96*sd), y-small.bar, unlist(x +1.96*sd), y+small.bar)
}
fun.plot.param.error.bar <- function(df.selected, var.x, param, small.bar, col.vec){
segments( df.selected[[var.x]], df.selected[[param]] - 1.96*df.selected[[paste(
......@@ -664,8 +687,8 @@ mm <- model.matrix(list.output$terms, newdat)
mm2 <- mm
mm2[, !colnames(mm2) %in% c(param, param.var)] <- 0
pred <- (mm2 %*% list.output$lmer.summary$fixed.coeff.E)/mean(mm[, param])
pvar1 <- diag(mm2 %*% tcrossprod(list.output$vcov, mm))
pred <- (mm2 %*% list.output$lmer.summary$fixed.coeff.E)#/mean(mm[, param])
pvar1 <- diag(mm2 %*% tcrossprod(list.output$vcov, mm2))
dat <- data.frame(pred = pred,
plo = pred-1.96*sqrt(pvar1),
phi = pred+1.96*sqrt(pvar1))
......@@ -691,7 +714,7 @@ list.res.MAT <- list(logG = 0,
Tf = 0,
sumTnBn = 0,
sumTfBn = 0,
sumBn = 0)
sumBn = 0)
list.res.MAP <- list(logG = 0,
logD = 0,
......@@ -707,7 +730,7 @@ list.res.MAP <- list(logG = 0,
Tf = 0,
sumTnBn = 0,
sumTfBn = 0,
sumBn = 0)
sumBn = 0)
new.data.MAT <- expand.grid(list.res.MAT)
......@@ -717,8 +740,10 @@ return(list.new.data)
}
plot.fun.ci.param.var.climate <- function(list.output, trait, param.vec, var.vec,
n.length = 100, add.mean.effect = FALSE, ...){
plot.fun.ci.param.var.climate <- function(list.output, param.vec,
var.vec,
n.length = 100,
add.mean.effect = FALSE, ...){
list.new.data <- fun.generate.pred.dat(n.length = 100)
par(mfrow=c(length(var.vec), length(param.vec)))
for (var in var.vec){
......@@ -729,7 +754,7 @@ for (var in var.vec){
dat.pred <- fun.predict.ci.lmer(list.output, newdat,
param, param.var)
plot(newdat[[var]], dat.pred$pred, type = 'l',
ylim = c(-1,1), xlab = var, ylab = param, ...)
ylim = c(-0.4,0.4), xlab = var, ylab = param, ...)
lines(newdat[[var]], dat.pred$plo, lty = 2)
lines(newdat[[var]], dat.pred$phi, lty = 2)
abline(h = 0, lwd = 2)
......@@ -748,10 +773,32 @@ for (var in var.vec){
}
## # test
## fm.MAT.MAP <- readRDS('output/lmer/all.no.log/Wood.density/species/lmer.LOGLIN.ER.AD.Tf.MAT.MAP/results.nolog.all.rds')
## fun.ci.param.var.climate(fm.MAT.MAP, 'Wood.density', param = 'sumTfBn', var = 'MAT')
## BOLKER FUNCTION FROM RPUB
easyPredCI <- function(model,newdata,alpha=0.05) {
## baseline prediction, on the linear predictor (logit) scale:
pred0 <- predict(model,re.form=NA,newdata=newdata)
## fixed-effects model matrix for new data
X <- model.matrix(formula(model,fixed.only=TRUE)[-2],
newdata)
beta <- fixef(model) ## fixed-effects coefficients
V <- vcov(model) ## variance-covariance matrix of beta
pred.se <- sqrt(diag(X %*% V %*% t(X))) ## std errors of predictions
## inverse-link (logistic) function: could also use plogis()
linkinv <- model@resp$family$linkinv
## construct 95% Normal CIs on the link scale and
## transform back to the response (probability) scale:
crit <- -qnorm(alpha/2)
linkinv(cbind(lwr=pred0-crit*pred.se,
upr=pred0+crit*pred.se))
}
### plot with random effect on top
......@@ -764,13 +811,17 @@ plot.param.BLUP <- function(list.res,
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("logD", "Tf","sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs"),
param.names = c('Size', 'Direct trait',
'Compet', 'Compet effect x trait',
'Compet response x trait',
'Compet x trait dissimilarity'),
col.vec,
pch.vec){
m <- matrix(c(1:4,5,5), 2, 3)
layout(m, widths=c(4, 4, 1),
heights=c(4, 4))
pch.vec,
names.bio){
m <- matrix(c(1:4), 1, 4)
layout(m, widths=c(7, 4, 4, 3 )/10,
heights=c(4)/10)
for (i in traits){
## screen((1:length(traits))[traits == i])
list.temp <- list.res[[paste("all.no.log_", i ,
"_species_", model,
sep = '')]]$lmer.summary
......@@ -779,44 +830,54 @@ for (i in traits){
names(param.std) <- names(list.temp$fixed.coeff.E)
param.std <- param.std[param.vec]
param.BLUP <- list.temp$set.BLUP
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-0.45, 0.25))
abline(h = 0)
axis(1, 1:length(param.vec), labels = param.vec, las = 3)
fun.plot.error.bar(1:length(param.vec),
param.mean,
if(i == traits[1]) {
par(mar=c(6,13,3,0))
}else{
par(mar=c(6,0,3,0))
}
plot(param.mean, 1:length(param.vec),
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5)
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
if(i == traits[1]) axis(2, 1:length(param.vec), labels = param.names, las = 1, cex.axis = 2,
mgp = c(1.5,0.55,0))
fun.plot.error.bar.horiz(param.mean,
1:length(param.vec),
param.std)
for (j in 1:length(param.vec)){
points(rep(j,nrow(param.BLUP)), param.BLUP[,param.vec[j]],
points(param.BLUP[,param.vec[j]], rep(j,nrow(param.BLUP)),
pch = pch.vec[rownames(param.BLUP)],
col = col.vec[rownames(param.BLUP)],
cex= 2)
cex= 1)
}
}
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names(col.vec), col = col.vec, pch = pch.vec,
bty = 'n')
legend("center",legend = names.bio, col = col.vec, pch = pch.vec,
bty = 'n', cex = 1)
}
##
plot.param <- function(list.res,
model = 'lmer.LOGLIN.ER.AD.Tf',
traits = c('Wood.density', 'SLA', 'Max.height'),
param.vec = c("logD", "Tf","sumBn", "sumTnBn",
"sumTfBn", "sumTnTfBn.abs"),
param.names = c('Size', 'Direct trait',
'Compet', 'Compet effect x trait',
'Compet response x trait',
'Compet x trait dissimilarity'),
col.vec,
pch.vec){
m <- matrix(c(1:4), 2, 2)
layout(m, widths=c(4, 4),
heights=c(4, 4))
pch.vec,
names.bio){
m <- matrix(c(1:3), 1, 3)
layout(m, widths=c(7, 4, 4 )/10,
heights=c(4)/10)
for (i in traits){
## screen((1:length(traits))[traits == i])
list.temp <- list.res[[paste("all.no.log_", i ,
"_species_", model,
sep = '')]]$lmer.summary
......@@ -824,17 +885,73 @@ for (i in traits){
param.std <- list.temp$fixed.coeff.Std.Error
names(param.std) <- names(list.temp$fixed.coeff.E)
param.std <- param.std[param.vec]
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-0.45, 0.25))
abline(h = 0)
axis(1, 1:length(param.vec), labels = param.vec, las = 3)
fun.plot.error.bar(1:length(param.vec),
param.mean,
param.BLUP <- list.temp$set.BLUP
if(i == traits[1]) {
par(mar=c(6,13,3,0))
}else{
par(mar=c(6,0,3,0))
}
plot(param.mean, 1:length(param.vec),
yaxt = 'n', xlab = paste('Effect size '), ylab = NA,
xlim = c(-0.45, 0.25) , pch = 16 , cex = 2, cex.lab = 1.5)
mtext(i, side=3, cex =1.2)
box(lwd= 2)
abline(v = 0)
if(i == traits[1]) axis(2, 1:length(param.vec), labels = param.names, las = 1, cex.axis = 2,
mgp = c(1.5,0.55,0))
fun.plot.error.bar.horiz(param.mean,
1:length(param.vec),
param.std)
}
}
plot.param.BLUP2 <- function(list.res, model = 'lmer.LOGLIN.ER.AD.Tf',
traits = c('Wood.density', 'SLA', 'Max.height'), param.vec = c("logD",
"Tf","sumBn", "sumTnBn", "sumTfBn", "sumTnTfBn.abs"), col.vec,
pch.vec){
m <- matrix(c(1:4,5,5), 2, 3)
layout(m, widths=c(4, 4, 1),
heights=c(4, 4))
for (i in traits){
list.temp <-
list.res[[paste("all.no.log_", i , "_species_", model, sep =
'')]]$lmer.summary
param.mean2 <- list.temp$fixed.coeff.E[param.vec]
param.mean <- rep(0, length(param.mean2))
param.std <-
list.temp$fixed.coeff.Std.Error
names(param.std) <-
names(list.temp$fixed.coeff.E)
param.std <- param.std[param.vec]
param.BLUP <- list.temp$set.BLUP
plot(param.mean, xaxt = 'n', ylab =
'std parameters', xlab = NA, main = i, ylim =
range(t(param.BLUP[names(pch.vec[!is.na(pch.vec)]),param.vec ]) -
param.mean2))
abline(h = 0)
axis(1, 1:length(param.vec), labels =
param.vec, las = 3)
for (j in 1:length(param.vec)){
points(rep(j,nrow(param.BLUP)), param.BLUP[,param.vec[j]] -
param.mean2[j], pch = pch.vec[rownames(param.BLUP)], col =
col.vec[rownames(param.BLUP)], cex= 1)
}
}
par(mar=c(0, 0, 0, 0))
plot(1,1,type="n", axes=F, xlab="", ylab="")
legend("center",legend = names(col.vec)[!is.na(pch.vec)],
col = col.vec[!is.na(pch.vec)],
pch = pch.vec[!is.na(pch.vec)],
bty = 'n')
}
##
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.MAT.MAP.r.species",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+MAT+MAP+Tf.A_EV+Tf.A_EV:MAT+Tf.A_EV:MAP+Tf.A_D+Tf.A_D:MAT+Tf.A_D:MAP+Tf.C+Tf.C:MAT+Tf.C:MAP+logD+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn.A_EV+sumTfBn.A_EV:MAT+sumTfBn.A_EV:MAP+sumTfBn.A_D+sumTfBn.A_D:MAT+sumTfBn.A_D:MAP+sumTfBn.C+sumTfBn.C:MAT+sumTfBn.C:MAP+sumTnBn.A_EV+sumTnBn.A_EV:MAT+sumTnBn.A_EV:MAP+sumTnBn.A_D+sumTnBn.A_D:MAT+sumTnBn.A_D:MAP+sumTnBn.C+sumTnBn.C:MAT+sumTnBn.C:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.MAT.MAP.r.species",
var.BLUP = 'species.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+MAT+MAP+Tf.A_EV+Tf.A_EV:MAT+Tf.A_EV:MAP+Tf.A_D+Tf.A_D:MAT+Tf.A_D:MAP+Tf.C+Tf.C:MAT+Tf.C:MAP+logD+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn.A_EV+sumTfBn.A_EV:MAT+sumTfBn.A_EV:MAP+sumTfBn.A_D+sumTfBn.A_D:MAT+sumTfBn.A_D:MAP+sumTfBn.C+sumTfBn.C:MAT+sumTfBn.C:MAP+sumTnBn.A_EV+sumTnBn.A_EV:MAT+sumTnBn.A_EV:MAP+sumTnBn.A_D+sumTnBn.A_D:MAT+sumTnBn.A_D:MAP+sumTnBn.C+sumTnBn.C:MAT+sumTnBn.C:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+Tf.A_EV+Tf.A_D+Tf.C+logD+sumBn+sumTfBn.A_EV+sumTfBn.A_D+sumTfBn.C+sumTnBn.A_EV+sumTnBn.A_D+sumTnBn.C+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|biomes.id)+(sumBn-1|biomes.id)+(sumTfBn.A_EV-1|biomes.id)+(sumTfBn.A_D-1|biomes.id)+(sumTfBn.C-1|biomes.id)+(sumTnBn.A_EV-1|biomes.id)+(sumTnBn.A_D-1|biomes.id)+(sumTnBn.C-1|biomes.id)+(sumTnTfBn.abs-1|biomes.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.r.species",
var.BLUP = 'species.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+Tf.A_EV+Tf.A_D+Tf.C+logD+sumBn+sumTfBn.A_EV+sumTfBn.A_D+sumTfBn.C+sumTnBn.A_EV+sumTnBn.A_D+sumTnBn.C+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|plot.id)+Tf+Tf:MAT+Tf:MAP+MAT+MAP+logD+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.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.MAT.MAP.r.species",
var.BLUP = 'species.id',
lmer.formula.tree.id=formula("logG~1+(1|plot.id)+Tf+Tf:MAT+Tf:MAP+MAT+MAP+logD+sumBn+sumBn:MAT+sumBn:MAP+sumTfBn+sumTfBn:MAT+sumTfBn:MAP+sumTnBn+sumTnBn:MAT+sumTnBn:MAP+sumTnTfBn.abs+sumTnTfBn.abs:MAT+sumTnTfBn.abs:MAP+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)+(MAT-1|species.id)+(MAP-1|species.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.biomes.species",
var.BLUP = 'biomes.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.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.ecocode.species",
var.BLUP = 'ecocode.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(1|species.id) +(logD-1|species.id) +(sumBn-1|species.id)+(Tf-1|ecocode.id)+(sumBn-1|ecocode.id)+(sumTfBn-1|ecocode.id)+(sumTnBn-1|ecocode.id)+(sumTnTfBn.abs-1|ecocode.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.r.species",
var.BLUP = 'species.id',
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(1|species.id)+(logD-1|species.id)+(sumBn-1|species.id)"))
}
......
##################### function to process data install all unstallled packages
source("R/packages.R")
git.root <- function() {
system("git rev-parse --show-toplevel", intern=TRUE)
}
source.root <- function(filename) {
source(file.path(git.root(), filename), chdir=TRUE)
}
path.root <- git.root()
source.root("R/packages.R")
check_packages(c("reshape", "data.table", "doParallel", "data.table","Rcpp"))
......@@ -306,7 +316,8 @@ fun.CWM.traits.all.plot.census <- function(census, obs.id, plot, diam, sp,
##############################################
##### FUNCTIONS FOR XY DATA
library(Rcpp)
sourceCpp("R/process.data/georges.cpp")
sourceCpp(file.path(path.root,
"R/process.data/georges.cpp"))
## function to compute local BA based on Xy and compute CWM for each traits
fun.CWM.traits.xy <- function(i, obs.id, traits.list, sp.num, sp.names,
sp.length, xy.table, BA.a, Rlim){
......
################################
## function for test.tree.CWM
source("R/utils/plot.R")
## source("R/utils/plot.R")
load.processed.data <- function(path, file.name = "data.tree.tot.no.std.csv"){
require(data.table)
fname <- file.path(path, file.name )
if(file.exists(fname)){
## data.all.sample <- read.csv(file = fname,
## stringsAsFactors = FALSE, nrows = 5000)
## classes <- sapply(data.all.sample, class)
## classes[classes=='integer'] <- "character"
## classes[classes=='logical'] <- "character"
## nrows <- as.numeric(system(paste('wc -l < ',fname),
## intern = TRUE))
## data <- read.csv(file = fname,
## stringsAsFactors = FALSE,
## nrows = nrows,
## colClasses = classes)
## ## data <- read.csv(fname, stringsAsFactors = FALSE)
## ## print(classes)
data <- fread(fname, stringsAsFactors = FALSE)
return(data)
}else{return(NULL)}
......@@ -658,4 +644,80 @@ return(sss)
predict.lmer.ci <- function(var.y, data){
require(lme4)
fm1 <- lmer(formula = eval(parse(text = paste(var.y,' ~shannon*biomes +(1|set)'))),
data = data)
newdat <- expand.grid( 0,
shannon=seq(0, 3.5, length.out = 20)
, biomes= levels(data$biomes))
names(newdat) <- c(var.y, 'shannon', 'biomes')
mm <- model.matrix(terms(fm1),newdat)
newdat[[var.y]] <- mm %*% fixef(fm1)
pvar1 <- diag(mm %*% tcrossprod(vcov(fm1),mm))
tvar1 <- pvar1+VarCorr(fm1)$Subject[1]
newdat <- data.frame(
newdat
, plo = newdat[[var.y]]-2*sqrt(pvar1)
, phi = newdat[[var.y]]+2*sqrt(pvar1)
)
return(newdat)
}
plot.line.by <- function(level, DF, col.vec, var.y,
var.x = 'shannon', cat = 'biomes', ...){
lines(DF[[var.x]][DF[[cat]] == level],
DF[[var.y]][DF[[cat]] == level],
col = col.vec[level], ...)
}
plot.poly.by <- function(level, DF, col.vec,
var.x = 'shannon', cat = 'biomes', ...){
polygon(c(DF[[var.x]][DF[[cat]] == level], rev(DF[[var.x]][DF[[cat]] == level])),
c(DF$plo[DF[[cat]] == level], rev(DF$phi[DF[[cat]] == level])) ,
col = adjustcolor(col.vec[level], alpha.f = 0.5),
...)
}
pred.shannon.vs.x.ci.lines <- function(DF, var, ...){
pred.dat <- predict.lmer.ci(var.y= var, data = DF)
plot(DF$shannon, DF[[var]], cex = 0, ...)
lapply(levels(DF$biomes), plot.line.by, DF= pred.dat,
col.vec = fun.col.pch.biomes()$col.vec[levels(DF$biomes)],
var.y = var, lwd= 3)
lapply(levels(DF$biomes), plot.poly.by, DF= pred.dat,
col.vec = fun.col.pch.biomes()$col.vec[levels(DF$biomes)],