Skip to content

GitLab

Commit 8baf852a authored by Georges Kunstler's avatar Georges Kunstler
Browse files

lmer global running

parent 5d4ce8a5
Hide whitespace changes
Inline Side-by-side
R/analysis/lmer.run.nolog.all.R
View file @ 8baf852a
......@@ -26,12 +26,6 @@ run.model.for.set.one.trait <- function(model.file,fun.model, trait, set,type.fi
#=====================================================================
# Run lmer() (in package lme4) for one ecoregion, one trait, one model
#=====================================================================
model.files.lmer.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.R.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.ER.R")
model.files.lmer.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.AD.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.R")
model.files.lmer.Tf.1 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.E.Tf.R",
"R/analysis/model.lmer.all/model.lmer.LOGLIN.R.Tf.R",
......@@ -42,20 +36,15 @@ model.files.lmer.Tf.2 <- c("R/analysis/model.lmer.all/model.lmer.LOGLIN.nocomp.T
"R/analysis/model.lmer.all/model.lmer.LOGLIN.simplecomp.Tf.R")
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){
lmer.output <- lmer(formula=formula,data=df.lmer,REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=20000) ) )
control=lmerControl(optCtrl=list(maxfun=40000) ) )
print(summary(lmer.output))
saveRDS(lmer.output,file=file.path(path.out, "results.nolog.all.rds"))
}
run.lmer <- function (model.file, trait,
run.lmer <- function (model.file, trait,
min.obs=10, sample.size=NA,
type.filling) {
require(lme4)
......@@ -68,7 +57,7 @@ run.lmer <- function (model.file, trait,
dir.create(path.out, recursive=TRUE, showWarnings=FALSE)
cat("run lmer for model",model.file," for trait",
trait,"\n")
df.lmer <- load.and.prepare.data.for.lmer(trait,
df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
cat("Ok data with Nobs",nrow(df.lmer),
......@@ -88,14 +77,24 @@ output.dir.lmer <- function (model, trait, set,type.filling) {
#============================================================
# Function to prepare data for lmer
#============================================================
load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
load.and.prepare.data.for.lmer <- function(trait,
min.obs, sample.size, type.filling,
base.dir = "output/processed/"){
### load data
library(data.table)
## library(data.table)
data.tree.tot <- fread(file.path(base.dir,"data.all.csv"),
stringsAsFactors = FALSE)
## data.tree.tot <- fread(file.path(base.dir,"data.all.csv"),
## stringsAsFactors = FALSE)
data.all.sample <- read.csv(file=file.path(base.dir, "data.all.csv"),
stringsAsFactors=FALSE,nrows=1000)
classes <- sapply(data.all.sample,class)
classes[classes=='integer'] <- "numeric"
nrows <- as.numeric(system('wc -l < output/processed/data.all.csv',intern=TRUE))
data.tree.tot <- read.csv(file=file.path(base.dir, "data.all.csv"),
stringsAsFactors=FALSE,
nrows=nrows,
colClasses=classes)
fun.data.for.lmer(data.tree.tot,trait,type.filling=type.filling)
}
......@@ -108,15 +107,15 @@ data.tree <- subset(data.tree,subset=(!is.na(data.tree[["BA.G"]])) &
## 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)
data.tree[["D"]]>9.9)
## select species with no missing traits
select.temp <- (!is.na(data.tree[[abs.CWM.tntf]]) &
!is.na(data.tree[[BATOT]]))
data.tree <- data.tree[select.temp,]
# select obs abov min perc.neigh
# 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]]))
# select species with minimum obs
# 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])
......@@ -136,7 +135,7 @@ logD <- fun.center.and.standardized.var(log(data.tree[["D"]]))
species.id <- factor(data.tree[["sp"]])
tree.id <- factor(data.tree[["tree.id"]])
plot.id <- factor(data.tree[["plot"]])
set.id <- factor(data.tree[["set"]])
ecocode.id <- factor(data.tree[['ecocode']])
......
R/analysis/run.local.R
View file @ 8baf852a
......@@ -13,196 +13,9 @@ source("R/analysis/lmer.run.nolog.all.R")
type.filling=type.filling) # return a DF
lmer.output.s <- lmer(formula=model$lmer.formula.tree.id,data=df.lmer[sample(1:nrow(df.lmer),200000), ],
REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=10000) ) )
lmer.output <- lmer(formula=model$lmer.formula.tree.id,data=df.lmer,
REML = FALSE,
control=lmerControl(optCtrl=list(maxfun=10000) ) )
lmer.output.s@optinfo$conv
control=lmerControl(optCtrl=list(maxfun=30000) ) )
lmer.output@optinfo$conv
##### SCRIPT TO TEST BEFORE TO SEND ON CLUSTER
source("R/analysis/lmer.run.nolog.R")
sets <- c('France','Spain','Sweden','Swiss','BCI','Fushan','Luquillo','NVS','Japan','Canada','Paracou','Mbaiki','US')
library(doParallel)
cl <- makeCluster(7)
registerDoParallel(cl)
foreach(set=sets) %dopar% run.models.for.set.all.traits(set,model.files.lmer.Tf.1, run.lmer,type.filling='species',std=FALSE)
foreach(set=sets) %dopar% run.models.for.set.all.traits(set,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=FALSE)
foreach(set=sets) %dopar% run.models.for.set.all.traits(set,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=TRUE)
source("R/analysis/lmer.run.R")
run.models.for.set.all.traits('France',model.files.lmer.Tf.1[3],run.lmer, traits='Wood.density',type.filling='species')
run.models.for.set.all.traits('Spain',model.files.lmer.Tf.1[3],run.lmer, traits='Wood.density',type.filling='species')
run.models.for.set.all.traits('Sweden',model.files.lmer.Tf.1[3],run.lmer, traits='Wood.density',type.filling='species')
### LOOK AT DIFERENCES OF THE TWO TYPE OF CI
### LOAD FOR WOOD DENSITY
fun.plot.compare.CI <- function(set,ecoregion,trait='Wood.density',type.filling='species'){
ER <- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.ER.Tf", "results.rds"))
CI.b<- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.ER.Tf", "results.ci.rds"))
fixed.coeff.E <- fixef(ER)
fixed.coeff.Std.Error <- sqrt(diag(vcov(ER)))
CI.sd <- cbind(fixed.coeff.E - 1.96*fixed.coeff.Std.Error,
fixed.coeff.E + 1.96*fixed.coeff.Std.Error)
mp <- barplot(fixed.coeff.E,ylim=range(CI.sd,CI.b),main=paste(set,ecoregion,trait))
segments(mp,CI.sd[,1],mp,CI.sd[,2])
segments(mp+0.1,CI.b[,1],mp+0.1,CI.b[,2],lty=2,col="red")
}
get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
}
pdf("figs/France.compare.CI.WD.pdf")
ecoregions <- get.ecoregions.for.set('France')
par(mfrow=c(2,3))
for (ecoregion in ecoregions)
try(fun.plot.compare.CI(set='France',ecoregion,trait='Wood.density',type.filling='species'))
dev.off()
pdf("figs/Spain.compare.CI.WD.pdf")
ecoregions <- get.ecoregions.for.set('Spain')
par(mfrow=c(2,4))
for (ecoregion in ecoregions)
try(fun.plot.compare.CI(set='Spain',ecoregion,trait='Wood.density',type.filling='species'))
dev.off()
## run.models.for.set.all.traits('France',model.files.lmer.Tf.2, run.lmer,type.filling='species')
#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R",
trait='Wood.density',set='France',ecoregion='HI',
type.filling='species')
run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",
trait='SLA',set='France',ecoregion='HI',
type.filling='species')
run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R",
trait='SLA',set='France',ecoregion='HI',
type.filling='species')
run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",
trait='SLA',set='France',ecoregion='HI',
type.filling='species')
run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.AD.Tf.R",
trait='SLA',set='France',ecoregion='HI',
type.filling='species')
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.boxplot.breaks <- function(x,y,Nclass=15,...){
breakss <- seq(min(x),max(x),length=Nclass+1)
x.cut <- cut(x,breakss)
mid.points <- breakss[-(Nclass+1)]+(breakss[2]-breakss[1])/2
mat <- cbind(x,y)
data <- as.data.frame(mat)
colors.dens <- densCols(mat)
plot(x,y, col=colors.dens, pch=20,...)
boxplot(y~x.cut,at=mid.points,add=TRUE,names=NA)
}
seq.from.to.quantile <- function(x,length.out=20,probs=c(0.001,0.999)){
qq <- quantile(x,probs)
return(seq(from=qq[1],to=qq[2],length.out=length.out))
}
fun.plot.residual.plus.regression.lines <- function(df.lmer,res.fix.no,res.fix.simple,
ERcomp,ADcomp,trait,set,ecoregion,
type.filling){
par(mfrow=c(2,3))
## Effect /reponse
fun.boxplot.breaks(df.lmer$sumBn,res.fix.no,cex=0.1,
xlab="sum of basal area",
ylab="growth residual",
main="Effect/response model")
mtext(paste(trait,set,ecoregion,type.filling), side=3, line =1,outer=TRUE)
lines(seq.from.to.quantile(df.lmer$sumBn),
seq.from.to.quantile(df.lmer$sumBn)*fixef(ERcomp)['sumBn'],
col='red')
fun.boxplot.breaks(df.lmer$sumTnBn,res.fix.simple,cex=0.1,
xlab="sum of basal area x Tn",
ylab="growth residual",
main="Effect/response model")
lines(seq.from.to.quantile(df.lmer$sumTnBn),
seq.from.to.quantile(df.lmer$sumTnBn)*fixef(ERcomp)['sumTnBn'],col='red')
fun.boxplot.breaks(df.lmer$sumTfBn,res.fix.simple,cex=0.1,
xlab="sum of basal area x Tf",
ylab="growth residual",
main="Effect/response model")
lines(seq.from.to.quantile(df.lmer$sumTfBn),
seq.from.to.quantile(df.lmer$sumTfBn)*fixef(ERcomp)['sumTfBn'],col='red')
## Absolute distance model
fun.boxplot.breaks(df.lmer$sumBn,res.fix.no,cex=0.1,
xlab="sum of basal area",
ylab="growth residual",
main="Absolute distance model")
lines(seq.from.to.quantile(df.lmer$sumBn),
seq.from.to.quantile(df.lmer$sumBn)*fixef(ADcomp)['sumBn'],
col='red')
fun.boxplot.breaks(df.lmer$sumTnTfBn.abs,res.fix.simple,cex=0.1,
xlab="sum of basal area x absolute trait distance",
ylab="growth residual",
main="Absolute distance model")
lines(seq.from.to.quantile(df.lmer$sumTnTfBn.abs),
seq.from.to.quantile(df.lmer$sumTnTfBn.abs)*fixef(ADcomp)['sumTnTfBn.abs'],col='red')
}
fun.load.data.and.plot.residual.plus.regression.lines <- function(trait,
set,
ecoregion,
type.filling)
df.lmer <- load.and.prepare.data.for.lmer(trait,set,ecoregion,
min.obs=10, sample.size=NA,
type.filling) # return a DF
simple <- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.simplecomp.Tf", "results.rds"))
nocomp <- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.nocomp.Tf", "results.rds"))
ERcomp <- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.ER.Tf", "results.rds"))
ADcomp <- readRDS(file.path("output/lmer", set,ecoregion,trait,type.filling,"lmer.LOGLIN.AD.Tf", "results.rds"))
res.fix.simple <- fun.compute.resid(simple)
res.fix.no <- fun.compute.resid(nocomp)
fun.plot.residual.plus.regression.lines(df.lmer,res.fix.no,res.fix.simple,
ERcomp,ADcomp,trait,set,ecoregion,
type.filling)
}
### GLMR
source("R/analysis/glmer.run.R")
saveRDS(lmer.output,file="output/lmer/res.test.all.rds")
#### TODO START TO WORK ON A WAY TO CHECK PREDICTION AGAINST OBSERVED TO SEE TRAIT EFFECT
run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
run.glmer("R/analysis/model.glmer/model.glmer.LOGLIN.ER.AD.Tf.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
R/analysis/stan-fun.R
View file @ 8baf852a
......@@ -3,120 +3,111 @@ library(rstan) #http://www.mc-stan.org/
stan.model <- "
data {
int<lower=0> N.indiv;
int<lower=0> N.species;
int<lower=0> N.plot;
int<lower=0> N.ecocode;
int[N.indiv] species.id;
int[N.indiv] plot.id;
int[N.indiv] ecocode.id;
real[N.indiv] logG;
real[N.indiv] logD;
real[N.indiv] sumBn;
real[N.indiv] Tf;
real[N.indiv] sumTfBn;
real[N.indiv] sumTnBn;
int<lower=0> N_indiv;
int<lower=0> N_species;
int<lower=0> N_plot;
int<lower=0> N_ecocode;
int species_id[N_indiv];
int plot_id[N_indiv];
int ecocode_id[N_indiv];
real logG[N_indiv];
real logD[N_indiv];
real sumBn[N_indiv];
real Tf[N_indiv];
real sumTfBn[N_indiv];
real sumTnBn[N_indiv];
}
parameters{
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean.logD;
real<lower=-10,upper=10> mean.Tf;
real<lower=-10,upper=10> mean.sumBn;
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.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;
real<lower=-1,upper=1> intercept;
real<lower=-10,upper=10> mean_logD;
real<lower=-10,upper=10> mean_Tf;
real<lower=-10,upper=10> mean_sumBn;
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_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;
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;
}
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]]
+ 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]
;
}
}
model {
# constants for prior
real sigma0;
################################################################################
######################## growth model with STAN #############################
###############################################
########### Hierarchical parameters ########
### species random param
for (n in 1:N.species)
{
param.logD.sp[n] ~ dnorm(mean.logD,sigma.logD.sp);
intercept.sp[n] ~ dnorm(intercept,sigma.inter.sp);
}
param_logD_sp ~ normal(mean_logD,sigma_logD_sp);
intercept_sp ~ normal(0,sigma_inter_sp);
### plot random param
for (n in 1:N.plot)
{
intercept.plot[n] ~ dnorm(0,sigma.inter.plot)
}
intercept_plot ~ normal(0,sigma_inter_plot);
### ecocode random param
for (n in 1:N.ecocode)
{
param.Tf.ecocode[n] ~ dnorm(mean.Tf,sigma.Tf.ecocode);
param.sumBn.ecocode[n] ~ dnorm(mean.sumBn,sigma.sumBn.ecocode);
param.sumTfBn.ecocode[n] ~ dnorm(mean.sumTfBn,sigma.sumTfBn.ecocode);
param.sunTnBn.ecocode[n] ~ dnorm(meanTnBn,sigma.sumTnBn.ecocode);
intercept.ecocode[n] ~ dnorm(0,sigma.inter.ecocode);
}
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);
############ Likelihood ###################
for (i in 1:N.indiv) {
real theo.g;
theo.g <- 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]
;
logG ~ dnorm(theo.g,sigma);
}
###############################################
########### Non-hierarchical parameters ########
# constants for prior
real sigma0;
sigma0 <- pow(1.0E-2,-2);
########### Non-Hierarchical parameters ########
# constant for prior
sigma0 <- 10;
# slope and intercept
intercept ~ dnorm(0,sigma0);
mean.logD ~ dnorm(0,sigma0);
mean.Tf ~ dnorm(0,sigma0);
mean.sumBn ~ dnorm(0,sigma0);
mean.sumTfBn ~ dnorm(0,sigma0);
mean.sumTnBn ~ dnorm(0,sigma0);
intercept ~ normal(0,sigma0);
mean_logD ~ normal(0,sigma0);
mean_Tf ~ normal(0,sigma0);
mean_sumBn ~ normal(0,sigma0);
mean_sumTfBn ~ normal(0,sigma0);
mean_sumTnBn ~ normal(0,sigma0);
# sigma
sigma.inter.sp~ uniform(0,6);
sigma.inter.ecocode~ 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_inter_sp~ uniform(0,6);
sigma_inter_ecocode~ 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~ uniform(0,6);
} # End of the stan model
}
###############################################
############ Likelihood ###################
logG ~ normal(theo_g,sigma);
}
"
......@@ -149,31 +140,43 @@ df.lmer <- load.and.prepare.data.for.lmer(trait,
min.obs, sample.size,
type.filling=type.filling) # return a DF
stan.dat <- list(N.indiv=nrow(df.lmer),
N.species=nlevels(df.lmer$species.id),
N.plot=nlevels(df.lmer$plot.id),
N.ecocode=nlevels(df.lmer$ecocode.id),
species.id=unclass(factor(df.lmer$species.id)),
plot.id=unclass(factor(df.lmer$plot.id)),
species.id=unclass(factor(df.lmer$ecocode.id)),
logG=df.lmer$logG,
logD=df.lmer$logD,
Tf=df.lmer$Tf,
sumBn=df.lmer$sumBn,
sumTfBn=df.lmer$sumTfBn,
sumTnBn=df.lmer$sumTnBn)
df.lmer.t <- df.lmer
stan.dat <- list(N_indiv=nrow(df.lmer.t),
N_species=nlevels(df.lmer.t$species.id),
N_plot=nlevels(df.lmer.t$plot.id),
N_ecocode=nlevels(df.lmer.t$ecocode.id),
species_id=unclass(factor(df.lmer.t$species.id)),
plot_id=unclass(factor(df.lmer.t$plot.id)),
ecocode_id=unclass(factor(df.lmer.t$ecocode.id)),
logG=df.lmer.t$logG,
logD=df.lmer.t$logD,
Tf=df.lmer.t$Tf,
sumBn=df.lmer.t$sumBn,
sumTfBn=df.lmer.t$sumTfBn,
sumTnBn=df.lmer.t$sumTnBn)
# param to save
param.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',
'param_Tf_ecocode',
'param_sumBn_ecocode','param_sumTfBn_ecocode',
'param_sumTnBn_ecocode','intercept_ecocode')
#'intercept_sp','param_logD_sp',
#warmup = burn-in
stan.out <- stan(
model_code=stan.code,
data=stan.dat,
init=list(stan.inits.LOGLIN.NULL),
iter=20000,
warmup=10000,
chains=1)
#stan.out2 <- extract(stan.out,pars=c(names(sr.true.pars),"lp__"),inc_warmup=F)
#ncols <- dim(stan.out2)[3]
#stan.out3 <- matrix(stan.out2[,,-ncols],ncol=(ncols-1),dimnames=list(NULL,dimnames(stan.out2)$parameters[-ncols]))
system.time(stan.out <- stan(
model_code=stan.model,
data=stan.dat,
pars=param.to.save,
init='random',
iter=500,
warmup=100,
chains=3))
saveRDS(stan.out,file="output/stan.test.rds")
launch_all_lmer.nolog.bash 0 → 100644
View file @ 8baf852a
#!/bin/bash
export LD_LIBRARY_PATH=/usr/lib64/R/library
mkdir -p trait.workshop
for site in "'BCI'" "'Canada'" "'France'" "'Fushan'" "'NVS'" "'Paracou'" "'Spain'" "'US'" "'Swiss'" "'Sweden'" "'NSW'" "'Mbaiki'" "'Luquillo'" "'Japan'"; do
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species',std=TRUE);print('done')\"" > trait.workshop/species1$site.sh
qsub trait.workshop/species1$site.sh -l nodes=1:ppn=1 -N "lmerspecies1$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=TRUE);print('done')\"" > trait.workshop/species2$site.sh
qsub trait.workshop/species2$site.sh -l nodes=1:ppn=1 -N "lmerspecies2$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus',std=TRUE);print('done')\"" > trait.workshop/genus1$site.sh
qsub trait.workshop/genus1$site.sh -l nodes=1:ppn=1 -N "lmergenus1$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus',std=TRUE);print('done')\"" > trait.workshop/genus2$site.sh
qsub trait.workshop/genus2$site.sh -l nodes=1:ppn=1 -N "lmergenus2$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species',std=FALSE);print('done')\"" > trait.workshop/species1std$site.sh
qsub trait.workshop/species1std$site.sh -l nodes=1:ppn=1 -N "lmerspecies1std$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species',std=FALSE);print('done')\"" > trait.workshop/species2std$site.sh
qsub trait.workshop/species2std$site.sh -l nodes=1:ppn=1 -N "lmerspecies2std$site" -q opt32G -j oe
echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.nolog.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus',std=FALSE);print('done')\"" > trait.workshop/genus1std$site.sh
qsub trait.workshop/genus1std$site.sh -l nodes=1:ppn=1 -N "lmergenus1std$site" -q opt32G -j oe