Commit 2a0ea9ec authored by kunstler's avatar kunstler
Browse files

compare with jags

parent b5f8f7c4
......@@ -9,8 +9,6 @@ source('R/analysis/lmer.run.R')
model.files.jags.Tf.1 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.species.R",
"R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.species.fixed.biomes.R")
model.files.jags.Tf.2 <-
c("R/analysis/model.jags/model.LOGLIN.ER.AD.Tf.r.set.sepcies.tree.R")
......
......@@ -78,7 +78,8 @@ run.lmer <- function (model.file, trait,
data.type ='simple',
var.sample = NA,
select.biome = NA,
merge.biomes.TF = FALSE) {
merge.biomes.TF = FALSE,
sample.vec.TF = FALSE) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
......@@ -95,7 +96,8 @@ run.lmer <- function (model.file, trait,
ecocode.var. = ecocode.var,
var.sample. = var.sample,
select.biome. = select.biome,
merge.biomes.TF = merge.biomes.TF)
merge.biomes.TF = merge.biomes.TF,
sample.vec.TF. = sample.vec.TF)
# return a list of DF
cat("Ok data with Nobs", nrow(list.df.lmer[[1]]),
"\n")
......
......@@ -55,11 +55,6 @@ for (j in 1:N_plot)
for (j in 1:N_set)
{
intercept_set[j] ~ dnorm(0,tau_inter_set)
}
## set
for (j in 1:N_set)
{
param_Tf_set[j] ~ dnorm(0, tau_Tf_set)T(-5, 5)
param_sumBn_set[j] ~ dnorm(0, tau_sumBn_set)T(-5, 5)
param_sumTfBn_set[j] ~ dnorm(0 ,tau_sumTfBn_set)T(-5, 5)
......
##
qt(p= 0.025, df = 100000)
qnorm(0.025)
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.02,0.01,0.01,0.05),2,2)
# MC
res.mc <- mvrnorm(n = 10000, mean.v, Sigma)
ratio.mc <- res.mc[,1]/res.mc[,2]
hist(ratio.mc, breaks = 100)
lines(v = quantile(ratio.mc, probs = c(0.025, 0.975)))
## Fieller's theorem
##because lot of data approximate teh student by a normal
ta <- qnorm(0.025)
g <- mean.v[1]*mean.v[2] - ta*Sigma[1,2]
d <- mean.v[2]^2 - ta*Sigma[2,2]
b <- mean.v[1]^2 - ta*Sigma[1,1]
ci.b <- sqrt(g^2-d*b)
......@@ -54,20 +54,20 @@ samplesize=$1
# 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.intra.1[2], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRA2${trait}.sh
# qsub Rscript_temp/allINTRA2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRA2${trait}" -q opt32G -j oe
# # INTRA 2
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.intra.2[1], run.lmer,$trait, data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB${trait}.sh
qsub Rscript_temp/allINTRAB${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB${trait}" -q opt32G -j oe
# # # INTRA 2
# 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.intra.2[1], run.lmer,$trait, data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB${trait}.sh
# qsub Rscript_temp/allINTRAB${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB${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.intra.2[2], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB2${trait}.sh
qsub Rscript_temp/allINTRAB2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB2${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.intra.2[2], merge.biomes.TF = TRUE, run.lmer,$trait,data.type = 'intra');print('done')\"" > Rscript_temp/allINTRAB2${trait}.sh
# qsub Rscript_temp/allINTRAB2${trait}.sh -d ~/trait.competition.workshop -l nodes=1:ppn=1,mem=8gb -N "lmerall2all.INTRAB2${trait}" -q opt32G -j oe
# # # 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[2], run.lmer,$trait, data.type = 'simple');print('done')\"" > Rscript_temp/allECO2${trait}.sh
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
......
......@@ -7,15 +7,16 @@ source("R/analysis/lmer.run.R")
## ## save list of all output for NA
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.simple.set.rds',
models = c(model.files.lmer.Tf.0, model.files.lmer.Tf.1,
model.files.lmer.Tf.3, model.files.lmer.Tf.4),
traits = c("SLA", "Wood.density", "Max.height"))
## format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
## list.file.name = 'list.lmer.out.all.NA.simple.set.rds',
## models = c(model.files.lmer.Tf.0, model.files.lmer.Tf.1,
## model.files.lmer.Tf.2, model.files.lmer.Tf.3),
## traits = c("SLA", "Wood.density", "Max.height"))
format.all.output.lmer(file.name = "NA.wwf.results.nolog.all.rds",
list.file.name = 'list.lmer.out.all.NA.intra.set.rds',
models = c(model.files.lmer.Tf.2),
models = c(model.files.lmer.Tf.intra.1,
model.files.lmer.Tf.intra.2),
traits = c("SLA", "Wood.density", "Max.height"),
data.type = 'intra')
......
......@@ -7,9 +7,6 @@ df <- gather(data.europe, country, density, -age.c)
library(dplyr)
# plot age dist
library(ggplot2)
theme_simple <- function(){
......@@ -38,6 +35,3 @@ ggplot(df, aes(x = age.c, y = density)) +
theme_simple() + facet_wrap(~ country)+ xlab("age (yr.)")
dev.off()
x11()
ggplot(data.age, aes(x = age.max)) + geom_histogram(aes(y = ..density..), colour = "black", fill="grey") + geom_density(, col = 'red')+ scale_x_continuous(limits=c(1, 300))+ theme_simple() + ggtitle("France") + xlab("age (yr.)")
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