process glmer and lmer adn new residual plots

R/analysis/glmer.output-fun.R

+#### function to analyse glmer output
+
+
+library(lme4)
+
+
+read.glmer.output <- function(file.name){
+  tryCatch(readRDS(file.name), error=function(cond)return(NULL))   # Choose a return value in case of error
+}
+
+
+
+summarise.glmer.output <- function(x){
+ list( nobs = nobs(x),
+       R2m =Rsquared.glmm.lmer(x)$Marginal,
+       R2c =Rsquared.glmm.lmer(x)$Conditional,
+       AIC = AIC(x),
+       deviance = deviance(x),
+       effect.response.var=variance.fixed.glmm.lmer.effect.and.response(x),
+       fixed.coeff.E=fixef(x),
+       fixed.coeff.Std.Error=sqrt(diag(vcov(x))),
+       fixed.var=variance.fixed.glmm.lmer(x))
+}
+
+
+summarise.glmer.output.list <- function(f ){
+    output.glmer <- read.glmer.output(f)
+    if(!is.null(output.glmer)){
+	res <- list(files.details=files.details(f),
+		   glmer.summary=summarise.glmer.output( output.glmer))
+    }else{
+        res <- NULL
+    }
+    return(res)
+}
+
+
+
+files.details <- function(x){
+	s <- data.frame(t(strsplit(x, "/", fixed=TRUE)[[1]]), stringsAsFactors= FALSE)
+	names(s)  <- c("d1", "d2", "set", "ecocode", "trait", "filling", "model", "file" )
+	s[-(1:2)]
+}
+
+
+
+#### R squred functions
+
+# Function rsquared.glmm requires models to be input as a list (can include fixed-
+# effects only models,but not a good idea to mix models of class "mer" with models
+# of class "lme") FROM http://jslefche.wordpress.com/2013/03/13/r2-for-linear-mixed-effects-models/
+
+Rsquared.glmm.lmer <- function(i){
+{
+      # Get variance of fixed effects by multiplying coefficients by design matrix
+      VarF=var(as.vector(fixef(i) %*% t(i@pp$X)))
+      # Get variance of random effects by extracting variance components
+      VarRand=colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
+      # Calculate marginal R-squared
+      Rm=VarF/(VarF+VarRand+pi^2/3)
+      # Calculate conditional R-squared (fixed effects+random effects/total variance)
+      Rc=(VarF+VarRand)/(VarF+VarRand+pi^2/3)
+      # Bind R^2s into a matrix and return with AIC values
+      Rsquared.mat=data.frame(Class=class(i),Family=summary(i)$family,Marginal=Rm,
+                              Conditional=Rc,AIC=AIC(i)) }
+
+return(Rsquared.mat)
+}
+
+
+variance.fixed.glmm.lmer <- function(i){
+ # Get variance of for each fixed effects by multiplying coefficients by design matrix
+var.vec <- apply(fixef(i) * t(i@pp$X),MARGIN=1,var)
+ # Get variance of fixed effects by multiplying coefficients by design matrix
+VarF=var(as.vector(fixef(i) %*% t(i@pp$X)))
+ # Get variance of random effects by extracting variance components
+VarRand=colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
+var.vec <- var.vec/(VarF+VarRand+pi^2/3)
+names(var.vec) <- paste(names(var.vec),"VAR",sep=".")
+return(var.vec)
+}
+
+variance.fixed.glmm.lmer.effect.and.response <- function(i){
+if(sum(c("sumTfBn","sumTnBn") %in% names(fixef(i)))==2){
+# Get variance of for each fixed effects by multiplying coefficients by design matrix
+var.vec <- var(as.vector(fixef(i)[c("sumTfBn","sumTnBn")] %*% t(i@pp$X[,c("sumTfBn","sumTnBn")])))
+# Get variance of fixed effects by multiplying coefficients by design matrix
+VarF <- var(as.vector(fixef(i) %*% t(i@pp$X)))
+# Get variance of random effects by extracting variance components
+VarRand <- colSums(do.call(rbind,lapply(VarCorr(i),function(j) j[1])))
+# Get residual variance
+var.vec <- var.vec/(VarF+VarRand)
+}else{
+var.vec <- NA
+}    
+names(var.vec) <- paste("effect.response","VAR",sep=".")
+return(var.vec)
+}
+
+
+## function to turn lmer output from list to DF
+fun.format.in.data.frame <- function(list.res,names.param){
+dat.t <- data.frame(t(rep(NA,3*length(names.param))))
+names(dat.t) <-  c(names.param,paste(names.param,"Std.Error",sep=".")
+                   ,paste(names.param,"VAR",sep="."))
+dat.t[,match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
+    list.res$glmer.summary$fixed.coeff.E
+dat.t[,match(names(list.res$glmer.summary$fixed.coeff.E),names(dat.t))] <-
+    list.res$glmer.summary$fixed.coeff.Std.Error
+dat.t[,match(names(list.res$glmer.summary$fixed.var),names(dat.t))] <-
+    list.res$glmer.summary$fixed.var
+res <- data.frame(list.res$files.details,list.res$glmer.summary[1:6],dat.t)
+return(res)
+}
+
+
+
+
+

R/analysis/glmer.output.R

+#!/usr/bin/env Rscript
+
+source("R/analysis/glmer.output-fun.R")
+files <- list.files("output/glmer", recursive=TRUE,  full.names =TRUE, pattern = "rds")
+out <- lapply(files, summarise.glmer.output.list)
+names(out) <- lapply(lapply(files,files.details),function(x) paste(as.vector(x[-6]),collapse="_"))
+saveRDS(out,file='output/glmer.list.out.rds')
+names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$glmer.summary$fixed.coeff.E))))
+DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
+DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
+saveRDS(DF.results,file='output/glmer.DF.rds')

R/analysis/glmer.run.R

@@ -38,10 +38,11 @@ model.files.glmer.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.R",
                  "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.R")
 
 model.files.glmer.Tf.1 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.E.Tf.R",
-                 "R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R.",
+                 "R/analysis/model.glmer/model.glmer.LOGLIN.R.Tf.R",
                  "R/analysis/model.glmer/model.glmer.LOGLIN.ER.Tf.R")
 model.files.glmer.Tf.2 <- c("R/analysis/model.glmer/model.glmer.LOGLIN.nocomp.Tf.R",
                  "R/analysis/model.glmer/model.glmer.LOGLIN.AD.Tf.R",
+                 "R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.R",
                  "R/analysis/model.glmer/model.glmer.LOGLIN.simplecomp.Tf.R")
 
 
@@ -105,7 +106,7 @@ data.tree <- subset(data.tree,subset=(!is.na(data.tree[["dead"]])) &
                                      (!is.na(data.tree[["D"]])) )
 ## remove tree with zero
 data.tree <- subset(data.tree,subset=data.tree[["D"]]>9.9)             
-## select species with no missing traits 
+## select species with missing traits 
 data.tree <- data.tree[(!is.na(data.tree[[abs.CWM.tntf]]) &
                     !is.na(data.tree[[BATOT]])),]
 # select  species with minimum obs    

R/analysis/lmer.output.R

@@ -8,5 +8,5 @@ saveRDS(out,file='output/lmer.list.out.rds')
 names.param <- unique(unlist(lapply(out,function(list.res) names(list.res$lmer.summary$fixed.coeff.E))))
 DF.results <- do.call("rbind",lapply(out,fun.format.in.data.frame,names.param=names.param))
 DF.results$id <- paste(DF.results$set,DF.results$ecocode,sep=".")
-saveRDS(out,file='output/lmer.DF.rds')
+saveRDS(DF.results,file='output/lmer.DF.rds')
 

R/analysis/lmer.run.R

@@ -38,10 +38,11 @@ model.files.lmer.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",
                  "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R")
 
 model.files.lmer.Tf.1 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.E.Tf.R",
-                 "R/analysis/model.lmer/model.lmer.LOGLIN.R.Tf.R.",
+                 "R/analysis/model.lmer/model.lmer.LOGLIN.R.Tf.R",
                  "R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R")
 model.files.lmer.Tf.2 <- c("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.Tf.R",
                  "R/analysis/model.lmer/model.lmer.LOGLIN.AD.Tf.R",
+                 "R/analysis/model.lmer/model.lmer.LOGLIN.HD.Tf.R",
                  "R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.Tf.R")
 
 

R/analysis/model.glmer/model.glmer.LOGLIN.HD.Tf.R

+load.model <- function () {
+    list(name="glmer.LOGLIN.HD.Tf",
+         glmer.formula=formula("dead~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnTfBn.diff"))
+}         

R/analysis/model.lmer/model.lmer.LOGLIN.HD.Tf.R

+load.model <- function () {
+    list(name="lmer.LOGLIN.HD.Tf",
+         lmer.formula=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|plot.species.id)+sumBn+sumTnTfBn.diff"),
+         lmer.formula.tree.id=formula("logG~1+Tf+logD+(1+logD|species.id)+(1|tree.id)+(1|plot.species.id)+sumBn+sumTnTfBn.diff"))
+}         
+         
+         
+        

R/analysis/plots.resid.regression.lines.R

+##### SCRIPT TO TEST plots resuls with regression lines
+source("R/analysis/lmer.run.R")
+
+
+
+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),oma=c(0,0,2,0))
+## 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=0.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)
+dir.create("figs/plot.resid", recursive=TRUE, showWarnings=FALSE)
+png(paste("figs/plot.resid/",
+          paste(trait,set,ecoregion,type.filling,"residual.png",sep="."),sep=""),
+    width = 1500, height = 1000)
+fun.plot.residual.plus.regression.lines(df.lmer,res.fix.no,res.fix.simple,
+                                         ERcomp,ADcomp,trait,set,ecoregion,
+                                         type.filling)
+dev.off()
+}
+
+
+## function to get all set and ecoregion to plot
+
+get.ecoregions.for.set <- function(set, base.dir = "./output/processed/"){
+  sub(paste(base.dir,set,"/",sep=""),"",list.dirs(paste(base.dir,set,sep="")))[-1]
+}
+
+plot.residual.for.set.all.traits  <- function(set,
+               traits = c("SLA", "Wood.density","Max.height","Leaf.N","Seed.mass"),
+               type.fillings=c("species","genus") , ...){
+ecoregions <- get.ecoregions.for.set(set, base.dir = "./output/processed/")    
+  for(trait in traits){
+      for (ecoregion in ecoregions){
+          for (fill in type.fillings){
+          try(fun.load.data.and.plot.residual.plus.regression.lines(trait,set,
+                                                          ecoregion,
+                                                          type.filling=fill))
+         }
+     }
+ }     
+}
+
+
+for (set in c("BCI","Canada","France","Fushan","NVS","Paracou","Spain","US","Swiss","Sweden","NSW","Mbaiki","Luquillo","Japan")){
+    ## 
+plot.residual.for.set.all.traits(set)
+
+}    
+
+
+### TODO NEED TO PLOT ALL REGRESSION LINES OF ALL SET ON SAME FIGURE FOR EACH TRAIT AND SEE WHY NOT GOING IN SAME DIRECTION

R/analysis/run.local.R

@@ -2,55 +2,109 @@
 source("R/analysis/lmer.run.R")
 
 #### 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.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
-## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='genus')
-## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
-## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
-## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
-## run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.HD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.R",trait='SLA',set='France',ecoregion='HI',type.filling='genus')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.AD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.simplecomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.nocomp.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
+run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.HD.R",trait='SLA',set='France',ecoregion='HI',type.filling='species')
 run.lmer("R/analysis/model.lmer/model.lmer.LOGLIN.ER.Tf.R",trait='SLA',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')
 
-df.lmer <- load.and.prepare.data.for.lmer(trait='SLA',set='France',ecoregion='HI',
-                                           min.obs=10, sample.size=NA,
-                                         type.filling='species') # return a DF
 
-         
-## simple <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp/results.rds")
-## nocomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp/results.rds")
-## ERcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER/results.rds")
-## ADcomp <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD/results.rds")
-simple.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.simplecomp.Tf/results.rds")
-nocomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.nocomp.Tf/results.rds")
-ERcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.Tf/results.rds")
-ADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.AD.Tf/results.rds")
-ERADcomp.Tf <- readRDS("output/lmer/France/HI/SLA/fill/lmer.LOGLIN.ER.AD.Tf/results.rds")
+output.dir.lmer <- function (model, trait, set, ecoregion,type.filling) {
+  file.path("output/lmer", set,ecoregion,trait,type.filling,model)
+}
+
+
 
 
-anova(simple,nocomp,ERcomp,ADcomp,simple.Tf,nocomp.Tf,ERcomp.Tf,ADcomp.Tf,ERADcomp.Tf)
+fun.compute.resid <- function(i){
+return(fitted(i) +resid(i) -i@pp$X %*%fixef(i))
+}    
 
-res.fix.simple <- fitted(simple)+resid(simple) -simple@pp$X %*% fixef(simple)
-res.fix.no <- fitted(nocomp)+resid(nocomp) -nocomp@pp$X %*% fixef(nocomp)
-res.fix.ER <- fitted(ERcomp)+resid(ERcomp) -ERcomp@pp$X %*% fixef(ERcomp)
+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)
+}
 
-par(mfrow=c(2,2))
- plot(df.lmer$sumBn,res.fix.no,cex=0.1)
-lines(0:15,0:15*fixef(simple)['sumBn'],col='red')
- plot(df.lmer$sumTnBn,res.fix.simple,cex=0.1)
-lines((-10):10,(-10):10*fixef(ERcomp)['sumTnBn'],col='red')
- plot(df.lmer$sumTfBn,res.fix.simple,cex=0.1)
-lines((-10):10,(-10):10*fixef(ERcomp)['sumTfBn'],col='red')
- plot(df.lmer$sumTnTfBn.abs,res.fix.simple,cex=0.1)
-lines((-10):10,(-10):10*fixef(ADcomp)['sumTnTfBn.abs'],col='red')
+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"))
 
- plot(df.lmer$Tf,res.fix.ER,cex=0.1)
-lines((-10):10,(-10):10*fixef(ERcomp.Tf)['Tf'],col='red')
+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

docs/analysis.outline/update.progress/progress.md

@@ -22,12 +22,12 @@ more likely to coexist locally than similar species
 @macarthur_limiting_1967). One way to quantify ecological similarity
 between species is via traits, such as leaf, seed and wood
 characteristics [@westoby_plant_2002]. Traits influence many aspects
-of plant performance, including resource acquisition. Under the *competition-niche similarity hypothesis* higher trait dissimilarity should results in higher resource partitioning at
+of plant performance, including resource acquisition. Under the *competition-niche similarity hypothesis* higher trait dissimilarity should results in higher resources partitioning at
 local scale and less intense competition. This idea underlies numerous ecological analyses
 [@kraft_functional_2008; @cornwell_community_2009]. However this
 assumption has rarely been tested against field or experimental
-outcomes. This is surprising because it is well known that competitive
-interactions among vascular plants are more complex. For instance,
+outcomes (but see @uriarte_trait_2010). This is surprising because it is well known that competitive
+interactions among vascular plants are more complex than simple resources partitioning. For instance,
 most plant species compete for the same limiting resources (water,
 light and nutrients), which makes simple local resources partitioning
 unlikely. The ranking of competitive ability for these common
@@ -36,11 +36,10 @@ interaction. If ranking processes are dominant, competitive outcomes
 should be more closely related to the hierarchy (or the hierarchical
 distance) of relevant functional traits
 [@mayfield_opposing_2010; @kunstler_competitive_2012] than trait
-dissimilarity. Recent analysis of competitive interactions at local
-scale between individual trees (using growth analysis with local
-competition index) in mountain forests in the French Alps
-[@kunstler_competitive_2012], support this view that competition is
-more related to trait hierarchy than trait similarity.
+dissimilarity. A recent analysis [@uriarte_trait_2010] of competitive interactions at local
+scale between individual trees (using growth adn survival with local
+competition index) have shown that it was possible to test the link between competition and traits with simple competition index. Following this approach Kunstler et al. [@kunstler_competitive_2012] have applied this method in mountain forests in the French Alps and shown that competition is
+more related to trait hierarchy than trait similarity for tree growth in this system.
 
 # Objective
 Given the importance in the ecological literature of the idea that
@@ -68,14 +67,15 @@ G_{f,p,i,t} = G\textrm{max}_{f,p,i} \times s(D_{i,t}) \times g\left(\sum_{n=1}^{
 where:
 
 - $G_{f,p,i,t}$ is the growth (diameter or basal area growth) of
-an individual $i$ from species $f$ growing in plot plot $p$ in census $t$,
+an individual $i$ from species $f$ growing in plot $p$ in census $t$,
 - $D_i$ is the diameter of the individual $i$,
-- $B_{n}$ is the the basal area of neighborhood tree of species $n$,
+- $B_{n}$ is the the basal area of neighborhood trees of species $n$,
 - $G\textrm{max}_{f,p,i}$ is the maximum growth rate of the focal species $f$ on the plot $p$ for the individual $i$,
 - $s$ and $g$ are functions representing the size and the competition effect respectively, and
 - $\lambda_{n,f}$ is a parameter representing the growth reduction for a
 unit of neighborhood basal area increase of species $n$ on species
-$f$.
+$f$,
+- $N_p$ is the number of species on plot $p$.
 
 
 
@@ -97,7 +97,7 @@ The central part of the analysis involves comparing alternative models for $\lam
 \lambda_{n,f} = a +b \times (t_{n} - t_{f}).
 \end{equation}
 
-The logic behind the hierarchical trait distance model, can be understand through a decomposition of competition in competitive effect and competitive response. The hierarchical trait distance model occurs when the traits conferring a high competitive effect also confer a high competition tolerance[^compreponse]. During the first day of the workshop we discussed the possibility of including a model with separate
+The logic behind the hierarchical trait distance model, can be understood through a decomposition of competition in competitive effect and competitive response. The hierarchical trait distance model occurs when the traits conferring a high competitive effect also confer a high competition tolerance[^compreponse]. During the first day of the workshop we discussed the possibility of including a model with separate
 links of traits with competitive effect and competitive response.  This
 model is connected to several papers by Goldberg *et al.*, where
  competition is framed in term of effect and response and their links to
@@ -135,16 +135,16 @@ functions for competitive response and effect respectively. A series of models w
 
 3. $\lambda$ is influenced by traits of the neighborhood and
    focal species (effect-response model):
-\begin{equation}  \label{response_effect_trait}
+\begin{equation}  \label{RE}
 \lambda_{n,f} = a +b \times t_{f} +c \times t_{n}.
 \end{equation}
 
 The trait hierarchical distance model eq. \ref{hier_dist_trait} is a
 sub-case of the effect and response model
-eq. \ref{response_effect_trait} where $b=-c$.
+eq. \ref{RE} \ref{hier_dist_trait} where $b=-c$.
 
 During the workshop David Coomes described how to express the model as a function of the community weighted mean of the trait of the neighborhood trees. For the
-most complex model eq. \ref{response_effect_trait} this gives:
+most complex model eq. \ref{RE} TODO SOLVE ISSUE WITH EQUATION REF this gives:
 
 \begin{equation}
 \sum_{n=1}^{N_p} \lambda_{n,f} \times B_n = \sum_{n=1}^{N_p} (a +b \times t_{f}
@@ -160,7 +160,7 @@ relative basal area abundance of species $n$, $B_n/B_\textrm{tot}$).
 
 Subsequent to the workshop, and in the material I presented at Ecotas13[^ecotas], I decided to
 compare the absolute trait distance model eq. \ref{abs_dist_trait} and the
-effect-response model eq. \ref{response_effect_trait}.
+effect-response model eq. \ref{RE}.
 
 [^ecotas]: Joint conference of the Ecological Society of Australia and the New Zealand Ecological Society, Nov 2013.
 
@@ -237,7 +237,7 @@ table \ref{table-perc} gives the percentage of the data for which at
 least 90% of neighborhood is covered with species or genus level
 trait. Paracou and M'Baiki are the only two sites with very low
 coverage (this because of missing traits but also because of missing
-taxonomic identification).
+taxonomic identification). TODO NEED TO CHECK WITH GHISLAIN NEW SPECIES LIST.
 
 ## Fitting of a mixed linear model
 
@@ -336,11 +336,55 @@ Most of the effect-response models fitted show a competitive effect (negative va
 - Explore if the decrease in the link between trait and competition at high MAP is related in a change in the packing of trait space in this communities.
 - Explore the possibility that trait effect may be different for evergreen/deciduous species (leaf traits) or angiosperm/conifer species (wood density). This could be done by fitting different parameters for the trait of evergreen deciduous and conifer in the effect-response model. This is not really possible for the absolute distance model.
 - Use an alternative way of dividing the NFI data than the ecoregion (class of MAP and MAT?).
-- Try to run a global analysis with all data (memory limit issue to solve).
+- Try to run a global analysis with all data (memory limit issue to solve). and run with global standarization.
 <!-- - build plots based on residual to look at the traits effects on -->
 <!--   competition. -->
 
 
+\newpage
+
+# Comments received
+
+@ **Nathan Swenson**: . It appears you are assessing growth/performance using diameter changes. This approach is very flawed from my perspective as it isn't a great metric of actual biomass/growth. Of course it is approximate, but the amount of biomass accumulated is also very dependent on the wood density. Simply put, the diameter increment increase of 1cm in a pioneer tree in our Puerto Rico dataset is not the same as a 1cm increase in a shade tolerant tree. The shade tolerant tree has likely put on much more biomass. One would also want to know height/dbh allometries and not knowing that will also propagate error. However, those are often not available data, but wood density is and I'm worried that you/we may be mischaracterizing performance. *\color{red} I'm using basal area growth and not diameter growth. Need to clarify that. Include use wood density, basal area growth and simple allomtric function is possible to estimate biomass growth. But using the same allometric function may lead to more error than the inclusion of wood density?*
+
+@ **Sylvie Gourlet-Fleury**:
+- Figure 5 6 and 7 unclear need to add a theorethical figure to explain that.
+
+- How are a, b and c related? Does "a" take most of the competition index effect? Is there any tendency revealed by the links between the three coeffs? And by the way, I wonder wether it could not be possible to constrain competition to be splitted between the focal tree and the neighbours (by species) only (a=0, only b and c non nul): would it make sense? This would avoid Btot absorb the greatest part of the effect? *\color{red} The intercept a account for the fact that a part of competition can be unrelated to the traits variation. Assuming $a=0$ would mean that competition effect for instance vary directly with trait value not sure this is realistic.*
+
+- 2) my skepticism regarding the possibility to estimate the number of parameters of your model (9) properly, along with the possibility to interpret them properly. It seems to me that you should keep the competiiton term as simple as possible (and for example $a=0$). *\color{red} There is no convergence problem in the estimation.*
+
+@ ** Mark Westoby**: 
+It's interesting that standardizing traits per data set doesn't produce very different results from standardizing globally. But I suppose it's true that for the four traits used here there may not be very strong climate-zone shifts in the mean, so long as one stays within forests.
+
+On the results: I have to say that I don't understand how the absolute distance models can be selected as the "best model" so frequently. Basically AD shows near-zero effect size (except for max height) as shown in Fig 4. Given that the effect size for effect-response is pretty small (R2 0.01 or so) and has more free parameters I can understand how AD can be selected as better than ER. But surely the no-effect models should then also be superior to the AD models. *\color{red} I think that I should maybe add a model with inter vs. intra competition to see if better. Need to test a HD model as well! sent on cluster*
+
+Anyhow, at this stage it is the effect-size results that make sense to me more so than the selected-best-model results.
+
+Further, the most striking result overall is how very small are the trait influences as they are currently being measured. So I totally agree that a top priority should be to express trait effect size relative to what is explained by total basal area of neighbours unweighted by traits. So far as I understand, the reference quantity here would be the difference in R2 between a model with just the Gmax and size effects, and a model that included all the lambda effects but without trait influence on either focal or neighbour side. And actually to generalize that point a bit, it seems to me that having some picture how much of the total growth variation arises from Gmax versus from individual size versus from competitive influence, might be helpful background to the narrative we want to build?  *\color{red} The best solution would be to compare the R^2 of a model with no trqait effect on competition vs a model with one $\lambda$ for each pair of species and then see what propotion of this variation is explained by traits with the R^2 of the trait model. But the full $\lambda$ model is probably difficult to estimate. An alternative is to fit a model with randome effect for $\lambda$. A last alternative is to compare the improvement in R^2 du to trait to the improvement of R^2 due to traits effect on competition. ONE SOLUTION THAT I STARTED TO IMPLEMENT IS TO LOOK AT RATIO OF VARIAnCE EXPLAINED BY TRAIT COMP EFFET OVER COMP ONLY EFFECT. THIS WITHIN THE SAME MODEL AND NOT COMPARING SEVERAL MODEL!!*
+
+TRY TO COMPUTE ${R^2_{trait} - R^2_{compet}} \over {R^2_{compet} - R^2_{nocompet}}$. OR do a variance decomposition. Fit a model with random lambda?
+
+For me, non-linear and multi-trait models would not be a high priority. By this I mean that I can envisage providing them in supplementary materials somewhere just to show that they were done, but I find it hard to imagine building the main narrative of a paper around them. *\color{red} OK.*
+
+Maybe there is a real biological question in the multi-trait matter somewhat along the following lines: suppose that max height and wood density both have some explanatory power in a particular instance -- is that actually the same component of variation being explained by the two traits, or are the effects additive? *\color{red} and if not does including two traits explain more?*
+
+I do somewhat understand how, if traits were allowed to be predictors of Gmax, this might substantially affect the apparent effects of traits on lambda, so that does seem to me worth investigating. The strong interaction between fits of different parameters is illustrated, is it not, in Fig 7? -- if I understand that correctly, that is the same response to total basal area parameter ("a") fitted in the four panels, the difference being only which traits are allowed to have effects on the other parameters b and c? *\color{red} tested no major effect.*
+
+In Figs 6 and 7, as I have expressed to you before, I have a generalized disquiet about the parameter-fitting process. There is considerable scatter in outcome between the different datasets. That in itself doesn't bother me -- it makes sense that there should be very considerable noise in these datasets, and one can think of different datasets as subsamples from some larger universe, and look at the broad features of the clouds of points in Figs 6 and 7. But at the same time, the model is attributing very tight confidence intervals to the slope parameters for each dataset and trait. If taken literally, that means that almost each pair of datasets has meaningfully different slope parameter for the trait response. But it's hard for me to imagine any credible biological interpretation of that. I feel that somehow these models are attributing confidence intervals that are much too tight, and this is a sentiment I have felt when looking at other Bayesian-fit papers also. But my understanding of Bayesian methods is not good enough to be able to guess in what way they are overestimating the precision with which we know things. *\color{red} need to discuss that! Could be good to discuss that with a statistition. Why are the sd error *
+
+Anyhow in summary, I would suggest as top priority a partitioning or quantification of the different effect sizes relative to each other (and then the sets of alternative models could be laid out in a way that corresponded to that).
+
+*\color{red} What about evergreen/decidious?*
+
+@ ** Francis Hui**; 
+
+- Try AIC and BIC, and see if the give different results. More generally, information criteria are a puzzler for mixed models. You need to think about: 1) whether to use the log-likelihood conditional or unconditional on the random effects, and 2) how you count the number of parameters. Does each variance estimate count as a parameter? There's no correct answer to both problems. Have a look at http://biomet.oxfordjournals.org/content/97/4/773.short and see if/how it can be applied to your work. 
+
+- If you want to implementing a model with multiple traits, it is also possible to do these in a linear mixed model framework rather than having to go Bayesian. Check out the glmmlasso package in R. It implements GLMMs with a lasso penalty. The lasso is a well know penalty based on the sum of the absolute value of the coefficients (L1 norm). One can think of as the frequentist version of putting a Laplace prior on the coefficients. Maximization of a penalized log-likelihood can shrink coefficients to exactly 0, thereby doing simultaneous estimation and variable selection.
+
+* NEED TO LOOK AT THAT*
+
 \newpage
 
 
@@ -701,7 +745,7 @@ Developing the multiplicative model gives
 (a +b \times t_{f}) \times (c +d \times t_{n}) = ac+bc \times t_f +ad \times t_n +bd \times t_f \times t_n
 \end{equation}
 
-This equation bears some similarity to the additive model plus interaction Equation \label{add-inter} - which is an extension of the effect/response model presented above (equation \label{response_effect_trait}) - which include an interaction between the traits $t_n$ and $t_f$ is:
+This equation bears some similarity to the additive model plus interaction Equation \label{add-inter} - which is an extension of the effect/response model presented above (equation \label{RE}) - which include an interaction between the traits $t_n$ and $t_f$ is:
 
 \begin{equation} \label{add-inter}
 \lambda_{n,f} = a' +b' \times t_{f} +c' \times t_{n}+d' \times t_{n} \times t_{f}
@@ -713,7 +757,7 @@ The two models are equal when:
 a'=ac \mspace{3mu} ;\mspace{3mu} b'=bc\mspace{3mu} ;\mspace{3mu} c'=ad \mspace{5mu} and \mspace{5mu} d'=bd
 \end{equation}
 
-The multiplicative model is more constraining than the additive model plus interaction. In other word the additive model with interaction can be fitted to any multiplicative model but the inverse is not true (This would requires adding an interaction in the multiplicative model). For instance, it is not possible to match the hierarchical distance because if $b'$ and $d' \neq 0$ then $d' \neq 0$ as well. More generally, if parameters $a$, $b$ , $c$ and $d$ vary between [-max.r, max.r] then $d'>b'*c'/(max.r^2)$ (or  $d'<b'*c'/(-max.r^2)$). Thus it is not possible to have a strong traits effect on response and effect and no interaction. From first principle I think it is difficult to decide which model (equation \label{response_effect_trait} or equation \label{multi-er}) is the most likely.
+The multiplicative model is more constraining than the additive model plus interaction. In other word the additive model with interaction can be fitted to any multiplicative model but the inverse is not true (This would requires adding an interaction in the multiplicative model). For instance, it is not possible to match the hierarchical distance because if $b'$ and $d' \neq 0$ then $d' \neq 0$ as well. More generally, if parameters $a$, $b$ , $c$ and $d$ vary between [-max.r, max.r] then $d'>b'*c'/(max.r^2)$ (or  $d'<b'*c'/(-max.r^2)$). Thus it is not possible to have a strong traits effect on response and effect and no interaction. From first principle I think it is difficult to decide which model (equation \label{RE} or equation \label{multi-er}) is the most likely.
 
 
 \pagebreak

launch_all_glmer.bash

@@ -6,16 +6,16 @@ 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/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies1$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='species');print('done')\"" > trait.workshop/glmerspecies1$site.sh
 	qsub trait.workshop/glmerspecies1$site.sh -l nodes=1:ppn=1 -N "glmerspecies1$site" -q opt32G -j oe
 
-  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='species');\n\"" > trait.workshop/glmerspecies2$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='species');print('done')\"" > trait.workshop/glmerspecies2$site.sh
 	qsub trait.workshop/glmerspecies2$site.sh -l nodes=1:ppn=1 -N "glmerspecies2$site" -q opt32G -j oe
 
-  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus1$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1, run.glmer,type.filling='genus');print('done')\"" > trait.workshop/glmergenus1$site.sh
 	qsub trait.workshop/glmergenus1$site.sh -l nodes=1:ppn=1 -N "glmergenus1$site" -q opt32G -j oe
 
-  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='genus');\n\"" > trait.workshop/glmergenus2$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2, run.glmer,type.filling='genus');print('done')\"" > trait.workshop/glmergenus2$site.sh
 	qsub trait.workshop/glmergenus2$site.sh -l nodes=1:ppn=1 -N "glmergenus2$site" -q opt32G -j oe
 
 done

launch_all_glmer_HD.bash

+#!/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/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2[3], run.glmer,type.filling='species');print('done')\"" > trait.workshop/glmerspecies1$site.sh
+	qsub trait.workshop/glmerspecies1$site.sh -l nodes=1:ppn=1 -N "glmerspecies1$site" -q opt32G -j oe
+
+
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.2[3], run.glmer,type.filling='genus');print('done')\"" > trait.workshop/glmergenus1$site.sh
+	qsub trait.workshop/glmergenus1$site.sh -l nodes=1:ppn=1 -N "glmergenus1$site" -q opt32G -j oe
+
+
+done
+

launch_all_glmer_R.bash

+#!/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/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1[2], run.glmer,type.filling='species');print('done')\"" > trait.workshop/glmerspecies1$site.sh
+	qsub trait.workshop/glmerspecies1$site.sh -l nodes=1:ppn=1 -N "glmerspecies1$site" -q opt32G -j oe
+
+
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/glmer.run.R');run.models.for.set.all.traits($site,model.files.glmer.Tf.1[2], run.glmer,type.filling='genus');print('done')\"" > trait.workshop/glmergenus1$site.sh
+	qsub trait.workshop/glmergenus1$site.sh -l nodes=1:ppn=1 -N "glmergenus1$site" -q opt32G -j oe
+
+
+done
+

launch_all_lmer.bash

@@ -6,16 +6,16 @@ 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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species');\n\"" > trait.workshop/species1$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='species');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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species');\n\"" > trait.workshop/species2$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='species');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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus1$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1, run.lmer,type.filling='genus');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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus');\n\"" > trait.workshop/genus2$site.sh
+  	echo "/usr/local/R/R-3.0.1/bin/Rscript -e \"source('R/analysis/lmer.run.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2, run.lmer,type.filling='genus');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
 
 done

launch_all_lmer_HD.bash

+#!/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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2[3], run.lmer,type.filling='species');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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.2[3], run.lmer,type.filling='genus');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
+
+done

launch_all_lmer_R.bash

+#!/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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1[2], run.lmer,type.filling='species');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.R');run.models.for.set.all.traits($site,model.files.lmer.Tf.1[2], run.lmer,type.filling='genus');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
+
+done