# HG changeset patch
# User mnhn65mo
# Date 1533299302 14400
# Node ID 5b126f770671fffff3ac493e96e2c136742c0f3f

Uploaded

diff -r 000000000000 -r 5b126f770671 ab_index.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ab_index.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,17 @@
+#!/usr/bin/env Rscript
+#library("getopt")
+#library(devtools)
+#library(RegionalGAM)
+
+args = commandArgs(trailingOnly=TRUE)
+source(args[1])
+
+
+tryCatch({input = read.table(args[2], header=TRUE,sep=" ")},finally={input = read.table(args[2], header=TRUE,sep=",")})
+pheno = read.table(args[3], header=TRUE,sep="	")
+
+if("TREND" %in% colnames(input)){
+    input <- input[input$TREND==1,c("SPECIES","SITE","YEAR","MONTH","DAY","COUNT")]
+}
+data.index <- abundance_index(input, pheno)
+write.table(data.index, file="data.index", row.names=FALSE, sep="	")
diff -r 000000000000 -r 5b126f770671 ab_index_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/ab_index_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,39 @@
+<tool id="ab_index" name="Abundance index" version="0.1">
+    <description>computation across species, sites and years</description>
+    <requirements>
+        <requirement type="package" version="3.4">r</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/ab_index.R' '$__tool_directory__/dennis_gam_initial_functions.R' '$input1' '$input2' '$output' ]]>
+    </command>
+    <inputs>
+        <param format="tabular,csv" name="input1" type="data" label="Count file" help="The file must contain the SPECIES, SITE, YEAR, MONTH, DAY and COUNT columns. CSV file separator ',' or tab."/>
+        <param format="tabular" name="input2" type="data" label="Flight curve output"/>
+    </inputs>
+    <outputs>
+        <data format="tabular" name="output" from_work_dir="data.index" />
+    </outputs>
+    <tests>
+        <test> <!-- FAILED -->
+            <param name="input1" value="gatekeeper_CM"/>
+            <param name="input2" value="flight_curve_result.tabular"/>
+            <output name="output" value="data_index.tabular"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+
+.. class:: infomark 
+
+==========================
+Compute abundance index
+==========================
+
+This tool is an implementation of the abundance_index function from `RegionalGAM package:  https://github.com/RetoSchmucki/regionalGAM/
+
+This function compute the Abundance Index across sites and years from your dataset and the regional flight curve.
+
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+    </citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 autocorr-res_acf.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/autocorr-res_acf.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,10 @@
+#!/usr/bin/env Rscript
+library(nlme)
+library(MASS)
+
+args = commandArgs(trailingOnly=TRUE)
+load(args[1])
+
+png('output-acf.png')
+graph<-acf(residuals(mod,type="normalized"))
+invisible(dev.off())
diff -r 000000000000 -r 5b126f770671 autocorr-res_acf_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/autocorr-res_acf_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,39 @@
+<tool id="autocor-acf" name="Autocorrelation test" version="0.1">
+    <description>check for temporal autocorrelation in the residuals</description>
+    <requirements>
+        <requirement type="package" version="3.1_131">r-nlme</requirement>
+        <requirement type="package" version="7.3_48">r-mass</requirement>
+        <requirement type="package" version="0.9.10">xorg-libxrender</requirement>
+        <requirement type="package" version="1.1.2">xorg-libsm</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/autocorr-res_acf.R' '$input1' '$output' ]]>
+    </command>
+    <inputs>
+        <param format="rdata" name="input1" type="data" label="gls model" help="Rdata file"/>
+    </inputs>
+    <outputs>
+        <data format="png" name="output" from_work_dir="output-acf.png" />
+    </outputs>
+    <tests>
+        <test>
+            <param name="input1" value="mod1"/>
+            <param output="output" value="res.png"/>
+        </test>
+    </tests>
+
+    <help><![CDATA[
+
+.. class:: infomark 
+
+==========================
+Model residuals autocorrelation check 
+==========================
+
+This tool is an implementation of the autocorr-res_acf function from `RegionalGAM package:  https://github.com/RetoSchmucki/regionalGAM/
+This function uses a simple linear model and explore for temporal autocorrelation that we will account in the final model.
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+</citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 dennis_gam_initial_functions.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/dennis_gam_initial_functions.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,638 @@
+### R-Script Adapted from script provided by the CEH, UK BY: Reto Schmucki [ reto.schmucki@mail.mcgill.ca]
+### DATE: 14 July 2014 function to run two stage model in DENNIS et al. 2013
+
+
+.onAttach <- function(libname, pkgname) {
+  packageStartupMessage(" The regionalGAM package that is no longer maintained, \n use the new rbms package instead. \n
+   devtools::install_github(\"RetoSchmucki/rbms\", force=TRUE)")
+}
+
+
+#' year_day_func Function
+#' This function generate the full sequence of days, months and include the observation to that file.
+#' @param sp_data A data.frame with your observation.
+#' @keywords year days
+#' @export
+#' @author Reto Schmucki
+#' @examples
+#' year_day_func()
+
+
+# FUNCTIONS
+
+year_day_func = function(sp_data) {
+
+    year <- unique(sp_data$YEAR)
+
+    origin.d <- paste(year, "01-01", sep = "-")
+    if ((year%%4 == 0) & ((year%%100 != 0) | (year%%400 == 0))) {
+        nday <- 366
+    } else {
+        nday <- 365
+    }
+
+    date.serie <- as.POSIXlt(seq(as.Date(origin.d), length = nday, by = "day"), format = "%Y-%m-%d")
+
+    dayno <- as.numeric(julian(date.serie, origin = as.Date(origin.d)) + 1)
+    month <- as.numeric(strftime(date.serie, format = "%m"))
+    week <- as.numeric(strftime(date.serie, format = "%W"))
+    week_day <- as.numeric(strftime(date.serie, format = "%u"))
+    day <- as.numeric(strftime(date.serie, format = "%d"))
+
+    site_list <- sp_data[!duplicated(sp_data$SITE), c("SITE")]
+
+    all_day_site <- data.frame(SPECIES = sp_data$SPECIES[1], SITE = rep(site_list, rep(nday, length(site_list))),
+        YEAR = sp_data$YEAR[1], MONTH = month, WEEK = week, DAY = day, DAY_WEEK = week_day, DAYNO = dayno,
+        COUNT = NA)
+
+    count_index <- match(paste(sp_data$SITE, sp_data$DAYNO, sep = "_"), paste(all_day_site$SITE, all_day_site$DAYNO,
+        sep = "_"))
+    all_day_site$COUNT[count_index] <- sp_data$COUNT
+    site_count_length <- aggregate(sp_data$COUNT, by = list(sp_data$SITE), function(x) list(1:length(x)))
+    names(site_count_length$x) <- as.character(site_count_length$Group.1)
+    site_countno <- utils::stack(site_count_length$x)
+    all_day_site$COUNTNO <- NA
+    all_day_site$COUNTNO[count_index] <- site_countno$values  # add count number to ease extraction of single count
+
+    # Add zero to close observation season two weeks before and after the first and last
+    first_obs <- min(all_day_site$DAYNO[!is.na(all_day_site$COUNT)])
+    last_obs <- max(all_day_site$DAYNO[!is.na(all_day_site$COUNT)])
+
+    closing_season <- c((first_obs - 11):(first_obs - 7), (last_obs + 7):(last_obs + 11))
+
+    # If closing season is before day 1 or day 365, simply set the first and last 5 days to 0
+    if (min(closing_season) < 1)
+        closing_season[1:5] <- c(1:5)
+    if (max(closing_season) > nday)
+        closing_season[6:10] <- c((nday - 4):nday)
+
+    all_day_site$COUNT[all_day_site$DAYNO %in% closing_season] <- 0
+    all_day_site$ANCHOR <- 0
+    all_day_site$ANCHOR[all_day_site$DAYNO %in% closing_season] <- 1
+
+    all_day_site <- all_day_site[order(all_day_site$SITE, all_day_site$DAYNO), ]
+
+    return(all_day_site)
+}
+
+
+#' trap_area Function
+#'
+#' This function compute the area under the curve using the trapezoid method.
+#' @param x A vector or a two columns matrix.
+#' @param y A vector, Default is NULL
+#' @keywords trapezoid
+#' @export
+#' @examples
+#' trap_area()
+
+
+trap_area = function(x, y = NULL) {
+    # If y is null and x has multiple columns then set y to x[,2] and x to x[,1]
+    if (is.null(y)) {
+        if (length(dim(x)) == 2) {
+            y = x[, 2]
+            x = x[, 1]
+        } else {
+            stop("ERROR: need to either specifiy both x and y or supply a two column data.frame/matrix to x")
+        }
+    }
+
+    # Check x and y are same length
+    if (length(x) != length(y)) {
+        stop("ERROR: x and y need to be the same length")
+    }
+
+    # Need to exclude any pairs that are NA for either x or y
+    rm_inds = which(is.na(x) | is.na(y))
+    if (length(rm_inds) > 0) {
+        x = x[-rm_inds]
+        y = y[-rm_inds]
+    }
+
+    # Determine values of trapezoids under curve Get inds
+    inds = 1:(length(x) - 1)
+    # Determine area using trapezoidal rule Area = ( (b1 + b2)/2 ) * h where b1 and b2 are lengths of bases
+    # (the parallel sides) and h is the height (the perpendicular distance between two bases)
+    areas = ((y[inds] + y[inds + 1])/2) * diff(x)
+
+    # total area is sum of all trapezoid areas
+    tot_area = sum(areas)
+
+    # Return total area
+    return(tot_area)
+}
+
+
+#' trap_index Function
+#'
+#' This function compute the area under the curve (Abundance Index) across species, sites and years
+#' @param sp_data A data.frame containing species count data generated from the year_day_func()
+#' @param y A vector, Default is NULL
+#' @keywords Abundance index
+#' @export
+#' @examples
+#' trap_index()
+
+
+
+trap_index = function(sp_data, data_col = "IMP", time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR")) {
+
+    # Build output data.frame
+    out_obj = unique(sp_data[, by_col])
+    # Set row.names to be equal to collapsing of output rows (will be unique, you need them to make uploading
+    # values back to data.frame will be easier)
+    row.names(out_obj) = apply(out_obj, 1, paste, collapse = "_")
+
+    # Using this row.names from out_obj above as index in by function to loop through values all unique combs
+    # of by_cols and fit trap_area to data
+    ind_dat = by(sp_data[, c(time_col, data_col)], apply(sp_data[, by_col], 1, paste, collapse = "_"), trap_area)
+
+    # Add this data to output object
+    out_obj[names(ind_dat), "SINDEX"] = round(ind_dat/7, 1)
+
+    # Set row.names to defaults
+    row.names(out_obj) = NULL
+
+    # Return output object
+    return(out_obj)
+}
+
+
+#' flight_curve Function
+#' This function compute the flight curve across and years
+#' @param your_dataset A data.frame containing original species count data. The data format is a csv (comma "," separated) with 6 columns with headers, namely "species","transect_id","visit_year","visit_month","visit_day","sp_count"
+#' @keywords standardize flight curve
+#' @export
+#' @examples
+#' flight_curve()
+
+
+flight_curve <- function(your_dataset) {
+
+    if("mgcv" %in% installed.packages() == "FALSE") {
+        print("mgcv package is not installed.")
+        x <- readline("Do you want to install it? Y/N")
+        if (x == 'Y') { 
+            install.packages("mgcv")
+        }
+        if (x == 'N') {
+            stop("flight curve can not be computed without the mgcv package, sorry")
+        }
+    }
+    your_dataset$DAYNO <- strptime(paste(your_dataset$DAY, your_dataset$MONTH,
+        your_dataset$YEAR, sep = "/"), "%d/%m/%Y")$yday + 1
+    dataset <- your_dataset[, c("SPECIES", "SITE", "YEAR", "MONTH",
+        "DAY", "DAYNO", "COUNT")]
+    sample_year <- unique(dataset$YEAR)
+    sample_year <- sample_year[order(sample_year)]
+    if (length(sample_year) >1 ) {
+        for (y in sample_year) {
+            dataset_y <- dataset[dataset$YEAR == y, ]
+            nsite <- length(unique(dataset_y$SITE))
+            # Determine missing days and add to dataset
+            sp_data_all <- year_day_func(dataset_y)
+            if (nsite > 200) {
+                sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
+                    200, replace = F)]), ]
+                sp_data_all <- sp_data_all
+            }
+            sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
+            print(paste("Fitting the GAM for",as.character(sp_data_all$SPECIES[1]),"and year",y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time()))
+            if(length(unique(sp_data_all$SITE))>1){
+                gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
+                    data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+            } 
+            else {
+                gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr")  -1,
+                    data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+            }
+            # Give a second try if the GAM does not converge the first time
+            if (class(gam_obj_site)[1] == "try-error") {
+            # Determine missing days and add to dataset
+                sp_data_all <- year_day_func(dataset_y)
+                if (nsite > 200) {
+                   sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
+                    200, replace = F)]), ]        
+                } 
+                else {
+                    sp_data_all <- sp_data_all
+                }
+                sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
+                print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time(),"second try"))
+                if(length(unique(sp_data_all$SITE))>1){
+                    gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
+                       data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+                } 
+                else {
+                    gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr")  -1,
+                        data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+                }
+                if (class(gam_obj_site)[1] == "try-error") {
+                    print(paste("Error in fitting the flight period for",sp_data_all$SPECIES[1],"at year", y,"no convergence after two trial"))
+                    sp_data_all[, "FITTED"] <- NA
+                    sp_data_all[, "COUNT_IMPUTED"] <- NA
+                    sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                    sp_data_all[, "NM"] <- NA
+                } 
+                else {
+                    # Generate a list of values for all days from the additive model and use
+                    # these value to fill the missing observations
+                    sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
+                        c("trimDAYNO", "SITE")], type = "response")
+                    # force zeros at the beginning end end of the year
+                    sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
+                    sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
+                    # if infinite number in predict replace with NA.
+                    if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
+                        print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
+                        sp_data_all[, "FITTED"] <- NA
+                        sp_data_all[, "COUNT_IMPUTED"] <- NA
+                        sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                        sp_data_all[, "NM"] <- NA    
+                    }
+                    else {
+                    sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
+                    sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
+                    # Define the flight curve from the fitted values and append them over
+                    # years (this is one flight curve per year for all site)
+                    site_sums <- aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
+                        FUN = sum)
+                    # Rename sum column
+                    names(site_sums)[names(site_sums) == "x"] <- "SITE_YR_FSUM"
+                    # Add data to sp_data data.frame (ensure merge does not sort the data!)
+                    sp_data_all = merge(sp_data_all, site_sums, by <- c("SITE"),
+                        all = TRUE, sort = FALSE)
+                    # Calculate normalized values
+                    sp_data_all[, "NM"] <- sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
+                    }
+                }
+            } 
+            else {
+                # Generate a list of values for all days from the additive model and use
+                # these value to fill the missing observations
+                sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
+                    c("trimDAYNO", "SITE")], type = "response")
+                # force zeros at the beginning end end of the year
+                sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
+                sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
+                # if infinite number in predict replace with NA.
+                if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
+                    print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
+                    sp_data_all[, "FITTED"] <- NA
+                    sp_data_all[, "COUNT_IMPUTED"] <- NA
+                    sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                    sp_data_all[, "NM"] <- NA    
+                }
+                else {
+                sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
+                sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
+                # Define the flight curve from the fitted values and append them over
+                # years (this is one flight curve per year for all site)
+                site_sums = aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
+                    FUN = sum)
+                # Rename sum column
+                names(site_sums)[names(site_sums) == "x"] = "SITE_YR_FSUM"
+                # Add data to sp_data data.frame (ensure merge does not sort the data!)
+                sp_data_all = merge(sp_data_all, site_sums, by = c("SITE"), all = TRUE,
+                    sort = FALSE)
+                # Calculate normalized values
+                sp_data_all[, "NM"] = sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
+                }
+            }
+            sp_data_filled <- sp_data_all
+            flight_curve <- data.frame(species = sp_data_filled$SPECIES, year = sp_data_filled$YEAR,
+                week = sp_data_filled$WEEK, DAYNO = sp_data_filled$DAYNO, DAYNO_adj = sp_data_filled$trimDAYNO,
+                nm = sp_data_filled$NM)[!duplicated(paste(sp_data_filled$YEAR,
+                sp_data_filled$DAYNO, sep = "_")), ]
+            flight_curve <- flight_curve[order(flight_curve$DAYNO), ]
+            # bind if exist else create
+            if (is.na(flight_curve$nm[1]))  next()
+            if ("flight_pheno" %in% ls()) {
+                flight_pheno <- rbind(flight_pheno, flight_curve)
+            } 
+            else {
+            flight_pheno <- flight_curve
+            }
+        }  # end of year loop
+    } 
+    else {
+        y <- unique(dataset$YEAR)
+        dataset_y <- dataset[dataset$YEAR == y, ]
+        nsite <- length(unique(dataset_y$SITE))
+        # Determine missing days and add to dataset
+        sp_data_all <- year_day_func(dataset_y)
+        if (nsite > 200) {
+            sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
+            200, replace = F)]), ]
+        } 
+        else {
+            sp_data_all <- sp_data_all
+        }
+        sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
+        print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,":",Sys.time()))
+        if(length(unique(sp_data_all$SITE))>1){
+            gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) -1,
+            data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+        } 
+        else {
+            gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr")  -1,
+            data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+        }
+        # Give a second try if the GAM does not converge the first time
+        if (class(gam_obj_site)[1] == "try-error") {
+            # Determine missing days and add to dataset
+            sp_data_all <- year_day_func(dataset_y)
+            if (nsite > 200) {
+                sp_data_all <- sp_data_all[as.character(sp_data_all$SITE) %in% as.character(unique(dataset_y$SITE)[sample(1:nsite,
+                200, replace = F)]), ]
+            } 
+            else {
+                sp_data_all <- sp_data_all
+            }
+            sp_data_all$trimDAYNO <- sp_data_all$DAYNO - min(sp_data_all$DAYNO) + 1
+            print(paste("Fitting the GAM for",sp_data_all$SPECIES[1],"at year", y,"with",length(unique(sp_data_all$SITE)),"sites :",Sys.time(),"second try"))
+            if(length(unique(sp_data_all$SITE))>1){
+                gam_obj_site <- try(mgcv::bam(COUNT ~ s(trimDAYNO, bs = "cr") + as.factor(SITE) - 1,
+                data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+            } 
+            else {
+                gam_obj_site <- try(mgcv::gam(COUNT ~ s(trimDAYNO, bs = "cr")  -1,
+                data = sp_data_all, family = poisson(link = "log")), silent = TRUE)
+            }
+            if (class(gam_obj_site)[1] == "try-error") {
+                print(paste("Error in fitting the flight period for",sp_data_all$SPECIES[1],"at year", y,"no convergence after two trial"))
+                sp_data_all[, "FITTED"] <- NA
+                sp_data_all[, "COUNT_IMPUTED"] <- NA
+                sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                sp_data_all[, "NM"] <- NA
+            } 
+            else {
+                # Generate a list of values for all days from the additive model and use
+                # these value to fill the missing observations
+                sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
+                c("trimDAYNO", "SITE")], type = "response")
+                # force zeros at the beginning end end of the year
+                sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
+                sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
+                # if infinite number in predict replace with NA.
+                if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
+                    print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
+                    sp_data_all[, "FITTED"] <- NA
+                    sp_data_all[, "COUNT_IMPUTED"] <- NA
+                    sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                    sp_data_all[, "NM"] <- NA    
+                }
+                else {
+                    sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
+                    sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
+                    # Define the flight curve from the fitted values and append them over
+                    # years (this is one flight curve per year for all site)
+                    site_sums <- aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
+                    FUN = sum)
+                    # Rename sum column
+                    names(site_sums)[names(site_sums) == "x"] <- "SITE_YR_FSUM"
+                    # Add data to sp_data data.frame (ensure merge does not sort the data!)
+                    sp_data_all = merge(sp_data_all, site_sums, by <- c("SITE"),
+                    all = TRUE, sort = FALSE)
+                    # Calculate normalized values
+                    sp_data_all[, "NM"] <- sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
+                }          
+            }
+        } 
+        else {
+            # Generate a list of values for all days from the additive model and use
+            # these value to fill the missing observations
+            sp_data_all[, "FITTED"] <- mgcv::predict.gam(gam_obj_site, newdata = sp_data_all[,
+            c("trimDAYNO", "SITE")], type = "response")
+            # force zeros at the beginning end end of the year
+            sp_data_all[sp_data_all$trimDAYNO < 60, "FITTED"] <- 0
+            sp_data_all[sp_data_all$trimDAYNO > 305, "FITTED"] <- 0
+            # if infinite number in predict replace with NA.
+            if(sum(is.infinite(sp_data_all[, "FITTED"]))>0){
+                print(paste("Error in the flight period for",sp_data_all$SPECIES[1],"at year", y,"weird predicted values"))
+                sp_data_all[, "FITTED"] <- NA
+                sp_data_all[, "COUNT_IMPUTED"] <- NA
+                sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- NA
+                sp_data_all[, "NM"] <- NA    
+            }
+            else {
+                sp_data_all[, "COUNT_IMPUTED"] <- sp_data_all$COUNT
+                sp_data_all[is.na(sp_data_all$COUNT), "COUNT_IMPUTED"] <- sp_data_all$FITTED[is.na(sp_data_all$COUNT)]
+                # Define the flight curve from the fitted values and append them over
+                # years (this is one flight curve per year for all site)
+                site_sums = aggregate(sp_data_all$FITTED, by = list(SITE = sp_data_all$SITE),
+                FUN = sum)
+                # Rename sum column
+                names(site_sums)[names(site_sums) == "x"] = "SITE_YR_FSUM"
+                # Add data to sp_data data.frame (ensure merge does not sort the data!)
+                sp_data_all = merge(sp_data_all, site_sums, by = c("SITE"), all = TRUE,
+                sort = FALSE)
+                # Calculate normalized values
+                sp_data_all[, "NM"] = sp_data_all$FITTED/sp_data_all$SITE_YR_FSUM
+            }
+        }
+        sp_data_filled <- sp_data_all
+        flight_curve <- data.frame(species = sp_data_filled$SPECIES, year = sp_data_filled$YEAR,
+        week = sp_data_filled$WEEK, DAYNO = sp_data_filled$DAYNO, DAYNO_adj = sp_data_filled$trimDAYNO,
+        nm = sp_data_filled$NM)[!duplicated(paste(sp_data_filled$YEAR,
+        sp_data_filled$DAYNO, sep = "_")), ]
+        flight_curve <- flight_curve[order(flight_curve$DAYNO), ]
+        if (is.na(flight_curve$nm[1])){
+            flight_pheno <- data.frame()
+        }
+        else {
+            # bind if exist else create
+            if ("flight_pheno" %in% ls()) {
+                flight_pheno <- rbind(flight_pheno, flight_curve)
+            } 
+            else {
+                flight_pheno <- flight_curve
+            }
+        }
+    }
+    return(flight_pheno)
+}
+
+
+#' abundance_index Function
+#'
+#' This function compute the Abundance Index across sites and years from your dataset and the regional flight curve
+#' @param your_dataset A data.frame containing original species count data. The data format is a csv (comma "," separated) with 6 columns with headers, namely	"species","transect_id","visit_year","visit_month","visit_day","sp_count"
+#' @param flight_pheno A data.frame for the regional flight curve computed with the function flight_curve
+#' @keywords standardize flight curve
+#' @export
+#' @examples
+#' abundance_index()
+
+abundance_index <- function(your_dataset,flight_pheno) {
+
+your_dataset$DAYNO <- strptime(paste(your_dataset$DAY, your_dataset$MONTH,
+    your_dataset$YEAR, sep = "/"), "%d/%m/%Y")$yday + 1
+
+dataset <- your_dataset[, c("SPECIES", "SITE", "YEAR", "MONTH",
+    "DAY", "DAYNO", "COUNT")]
+
+sample_year <- unique(dataset$YEAR)
+sample_year <- sample_year[order(sample_year)]
+
+
+if (length(sample_year)>1){
+
+for (y in sample_year) {
+
+    year_pheno <- flight_pheno[flight_pheno$year == y, ]
+
+    dataset_y <- dataset[dataset$YEAR == y, ]
+
+    sp_data_site <- year_day_func(dataset_y)
+    sp_data_site$trimDAYNO <- sp_data_site$DAYNO - min(sp_data_site$DAYNO) + 1
+
+    sp_data_site <- merge(sp_data_site, year_pheno[, c("DAYNO", "nm")],
+        by = c("DAYNO"), all.x = TRUE, sort = FALSE)
+
+    # compute proportion of the flight curve sampled due to missing visits
+    pro_missing_count <- data.frame(SITE = sp_data_site$SITE, WEEK = sp_data_site$WEEK,
+        NM = sp_data_site$nm, COUNT = sp_data_site$COUNT, ANCHOR = sp_data_site$ANCHOR)
+    pro_missing_count$site_week <- paste(pro_missing_count$SITE, pro_missing_count$WEEK,
+        sep = "_")
+    siteweeknocount <- aggregate(pro_missing_count$COUNT, by = list(pro_missing_count$site_week),
+        function(x) sum(!is.na(x)) == 0)
+    pro_missing_count <- pro_missing_count[pro_missing_count$site_week %in%
+        siteweeknocount$Group.1[siteweeknocount$x == TRUE], ]
+    pro_count_agg <- aggregate(pro_missing_count$NM, by = list(pro_missing_count$SITE),
+        function(x) 1 - sum(x, na.rm = T))
+    names(pro_count_agg) <- c("SITE", "PROP_PHENO_SAMPLED")
+
+    # remove samples outside the monitoring window
+    sp_data_site$COUNT[sp_data_site$nm==0] <- NA
+
+    # Compute the regional GAM index
+
+    if(length(unique(sp_data_site$SITE))>1){
+        glm_obj_site <- glm(COUNT ~ factor(SITE) + offset(log(nm)) - 1, data = sp_data_site,
+            family = quasipoisson(link = "log"), control = list(maxit = 100))
+    } else {
+       glm_obj_site <- glm(COUNT ~ offset(log(nm)) - 1, data = sp_data_site,
+           family = quasipoisson(link = "log"), control = list(maxit = 100))
+    }
+
+    sp_data_site[, "FITTED"] <- predict.glm(glm_obj_site, newdata = sp_data_site,
+        type = "response")
+    sp_data_site[, "COUNT_IMPUTED"] <- sp_data_site$COUNT
+    sp_data_site[is.na(sp_data_site$COUNT), "COUNT_IMPUTED"] <- sp_data_site$FITTED[is.na(sp_data_site$COUNT)]
+
+   ## add fitted value for missing mid-week data
+    sp_data_site <- sp_data_site[!paste(sp_data_site$DAY_WEEK, sp_data_site$COUNT) %in%
+        c("1 NA", "2 NA", "3 NA", "5 NA", "6 NA", "7 NA"), ]
+
+    ## remove all added mid-week values for weeks with real count
+    ## (observation)
+    sp_data_site$site_week <- paste(sp_data_site$SITE, sp_data_site$WEEK,
+        sep = "_")
+    siteweekcount <- aggregate(sp_data_site$COUNT, by = list(sp_data_site$site_week),
+        function(x) sum(!is.na(x)) > 0)
+    sp_data_site <- sp_data_site[!(is.na(sp_data_site$COUNT) & (sp_data_site$site_week %in%
+        siteweekcount$Group.1[siteweekcount$x == TRUE])), names(sp_data_site) !=
+        "site_week"]
+
+    ## Compute the regional GAM index
+    print(paste("Compute index for",sp_data_site$SPECIES[1],"at year", y,"for",length(unique(sp_data_site$SITE)),"sites:",Sys.time()))
+    regional_gam_index <- trap_index(sp_data_site, data_col = "COUNT_IMPUTED",
+        time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR"))
+
+    cumu_index <- merge(regional_gam_index, pro_count_agg, by = c("SITE"),
+        all.x = TRUE, sort = FALSE)
+    names(cumu_index) <- c("SITE", "SPECIES", "YEAR", "regional_gam", "prop_pheno_sampled")
+
+    cumu_index <- cumu_index[order(cumu_index$SITE), ]
+
+    # bind if exist else create
+    if ("cumullated_indices" %in% ls()) {
+        cumullated_indices <- rbind(cumullated_indices, cumu_index)
+    } else {
+        cumullated_indices <- cumu_index
+    }
+
+}  # end of year loop
+
+} else {
+
+    y <- unique(dataset$YEAR)
+    year_pheno <- flight_pheno[flight_pheno$year == y, ]
+
+    dataset_y <- dataset[dataset$YEAR == y, ]
+
+    sp_data_site <- year_day_func(dataset_y)
+    sp_data_site$trimDAYNO <- sp_data_site$DAYNO - min(sp_data_site$DAYNO) + 1
+
+    sp_data_site <- merge(sp_data_site, year_pheno[, c("DAYNO", "nm")],
+        by = c("DAYNO"), all.x = TRUE, sort = FALSE)
+
+    # compute proportion of the flight curve sampled due to missing visits
+    pro_missing_count <- data.frame(SITE = sp_data_site$SITE, WEEK = sp_data_site$WEEK,
+        NM = sp_data_site$nm, COUNT = sp_data_site$COUNT, ANCHOR = sp_data_site$ANCHOR)
+    pro_missing_count$site_week <- paste(pro_missing_count$SITE, pro_missing_count$WEEK,
+        sep = "_")
+    siteweeknocount <- aggregate(pro_missing_count$COUNT, by = list(pro_missing_count$site_week),
+        function(x) sum(!is.na(x)) == 0)
+    pro_missing_count <- pro_missing_count[pro_missing_count$site_week %in%
+        siteweeknocount$Group.1[siteweeknocount$x == TRUE], ]
+    pro_count_agg <- aggregate(pro_missing_count$NM, by = list(pro_missing_count$SITE),
+        function(x) 1 - sum(x, na.rm = T))
+    names(pro_count_agg) <- c("SITE", "PROP_PHENO_SAMPLED")
+
+    # remove samples outside the monitoring window
+    sp_data_site$COUNT[sp_data_site$nm==0] <- NA
+
+    # Compute the regional GAM index
+    if(length(unique(sp_data_site$SITE))>1){
+        glm_obj_site <- glm(COUNT ~ factor(SITE) + offset(log(nm)) - 1, data = sp_data_site,
+            family = quasipoisson(link = "log"), control = list(maxit = 100))
+    } else {
+       glm_obj_site <- glm(COUNT ~  offset(log(nm)) - 1, data = sp_data_site,
+           family = quasipoisson(link = "log"), control = list(maxit = 100))
+    }
+
+    sp_data_site[, "FITTED"] <- predict.glm(glm_obj_site, newdata = sp_data_site,
+        type = "response")
+    sp_data_site[, "COUNT_IMPUTED"] <- sp_data_site$COUNT
+    sp_data_site[is.na(sp_data_site$COUNT), "COUNT_IMPUTED"] <- sp_data_site$FITTED[is.na(sp_data_site$COUNT)]
+
+    # add fitted value for missing mid-week data
+    sp_data_site <- sp_data_site[!paste(sp_data_site$DAY_WEEK, sp_data_site$COUNT) %in%
+        c("1 NA", "2 NA", "3 NA", "5 NA", "6 NA", "7 NA"), ]
+
+    # remove all added mid-week values for weeks with real count
+    # (observation)
+    sp_data_site$site_week <- paste(sp_data_site$SITE, sp_data_site$WEEK,
+        sep = "_")
+    siteweekcount <- aggregate(sp_data_site$COUNT, by = list(sp_data_site$site_week),
+        function(x) sum(!is.na(x)) > 0)
+    sp_data_site <- sp_data_site[!(is.na(sp_data_site$COUNT) & (sp_data_site$site_week %in%
+        siteweekcount$Group.1[siteweekcount$x == TRUE])), names(sp_data_site) !=
+        "site_week"]
+
+    # Compute the regional gam index
+    print(paste("Compute index for",sp_data_site$SPECIES[1],"at year", y,"for",length(unique(sp_data_site$SITE)),"sites:",Sys.time()))
+    regional_gam_index <- trap_index(sp_data_site, data_col = "COUNT_IMPUTED",
+        time_col = "DAYNO", by_col = c("SPECIES", "SITE", "YEAR"))
+
+    cumu_index <- merge(regional_gam_index, pro_count_agg, by = c("SITE"),
+        all.x = TRUE, sort = FALSE)
+    names(cumu_index) <- c("SITE", "SPECIES", "YEAR", "regional_gam", "prop_pheno_sampled")
+
+    cumu_index <- cumu_index[order(cumu_index$SITE), ]
+
+    # bind if exist else create
+    if ("cumullated_indices" %in% ls()) {
+        cumullated_indices <- rbind(cumullated_indices, cumu_index)
+    } else {
+        cumullated_indices <- cumu_index
+    }
+
+}
+
+return(cumullated_indices)
+
+}
diff -r 000000000000 -r 5b126f770671 flight_curve.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flight_curve.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,12 @@
+#!/usr/bin/env Rscript
+#library('getopt')
+#library(devtools)
+
+args = commandArgs(trailingOnly=TRUE)
+source(args[1]) #TODO replace by library(regionalGAM) if available as official package from bioconda
+
+tryCatch({input = read.table(args[2], header=TRUE,sep="	")},finally={input = read.table(args[2], header=TRUE,sep=",")})
+dataset1 <- input[,c("SPECIES","SITE","YEAR","MONTH","DAY","COUNT")]
+pheno <- flight_curve(dataset1)
+
+write.table(pheno, file="pheno", row.names=FALSE, sep="	")
diff -r 000000000000 -r 5b126f770671 flight_curve_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/flight_curve_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,36 @@
+<tool id="flight-curve" name="Flight curve">
+<description>compute the regional expected pattern of abundance</description>
+    <requirements>
+        <requirement type="package" version="3.4">r</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/flight_curve.R' '$__tool_directory__/dennis_gam_initial_functions.R' '$input' '$output' ]]>
+    </command>
+    <inputs>
+        <param format="tabular,csv" name="input" type="data" label="Count file" help="The file must contain the SPECIES, SITE, YEAR, MONTH, DAY and COUNT columns. CSV file separator ',' or tab."/>
+    </inputs>
+    <outputs>
+        <data format="tabular" name="output" from_work_dir="pheno" />
+    </outputs>
+    <tests>
+        <test> <!--FAILED-->
+            <param name="input" value="gatekeeper_CM"/>
+            <output name="output" value="flight_curve_result.tabular"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+    
+.. class:: infomark 
+
+==========================
+Regional phenology
+==========================
+
+This tool is an implementation of the flight_curve function `RegionalGAM package:  https://github.com/RetoSchmucki/regionalGAM/
+This function computes the annual phenology on a specific region.
+
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+    </citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 glmmpql.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/glmmpql.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,18 @@
+#!/usr/bin/env Rscript
+library(nlme)
+library(MASS)
+
+args = commandArgs(trailingOnly=TRUE)
+input = read.table(args[1], header=TRUE,sep="	") #input=data.index =ab_index-output
+
+glmm.mod_fullyear <- glmmPQL(regional_gam~ as.factor(YEAR)-1,data=input,family=quasipoisson,random=~1|SITE, correlation = corAR1(form = ~ YEAR | SITE),verbose = FALSE)
+
+col.index <- as.numeric(glmm.mod_fullyear$coefficients$fixed)
+year <- unique(input$YEAR)
+
+write.table(col.index, file="output-glmmpql", row.names=FALSE, sep=" ")
+#write.table(col.index, file="output-glmmpql", row.names=FALSE)
+
+png('output-plot.png')
+plot(year,col.index,type='o', xlab="year",ylab="collated index")
+invisible(dev.off())
diff -r 000000000000 -r 5b126f770671 glmmpql_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/glmmpql_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,40 @@
+<tool id="glmmpql" name="Expected temporal trend" version="0.1">
+    <description>of species abundance</description>
+    <requirements>
+        <requirement type="package" version="3.1_131">r-nlme</requirement>
+        <requirement type="package" version="7.3_48">r-mass</requirement>
+        <requirement type="package" version="0.9.10">xorg-libxrender</requirement>
+        <requirement type="package" version="1.2.2">xorg-libsm</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/glmmpql.R' '$input1' '$output' '$output2' ]]>
+    </command>
+    <inputs>
+        <param format="tabular" name="input1" type="data" label="Tabular file generated by the ab_index tool"/>
+    </inputs>
+    <outputs>
+        <data format="png" name="output" from_work_dir="output-plot.png" />
+        <data format="tabular" name="output2" from_work_dir="output-glmmpql" />
+    </outputs>
+    <tests>
+        <test>
+            <param name="input1" value="data_index.tabular"/>
+            <param name="output" value="glmmpql-plot.png"/>
+            <output name="output2" value="glmmpql.tabular"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+    
+.. class:: infomark 
+
+==========================
+Temporal trends of species abundance
+==========================
+
+This tool is an implementation of the glmmpql function from `RegionalGAM package:  https://github.com/RetoSchmucki/regionalGAM/
+This function computes a collated index for each year and estimates the temporal trend.
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+    </citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 gls-adjusted.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gls-adjusted.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,16 @@
+#!/usr/bin/env Rscript
+library(nlme)
+library(MASS)
+
+args = commandArgs(trailingOnly=TRUE)
+input1 = read.table(args[1], header=TRUE) #input1=col.index =glmmpql-output
+input2 = read.table(args[2], header=TRUE,sep="	") #input2=data.index =abundance_index-output
+
+input1<-as.matrix(input1)
+
+year <- unique(input2$YEAR)
+mod <- gls(input1~year, correlation = corARMA(p=2))
+summary<-summary(mod)
+
+save(mod, file = "mod_adjusted.rda")
+capture.output(summary, file="mod_adjusted-summary.txt")
diff -r 000000000000 -r 5b126f770671 gls-adjusted_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gls-adjusted_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,41 @@
+<tool id="gls-adjusted" name="Linear regression ajusted" version="0.1">
+    <description>for autocorrelation in the residuals</description>
+    <requirements>
+        <requirement type="package" version="3.1_131">r-nlme</requirement>
+        <requirement type="package" version="7.3_48">r-mass</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/gls-adjusted.R' '$input1' '$input2' '$output' '$output2' ]]>
+    </command>
+    <inputs>
+        <param format="tabular" name="input1" type="data" label="File generated by the glmmpql/Expected temporal trend tool"/>
+        <param format="tabular" name="input2" type="data" label="File generated by the ab_index tool"/>
+    </inputs>
+    <outputs>
+        <data format="txt" name="output" from_work_dir="mod_adjusted-summary.txt" />
+        <data format="rdata" name="output2" from_work_dir="mod_adjusted.rda" />
+    </outputs>
+    <tests>
+        <test> <!--FAILED-->
+            <param name="input1" value="glmmpql.tabular"/>
+            <param name="input2" value="data_index.tabular"/>
+            <output name="output" value="gls-adju_sum.txt"/>
+            <output name="output2" value="gls-adju_rda.rda"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+    
+.. class:: infomark 
+
+==========================
+STOC eps for specialists groups
+==========================
+
+This tool is an implementation of the gls_adjusted function from `RegionalGAM package: https://github.com/RetoSchmucki/regionalGAM/
+This function adjust the model to account for autocorrelation in the residuals.
+
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+    </citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 gls.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gls.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,16 @@
+#!/usr/bin/env Rscript
+library(nlme)
+library(MASS)
+
+args = commandArgs(trailingOnly=TRUE)
+input1 = read.table(args[1], header=TRUE) #input1=col.index =glmmpql-output
+input2 = read.table(args[2], header=TRUE,sep="	") #input2=data.index =abundance_index-output
+
+input1<-as.matrix(input1)
+
+year <- unique(input2$YEAR)
+mod <- gls(input1~year)
+summary<-summary(mod)
+
+save(mod, file = "mod.rda")
+capture.output(summary, file="mod-summary.txt")
diff -r 000000000000 -r 5b126f770671 gls_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/gls_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,41 @@
+<tool id="gls" name="Model temporal trend" version="0.1">
+    <description>with a simple linear regression</description>
+    <requirements>
+        <requirement type="package" version="3.1_131">r-nlme</requirement>
+        <requirement type="package" version="7.3_48">r-mass</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/gls.R' '$input1' '$input2' '$output' '$output2' ]]>
+    </command>
+    <inputs>
+        <param name="input1" type="data" format="tabular" label="File generated by the glmmpql/Expected temporal trend tool"/>
+        <param name="input2" type="data" format="tabular" label="File generated by the outil ab_index tool"/>
+    </inputs>
+    <outputs>
+        <data format="txt" name="output" from_work_dir="mod-summary.txt" />
+        <data format="rdata" name="output2" from_work_dir="mod.rda" />
+    </outputs>
+    <tests>
+        <test> <!--FAILED for binary-->
+            <param name="input1" value="glmmpql.tabular"/>
+            <param name="input2" value="data_index.tabular"/>
+            <output name="output" value="gls_sum.txt"/>
+            <output name="output2" value="gls_rda.rda"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+
+.. class:: infomark 
+
+==========================
+Model temporal trend
+==========================
+
+This tool is an implementation of the gls function from `RegionalGAM package:  https://github.com/RetoSchmucki/regionalGAM/
+This function model temporal trend with a simple linear regression.
+
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+</citations>
+</tool>
diff -r 000000000000 -r 5b126f770671 plot_trend.R
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/plot_trend.R	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,17 @@
+#!/usr/bin/env Rscript
+library(nlme)
+library(MASS)
+
+args = commandArgs(trailingOnly=TRUE)
+input = read.table(args[1], header=TRUE,sep="	") #input=data.index =ab_index-output
+load(args[2])
+
+glmm.mod_fullyear <- glmmPQL(regional_gam~ as.factor(YEAR)-1,data=input,family=quasipoisson,random=~1|SITE, correlation = corAR1(form = ~ YEAR | SITE),verbose = FALSE)
+
+col.index <- as.numeric(glmm.mod_fullyear$coefficients$fixed)
+year <- unique(input$YEAR)
+
+png('output-plot-trend.png')
+plot(year,col.index, type='o',xlab="year",ylab="collated index")
+abline(mod,lty=2,col="red")
+dev.off()
diff -r 000000000000 -r 5b126f770671 plot_trend_en.xml
--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/plot_trend_en.xml	Fri Aug 03 08:28:22 2018 -0400
@@ -0,0 +1,40 @@
+<tool id="plot_trend" name="Plot abundance" version="0.1">
+    <description>with trend line</description>
+    <requirements>
+        <requirement type="package" version="3.1_131">r-nlme</requirement>
+        <requirement type="package" version="7.3_48">r-mass</requirement>
+        <requirement type="package" version="0.9.10">xorg-libxrender</requirement>
+        <requirement type="package" version="1.1.2">xorg-libsm</requirement>
+    </requirements>
+    <command detect_errors="exit_code"><![CDATA[
+        Rscript '$__tool_directory__/plot_trend.R' '$input1' '$input2' '$output' ]]>
+    </command>
+    <inputs>
+        <param format="tabular" name="input1" type="data" label="File generated by the ab_index tool"/>
+        <param format="rdata" name="input2" type="data" label="gls model" help="Rdata File"/>
+    </inputs>
+    <outputs>
+        <data format="png" name="output" from_work_dir="output-plot-trend.png" />
+    </outputs>
+    <tests>
+        <test>
+            <param name="input1" value="data_index.tabular"/>
+            <param name="input2" value="gls-adju_rda.rda"/>
+            <output name="output" value="trend.png"/>
+        </test>
+    </tests>
+    <help><![CDATA[
+
+.. class:: infomark 
+
+==========================
+Plot abundance trend 
+==========================
+
+This tool is an implementation of the plot_trend function from `RegionalGAM package: https://github.com/RetoSchmucki/regionalGAM/
+This function plot abundance and add a trend line.
+    ]]></help>
+    <citations>
+        <citation type="doi">10.1111/1365-2664.12561</citation>  
+    </citations>
+</tool>