Skip to content

GitLab

Commit 16b13c19 authored by Kunstler Georges's avatar Kunstler Georges
Browse files

new ecocode

parent c7f885b5
Hide whitespace changes
Inline Side-by-side
R/analysis/lmer.output-fun.R
View file @ 16b13c19
......@@ -248,7 +248,7 @@ seq.BATOT <- seq(from = 0,
return(list(trait.quant, seq.BATOT))
}
fun.generate.pred.param.dat <- function(type, list.sd, Tf.low, Tf.high, seq.sumBn, N.pred = 100){
fun.generate.pred.param.dat <- function(list.sd, Tf.low, Tf.high, seq.sumBn, N.pred = 100){
Tf.mean <- 0
D.mean <- 0
sd_sumBn <- list.sd$sd.sumBn
......@@ -269,6 +269,59 @@ fun.generate.pred.param.dat <- function(type, list.sd, Tf.low, Tf.high, seq.sum
return(df)
}
fun.generate.pred.param.kikj.dat <- function(list.sd, Tf.low,
Tf.high, seq.sumBn,
N.pred = 100){
Tf.mean <- 0
D.mean <- 0
sd_sumBn.intra <- list.sd$sd.sumBn.intra
sd_sumBn.inter <- list.sd$sd.sumBn.inter
sd_sumTfBn <- list.sd$sd.sumTfBn
sd_sumTnBn <- list.sd$sd.sumTnBn
sd_sumTnTfBn.abs <- list.sd$sd.sumTnTfBn.abs
mean.sumBn.intra <- 1#mean(seq.sumBn)
mean.sumBn.inter <- 1#mean(seq.sumBn)
print(mean.sumBn)
seq.Tf <- seq(from = Tf.low, to = Tf.high, length.out = N.pred)
df <- data.frame('logG' = rep(0 , N.pred),
'logD' = rep(D.mean, N.pred),
'Tf' =(seq.Tf-Tf.low),
'sumBn.intra' = rep(mean.sumBn.intra, N.pred)/sd_sumBn,
'sumBn.inter' = rep(mean.sumBn.inter, N.pred)/sd_sumBn,
'sumTfBn' = (seq.Tf-Tf.low)*mean.sumBn.intra/sd_sumTfBn,
'sumTnBn' = seq.Tf*mean.sumBn.intra/sd_sumTnBn,
'sumTnTfBn.abs' = abs(seq.Tf-Tf.low)*
mean.sumBn.inter/sd_sumTnTfBn.abs)
return(df)
}
fun.generate.pred.param.rho.dat <- function(list.sd, Tf.low,
Tf.high, seq.sumBn,
N.pred = 100){
Tf.mean <- 0
D.mean <- 0
sd_sumBn.intra <- list.sd$sd.sumBn.intra
sd_sumBn.inter <- list.sd$sd.sumBn.inter
sd_sumTfBn <- list.sd$sd.sumTfBn
sd_sumTnBn <- list.sd$sd.sumTnBn
sd_sumTnTfBn.abs <- list.sd$sd.sumTnTfBn.abs
mean.sumBn.intra <- 1#mean(seq.sumBn)
mean.sumBn.inter <- 1#mean(seq.sumBn)
print(mean.sumBn)
seq.Tf <- seq(from = Tf.low, to = Tf.high, length.out = N.pred)
df <- data.frame('logG' = rep(0 , N.pred),
'logD' = rep(D.mean, N.pred),
'Tf' =seq.Tf,
'sumBn.intra' = -rep(mean.sumBn.intra, N.pred)/sd_sumBn,
'sumBn.inter' = rep(mean.sumBn.inter, N.pred)/sd_sumBn,
'sumTfBn' = seq.Tf*mean.sumBn.intra/sd_sumTfBn,
'sumTnBn' = seq.Tf*mean.sumBn.intra/sd_sumTnBn,
'sumTnTfBn.abs' = abs(seq.Tf-Tf.low)*
mean.sumBn.inter/sd_sumTnTfBn.abs)
return(df)
}
fun.generate.pred.dat <- function(type, list.sd, Tf.low, Tf.high, seq.sumBn){
Tf.mean <- 0
......@@ -400,7 +453,8 @@ easyPredCI <- function(list.res, newdata, alpha=0.05) {
## BOLKER FUNCTION FROM RPUB
easyPredCI.param <- function(list.res, type, newdata, alpha=0.05) {
if (! type %in% c('maxG', 'alphar', 'alphae', 'alphal')) stop ('error in type')
if (! type %in% c('maxG', 'alphar', 'alphae', 'alphal'))
stop ('error in type')
beta <- list.res$lmer.summary$fixed.coeff.E
V <- list.res$vcov
form <- formula(list.res$terms)
......@@ -428,6 +482,35 @@ easyPredCI.param <- function(list.res, type, newdata, alpha=0.05) {
}
## BOLKER FUNCTION FROM RPUB
easyPredCI.stabl <- function(list.res, type, newdata, alpha=0.05) {
if (! type %in% c('kikj', 'rho'))
stop ('error in type')
beta <- list.res$lmer.summary$fixed.coeff.E
V <- list.res$vcov
form <- formula(list.res$terms)
## fixed-effects model matrix for new data
X <- model.matrix(form,
newdata)
sel.keep <- switch(type ,
kikj = c(2,6),
rho = c(4,5,8))
X[, -sel.keep] <- 0
pred <- X %*% beta
pred.se <- sqrt(diag(X %*% V %*% t(X))) ## std errors of predictions
## inverse-link (logistic) function: could also use plogis()
## construct 95% Normal CIs on the link scale and
## transform back to the response (probability) scale:
crit <- -qnorm(alpha/2)
if (type %in% c('alphar', 'alphae', 'alphal')) inter <- beta[4]
if (type == 'maxG') inter <- beta[1]
cbind(newdata,
pred = pred - inter,
lwr = pred-crit*pred.se - inter,
upr = pred+crit*pred.se - inter)
}
## predict for one traits
predict.for.one.traits <- function(trait, type,
dir.root, list.res,
......@@ -448,7 +531,7 @@ predict.for.one.traits.3params <- function(trait,
dir.root, list.res,
alpha = 0.05){
list.var <- get.predict.var.scaled(trait, dir.root)
new.data <- fun.generate.pred.param.dat(type, list.sd = list.res$list.sd,
new.data <- fun.generate.pred.param.dat(list.sd = list.res$list.sd,
Tf.low = list.var[[1]][['ql']],
Tf.high = list.var[[1]][['qh']],
seq.sumBn = list.var[[2]])
......@@ -483,6 +566,77 @@ return(rbind(pred.res.max, pred.res.ar, pred.res.ae, pred.res.al))
}
predict.for.one.traits.stabl.intra <- function(trait,
dir.root, list.res,
alpha = 0.05){
list.var <- get.predict.var.scaled(trait, dir.root)
new.data.rho <- fun.generate.pred.param.rho.dat( list.sd = list.res$list.sd,
Tf.low = list.var[[1]][['ql']],
Tf.high = list.var[[1]][['qh']],
seq.sumBn = list.var[[2]])
new.data.kikj <- fun.generate.pred.param.kikj.dat( list.sd = list.res$list.sd,
Tf.low = list.var[[1]][['ql']],
Tf.high = list.var[[1]][['qh']],
seq.sumBn = list.var[[2]])
pred.res.rho <- easyPredCI.stabl(list.res, type = 'rho', new.data.rho, alpha)
pred.res.rho$Tf <- seq(from = list.var[[1]][['ql.o']],
to = list.var[[1]][['qh.o']],
length.out = 100)
pred.res.rho$param.type <- 'rho'
pred.res.rho$pred <- exp(pred.res.rho$pred)
pred.res.rho$lwr <- exp(pred.res.rho$lwr)
pred.res.rho$upr <- exp(pred.res.rho$upr)
pred.res.kikj <- easyPredCI.stabl(list.res, type = 'kikj', new.data.kikj, alpha)
pred.res.kikj$Tf <- seq(from = list.var[[1]][['ql.o']],
to = list.var[[1]][['qh.o']],
length.out = 100)
pred.res.kikj$param.type <- 'kikj'
pred.res.kikj$pred <- exp(pred.res.kikj$pred)
pred.res.kikj$lwr <- exp(pred.res.kikj$lwr)
pred.res.kikj$upr <- exp(pred.res.kikj$upr)
return(rbind(pred.res.rho, pred.res.kikj))
}
predict.for.one.traits.stabl <- function(trait,
dir.root, list.res,
alpha = 0.05){
list.var <- get.predict.var.scaled(trait, dir.root)
new.data.rho <- fun.generate.pred.param.dat(list.sd = list.res$list.sd,
Tf.low = list.var[[1]][['ql']],
Tf.high = list.var[[1]][['qh']],
seq.sumBn = list.var[[2]])
new.data.kikj <- fun.generate.pred.param.kiky.dat( list.sd = list.res$list.sd,
Tf.low = list.var[[1]][['ql']],
Tf.high = list.var[[1]][['qh']],
seq.sumBn = list.var[[2]])
pred.res.rho <- easyPredCI.param(list.res, type = 'alphal', new.data.rho, alpha)
pred.res.rho$Tf <- seq(from = list.var[[1]][['ql.o']],
to = list.var[[1]][['qh.o']],
length.out = 100)
pred.res.rho$param.type <- 'alphal'
pred.res.rho$pred <- exp(pred.res.rho$pred)
pred.res.rho$lwr <- exp(pred.res.rho$lwr)
pred.res.rho$upr <- exp(pred.res.rho$upr)
pred.res.rho <- easyPredCI.stabl(list.res, type = 'rho', new.data.rho, alpha)
pred.res.rho$Tf <- seq(from = list.var[[1]][['ql.o']],
to = list.var[[1]][['qh.o']],
length.out = 100)
pred.res.rho$param.type <- 'rho'
pred.res.kikj <- easyPredCI.stabl(list.res, type = 'kikj', new.data.kikj, alpha)
pred.res.kikj$Tf <- seq(from = list.var[[1]][['ql.o']],
to = list.var[[1]][['qh.o']],
length.out = 100)
pred.res.kikj$param.type <- 'kikj'
pred.res.kikj$pred <- exp(pred.res.kikj$pred)
pred.res.kikj$lwr <- exp(pred.res.kikj$lwr)
pred.res.kikj$upr <- exp(pred.res.kikj$upr)
return(rbind(pred.res.rho, pred.res.kikj))
}
# turn pred log(BA.G) in BA.G
inv.link.BA.G <- function(pred, sd.BA.G, mean.BA.G, min.BA.G = 40){
return(exp(pred * sd.BA.G + mean.BA.G) - min.BA.G -1)
......@@ -1365,3 +1519,16 @@ fun.plot.param.tf(df = filter(data.param,
'Competitive effect ',
(alpha[e] %*% t[c]))))
}
# compute Ci of ratio based on Fieller's methods
ratio.ci <- function(mean.v, Sigma, alpha = 0.05){
ta <- -qnorm(alpha/2)
theta <- mean.v[1]/mean.v[2]
k <- ta^2*Sigma[2,2]/(theta)^2
mm <- theta+(k/(1-k))*(theta+Sigma[1,2]/Sigma[2,2])
ci <- ta/(mean.v[2]*(1-k))*sqrt(Sigma[1,1]+2*theta^2*Sigma[2,2] -
k*(Sigma[1,1]- Sigma[1,2]^2/Sigma[2,2]))
return(c(mean = theta, q.l = mm-ci, q.h = mm+ci))
}
R/analysis/lmer.run.R
View file @ 16b13c19
......@@ -337,7 +337,7 @@ plot.id <- as.character(factor(data.tree[["plot"]]))
tree.id <- as.character(factor(data.tree[["tree.id"]]))
set.id <- as.character(factor(data.tree[["set"]]))
ecocode.id <- as.character(factor(data.tree[[ecocode.var]]))
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
biomes.id <- factor(data.tree[['biomes']])
#= multiply CWMs by BATOT
......@@ -377,7 +377,7 @@ plot.id <- as.character(factor(data.tree[["plot"]]))
tree.id <- as.character(factor(data.tree[["tree.id"]]))
set.id <- as.character(factor(data.tree[["set"]]))
ecocode.id <- as.character(factor(data.tree[[ecocode.var]]))
ecocode.id <- as.character(factor(paste(data.tree[[ecocode.var]], data.tree[["set"]])))
biomes.id <- factor(data.tree[['biomes']])
#= multiply CWMs by BATOT
......
R/analysis/ratio.CI.R
View file @ 16b13c19
##
library(MASS)
mvrnorm <-
function(n = 1, mu, Sigma, tol=1e-6, empirical = FALSE, EISPACK = FALSE)
{
p <- length(mu)
if(!all(dim(Sigma) == c(p,p))) stop("incompatible arguments")
if(EISPACK) stop("'EISPACK' is no longer supported by R", domain = NA)
eS <- eigen(Sigma, symmetric = TRUE)
ev <- eS$values
if(!all(ev >= -tol*abs(ev[1L]))) stop("'Sigma' is not positive definite")
X <- matrix(rnorm(p * n), n)
if(empirical) {
X <- scale(X, TRUE, FALSE) # remove means
X <- X %*% svd(X, nu = 0)$v # rotate to PCs
X <- scale(X, FALSE, TRUE) # rescale PCs to unit variance
}
X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*% t(X)
nm <- names(mu)
if(is.null(nm) && !is.null(dn <- dimnames(Sigma))) nm <- dn[[1L]]
dimnames(X) <- list(nm, NULL)
if(n == 1) drop(X) else t(X)
}
## define param
mean.v <- c(0.15, 0.25)
Sigma <- matrix(c(0.002,0.0001,0.0001,0.005),2,2)
# MC
res.mc <- mvrnorm(n = 10000, mean.v, Sigma)
ratio.mc <- res.mc[,1]/res.mc[,2]
par(mfrow = c(2,2))
hist(res.mc[,1], breaks = 100)
abline(v = quantile(res.mc[, 1], probs = c(0.025, 0.975)))
hist(res.mc[,2], breaks = 100)
abline(v = quantile(res.mc[, 2], probs = c(0.025, 0.975)))
hist(ratio.mc, breaks = 1000, xlim = c(0, 1.5))
abline(v = quantile(ratio.mc, probs = c(0.025, 0.975)))
abline(v = mean(res.mc[, 1])/mean(res.mc[, 2]), col = 'red')
# compute Ci of ratio based on Fieller's methods
ratio.ci <- function(mean.v, Sigma, alpha = 0.05){
......@@ -48,11 +10,5 @@ k <- ta^2*Sigma[2,2]/(theta)^2
mm <- theta+(k/(1-k))*(theta+Sigma[1,2]/Sigma[2,2])
ci <- ta/(mean.v[2]*(1-k))*sqrt(Sigma[1,1]+2*theta^2*Sigma[2,2] -
k*(Sigma[1,1]- Sigma[1,2]^2/Sigma[2,2]))
return(c(q.l = mm-ci, q.h = mm+ci))
}
abline(v = ratio.ci(mean.v, Sigma), col = 'green')
quantile(ratio.mc, probs = c(0.025, 0.975))
ratio.ci(mean.v, Sigma)
launch.cluster/launch_all_lmer.bash
View file @ 16b13c19
......@@ -64,8 +64,8 @@ samplesize=$1
# # # ecocode 3
# echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO${trait}.sh
# qsub Rscript_temp/allECO${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[1], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO${trait}.sh
qsub Rscript_temp/allECO${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO${trait}" -q opt32G -j oe
echo "/usr/local/R/R-3.1.1/bin/Rscript -e \"source('R/analysis/lmer.run.R'); run.multiple.model.for.set.one.trait(model.files.lmer.Tf.2[2], run.lmer,merge.biomes.TF = TRUE,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO2${trait}.sh
qsub Rscript_temp/allECO2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.ECO2${trait}" -q opt32G -j oe
......
scripts/analysis/results.fig.R
View file @ 16b13c19
......@@ -76,6 +76,8 @@ plot.param.mean.and.biomes.fixed(list.all.results.intra , data.type = "intra",
add.param.descrip.TF =2)
dev.off()
pdf('figs/figres12.TP.pdf', height = 14, width = 16)
plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
models = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species',
......@@ -97,9 +99,53 @@ plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
dev.off()
pdf('figs/figres12.TP.intra.pdf', height = 14, width = 16)
plot.param.mean.and.biomes.fixed(list.all.results.intra , data.type = "intra",
models = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.species',
'lmer.LOGLIN.ER.AD.Tf.MAT.MAP.intra.r.set.fixed.biomes.species'),
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn.intra", "sumBn.inter", "Tf"),
param.print = 1:6,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha['0 intra'])),
expression('Trait indep'(alpha['0 inter'])),
expression("Direct trait "(m[1])),
expression("Size "(gamma %*% log('D'))) ) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.28, 0.31),
add.param.descrip.TF =2)
dev.off()
pdf('figs/figres12.ecocode.pdf', height = 14, width = 16)
plot.param.mean.and.biomes.fixed(list.all.results.set , data.type = "simple",
models = c('lmer.LOGLIN.ER.AD.Tf.r.ecocode.species',
'lmer.LOGLIN.ER.AD.Tf.r.ecocode.fixed.biomes.species'),
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
param.print = 1:5,
param.names = c(expression('Trait sim '(alpha['s'])),
expression('Tolerance '(alpha['t'])),
expression('Effect '(alpha['e'])),
expression('Trait indep'(alpha[0])),
expression("Direct trait "(m[1])),
expression("Size "(gamma %*% log('D'))) ) ,
col.vec = fun.col.pch.biomes()$col.vec,
pch.vec = fun.col.pch.biomes()$pch.vec,
names.bio = names.biomes ,
xlim = c(-0.28, 0.27))
dev.off()
pdf('figs/figres1.ecocode.pdf', height = 14, width = 16)
plot.param(list.all.results.set , data.type = 'simple',
model = c('lmer.LOGLIN.ER.AD.Tf.MAT.MAP.r.set.species'),
model = c('lmer.LOGLIN.ER.AD.Tf.r.ecocode.species'),
traits = c('Wood.density' , 'SLA', 'Max.height'),
param.vec = c("sumTnTfBn.abs", "sumTfBn","sumTnBn",
"sumBn", "Tf"),
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment