R/detrend.series.R

In AndyBunn/dplR: Dendrochronology Program Library in R

`detrend.series` <-
  function(y, y.name = "", make.plot = TRUE,
           method = c("Spline", "ModNegExp", "Mean", "Ar", "Friedman",
                      "ModHugershoff", "AgeDepSpline"),
           nyrs = NULL, f = 0.5, pos.slope = FALSE,
           constrain.nls = c("never", "when.fail", "always"),
           verbose = FALSE, return.info = FALSE,
           wt, span = "cv", bass = 0, difference = FALSE)
  {
    check.flags(make.plot, pos.slope, verbose, return.info)
    dirtyDog <- FALSE # this will be used to warn the user in any fits are <0
    
    if (length(y.name) == 0) {
      y.name2 <- ""
    } else {
      y.name2 <- as.character(y.name)[1]
      stopifnot(Encoding(y.name2) != "bytes")
    }
    known.methods <- c("Spline", "ModNegExp", "Mean", "Ar", "Friedman",
                       "ModHugershoff", "AgeDepSpline")
    constrain2 <- match.arg(constrain.nls)
    method2 <- match.arg(arg = method,
                         choices = known.methods,
                         several.ok = TRUE)
    wt.missing <- missing(wt)
    wt.description <- NULL
    if (verbose) {
      widthOpt <- getOption("width")
      indentSize <- 1
      indent <- function(x) {
        paste0(paste0(rep.int(" ", indentSize), collapse = ""), x)
      }
      sepLine <-
        indent(paste0(rep.int("~", max(1, widthOpt - 2 * indentSize)),
                      collapse = ""))
      cat(sepLine,
          gettext("Verbose output: ", domain="R-dplR"), y.name2,
          sep = "\n")
      wt.description <- if (wt.missing) "default" else deparse(wt)
      opts <- c("make.plot" = make.plot,
                "method(s)" = deparse(method2),
                "nyrs" = if (is.null(nyrs)) "NULL" else nyrs,
                "f" = f,
                "pos.slope" = pos.slope,
                "constrain.nls" = constrain2,
                "verbose" = verbose,
                "return.info" = return.info,
                "wt" = wt.description,
                "span" = span,
                "bass" = bass,
                "difference" = difference)
      optNames <- names(opts)
      optChar <- c(gettext("Options", domain="R-dplR"),
                   paste(str_pad(optNames,
                                 width = max(nchar(optNames)),
                                 side = "right"),
                         opts, sep = "  "))
      cat(sepLine, indent(optChar), sep = "\n")
    }
    
    ## Remove NA from the data (they will be reinserted later)
    good.y <- which(!is.na(y))
    if(length(good.y) == 0) {
      stop("all values are 'NA'")
    } else if(any(diff(good.y) != 1)) {
      stop("'NA's are not allowed in the middle of the series")
    }
    y2 <- y[good.y]
    nY2 <- length(y2)
    
    ## Recode any zero values to 0.001
    if (verbose || return.info) {
      years <- names(y2)
      if (is.null(years)) {
        years <- good.y
      }
      zeroFun <- function(x) list(zero.years = years[is.finite(x) & x == 0])
      nFun <- function(x) list(n.zeros = length(x[[1]]))
      zero.years.data <- zeroFun(y2)
      n.zeros.data <- nFun(zero.years.data)
      dataStats <- c(n.zeros.data, zero.years.data)
      if (verbose) {
        cat("", sepLine, sep = "\n")
        if (n.zeros.data[[1]] > 0){
          if (is.character(years)) {
            cat(indent(gettext("Zero years in input series:\n",
                               domain="R-dplR")))
          } else {
            cat(indent(gettext("Zero indices in input series:\n",
                               domain="R-dplR")))
          }
          cat(indent(paste(zero.years.data[[1]], collapse = " ")),
              "\n", sep = "")
        } else {
          cat(indent(gettext("No zeros in input series.\n",
                             domain="R-dplR")))
        }
      }
    }
    y2[y2 == 0] <- 0.001
    
    resids <- list()
    curves <- list()
    modelStats <- list()
    ################################################################################    
    ################################################################################    
    # Ok. Let's start the methods
    
    ################################################################################    
    if("ModNegExp" %in% method2){
      ## Nec or lm
      nec.func <- function(Y, constrain) {
        nY <- length(Y)
        a <- mean(Y[seq_len(max(1, floor(nY * 0.1)))])
        b <- -0.01
        k <- mean(Y[floor(nY * 0.9):nY])
        nlsForm <- Y ~ I(a * exp(b * seq_along(Y)) + k)
        nlsStart <- list(a=a, b=b, k=k)
        checked <- FALSE
        constrained <- FALSE
        ## Note: nls() may signal an error
        if (constrain == "never") {
          nec <- nls(formula = nlsForm, start = nlsStart)
        } else if (constrain == "always") {
          nec <- nls(formula = nlsForm, start = nlsStart,
                     lower = c(a=0, b=-Inf, k=0),
                     upper = c(a=Inf, b=0, k=Inf),
                     algorithm = "port")
          constrained <- TRUE
        } else {
          nec <- nls(formula = nlsForm, start = nlsStart)
          coefs <- coef(nec)
          if (coefs[1] <= 0 || coefs[2] >= 0) {
            stop()
          }
          fits <- predict(nec)
          if (fits[nY] > 0) {
            checked <- TRUE
          } else {
            nec <- nls(formula = nlsForm, start = nlsStart,
                       lower = c(a=0, b=-Inf, k=0),
                       upper = c(a=Inf, b=0, k=Inf),
                       algorithm = "port")
            constrained <- TRUE
          }
        }
        if (!checked) {
          coefs <- coef(nec)
          if (coefs[1] <= 0 || coefs[2] >= 0) {
            stop()
          }
          fits <- predict(nec)
          if (fits[nY] <= 0) {
            ## This error is a special case that needs to be
            ## detected (if only for giving a warning).  Any
            ## smarter way to implement this?
            return(NULL)
          }
        }
        tmpFormula <- nlsForm
        formEnv <- new.env(parent = environment(detrend.series))
        formEnv[["Y"]] <- Y
        formEnv[["a"]] <- coefs["a"]
        formEnv[["b"]] <- coefs["b"]
        formEnv[["k"]] <- coefs["k"]
        environment(tmpFormula) <- formEnv
        structure(fits, constrained = constrained,
                  formula = tmpFormula, summary = summary(nec))
      }
      ModNegExp <- try(nec.func(y2, constrain2), silent=TRUE)
      mneNotPositive <- is.null(ModNegExp)
      
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(gettext("Detrend by ModNegExp.\n", domain = "R-dplR")))
        cat(indent(gettext("Trying to fit nls model...\n",
                           domain = "R-dplR")))
      }
      if (mneNotPositive || inherits(ModNegExp,"try-error")) {
        if (verbose) {
          cat(indent(gettext("nls failed... fitting linear model...",
                             domain = "R-dplR")))
        }
        ## Straight line via linear regression
        if (mneNotPositive) {
          dirtyDog <- TRUE
          msg <- gettext("Fits from method==\'ModNegExp\' are not all positive. \n  See constrain.nls argument in detrend.series. \n  ARSTAN would tell you to plot that dirty dog at this point.\n  Proceed with caution.",
                         domain = "R-dplR")
          if(y.name2==""){
            msg2 <- gettext(msg, domain = "R-dplR")
          }
          else {
            msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                      gettext(msg, domain = "R-dplR"))
          }
          warning(msg2)
          if (verbose) {
            cat(sepLine, indent(msg), sepLine, sep = "\n")
          }
        }
        x <- seq_len(nY2)
        lm1 <- lm(y2 ~ x)
        coefs <- coef(lm1)
        xIdx <- names(coefs) == "x"
        coefs <- c(coefs[!xIdx], coefs[xIdx])
        if (verbose) {
          cat(indent(c(gettext("Linear model fit", domain = "R-dplR"),
                       gettextf("Intercept: %s", format(coefs[1]),
                                domain = "R-dplR"),
                       gettextf("Slope: %s", format(coefs[2]),
                                domain = "R-dplR"))),
              sep = "\n")
        }
        if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
          tm <- cbind(1, x)
          ModNegExp <- drop(tm %*% coefs)
          useMean <- !isTRUE(ModNegExp[1] > 0 &&
                               ModNegExp[nY2] > 0)
          if (useMean) {
            dirtyDog <- TRUE
            msg <- gettext("Linear fit (backup of method==\'ModNegExp\') is not all positive. \n  Proceed with caution. \n  ARSTAN would tell you to plot that dirty dog at this point.",
                           domain = "R-dplR")
            if(y.name2==""){
              msg2 <- gettext(msg, domain = "R-dplR")
            }
            else {
              msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                        gettext(msg, domain = "R-dplR"))
            }
            warning(msg2)
            if (verbose) {
              cat(sepLine, indent(msg), sepLine, sep = "\n")
            }
          }
        } else {
          useMean <- TRUE
        }
        if (useMean) {
          theMean <- mean(y2)
          if (verbose) {
            cat(indent(c(gettext("lm has a positive slope",
                                 "pos.slope = FALSE",
                                 "Detrend by mean.",
                                 domain = "R-dplR"),
                         gettextf("Mean = %s", format(theMean),
                                  domain = "R-dplR"))),
                sep = "\n")
          }
          ModNegExp <- rep.int(theMean, nY2)
          mneStats <- list(method = "Mean", mean = theMean)
        } else {
          mneStats <- list(method = "Line", coefs = coef(summary(lm1)))
        }
      } else if (verbose || return.info) {
        mneSummary <- attr(ModNegExp, "summary")
        mneCoefs <- mneSummary[["coefficients"]]
        mneCoefsE <- mneCoefs[, 1]
        if (verbose) {
          cat(indent(c(gettext("nls coefs", domain = "R-dplR"),
                       paste0(names(mneCoefsE), ": ",
                              format(mneCoefsE)))),
              sep = "\n")
        }
        mneStats <- list(method = "NegativeExponential",
                         is.constrained = attr(ModNegExp, "constrained"),
                         formula = attr(ModNegExp, "formula"),
                         coefs = mneCoefs)
      } else {
        mneStats <- NULL
      }
      if(difference){ resids$ModNegExp <- y2 - ModNegExp }
      else{ resids$ModNegExp <- y2 / ModNegExp }
      curves$ModNegExp <- ModNegExp
      modelStats$ModNegExp <- mneStats
      do.mne <- TRUE
    } else {
      do.mne <- FALSE
    }
    ################################################################################    
    if("ModHugershoff" %in% method2){
      ## hug or lm
      hug.func <- function(Y, constrain) {
        nY <- length(Y)
        a <- mean(Y[floor(nY * 0.9):nY])
        b <- 1
        g <- 0.1
        d <- mean(Y[floor(nY * 0.9):nY])
        nlsForm <- Y ~ I(a*seq_along(Y)^b*exp(-g*seq_along(Y))+d)
        nlsStart <- list(a=a, b=b, g=g, d=d)
        checked <- FALSE
        constrained <- FALSE
        ## Note: nls() may signal an error
        if (constrain == "never") {
          hug <- nls(formula = nlsForm, start = nlsStart)
        } else if (constrain == "always") {
          hug <- nls(formula = nlsForm, start = nlsStart,
                     lower = c(a=0, b=-Inf, g=0, d=0),
                     upper = c(a=Inf, b=0, g=Inf, d=Inf),
                     algorithm = "port")
          constrained <- TRUE
        } else {
          hug <- nls(formula = nlsForm, start = nlsStart)
          coefs <- coef(hug)
          if (coefs[1] <= 0 || coefs[2] <= 0) {
            stop()
          }
          fits <- predict(hug)
          if (fits[nY] > 0) {
            checked <- TRUE
          } else {
            hug <- nls(formula = nlsForm, start = nlsStart,
                       lower = c(a=0, b=-Inf, g=0, d=0),
                       upper = c(a=Inf, b=0, g=Inf, d=Inf),
                       algorithm = "port")
            constrained <- TRUE
          }
        }
        if (!checked) {
          coefs <- coef(hug)
          if (coefs[1] <= 0 || coefs[2] <= 0) {
            stop()
          }
          fits <- predict(hug)
          if (fits[nY] <= 0) {
            ## This error is a special case that needs to be
            ## detected (if only for giving a warning).  Any
            ## smarter way to implement this?
            return(NULL)
          }
        }
        tmpFormula <- nlsForm
        formEnv <- new.env(parent = environment(detrend.series))
        formEnv[["Y"]] <- Y
        formEnv[["a"]] <- coefs["a"]
        formEnv[["b"]] <- coefs["b"]
        formEnv[["g"]] <- coefs["g"]
        formEnv[["d"]] <- coefs["d"]
        environment(tmpFormula) <- formEnv
        structure(fits, constrained = constrained,
                  formula = tmpFormula, summary = summary(hug))
      }
      ModHugershoff <- try(hug.func(y2, constrain2), silent=TRUE)
      hugNotPositive <- is.null(ModHugershoff)
      
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(gettext("Detrend by ModHugershoff.\n", domain = "R-dplR")))
        cat(indent(gettext("Trying to fit nls model...\n",
                           domain = "R-dplR")))
      }
      if (hugNotPositive || inherits(ModHugershoff,"try-error")) {
        if (verbose) {
          cat(indent(gettext("nls failed... fitting linear model...",
                             domain = "R-dplR")))
        }
        ## Straight line via linear regression
        if (hugNotPositive) {
          dirtyDog <- TRUE
          msg <- gettext("Fits from method==\'ModHugershoff\' are not all positive. \n  See constrain.nls argument in detrend.series. \n  ARSTAN would tell you to plot that dirty dog at this point.\n  Proceed with caution.",
                         domain = "R-dplR")
          if(y.name2==""){
            msg2 <- gettext(msg, domain = "R-dplR")
          }
          else {
            msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                      gettext(msg, domain = "R-dplR"))
          }
          warning(msg2)
          if (verbose) {
            cat(sepLine, indent(msg), sepLine, sep = "\n")
          }
        }
        x <- seq_len(nY2)
        lm1 <- lm(y2 ~ x)
        coefs <- coef(lm1)
        xIdx <- names(coefs) == "x"
        coefs <- c(coefs[!xIdx], coefs[xIdx])
        if (verbose) {
          cat(indent(c(gettext("Linear model fit", domain = "R-dplR"),
                       gettextf("Intercept: %s", format(coefs[1]),
                                domain = "R-dplR"),
                       gettextf("Slope: %s", format(coefs[2]),
                                domain = "R-dplR"))),
              sep = "\n")
        }
        if (all(is.finite(coefs)) && (coefs[2] <= 0 || pos.slope)) {
          tm <- cbind(1, x)
          ModHugershoff <- drop(tm %*% coefs)
          useMean <- !isTRUE(ModHugershoff[1] > 0 &&
                               ModHugershoff[nY2] > 0)
          if (useMean) {
            dirtyDog <- TRUE
            msg <- gettext("Linear fit (backup of method==\'ModHugershoff\') is not all positive. \n  ARSTAN would tell you to plot that dirty dog at this point. \n  Proceed with caution.",
                           domain = "R-dplR")
            if(y.name2==""){
              msg2 <- gettext(msg, domain = "R-dplR")
            }
            else {
              msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                        gettext(msg, domain = "R-dplR"))
            }
            warning(msg2)
            if (verbose) {
              cat(sepLine, indent(msg), sepLine, sep = "\n")
            }
            
          }
        } else {
          useMean <- TRUE
        }
        if (useMean) {
          theMean <- mean(y2)
          if (verbose) {
            cat(indent(c(gettext("lm has a positive slope",
                                 "pos.slope = FALSE",
                                 "Detrend by mean.",
                                 domain = "R-dplR"),
                         gettextf("Mean = %s", format(theMean),
                                  domain = "R-dplR"))),
                sep = "\n")
          }
          ModHugershoff <- rep.int(theMean, nY2)
          hugStats <- list(method = "Mean", mean = theMean)
        } else {
          hugStats <- list(method = "Line", coefs = coef(summary(lm1)))
        }
      } else if (verbose || return.info) {
        hugSummary <- attr(ModHugershoff, "summary")
        hugCoefs <- hugSummary[["coefficients"]]
        hugCoefsE <- hugCoefs[, 1]
        if (verbose) {
          cat(indent(c(gettext("nls coefs", domain = "R-dplR"),
                       paste0(names(hugCoefsE), ": ",
                              format(hugCoefsE)))),
              sep = "\n")
        }
        hugStats <- list(method = "Hugershoff",
                         is.constrained = attr(ModHugershoff, "constrained"),
                         formula = attr(ModHugershoff, "formula"),
                         coefs = hugCoefs)
      } else {
        hugStats <- NULL
      }
      if(difference){ resids$ModHugershoff <- y2 - ModHugershoff }
      else{ resids$ModHugershoff <- y2 / ModHugershoff }
      curves$ModHugershoff <- ModHugershoff
      modelStats$ModHugershoff <- hugStats
      do.hug <- TRUE
    } else {
      do.hug <- FALSE
    }
    ################################################################################    
    if("AgeDepSpline" %in% method2){
      ## Age dep smoothing spline with nyrs (50 default) as the init stiffness
      ## are NULL
      if(is.null(nyrs))
        nyrs2 <- 50
      else
        nyrs2 <- nyrs
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(c(gettext(c("Detrend by age-dependent spline.",
                               "Spline parameters"), domain = "R-dplR"),
                     paste0("nyrs = ", nyrs2, ", pos.slope = ", pos.slope))),
            sep = "\n")
      }
      AgeDepSpline <- ads(y=y2, nyrs0=nyrs2, pos.slope = pos.slope)
      if (any(AgeDepSpline <= 0)) {
        dirtyDog <- TRUE
        msg <- "Fits from method==\'AgeDepSpline\' are not all positive. \n  This is extremely rare. Series will be detrended with method==\'Mean\'. \n  This might not be what you want. \n  ARSTAN would tell you to plot that dirty dog at this point. \n  Proceed with caution."
        if(y.name2==""){
          msg2 <- gettext(msg, domain = "R-dplR")
        }
        else {
          msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                    gettext(msg, domain = "R-dplR"))
        }
        warning(msg2)
        if (verbose) {
          cat(sepLine, indent(msg), sepLine, sep = "\n")
        }
        theMean <- mean(y2)
        AgeDepSpline <- rep.int(theMean, nY2)
        AgeDepSplineStats <- list(method = "Mean", mean = theMean)
      } else {
        AgeDepSplineStats <- list(method = "Age-Dep Spline", nyrs = nyrs2, pos.slope=pos.slope)
      }
      if(difference){ resids$AgeDepSpline <- y2 - AgeDepSpline }
      else{ resids$AgeDepSpline <- y2 / AgeDepSpline }
      curves$AgeDepSpline <- AgeDepSpline
      modelStats$AgeDepSpline <- AgeDepSplineStats
      
      do.ads <- TRUE
    } else {
      do.ads <- FALSE
    }
    ################################################################################    
    if("Spline" %in% method2){
      ## Smoothing spline
      ## "n-year spline" as the spline whose frequency response is
      ## 50%, or 0.50, at a wavelength of 67%n years if nyrs and f
      ## are NULL
      if(is.null(nyrs))
        nyrs2 <- floor(nY2 * 0.67)
      else
        nyrs2 <- nyrs
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(c(gettext(c("Detrend by spline.",
                               "Spline parameters"), domain = "R-dplR"),
                     paste0("nyrs = ", nyrs2, ", f = ", f))),
            sep = "\n")
      }
      #Spline <- ffcsaps(y=y2, x=seq_len(nY2), nyrs=nyrs2, f=f)
      Spline <- caps(y=y2, nyrs=nyrs2, f=f)
      if (any(Spline <= 0)) {
        dirtyDog <- TRUE
        msg <- "Fits from method==\'Spline\' are not all positive. \n  Series will be detrended with method==\'Mean\'. \n  This might not be what you want. \n  ARSTAN would tell you to plot that dirty dog at this point. \n  Proceed with caution."
        if(y.name2==""){
          msg2 <- gettext(msg, domain = "R-dplR")
        }
        else {
          msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                    gettext(msg, domain = "R-dplR"))
        }
        warning(msg2)
        if (verbose) {
          cat(sepLine, indent(msg), sepLine, sep = "\n")
        }
        theMean <- mean(y2)
        Spline <- rep.int(theMean, nY2)
        splineStats <- list(method = "Mean", mean = theMean)
      } else {
        splineStats <- list(method = "Spline", nyrs = nyrs2, f = f)
      }
      if(difference){ resids$Spline <- y2 - Spline }
      else{ resids$Spline <- y2 / Spline }
      curves$Spline <- Spline
      modelStats$Spline <- splineStats
      
      do.spline <- TRUE
    } else {
      do.spline <- FALSE
    }
    ################################################################################    
    if("Mean" %in% method2){
      ## Fit a horiz line
      theMean <- mean(y2)
      Mean <- rep.int(theMean, nY2)
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(c(gettext("Detrend by mean.", domain = "R-dplR"),
                     paste("Mean = ", format(theMean)))),
            sep = "\n")
      }
      meanStats <- list(method = "Mean", mean = theMean)
      if(difference){ resids$Mean <- y2 - Mean }
      else{ resids$Mean <- y2 / Mean }
      curves$Mean <- Mean
      modelStats$Mean <- meanStats
      do.mean <- TRUE
    } else {
      do.mean <- FALSE
    }
    ################################################################################    
    if("Ar" %in% method2){
      ## Fit an ar model - aka prewhiten
      Ar <- ar.func(y2, model = TRUE)
      arModel <- attr(Ar, "model")
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(gettext("Detrend by prewhitening.", domain = "R-dplR")))
        print(arModel)
      }
      arStats <- list(method = "Ar", order = arModel[["order"]],
                      ar = arModel[["ar"]])
      # This will propagate NA to rwi as a result of detrending.
      # Other methods don't. Problem when interacting with other
      # methods?
      # Also, this can (and does!) produce negative RWI values.
      # See example using CAM011. Thus:
      if (any(Ar <= 0, na.rm = TRUE)) {
        dirtyDog <- TRUE
        msg <- "Fits from method==\'Ar\' are not all positive. \n  Setting values <0 to 0 before rescaling.  \n  This might not be what you want. \n  ARSTAN would tell you to plot that dirty dog at this point. \n  Proceed with caution."
        if(y.name2==""){
          msg2 <- gettext(msg, domain = "R-dplR")
        }
        else {
          msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                    gettext(msg, domain = "R-dplR"))
        }
        warning(msg2)
        if (verbose) {
          cat(sepLine, indent(msg), sepLine, sep = "\n")
        }
        
        Ar[Ar<0] <- 0
      }
      if(difference){ resids$Ar <- Ar - mean(Ar,na.rm=TRUE) }
      else{ resids$Ar <- Ar / mean(Ar,na.rm=TRUE) }
      curves$Ar <- mean(Ar,na.rm=TRUE)
      modelStats$Ar <- arStats
      do.ar <- TRUE
    } else {
      do.ar <- FALSE
    }
    ################################################################################    
    if ("Friedman" %in% method2) {
      if (is.null(wt.description)) {
        wt.description <- if (wt.missing) "default" else deparse(wt)
      }
      if (verbose) {
        cat("", sepLine, sep = "\n")
        cat(indent(c(gettext(c("Detrend by Friedman's super smoother.",
                               "Smoother parameters"), domain = "R-dplR"),
                     paste0("span = ", span, ", bass = ", bass),
                     paste0("wt = ", wt.description))),
            sep = "\n")
      }
      if (wt.missing) {
        Friedman <- supsmu(x = seq_len(nY2), y = y2, span = span,
                           periodic = FALSE, bass = bass)[["y"]]
      } else {
        Friedman <- supsmu(x = seq_len(nY2), y = y2, wt = wt, span = span,
                           periodic = FALSE, bass = bass)[["y"]]
      }
      if (any(Friedman <= 0)) {
        dirtyDog <- TRUE
        msg <- "Fits from method==\'Friedman\' are not all positive. \n  Series will be detrended with method==\'Mean\'. \n  This might not be what you want. \n  ARSTAN would tell you to plot that dirty dog at this point. \n  Proceed with caution."
        if(y.name2==""){
          msg2 <- gettext(msg, domain = "R-dplR")
        }
        else {
          msg2 <- c(gettextf("In raw series %s: ", y.name2, domain = "R-dplR"),
                    gettext(msg, domain = "R-dplR"))
        }
        warning(msg2)
        if (verbose) {
          cat(sepLine, indent(msg), sepLine, sep = "\n")
        }
        
        
        theMean <- mean(y2)
        Friedman <- rep.int(theMean, nY2)
        friedmanStats <- list(method = "Mean", mean = theMean)
      } else {
        friedmanStats <- list(method = "Friedman", wt = wt.description, span = span, bass = bass)
      }
      
      if(difference){ resids$Friedman <- y2 - Friedman }
      else{ resids$Friedman <- y2 / Friedman }
      curves$Friedman <- Friedman
      modelStats$Friedman <-
        list(method = "Friedman",
             wt = if (wt.missing) "default" else wt,
             span = span, bass = bass)
      do.friedman <- TRUE
    } else {
      do.friedman <- FALSE
    }
    ################################################################################    
    ################################################################################    
    
    resids <- data.frame(resids)
    curves <- data.frame(curves)
    if (verbose || return.info) {
      zero.years <- lapply(resids, zeroFun)
      n.zeros <- lapply(zero.years, nFun)
      modelStats <- mapply(c, modelStats, n.zeros, zero.years,
                           SIMPLIFY = FALSE)
      if (verbose) {
        n.zeros2 <- unlist(n.zeros, use.names = FALSE)
        zeroFlag <- n.zeros2 > 0
        methodNames <- names(modelStats)
        if (any(zeroFlag)) {
          cat("", sepLine, sep = "\n")
          for (i in which(zeroFlag)) {
            if (is.character(years)) {
              cat(indent(gettextf("Zero years in %s series:\n",
                                  methodNames[i], domain="R-dplR")))
            } else {
              cat(indent(gettextf("Zero indices in %s series:\n",
                                  methodNames[i], domain="R-dplR")))
            }
            cat(indent(paste(zero.years[[i]][[1]], collapse = " ")),
                "\n", sep = "")
          }
        }
      }
    }
    
    if(make.plot){
      cols <- c("#24492e","#015b58","#2c6184","#59629b","#89689d","#ba7999","#e69b99")
      op <- par(no.readonly=TRUE)
      on.exit(par(op))
      n.methods <- ncol(resids)
      par(mar=c(2.1, 2.1, 2.1, 2.1), mgp=c(1.1, 0.1, 0),
          tcl=0.5, xaxs="i")
      if (n.methods > 4) {
        par(cex.main = min(1, par("cex.main")))
      }
      mat <- switch(n.methods,
                    matrix(c(1,2), nrow=2, ncol=1, byrow=TRUE),
                    matrix(c(1,1,2,3), nrow=2, ncol=2, byrow=TRUE),
                    matrix(c(1,2,3,4), nrow=2, ncol=2, byrow=TRUE),
                    matrix(c(1,1,2,3,4,5), nrow=3, ncol=2, byrow=TRUE),
                    matrix(c(1,1,1,2,3,4,5,6,0), nrow=3, ncol=3, byrow=TRUE),
                    matrix(c(1,1,1,2,3,4,5,6,7), nrow=3, ncol=3, byrow=TRUE),
                    matrix(c(1,2,3,4,5,6,7,8), nrow=4, ncol=2, byrow=TRUE))
      
      layout(mat,
             widths=rep.int(0.5, ncol(mat)),
             heights=rep.int(1, nrow(mat)))
      
      # 1
      plot(y2, type="l", ylab="mm", col = "grey",
           xlab=gettext("Age (Yrs)", domain="R-dplR"),
           main=gettextf("Raw Series %s", y.name2, domain="R-dplR"))
      if(do.spline) lines(Spline, col=cols[1], lwd=2)
      if(do.mne) lines(ModNegExp, col=cols[2], lwd=2)
      if(do.mean) lines(Mean, col=cols[3], lwd=2)
      if(do.friedman) lines(Friedman, col=cols[5], lwd=2)
      if(do.hug) lines(ModHugershoff, col=cols[6], lwd=2)
      if(do.ads) lines(AgeDepSpline, col=cols[7], lwd=2)
      
      # 1
      if(do.spline){
        plot(resids$Spline, type="l", col=cols[1],
             main=gettext("Spline", domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
      }
      # 2
      if(do.mne){
        plot(resids$ModNegExp, type="l", col=cols[2],
             main=gettext("Neg. Exp. Curve or Straight Line",
                          domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        
      }
      # 3
      if(do.mean){
        plot(resids$Mean, type="l", col=cols[3],
             main=gettext("Horizontal Line (Mean)", domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        
      }
      # 4
      if(do.ar){
        plot(resids$Ar, type="l", col=cols[4],
             main=gettextf("Ar", domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        mtext(text="(Not plotted with raw series)",side=3,line=-1,cex=0.75)
      }
      # 5
      if (do.friedman) {
        plot(resids$Friedman, type="l", col=cols[5],
             main=gettext("Friedman's Super Smoother", domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        
      }
      # 6
      if(do.hug){
        plot(resids$ModHugershoff, type="l", col=cols[6],
             main=gettext("Hugershoff or Straight Line",
                          domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        
      }
      # 7
      if(do.ads){
        plot(resids$AgeDepSpline, type="l", col=cols[7],
             main=gettext("Age Dep Spline",
                          domain="R-dplR"),
             xlab=gettext("Age (Yrs)", domain="R-dplR"),
             ylab=gettext("RWI", domain="R-dplR"))
        if(difference){ abline(h=0) }
        else{ abline(h=1) }
        
      }
    }
      # Done get output
      resids2 <- matrix(NA, ncol=ncol(resids), nrow=length(y))
      resids2 <- data.frame(resids2)
      names(resids2) <- names(resids)
      if(!is.null(names(y))) row.names(resids2) <- names(y)
      resids2[good.y, ] <- resids
      
      curves2 <- matrix(NA, ncol=ncol(curves), nrow=length(y))
      curves2 <- data.frame(curves2)
      names(curves2) <- names(curves)
      if(!is.null(names(y))) row.names(curves2) <- names(y)
      curves2[good.y, ] <- curves
      ## Reorder columns of output to match the order of the argument
      ## "method".
      resids2 <- resids2[, method2]
      curves2 <- curves2[, method2]
      ## Make sure names (years) are included if there is only one method
      if(!is.data.frame(resids2)) names(resids2) <- names(y)
      if (return.info) {
        list(series = resids2,
             curves = curves2,
             model.info = modelStats[method2],
             data.info = dataStats,
             dirtyDog = dirtyDog)
      } else {
        resids2
      }
    }