R/false_ring_proportion.R

In iadf: Analysis of Intra Annual Density Fluctuations

globalVariables(".")

#frp----------------------------------------------------------------------------
#' @title false ring proportion
#' @description Calculate the false ring proportion of a set of
#'   binary false ring assignments.
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked binary with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @return a data frame
#' @seealso \code{\link{afrp}}
#' @export
frp <- function(iadf) {

  if (!is.data.frame(iadf) && !is.matrix(iadf)) {
    stop("x has to be a data.frame or matrix")
  }

  if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
    stop("the data frame must only contain 0, 1, and NA values")
  }

  n <- ncol(iadf) - rowSums(is.na(iadf), na.rm = TRUE)

  if (any(n == 0)) {
    warning("there are no trees in year: ",
            paste(names(n[n == 0]), collapse = ", "))
  }

  iadfp <- rowMeans(iadf, na.rm = TRUE)

  out <- data.frame(year = as.numeric(names(iadfp)), frp = iadfp)
  return(out)
}

# afrp  ------------------------------------------------------------------------
#' @title adjusted false ring proportion
#' @description Calculate the adjusted false ring proportion, as
#'   suggested by Osborn et. al. (1997), of a set of binary false ring
#'   assignments.
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @references Osborn TJ, Briffa KR and Jones PD (1997) Adjusting variance for
#'   sample-size in tree-ring chronologies and other regional mean time-series.
#'   Dendrochronologia 15, 89-99.
#' @return a data frame
#' @seealso \code{\link{frp}}
#' @export
afrp <- function(iadf) {

  if (!is.data.frame(iadf) && !is.matrix(iadf)) {
    stop("x has to be a data.frame or matrix")
  }

  if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
    stop("the data frame must only contain 0, 1, and NA values")
  }

  n <- ncol(iadf) - rowSums(is.na(iadf), na.rm = TRUE)

  if (any(n == 0)){
    warning("there are no trees in year: ",
            paste(names(n[n == 0]), collapse = ", "))
  }

  iadfp <- rowMeans(iadf, na.rm = TRUE)

  aiadfp <- iadfp * sqrt(n)

  out <- data.frame(year = as.numeric(names(aiadfp)), afrp = aiadfp)
  return(out)
}

# novak_freq -------------------------------------------------------------------
#' @title iadf frequency per cambial age
#' @description Calculate the frequency per cambial age as suggested
#'  by Novak et al. (2013).
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked binary with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @param po a data frame with pith offsets with series names in the first and
#'   pith offset as number of rings in the second column
#' @references
#'   Novak, Klemen and Sánchez, Miguel Angel Saz and Čufar, Katarina and
#'   Raventós, Josep and de Luis, Martin. Age, climate and intra-annual density
#'   fluctuations in in Spain, IAWA Journal, 34, 459-474 (2013),
#'   doi:10.1163/22941932-00000037
#' @return a data frame
#' @export
#' @seealso \code{\link{novak_weibull}},
#'   \code{\link{novak_index}}
#' @examples
#' data('example_iadf')
#' model <- novak_weibull(novak_freq(example_iadf), 15)
#' novak_index(example_iadf, model)

novak_freq <- function(iadf, po = NULL){

  if (!is.data.frame(iadf) && !is.matrix(iadf)) {
    stop("iadf has to be a data.frame or matrix")
  }

  if (!is.null(po)) {
    if (!is.data.frame(po) || ncol(po) != 2) {
      stop("po has to be a data.frame with ncol == 2")
    }
  }

  if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
    stop("iadf must only contain 0, 1, and NA values")
  }

  aligned <- to_cambial_age(iadf, po)
  freq <- rowMeans(aligned, na.rm = TRUE)
  n <- rowSums(!is.na(aligned))
  out <- data.frame(cambial.age = seq_along(freq), freq = freq,
                    sample.depth = n)
  class(out) <- c("data.frame", "novak.freq")
  return(out)
}


# novak weibull find start values ----------------------------------------------
#' @title novak_weibull_find_start
#' @description Find good start values manually in case
#' \code{\link{novak_weibull}}
#'   returns an error caused by insufficient default starting values.
#' @param novak_freq_object A novak_freq_object as obtained from
#'   \code{\link[iadf]{novak_freq}}
#' @param min.n minimum number of samples within each cambial age to be included
#' in model estimation
#' @param max_a maximum value of manipulate slider for parameter a
#' @param max_b maximum value of manipulate slider for parameter b
#' @param max_c maximum value of manipulate slider for parameter c
#'
#' @return a list which can be used as input argument 'start' in
#' \code{\link{novak_weibull}}
#' @export

novak_weibull_find_start <- function(novak_freq_object, min.n = 15, max_a = 10,
                                     max_b = 3, max_c = 30) {
  start <- list(a = 4, b = 0.33, c = 15.5)
  observe <- FALSE

  a0 <- c(); b0 <- c(); c0 <- c(); x <- c(); return_value <- c()

  manipulate::manipulate({
    plot(novak_freq_object[["freq"]] ~ novak_freq_object[["cambial.age"]],
         pch = 16, col = ifelse(novak_freq_object[["sample.depth"]] >= min.n,
                                "black", "gray80"),
         xlab = "cambial age", ylab = "frequency")
    a <- a0;  b <- b0;  c <- c0
    curve(a * b * c * x^2 * exp(-a * x^b), add = TRUE)
    start <<- list(a = a, b = b, c = c)
    if (return_value) {
      observe <<- TRUE
    }
  },
  a0 = manipulate::slider(0, max_a, step = 0.0001, initial = start$a),
  b0 = manipulate::slider(0, max_b, step = 0.0001, initial = start$b),
  c0 = manipulate::slider(0, max_c, step = 0.0001, initial = start$c),
  return_value = manipulate::button("Return value")
  )

  repeat{
    if (observe) {
      return(start)
    }
  }
}

# novak_weibull ----------------------------------------------------------------
#' @title novak_weibull
#' @description Fit a Weibull function for the calculation of
#'   age corrected IADF frequencies according to Novak et al. (2013).
#' @param novak_freq_object A novak_freq_object as obtained from
#'   \code{\link[iadf]{novak_freq}}
#' @param min.n minimum number of samples within each cambial age to be included
#' in model estimation
#' @param start set custom start values - default to
#' \code{list(a = 4, b = 0.33, c = 15.5)}
#' @param max.iter maximum iterations for internally used
#' \code{\link[stats]{nls}}
#' @param make.plot logical
#' @param ... additional plotting arguments
#' @references
#'   Novak, Klemen and Sánchez, Miguel Angel Saz and Čufar, Katarina and
#'   Raventós, Josep and de Luis, Martin. Age, climate and intra-annual density
#'   fluctuations in in Spain, IAWA Journal, 34, 459-474 (2013),
#'   doi:10.1163/22941932-00000037
#' @return a model object of class "nls"
#' @export
#' @seealso \code{\link{novak_freq}},
#'   \code{\link{novak_index}}
#' @examples
#' data('example_iadf')
#' model <- novak_weibull(novak_freq(example_iadf), 15)
#' novak_index(example_iadf, model)
novak_weibull <- function(novak_freq_object, min.n = 15, start = NULL,
                          max.iter = 500, make.plot = TRUE, ...) {

  if (!any(class(novak_freq_object) == "novak.freq")) {
    stop("input must be derived from novak_freq()")
  }

  tmp <- novak_freq_object[novak_freq_object[["sample.depth"]] >= min.n &
                             rowSums(is.na(novak_freq_object)) < 3, ]

  cambial_age <- tmp[["cambial.age"]]
  freq <- tmp[["freq"]]

  if (is.null(start)) {
    start <- list(a = 4, b = 0.33, c = 15.5)
  }

  try({
    wbu <- stats::nls(freq ~ a * b * c * cambial_age^2 * exp(-a * cambial_age^b),
                      start = start, control = list(maxiter = max.iter))

    if (make.plot) {
      ndata <- data.frame(cambial_age = seq(0, max(cambial_age, na.rm = TRUE),
                                            by = 0.1))
      wpred <- predict(wbu, newdata = ndata)
      plot(novak_freq_object[["freq"]] ~ novak_freq_object[["cambial.age"]],
           pch = 16, col = ifelse(novak_freq_object[["sample.depth"]] >= min.n,
                                  "black", "gray80"),
           xlab = "cambial age", ylab = "frequency", ...)
      lines(ndata$cambial_age, wpred, col = "red")
    }

    return(wbu)
  })
}


# novak_index ------------------------------------------------------------------
#' @title novak_index
#' @description  Calculation of age corrected IADF frequencies according
#'   to Novak et al. (2013).
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked binary with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @param model a model, output of either  \code{\link[iadf]{novak_weibull}}
#' @return a data frame
#' @param po an optional data frame of pith offsets with series names in the
#'   first and pith offsets in the second column
#' @param method method for the RCS detrending, 'quotient' or 'difference'
#' @references
#'   Novak, Klemen and Sánchez, Miguel Angel Saz and Čufar, Katarina and
#'   Raventós, Josep and de Luis, Martin.  Age, climate and intra-annual density
#'   fluctuations in in Spain, IAWA Journal, 34, 459-474 (2013),
#'   doi:10.1163/22941932-00000037
#' @export
#' @seealso \code{\link{novak_freq}}, \code{\link{novak_weibull}}
#' @examples
#' data('example_iadf')
#' model <- novak_weibull(novak_freq(example_iadf), 15)
#' novak_index(example_iadf, model)
novak_index <- function(iadf, model, po = NULL, method = "difference"){

  if (is.null(po)) {
    po <- data.frame(series = names(iadf), po = 1)
  }

  if (!method %in% c("quotient", "difference")) {
    stop("please choose either 'quotient' or 'difference' as method")
  }

  if (!is.data.frame(iadf)) {
    stop("iadf must be of class data.frame")
  }

  if (!(is.data.frame(po))) {
    stop("po must be of class data.frame")
  }

  if (!setequal(po[, 1], names(iadf))) {
    stop("series names in po are not the same as provided in iadf")
  }

  longest <- max(series_length(iadf) + po[, 2])
  rc <- stats::predict(model, newdata = list(cambial_age = seq_len(longest)))
  detrended <- detrend_given_rc(iadf, rc, po, method = method)
  tmp <- rowMeans(detrended, na.rm = TRUE)
  out <- data.frame(year = as.numeric(names(tmp)), index = as.vector(tmp))
  return(out)
}


# campelo_freq -----------------------------------------------------------------
#' @title iadf frequency per ring width class
#' @description Calculate the frequency per ring width class as suggested
#'  by Campelo (2015).
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked binary with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @param rwl data frame containing ring widths with years in rows and
#'   series in columns
#' @param n number of ring width classes
#' @return a data frame
#' @references
#'   Campelo, F., Vieira, J., Battipaglia, G. et al. Which matters most for the
#'   formation of intra-annual density fluctuations in Pinus pinaster: age or
#'   size? Trees (2015) 29: 237. doi:10.1007/s00468-014-1108-9
#' @export
#' @seealso  \code{\link{campelo_chapman}}, \code{\link{campelo_index}}
#' @importFrom rlang .data
#' @importFrom dplyr %>%
#' @examples
#' data('example_iadf')
#' data('example_rwl')
#' model <- campelo_chapman(campelo_freq(example_iadf, example_rwl))
#' campelo_index(example_iadf, example_rwl, model)
campelo_freq <- function(iadf, rwl, n = 20){
  if (utils::packageVersion("dplyr") > "0.5.0") {
    if (!is.data.frame(iadf)) {
      stop("iadf has to be a data.frame")
    }

    if (!is.data.frame(rwl)) {
      stop("rwl has to be a data.frame")
    }

    if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
      stop("iadf must only contain 0, 1, and NA values")
    }

    nclass <- rlang::quo(n)

    iadf_tidy <- iadf %>%
      tibble::rownames_to_column("year") %>%
      tibble::as_tibble() %>%
      dplyr::mutate(year = as.integer(.data$year)) %>%
      tidyr::gather("series", "iadf", -.data$year) %>%
      dplyr::filter(!is.na(.data$iadf))

    rwl_tidy <- tidy_rwl(rwl) %>%
      dplyr::mutate(class = cut(.data$rwl, !!nclass))

    tmp <- dplyr::left_join(iadf_tidy, rwl_tidy, by = c("year", "series"))

    out <- tmp %>%
      dplyr::group_by(.data$class) %>%
      dplyr::summarise(freq = mean(.data$iadf, na.rm = TRUE),
                       class.mean.rwl = mean(.data$rwl, na.rm = TRUE),
                       sample.depth = length(.data$iadf)) %>%
      as.data.frame

    class(out) <- c("data.frame", "campelo.freq")

    return(out)

  } else {# keep code for dplyr < 0.6 for backwards compatibility

    if (!is.data.frame(iadf)) {
      stop("iadf has to be a data.frame")
    }

    if (!is.data.frame(rwl)) {
      stop("rwl has to be a data.frame")
    }

    if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
      stop("iadf must only contain 0, 1, and NA values")
    }

    iadf_tidy <- iadf %>%
      tibble::rownames_to_column("year") %>%
      tibble::as_tibble() %>%
      dplyr::mutate_(.dots = setNames(
        list(lazyeval::interp(~as.numeric(var),
                              var = as.name("year"))), "year")) %>%
      tidyr::gather_("series", "iadf",
                     dplyr::select_vars_(names(.), names(.),
                                         exclude = "year")) %>%
      dplyr::filter_(lazyeval::interp(~(!is.na(nam)), nam = as.name("iadf")))

    rwl_tidy <- tidy_rwl(rwl) %>%
      dplyr::mutate_(.dots = setNames(
        list(lazyeval::interp(~cut(r, n),
                              r = as.name("rwl"))), "class"))

    tmp <- dplyr::left_join(iadf_tidy, rwl_tidy, by = c("year", "series"))

    freq <- tapply(tmp[["iadf"]], tmp[["class"]],
                   function(y) mean(y, na.rm = TRUE))
    class_mean <- tapply(tmp[["rwl"]], tmp[["class"]],
                         function(y) mean(y, na.rm = TRUE))
    sample_depth <- tapply(tmp[["iadf"]], tmp[["class"]],
                           function(y) length(na.omit(y)))
    out <- data.frame(class = factor(names(freq), levels = names(freq)),
                      freq = as.vector(freq),
                      class.mean.rwl = as.vector(class_mean),
                      sample.depth = as.vector(sample_depth))
    class(out) <- c("data.frame", "campelo.freq")
    return(out)
  }
}

# campelo chapman find start values --------------------------------------------
#' @title campelo_chapman_find_start
#' @description Find good start values manually in case
#' \code{\link{campelo_chapman}}
#'   returns an error caused by insufficient default starting values.
#' @param campelo_freq_object a campelo frequency object,
#' output of \code{\link[iadf]{campelo_freq}}
#' @param min.n minimum number of samples within each group to be included in
#'   model estimation
#' @param max_a maximum value of manipulate slider for parameter a
#' @param max_b maximum value of manipulate slider for parameter b
#' @param max_c maximum value of manipulate slider for parameter c
#'
#' @return a list which can be used as input argument 'start' in
#' \code{\link{campelo_chapman}}
#' @export
campelo_chapman_find_start <- function(campelo_freq_object, min.n = 15,
                                       max_a = 3, max_b = 1, max_c = 17) {

  tmp <- campelo_freq_object[campelo_freq_object[["sample.depth"]] >= min.n &
                               rowSums(is.na(campelo_freq_object)) < 3, ]

  ring_width <- tmp[["class.mean.rwl"]]
  freq <- tmp[["freq"]]

  start <- list(a = max_a / 2, b = max_b / 2, c = max_c / 2)
  observe <- FALSE

  a0 <- c(); b0 <- c(); c0 <- c(); x <- c(); return_value <- c()

  manipulate::manipulate({
    plot(campelo_freq_object[["freq"]] ~ campelo_freq_object[["class.mean.rwl"]],
         pch = 16, col = ifelse(campelo_freq_object[["sample.depth"]] >= min.n,
                                "black", "gray80"),
         xlab = "mean class ring width", ylab = "frequency")
    a <- a0;  b <- b0;  c <- c0
    curve(a * (1 - exp(-b * x))^c, add = TRUE)
    start <<- list(a = a, b = b, c = c)
    if (return_value) {
      observe <<- TRUE
    }
  },
  a0 = manipulate::slider(0, max_a, step = 0.0001, initial = start$a),
  b0 = manipulate::slider(0, max_b, step = 0.0001, initial = start$b),
  c0 = manipulate::slider(0, max_c, step = 0.0001, initial = start$c),
  return_value = manipulate::button("Return value")
  )

  repeat{
    if (observe){
      return(start)
    }
  }
}

# campelo_chapman --------------------------------------------------------------
#' @title campelo_chapman
#' @description Chapman model fitting to size classes for the calculation of
#' size corrected IADF frequencies according to Campelo et al. (2015).
#' @param campelo_freq_object a campelo frequency object,
#' output of \code{\link[iadf]{campelo_freq}}
#' @param min.n minimum number of samples within each group to be included in
#'   model estimation
#' @param start set custom start values - default to
#' \code{list(a = 0.8, b = 0.03, c = 12.5)}
#' @param max.iter maximum iterations for internally used
#' \code{\link[stats]{nls}}
#' @param make.plot logical
#' @param ... additional plotting arguments
#' @return a model object of class "nls"
#' @export
#' @seealso \code{\link{campelo_freq}}, \code{\link{campelo_index}}
#' @references
#'   Campelo, F., Vieira, J., Battipaglia, G. et al. Which matters most for the
#'   formation of intra-annual density fluctuations in Pinus pinaster: age or
#'   size? Trees (2015) 29: 237. doi:10.1007/s00468-014-1108-9
#' @examples
#' data('example_iadf')
#' data('example_rwl')
#' model <- campelo_chapman(campelo_freq(example_iadf, example_rwl))
#' campelo_index(example_iadf, example_rwl, model)
campelo_chapman <- function(campelo_freq_object, min.n = 15, start = NULL,
                            make.plot = TRUE, max.iter = 500, ...){

  if (!any(class(campelo_freq_object) == "campelo.freq")) {
    stop("input must be derived from campelo_freq()")
  }

  tmp <- campelo_freq_object[campelo_freq_object[["sample.depth"]] >= min.n &
                               rowSums(is.na(campelo_freq_object)) < 3, ]

  ring_width <- tmp[["class.mean.rwl"]]
  freq <- tmp[["freq"]]

  if (is.null(start)) {
    start <- list(a = 0.8, b = 0.03, c = 12.5)
  }

  try({
    chapm <- stats::nls(freq ~ a * (1 - exp(-b * ring_width))^c, start = start,
                        control = list(maxiter = max.iter))

    if (make.plot){
      ndata <- data.frame(ring_width = seq(0, max(ring_width, na.rm = TRUE),
                                           by = 0.1))
      cpred <- predict(chapm, newdata = ndata)
      plot(campelo_freq_object[["freq"]] ~ campelo_freq_object[["class.mean.rwl"]],
           pch = 16, col = ifelse(campelo_freq_object[["sample.depth"]] >= min.n,
                                  "black", "gray80"),
           xlab = "mean class ring width", ylab = "frequency", ...)
      lines(ndata$ring_width, cpred, col = "red")
    }

    return(chapm)
  })
}


# campelo_index ----------------------------------------------------------------
#' @title campelo_index
#' @description Calculation of size corrected IADF frequencies according
#'   to Campelo et al. (2015)
#' @param iadf A data frame with numeric columns representing individual series
#'   and years as rownames where years with IADF are marked binary with 1,
#'   those without with 0, years not covered by the series are set to NA.
#' @param rwl a rwl/data.frame object
#' @param model a chapman model, output of \code{\link[iadf]{campelo_chapman}}
#' @return a data frame
#' @export
#' @seealso \code{\link{campelo_freq}}, \code{\link{campelo_chapman}}
#' @importFrom rlang .data
#' @importFrom dplyr %>%
#' @references
#'   Campelo, F., Vieira, J., Battipaglia, G. et al. Which matters most for the
#'   formation of intra-annual density fluctuations in Pinus pinaster: age or
#'   size? Trees (2015) 29: 237. doi:10.1007/s00468-014-1108-9
#' @examples
#' data('example_iadf')
#' data('example_rwl')
#' model <- campelo_chapman(campelo_freq(example_iadf, example_rwl))
#' campelo_index(example_iadf, example_rwl, model)
campelo_index <- function(iadf, rwl, model){

  if (utils::packageVersion("dplyr") > "0.5.0") {

    if (!is.data.frame(iadf)) {
      stop("iadf has to be a data.frame")
    }

    if (!is.data.frame(rwl)) {
      stop("rwl has to be a data.frame")
    }

    if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
      stop("iadf must only contain 0, 1, and NA values")
    }

    iadf_tidy <- iadf %>%
      tibble::rownames_to_column("year") %>%
      tibble::as_tibble() %>%
      dplyr::mutate(year = as.numeric(.data$year)) %>%
      tidyr::gather("series", "iadf", -.data$year) %>%
      dplyr::filter(!is.na(.data$iadf))

    rwl_tidy <- tidy_rwl(rwl)

    tmp <- dplyr::left_join(iadf_tidy, rwl_tidy, by = c("year", "series")) %>%
      dplyr::mutate(prediction =
                      stats::predict(model,
                                     newdata = list(ring_width = .data$rwl))) %>%
      dplyr::mutate(index = .data$iadf - .data$prediction) %>%
      dplyr::select("year", "series", "index") %>%
      tidyr::spread("series", "index")

    mean_index <- cbind(tmp["year"], rowMeans(dplyr::select_(tmp, "-year"),
                                              na.rm = TRUE))
    names(mean_index) <- c("year", "index")

    spline <- dplR::ffcsaps(mean_index[["index"]])
    detrend_index <- dplyr::mutate(mean_index,
                                   detrended_index = .data$index - spline)

    return(detrend_index)

  } else {#keep code for dplyr < 0.6 for backwards compatibility

    if (!is.data.frame(iadf)) {
      stop("iadf has to be a data.frame")
    }

    if (!is.data.frame(rwl)) {
      stop("rwl has to be a data.frame")
    }

    if (!all(as.matrix(iadf) %in% c(0, 1, NA))) {
      stop("iadf must only contain 0, 1, and NA values")
    }

    iadf_tidy <- iadf %>%
      tibble::rownames_to_column("year") %>%
      tibble::as_tibble() %>%
      dplyr::mutate_(.dots = setNames(
        list(lazyeval::interp(~as.numeric(var),
                              var = as.name("year"))), "year")) %>%
      tidyr::gather_("series", "iadf",
                     dplyr::select_vars_(names(.), names(.),
                                         exclude = "year")) %>%
      dplyr::filter_(lazyeval::interp(~(!is.na(nam)), nam = as.name("iadf")))

    rwl_tidy <- tidy_rwl(rwl)

    tmp <- dplyr::left_join(iadf_tidy, rwl_tidy, by = c("year", "series")) %>%
      dplyr::mutate_(.dots = setNames(list(
        lazyeval::interp(~predict(md, newdata = list(ring_width = .[["rwl"]])),
                         md = as.name("model"))), "prediction")) %>%
      dplyr::mutate_(.dots = setNames(list(
        lazyeval::interp(~ia - pr, ia = as.name("iadf"),
                         pr = as.name("prediction"))), "index")) %>%
      dplyr::select_("year", "series", "index") %>%
      tidyr::spread_("series", "index")

    mean_index <- cbind(tmp["year"], rowMeans(dplyr::select_(tmp, "-year"),
                                              na.rm = TRUE))
    names(mean_index) <- c("year", "index")

    spline <- dplR::ffcsaps(mean_index[["index"]])
    detrend_index <- dplyr::mutate_(mean_index, .dots = setNames(list(
      lazyeval::interp(~ i - s, i = as.name("index"), s = as.name("spline"))),
      "detrended_index"))

    return(detrend_index)

  }
}