Commit 2454c63b authored by Georges Kunstler's avatar Georges Kunstler
Browse files

progress on cat models

parent c5de53b5
#######################################
#######################################
###### EXPLORE DATA SET BEFORE ANALYSIS
rm(list = ls())
source("R/utils/mem.R")
source("R/analysis/lmer.run.nolog.all.R")
filedir <- "output/processed"
traits <- c('Wood.density', "SLA", 'Leaf.N', 'Max.height')
for (t in traits){
print(t)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE)
print(1)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, cat.TF=TRUE)
print(2)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, fname = 'data.all.csv')
print(3)
load.and.save.data.for.lmer('Wood.density', data.table.TF = TRUE, cat.TF=TRUE, fname = 'data.all.csv')
print(4)
}
......@@ -7,17 +7,18 @@ library(lme4)
run.multiple.model.for.set.one.trait <- function(model.files,fun.model, trait,
type.filling, cat.TF, ...){
type.filling, fname = 'data.all.no.std.csv',
cat.TF = FALSE, ...){
for (m in model.files)
try(run.model.for.set.one.trait (m, fun.model,trait,
type.filling=type.filling, cat.TF = cat.TF, ...))
type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
run.model.for.set.one.trait <- function(model.file,fun.model, trait,
type.filling, cat.TF, ...){
type.filling, fname, cat.TF, ...){
fun.model <- match.fun(fun.model)
try(fun.model(model.file, trait, type.filling=type.filling, cat.TF = cat.TF, ...))
try(fun.model(model.file, trait, type.filling=type.filling, fname = fname, cat.TF = cat.TF, ...))
}
......@@ -62,7 +63,7 @@ fun.call.glmer.and.save <- function(formula,df.lmer,path.out){
run.glmer <- function (model.file, trait,
min.obs=5, sample.size=NA,
type.filling, fname = 'data.all.no.std.csv', cat.TF) {
type.filling, fname , cat.TF) {
require(lme4)
source(model.file, local = TRUE)
model <- load.model()
......
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.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[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
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,"results.nolog.all.rds"))
}
}
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="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.no.log.rds')
#!/usr/bin/env Rscript
source("R/analysis/lmer.output-fun.R")
source("R/analysis/lmer.run.nolog.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[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model.obj$name, trait, 'all.log',
type.filling=type.filling)
files <- c(files,file.path(pathout,"results.nolog.all.rds"))
}
}
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="_"))
### remove missing
out <- out[!unlist(lapply(out,FUN=function(x) is.null(x$lmer.summary)))]
saveRDS(out,file='output/list.lmer.out.all.rds')
......@@ -8,8 +8,9 @@ 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[-4])){
for(model in c(model.files.lmer.Tf.1[3], model.files.lmer.Tf.2[3],
model.files.lmer.Tf.3[3], model.files.lmer.Tf.4[3],
model.files.lmer.Tf.5,model.files.lmer.Tf.6)){
source(model, local = TRUE)
model.obj <- load.model()
pathout <- output.dir.lmer(model.obj$name, trait, 'all',
......
......@@ -21,6 +21,8 @@ list.all.results <- readRDS('output/list.lmer.out.no.log.all.temp.rds')
list.all.results <- readRDS('output/list.lmer.out.nolog.all.old.rds')
list.all.results <- readRDS('output/list.lmer.out.nolog.all.rds')
list.all.results <- readRDS('output/list.lmer.out.all.no.log.rds')
names.param <- unique(unlist(lapply(list.all.results,
function(x) names(x$lmer.summary$fixed.coeff.E))))
names.param2 <- names.param
......@@ -47,23 +49,40 @@ print(pander(res))
}
## plot the parameters
## add set random effect estim
pdf('figs/param.global.model.ER.AD.pdf')
par(mfrow = c(2,2))
for (i in c('Wood.density', 'SLA', 'Leaf.N', 'Max.height')){
DF.t <- DF.global[DF.global$trait == i, ]
param.vec <- c("logD", "Tf","sumBn", "sumTnBn","sumTfBn", "sumTnTfBn.abs")
list.temp <- list.all.results[[paste("all.no.log_", i ,
"_species_lmer.LOGLIN.ER.AD.Tf",
sep = '')]]$lmer.summary
param.mean <- list.temp$fixed.coeff.E[param.vec]
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
print(i)
param.vec <- c("Tf", "logD","sumBn", "sumTnBn","sumTfBn", "sumTnTfBn.abs")
param.vec.std <- paste(param.vec, 'Std.Error', sep = '.')
plot(unlist(DF.t[DF.t$model == "lmer.LOGLIN.ER.AD.Tf", param.vec]),
plot(param.mean,
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-1, 1))
abline(h = 0)
axis(1, 1:length(param.vec), labels = param.vec, las = 3)
fun.plot.error.bar(1:length(param.vec),
DF.t[DF.t$model == "lmer.LOGLIN.ER.AD.Tf", param.vec],
DF.t[DF.t$model == "lmer.LOGLIN.ER.AD.Tf", param.vec.std])
param.mean,
param.std)
for (j in 1:length(param.vec)){
points(rep(j,nrow(param.BLUP)), param.BLUP[,param.vec[j]],
pch = pch.vec[rownames(param.BLUP)],
col = col.vec[rownames(param.BLUP)])
}
}
legend("topright",legend = names(col.vec), col = col.vec, pch = pch.vec,
ncol =2, bty = 'n')
dev.off()
......@@ -72,7 +91,7 @@ par(mfrow = c(2,2))
for (i in c('Wood.density', 'SLA', 'Leaf.N', 'Max.height')){
DF.t <- DF.global[DF.global$trait == i, ]
print(i)
param.vec <- c("Tf", 'MAT.Tf', 'MAP.Tf', 'MAT', 'MAP', "logD","sumBn", 'MAT.sumBn', 'MAP.sumBn',
param.vec <- c('MAT', 'MAP', "logD", "Tf", 'MAT.Tf', 'MAP.Tf', "sumBn", 'MAT.sumBn', 'MAP.sumBn',
"sumTnBn", 'MAT.sumTnBn', 'MAP.sumTnBn', "sumTfBn", 'MAT.sumTfBn',
'MAP.sumTfBn', "sumTnTfBn.abs" , 'MAT.sumTnTfBn.abs',
'MAP.sumTnTfBn.abs')
......@@ -94,12 +113,12 @@ par(mfrow = c(2,2))
for (i in c('Wood.density', 'SLA', 'Leaf.N', 'Max.height')){
DF.t <- DF.global[DF.global$trait == i, ]
print(i)
param.vec <- c("Tf.A_EV", "Tf.A_D", 'Tf.C',"logD","sumBn",
param.vec <- c("logD","Tf.A_EV", "Tf.A_D", 'Tf.C',"sumBn",
"sumTnBn.A_EV", "sumTnBn.A_D", "sumTnBn.C",
"sumTfBn.A_EV", "sumTfBn.A_D","sumTfBn.C",
"sumTnTfBn.abs.A_EV" , "sumTnTfBn.abs.A_D", "sumTnTfBn.abs.C")
param.vec.std <- paste(param.vec, 'Std.Error', sep = '.')
mode.selected <- 'lmer.LOGLIN.ER.AD.Tf.norandom'
mode.selected <- 'lmer.LOGLIN.ER.AD.Tf.CAT.norandom'
plot(unlist(DF.t[DF.t$model == mode.selected, param.vec]),
xaxt = 'n', ylab = 'std parameters', xlab = NA,
main = i, ylim = c(-1, 1))
......
......@@ -57,11 +57,14 @@ model.files.lmer.Tf.4 <-
"R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.MAT.MAP.R")
model.files.lmer.Tf.5 <-
c("R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.norandom.R")
c("R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.norandom.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.MAT.MAP.norandom.R")
model.files.lmer.Tf.6 <-
c("R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.CAT.norandom.R")
c("R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.CAT.norandom.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.CAT.MAT.MAP.norandom.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.AD.Tf.CAT.R")
fun.call.lmer.and.save <- function(formula, df.lmer, path.out){
lmer.output <- lmer(formula = formula, data = df.lmer, REML = FALSE,
......@@ -87,7 +90,7 @@ run.lmer <- function (model.file, trait,
trait, "\n")
df.lmer <- load.and.prepare.data.for.lmer(trait,fname = fname, cat.TF = cat.TF)
# return a DF
cat("Ok data with Nobs", nrow(df.lmer),
"\n")
print(mem())
......@@ -177,8 +180,8 @@ select.temp <-(1:nrow(data.tree))[!is.na(data.tree[["BA.G"]]) &
data.tree[[perc.neigh.gs]] >
min.perc.neigh.gs &
!is.na(data.tree[[perc.neigh.gs]])]
## remove tree with NA
data.tree <- data.tree[select.temp ,]
# select species with minimum obs
......@@ -252,11 +255,9 @@ ecocode.id <- factor(data.tree[['ecocode']])
#get the three cwm per cat
vec.CWM.tn <- paste(CWM.tn, c('A_EV', 'A_D', 'C'), sep = '.')
vec.abs.CWM.tntf <- paste(abs.CWM.tntf, c('A_EV', 'A_D', 'C'), sep = '.')
vec.abs.CWM.tntf <- abs.CWM.tntf
sumTnTfBn.abs.A_EV<- data.tree[[vec.abs.CWM.tntf[1]]]
sumTnTfBn.abs.A_D <- data.tree[[vec.abs.CWM.tntf[2]]]
sumTnTfBn.abs.C <- data.tree[[vec.abs.CWM.tntf[3]]]
sumTnTfBn.abs<- data.tree[[vec.abs.CWM.tntf]]
sumTnBn.A_EV<- data.tree[[vec.CWM.tn[1]]]
sumTnBn.A_D <- data.tree[[vec.CWM.tn[2]]]
......@@ -267,9 +268,9 @@ sumTfBn.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *sumTfBn
sumTfBn.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *sumTfBn
sumTfBn.C <- as.numeric(data.tree[['SLA.cat']] == 3) *sumTfBn
sumBn <- data.tree[[BATOT]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
Tf.A_EV <- as.numeric(data.tree[['SLA.cat']] == 1) *data.tree[[tf]]
Tf.A_D <- as.numeric(data.tree[['SLA.cat']] == 2) *data.tree[[tf]]
Tf.C <- as.numeric(data.tree[['SLA.cat']] == 3) *data.tree[[tf]]
return(
data.frame(logG = logG,
logD = logD,
......@@ -280,12 +281,8 @@ data.frame(logG = logG,
set.id = set.id,
ecocode.id = ecocode.id,
plot.id = plot.id,
sumTnTfBn.abs.A_EV =
fun.center.and.standardized.var(sumTnTfBn.abs.A_EV),
sumTnTfBn.abs.A_D =
fun.center.and.standardized.var(sumTnTfBn.abs.A_D),
sumTnTfBn.abs.C =
fun.center.and.standardized.var(sumTnTfBn.abs.C),
sumTnTfBn.abs =
fun.center.and.standardized.var(sumTnTfBn.abs),
Tf.A_EV= fun.center.and.standardized.var(Tf.A_EV),
Tf.A_D= fun.center.and.standardized.var(Tf.A_D),
Tf.C= fun.center.and.standardized.var(Tf.C),
......
load.model <- function(){
list(
name = "LOGLIN.ER.AD.Tf",
## inits function
initial.values = function(){
list("intercept.trait"=runif(1,min=-1,max=1),
"slope.trait.e"=runif(1,min=-2,max=2),
"slope.trait.r"=runif(1,min=-2,max=2),
"mean.B"=runif(1,min=-0.5,max=1),
"sigma.B"=runif(1,min=0.0001,max=1.5),
"sigma"=runif(1,min=0.0001,max=4),
"sigma.alpha"=runif(1,min=0.0001,max=4),
"sigma.plot"=runif(1,min=0.0001,max=4))},
## parameters to save
parameters.to.save =c("alpha","B1",
"slope.trait.e","slope.trait.r","intercept.trait",
"sigma","sigma.plot","mean.B","sigma.B",
"sigma.alpha"),
bug =
"#######################################################
######################## growth model with JAGS #######
model {
############ Likelihood ###################
for (i in 1:N.indiv) {
logG[i] ~ dnorm(theo.g[i],tau)
theo.g[i] <- alpha[species.id[i]] + B1[species.id[i]]*logD[i] +B2[set.id[i]]*sumBn +random.plot[plot.species.id[i]]
}
################################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N.species)
{
B1[n] ~ dnorm(mean.B[1],tau.B[1])
alpha[n] ~ dnorm(0,tau.alpha)
for ( j in 1:Ncompet.sp)
{
param.compet[n,j] <- slope.trait.e*trait.compet[j] +slope.trait.r*trait.target[n] +intercept.trait
}
}
## plot effect
for (j in 1:N.plot.species)
{
random.plot[j] ~ dnorm(0,tau.plot)
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
tau0 <- 1.0E-4
# variance error
tau <- pow(sigma,-2)
sigma ~ dunif(0.00001,5)
tau.plot <- pow(sigma.plot,-2)
sigma.plot ~ dunif(0.00001,5)
## traits parameters
slope.trait.e ~ dnorm(0,tau0)
slope.trait.r ~ dnorm(0,tau0)
intercept.trait ~ dnorm(0,tau0)
## prior for hyper parameter mean and sd
for (l in 1:1)
{
mean.B[l]~ dnorm(0,tau0)
tau.B[l] <- pow(sigma.B[1],-2)
sigma.B[l] ~ dunif(0.000001,5)
}
tau.alpha <- pow(sigma.alpha,-2)
sigma.alpha ~ dunif(0.000001,5)
} # End of the jags model
")
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.MAT.MAP.norandom",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.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+(logD-1|species.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.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+(logD-1|species.id)+(Tf.A_EV-1|set.id)+(Tf.A_D-1|set.id)+(Tf.C-1|set.id)+(sumBn-1|set.id)+(sumTfBn.A_EV-1|set.id)+(sumTfBn.A_D-1|set.id)+(sumTfBn.C-1|set.id)+(sumTnBn.A_EV-1|set.id)+(sumTnBn.A_C-1|set.id)+(sumTnBn.C-1|set.id)+(sumTnTfBn.abs-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.CAT.norandom",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.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.A_EV+sumTnTfBn.abs.A_D+sumTnTfBn.abs.C+(logD-1|species.id)"))
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.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+(logD-1|species.id)"))
}
......
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.MAT.MAP.norandom",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(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+(logD-1|species.id)+(Tf-1|set.id)"))
}
load.model <- function () {
list(name="lmer.LOGLIN.ER.AD.Tf.norandom",
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumTfBn+sumBn+sumTnBn+sumTnTfBn.abs+(logD-1|species.id)"))
lmer.formula.tree.id=formula("logG~1+(1|set.id)+(1|species.id)+(1|plot.id)+Tf+logD+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(logD-1|species.id)"))
}
......
load.model <- function(){
list(
name = "LOGLIN.RE",
## parameters to save
parameters.to.save <- c('intercept','mean_logD','mean_Tf','mean_sumBn',
parameters.to.save <- c('intercept','mean_logD','mean_Tf',
'mean_sumBn',
'mean_sumTfBn','mean_sumTnBn','sigma_inter_sp',
'sigma_inter_ecocode','sigma_inter_plot','sigma_logD_sp',
'sigma_Tf_ecocode','sigma_sumBn_ecocode',
'sigma_sumTfBn_ecocode','sigma_sumTnBn_ecocode',
'sigma_inter_set','sigma_inter_plot','sigma_logD_sp',
'sigma_Tf_set','sigma_sumBn_set',
'sigma_sumTfBn_set','sigma_sumTnBn_set',
'sigma',
'param_Tf_ecocode',
'param_sumBn_ecocode','param_sumTfBn_ecocode',
'param_sumTnBn_ecocode','intercept_ecocode')
'param_Tf_set',
'param_sumBn_set','param_sumTfBn_set',
'param_sumTnBn_set','intercept_set')
,
stan = "
......@@ -21,11 +22,11 @@ data {
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_ecocode;
int<lower=0> N_set;
int species_id[N_indiv];
int plot_id[N_indiv];
int ecocode_id[N_indiv];
int set_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
......@@ -43,37 +44,37 @@ parameters{
real<lower=-10,upper=10> mean_sumTfBn;
real<lower=-10,upper=10> mean_sumTnBn;
real<lower=0.0001,upper=5> sigma_inter_sp;
real<lower=0.0001,upper=5> sigma_inter_ecocode;
real<lower=0.0001,upper=5> sigma_inter_set;
real<lower=0.0001,upper=5> sigma_inter_plot;
real<lower=0.0001,upper=5> sigma_logD_sp;
real<lower=0.0001,upper=5> sigma_Tf_ecocode;
real<lower=0.0001,upper=5> sigma_sumBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTfBn_ecocode;
real<lower=0.0001,upper=5> sigma_sumTnBn_ecocode;
real<lower=0.0001,upper=5> sigma_Tf_set;
real<lower=0.0001,upper=5> sigma_sumBn_set;
real<lower=0.0001,upper=5> sigma_sumTfBn_set;
real<lower=0.0001,upper=5> sigma_sumTnBn_set;
real<lower=0.0001,upper=5> sigma;
vector[N_species] param_logD_sp;
vector[N_species] intercept_sp;
vector[N_plot] intercept_plot;
vector[N_ecocode] param_Tf_ecocode;
vector[N_ecocode] param_sumBn_ecocode;
vector[N_ecocode] param_sumTfBn_ecocode;
vector[N_ecocode] param_sumTnBn_ecocode;
vector[N_ecocode] intercept_ecocode;
vector[N_set] param_Tf_set;
vector[N_set] param_sumBn_set;
vector[N_set] param_sumTfBn_set;
vector[N_set] param_sumTnBn_set;
vector[N_set] intercept_set;
}
transformed parameters {
# vector for prediction
vector[N_indiv] theo_g;
for (i in 1:N_indiv){
theo_g[i] <- intercept + intercept_sp[species_id[i]]
theo_g[i] <- intercept + intercept_sp[species_id[i]]
+ intercept_plot[plot_id[i]]
+ intercept_ecocode[ecocode_id[i]]
+ param_logD_sp[species_id[i]]*logD[i]
+ param_Tf_ecocode[ecocode_id[i]]*Tf[i]
+ param_sumBn_ecocode[ecocode_id[i]]*sumBn[i]
+ param_sumTfBn_ecocode[ecocode_id[i]]*sumTfBn[i]
+ param_sumTnBn_ecocode[ecocode_id[i]]*sumTnBn[i]
;
+ intercept_set[set_id[i]]
+ param_logD_sp[species_id[i]]*logD[i]
+ param_Tf_set[set_id[i]]*Tf[i]
+ param_sumBn_set[set_id[i]]*sumBn[i]
+ param_sumTfBn_set[set_id[i]]*sumTfBn[i]
+ param_sumTnBn_set[set_id[i]]*sumTnBn[i]
;
}
}
......@@ -91,12 +92,12 @@ model {
intercept_sp ~ normal(0,sigma_inter_sp);
### plot random param
intercept_plot ~ normal(0,sigma_inter_plot);
### ecocode random param
param_Tf_ecocode ~ normal(mean_Tf,sigma_Tf_ecocode);
param_sumBn_ecocode ~ normal(mean_sumBn,sigma_sumBn_ecocode);
param_sumTfBn_ecocode ~ normal(mean_sumTfBn,sigma_sumTfBn_ecocode);
param_sumTnBn_ecocode ~ normal(mean_sumTnBn,sigma_sumTnBn_ecocode);
intercept_ecocode ~ normal(0,sigma_inter_ecocode);
### set random param
param_Tf_set ~ normal(mean_Tf,sigma_Tf_set);
param_sumBn_set ~ normal(mean_sumBn,sigma_sumBn_set);
param_sumTfBn_set ~ normal(mean_sumTfBn,sigma_sumTfBn_set);
param_sumTnBn_set ~ normal(mean_sumTnBn,sigma_sumTnBn_set);
intercept_set ~ normal(0,sigma_inter_set);
###############################################
########### Non-Hierarchical parameters ########
......@@ -111,19 +112,19 @@ model {
mean_sumTnBn ~ normal(0,sigma0);
# sigma
sigma_inter_sp~ uniform(0,6);
sigma_inter_ecocode~ uniform(0,6);
sigma_inter_set~ uniform(0,6);
sigma_inter_plot~ uniform(0,6);
sigma_logD_sp~ uniform(0,6);
sigma_Tf_ecocode~ uniform(0,6);
sigma_sumBn_ecocode~ uniform(0,6);
sigma_sumTfBn_ecocode~ uniform(0,6);
sigma_sumTnBn_ecocode~ uniform(0,6);
sigma_Tf_set~ uniform(0,6);
sigma_sumBn_set~ uniform(0,6);
sigma_sumTfBn_set~ uniform(0,6);
sigma_sumTnBn_set~ uniform(0,6);
sigma~ uniform(0,6);
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
")
")
}
library(lme4)
packageVersion("lme4")
dat <- readRDS('output/processed/data.glmer.SLA.no.log.no.cat.rds')
dat$logyear2 <- dat$logyear/mean(dat$logyear)
hist(log(exp(dat$logyear)/10))
range(exp(dat$logyear))
table(dat$set[exp(dat$logyear)<1])
hist(dat$logyear/mean(dat$logyear))
form <-"dead~1+offset(logyear2)+Tf+logD+(logD|species.id)+(1|plot.id)+(1|set.id)+sumBn+sumTfBn+sumTnBn+sumTnTfBn.abs+(Tf-1|set.id)+(sumTfBn-1|set.id)"
## +(sumBn-1|set.id)+(sumTnBn-1|set.id)+(sumTnTfBn.abs-1|set.id)
glmer.test <- glmer(form,data = dat[sample(1:nrow(dat), 10000), ], family=binomial(link = 'cloglog'))
start.v <- getME(glmer.test,c("theta","beta"))
names(start.v) <- c('theta', 'fixef')
glmer.test2 <- update(glmer.test,start=start.v)
glmer.NM <- update(glmer.test,control=glmerControl(optimizer="Nelder_Mead"))