#### function to analyse lmer output


library(lme4)


read.lmer.output <- function(file.name){
  tryCatch(readRDS(file.name), error=function(cond)return(NULL))   # Choose a return value in case of error
}



summarise.lmer.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),
       conv=x@optinfo$conv,
       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.lmer.output.list <- function(f ){
    output.lmer <- read.lmer.output(f)
    if(!is.null(output.lmer)){
	res <- list(files.details=files.details(f),
		   lmer.summary=summarise.lmer.output( output.lmer))
    }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])))
      # Get residual variance
      VarResid <- attr(VarCorr(i),"sc")^2
      # Calculate marginal R-squared (fixed effects/total variance)
      Rm <- VarF/(VarF+VarRand+VarResid)
      # Calculate conditional R-squared (fixed effects+random effects/total variance)
      Rc <- (VarF+VarRand)/(VarF+VarRand+VarResid)
      # Bind R^2s into a matrix and return with AIC values
      Rsquared.mat <- data.frame(Class=class(i),Family="Gaussian",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])))
# Get residual variance
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
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
VarResid <- attr(VarCorr(i),"sc")^2
var.vec <- var.vec/(VarF+VarRand+VarResid)
}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$lmer.summary$fixed.coeff.E),names(dat.t))] <-
    list.res$lmer.summary$fixed.coeff.E
dat.t[,length(names.param)+match(names(list.res$lmer.summary$fixed.coeff.E),names(dat.t))] <-
    list.res$lmer.summary$fixed.coeff.Std.Error
dat.t[,match(names(list.res$lmer.summary$fixed.var),names(dat.t))] <-
    list.res$lmer.summary$fixed.var
res <- data.frame(list.res$files.details,list.res$lmer.summary[1:7],dat.t)
return(res)
}


#################################################################
#################################################################
##### FUNCTION TO analyse the results AIC effect size

## function to compute delat R2
fun.compute.criteria.diff <- function(i,DF.results,criteria.selected){
select.simple.compet <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i] &
                    DF.results$model=='lmer.LOGLIN.simplecomp.Tf'
select.no.compet <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i] &
                    DF.results$model=='lmer.LOGLIN.nocomp.Tf'
if(sum(select.simple.compet)==1){
diff.criteria.simple.compet <- DF.results[[criteria.selected]][i] - DF.results[[criteria.selected]][
                                                                                select.simple.compet]
}else{
diff.criteria.simple.compet <- NA
}
if(sum(select.no.compet)==1){
diff.criteria.no.compet <- DF.results[[criteria.selected]][i] - DF.results[[criteria.selected]][
                                                                              select.no.compet]
}else{
diff.criteria.no.compet <- NA
}

df.res <- data.frame(diff.criteria.simple.compet,diff.criteria.no.compet)
names(df.res) <- paste(criteria.selected,c('simplecomp','nocomp'),sep=".")

return(df.res)
}



fun.compute.delta.AIC <- function(i,DF.results){
select.model.trait.fill <- DF.results$id==DF.results$id[i] &
                    DF.results$trait==DF.results$trait[i] &
                    DF.results$filling==DF.results$filling[i]
if(sum(select.model.trait.fill)>0){
delta.AIC <- DF.results[['AIC']][i] - min(DF.results[['AIC']][select.model.trait.fill])
if (sum( DF.results$nobs[select.model.trait.fill]!=DF.results$nobs[i])>0)
    stop('no same number of observation')
}else{
delta.AIC <- NA
}
df.res <- data.frame(delta.AIC=delta.AIC)
return(df.res)
}

## function to compute ratio of variance explained by a trait variable over the variance explained by the BATOT
fun.ratio.var.fixed.effect <- function(DF.results){
mat.ratio <- DF.results[,c('sumTnTfBn.abs.VAR','sumTfBn.VAR','sumTnBn.VAR','effect.response.var')]/
    DF.results[,'sumBn.VAR']
names(mat.ratio) <- c('abs.dist','Response','Effect','Effect.Response')
return(mat.ratio)
}

### FUNCTION TO REPORT BEST MODEL PER ECOREGION AND TRAITS
fun.AIC <- function(id2.one,DF.results){
 models <- c('lmer.LOGLIN.nocomp.Tf', 'lmer.LOGLIN.simplecomp.Tf','lmer.LOGLIN.HD.Tf',
             'lmer.LOGLIN.E.Tf','lmer.LOGLIN.R.Tf','lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
 best <- as.vector(DF.results[DF.results$id2==id2.one,c('id2','trait','set','ecocode','filling','MAT','MAP','model')])[which.min(DF.results$AIC[DF.results$id2==id2.one]),]
 AIC.all <- as.vector(DF.results[DF.results$id2==id2.one,c('AIC')])
 names(AIC.all) <- as.vector(DF.results[DF.results$id2==id2.one,c('model')])
 AIC.all <- AIC.all[models]-min(AIC.all)
 res <- data.frame((best),t(AIC.all))
 names(res) <- c('id2','trait','set','ecocode','filling','MAT','MAP','best.model',paste('AIC',models,sep='.'))
 return(res)
}

fun.AICc <- function(id2.one,DF.results){
 models <- c('lmer.LOGLIN.nocomp.Tf', 'lmer.LOGLIN.simplecomp.Tf','lmer.LOGLIN.HD.Tf','lmer.LOGLIN.E.Tf','lmer.LOGLIN.R.Tf','lmer.LOGLIN.ER.Tf','lmer.LOGLIN.AD.Tf')
 Deviance.all <- DF.results[DF.results$id2==id2.one,'deviance']
 names(Deviance.all) <- DF.results[DF.results$id2==id2.one,'model']
 Deviance.all <- Deviance.all[models]
 nobs.all <- DF.results[DF.results$id2==id2.one,'nobs']
 names(nobs.all) <- DF.results[DF.results$id2==id2.one,'model']
 nobs.all <- nobs.all[models]
 n.param <- c(2,3,4,4,4,5,4)
 AICc <-  Deviance.all+2*n.param*(nobs.all)/(nobs.all-n.param-1)
 id2.n <- unique(DF.results[DF.results$id2==id2.one,c('id2')])
 res <- data.frame(id2.n,models[which.min(AICc)],t(AICc),row.names=NULL)
 names(res) <- c('id2','best.model',models)
 return(res)
}

#################################333

### function to get the data for a given model with criteria to select
fun.select.ecoregions.trait <- function(DF.results,trait.name,model.selected,
                                      nobs.min=1000,filling.selected="species",
                                      threshold.delta.AIC){
DF.results[DF.results$nobs>nobs.min &
           DF.results$filling==filling.selected &
           DF.results$trait==trait.name &
           DF.results$model %in% model.selected &
           DF.results$delta.AIC<threshold.delta.AIC,]
}

### function to get the data for a given model with criteria to select only site with competition
fun.select.ecoregions.trait.compet <- function(DF.results,trait.name,model.selected,
                                      nobs.min=1000,filling.selected="species",
                                      threshold.delta.AIC){
DF.results[DF.results$nobs>nobs.min &
           DF.results$filling==filling.selected &
           DF.results$trait==trait.name &
           DF.results$model %in% model.selected &
           DF.results$sumBn < 0.0 &
           ## DF.results$delta.AIC==0,]
           DF.results$delta.AIC<threshold.delta.AIC,]
}




#########################
##### FUNCTIONS FOR PLOTS
fun.plot.lmer.res.x.y <- function(model.selected,trait.name,DF.results,var.x,var.y,threshold.delta.AIC,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected,threshold.delta.AIC=10000000)
plot(df.selected[[var.x]],df.selected[[var.y]],...)
}

fun.plot.lmer.res.x.y.2 <- function(model.selected,trait.name,DF.results,var.x,var.y,col.vec,pch.AIC=TRUE,cex.AIC=TRUE,col.set=TRUE,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected,threshold.delta.AIC=10000000)
if(pch.AIC) {pch.vec <- c(1,16)[as.numeric(df.selected[['delta.AIC']]==0)+1]}else{pch.vec <- 1}
if(cex.AIC) {cex.vec <- c(1,1.5)[as.numeric(df.selected[['delta.AIC']]==0)+1]}else{cex.vec <- 1}
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
plot(df.selected[[var.x]],df.selected[[var.y]],
     pch=pch.vec,
     cex=cex.vec,
     col=col.vec2,...)
}

fun.plot.lmer.res.boxplot <- function(model.selected,trait.name,DF.results,var.y,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
boxplot(df.selected[[var.y]],...)
}

fun.plot.lmer.res.boxplot.compare.model <- function(model.selected,trait.name,DF.results,var.y,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
boxplot(df.selected[[var.y]]~df.selected[['model']],...)
## compute percentage of model beteer than null
print(tapply(df.selected[['AIC.simplecomp']]<0,
       INDEX=df.selected[['model']],
       FUN=sum)/tapply(df.selected[['AIC.simplecomp']]<0,
                       INDEX=df.selected[['model']],
                       FUN=length))
}




fun.plot.param.error.bar <- function(df.selected,var.x,param,small.bar,col.vec){
segments( df.selected[[var.x]],df.selected[[param]] - 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],
          df.selected[[var.x]],df.selected[[param]] + 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
segments( df.selected[[var.x]]-small.bar,df.selected[[param]] - 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],
          df.selected[[var.x]]+small.bar,df.selected[[param]] - 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
segments( df.selected[[var.x]]-small.bar,df.selected[[param]] + 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],
          df.selected[[var.x]]+small.bar,df.selected[[param]] + 1.96*df.selected[[paste(param,"Std.Error",sep=".")]],col=col.vec)
}

fun.plot.all.param.x.y.c <- function(model.selected,trait.name,DF.results,var.x,params,
                                     small.bar,threshold.delta.AIC,col.vec,col.set=TRUE,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected,threshold.delta.AIC=threshold.delta.AIC)
if(col.set) {col.vec2 <- col.vec[unclass(df.selected[['set']])]}else{col.vec2 <- 1}
ylim <- range(c(df.selected[[params[1]]] - 1.96*df.selected[[paste(params[1],"Std.Error",sep=".")]],
                df.selected[[params[1]]] + 1.96*df.selected[[paste(params[1],"Std.Error",sep=".")]]),na.rm=TRUE)
plot(df.selected[[var.x]],df.selected[[params[1]]],col=col.vec2,ylim=ylim,...)
fun.plot.param.error.bar(df.selected,var.x,param=params[1],small.bar,col=col.vec2)
}


fun.plot.all.param.boxplot <- function(model.selected,trait.name,DF.results,params,small.bar,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
 if(length(params)>1){
DF.t <- data.frame(param=rep(names(params),each=nrow(df.selected)),value=c(df.selected[[params[1]]],df.selected[[params[2]]]))
boxplot(DF.t[['value']]~DF.t[['param']],...)
}else{
boxplot(df.selected[[params[1]]],...)
}
}

fun.plot.all.param.er.diff.MAP <- function(model.selected,trait.name,DF.results,...){
df.selected <- fun.select.ecoregions.trait(DF.results,trait.name=trait.name,model.selected=model.selected)
plot(df.selected[['MAP']],df.selected[['sumTnBn']]-df.selected[['sumTfBn']],...)
}



fun.plot.panel.lmer.res.x.y <- function(models,traits,DF.results,var.x,var.y.l,express,...){
    ncols = length(traits)
    nrows = length(models)
    list.models <- as.list(names(models))
    names(list.models) <- rep('model',length(list.models))
    DF.results$set <- factor(DF.results$set)
    col.vec <- niceColors(n=nlevels(DF.results$set))
    par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
    for(i in 1:nrows)
        for(j in 1:ncols){
              fun.plot.lmer.res.x.y.2(models[i],traits[j],
                          DF.results,var.x,var.y=var.y.l[[i]],col.vec,...)
              abline(h=0,lty=3)
            if(i==1 )
                mtext(traits[j], side=3, line =1)
            if(i==nrows)
                mtext(var.x, side=1, line =4)
            if(j==1)
                mtext(paste('Effect size',list.models[i]), side=2, line =4,cex=0.9)
            if(i==nrows & j==ncols)
                legend('topright',legend=levels(DF.results$set),pch=16,
                       col=col.vec,bty='n',ncol=2)
        }
}


fun.plot.panel.lmer.res.boxplot <- function(models,traits,DF.results,var.y,express,...){
    ncols = length(traits)
    nrows = length(models)
    list.models <- as.list(names(models))
    names(list.models) <- rep('model',length(list.models))
    par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
    for(i in 1:nrows)
        for(j in 1:ncols){
              fun.plot.lmer.res.boxplot(models[i],traits[j],
                          DF.results,var.x,var.y,...)
              abline(h=0,lty=3)
            if(i==1 )
                mtext(traits[j], side=3, line =1)
            if(j==1)
                mtext(paste("Effect size", list.models[i],sep=" "), side=2, line =4,cex=0.9)
        }
}


fun.plot.panel.lmer.res.boxplot.compare <- function(models,traits,DF.results,var.y,express,...){
    ncols = length(traits)
    list.models <- as.list(names(models))
    names(list.models) <- rep('model',length(list.models))
    par(mfrow = c(1, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
    for(j in 1:ncols){
              fun.plot.lmer.res.boxplot.compare.model(models,traits[j],
                          DF.results,var.y,names=names(models),...)
              abline(h=0,lty=3)
                mtext(traits[j], side=3, line =1,cex=2)
            if(j==1 )
                mtext("Effect size", side=2, line =4,cex=1.5)
        }
}

fun.plot.panel.lmer.parameters.c <- function(models,traits,DF.results,var.x,list.params,threshold.delta.AIC,small.bar=10,...){
    ncols = length(traits)
    nrows = length(models)
    par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
    DF.results$set <- factor(DF.results$set)
    col.vec <- niceColors(n=nlevels(DF.results$set))

## ### TO COMPARE THE PARAMTERS WE NEED TO DIVIDE THE ABS>DIST per two
    ## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
    for(i in 1:nrows)
        for(j in 1:ncols){
              fun.plot.all.param.x.y.c(models[i],traits[j],DF.results,var.x,params=list.params[[i]],
                                       small.bar=small.bar,
                                       threshold.delta.AIC=threshold.delta.AIC,col.vec=col.vec,col.set=TRUE,...)
              abline(h=0,lty=3)
              if(length(list.params[[i]])>1)
                  legend("topright",names(list.params[[i]]),
                         pch=rep(1,length(list.params[[i]])),
                         col=1:length(list.params[[i]]),bty='n',cex=1)
            if(i==1 )
                mtext(traits[j], side=3, line =1)
            if(i==nrows)
                mtext(var.x, side=1, line =4)
            if(j==1)
                mtext(paste('param',names(models)[i]), side=2, line =4,cex=1)
            if(i==nrows & j==ncols)
                legend('topright',legend=levels(DF.results$set),pch=16,
                       col=col.vec,bty='n',ncol=2)
        }
}



fun.plot.panel.lmer.parameters.boxplot <- function(models,traits,DF.results,list.params,small.bar=10,...){
    ncols = length(traits)
    nrows = length(models)
    par(mfrow = c(nrows, ncols), mar = c(2,2,1,1), oma = c(4,4,4,1) )
### TO COMPARE THE PARAMTERS WE NEED TO DIVIDE THE ABS>DIST per two
    ## DF.results$sumTnTfBn.abs <- DF.results$sumTnTfBn.abs/2
    for(i in 1:nrows)
        for(j in 1:ncols){
              fun.plot.all.param.boxplot(models[i],traits[j],DF.results,params=list.params[[i]],small.bar=small.bar,...)
              abline(h=0,lty=3)
            if(i==1 )
                mtext(traits[j], side=3, line =1)
            if(j==1)
                mtext(paste('param',names(models)[i]), side=2, line =4,cex=1)
        }
}