R/boot_nls.R

In nlraa: Nonlinear Regression for Agricultural Applications

#' Bootstraping for nonlinear models 
#' 
#' @title Bootstrapping for nonlinear models 
#' @name boot_nls
#' @param object object of class \code{\link[stats]{nls}}
#' @param f function to be applied (and bootstrapped), default coef 
#' @param R number of bootstrap samples, default 999
#' @param psim simulation level for \code{\link{simulate_nls}}
#' @param resid.type either \dQuote{resample}, \dQuote{normal} or \dQuote{wild}.
#' @param data optional data argument (useful/needed when data are not in an available environment).
#' @param verbose logical (default TRUE) whether to print a message if model does not converge. 
#' @param ... additional arguments to be passed to function \code{\link[boot]{boot}}
#' @details The residuals can either be generated by resampling with replacement 
#' (default or non-parametric), from a normal distribution (parameteric) or by changing 
#' their signs (wild). This last one is called \dQuote{wild bootstrap}. 
#' There is more information in \code{\link{boot_lm}}.
#' @note at the moment, when the argument data is used, it is not possible to check that it matches the 
#' original data used to fit the model. It will also override the fetching of data.
#' @seealso \code{\link[car]{Boot}}
#' @export
#' @examples 
#' \donttest{
#' require(car)
#' data(barley, package = "nlraa")
#' ## Fit a linear-plateau
#' fit.nls <- nls(yield ~ SSlinp(NF, a, b, xs), data = barley)
#' 
#' ## Bootstrap coefficients by default
#' ## Keeping R small for simplicity, increase R for a more realistic use
#' fit.nls.bt <- boot_nls(fit.nls, R = 1e2)
#' ## Compute confidence intervals
#' confint(fit.nls.bt, type = "perc")
#' ## Visualize
#' hist(fit.nls.bt, 1, ci = "perc", main = "Intercept")
#' hist(fit.nls.bt, 2, ci = "perc", main = "linear term")
#' hist(fit.nls.bt, 3, ci = "perc", main = "xs break-point term")
#' }
#' 

boot_nls <- function(object, 
                    f = NULL, 
                    R = 999, 
                    psim = 2, 
                    resid.type = c("resample","normal","wild"),
                    data = NULL,
                    verbose = TRUE,
                    ...){
  ## I chose to write it in this way instead of UseMethod
  ## because I feel it is more efficient and results in less code
  ## Maybe I'm wrong and will change this in the future
  ## Error checking
  if(!inherits(object, "nls")) stop("object should be of class 'nls'")
  
  if(psim < 2) stop("psim should be 2 or higher")
  
  resid.type <- match.arg(resid.type)
  
  ## extract the data
  dat <- try(eval(getCall(object)$data), silent = TRUE)
  if(inherits(dat, "try-error") && missing(data)){
    stop("'data' argument is required. It is likely you are using boot_lm inside another function. 
          The data are not in an available environment.", call. = FALSE)    
  } 
  
  if(!missing(data)){
    dat <- data
  }

  if(missing(f)) f <- coef
  
  f0 <- f(object) ## To model the behavior of the orignial function
  NA.j <- 0
  
  boot_fun_resid <- function(data, fn, model, psim,...){
    
    ## Copy that will be later modified
    bdat <- data
    ## I need a hack to make the first iteration be t0 for boot
    if(identical(get(".k.boot.nls", envir = nlraa.nls.env), 0L)){
      ## This makes the assumption that the first time
      ## boot calls 'boot_fun_resid' it will generate 't0' 
      new.y <- fitted(model) + resid(model)
      assign(".k.boot.nls", 1L, envir = nlraa.nls.env)
    }else{
      ## Fitted values, which are simulated values if psim = 1
      new.y <- simulate_nls_one(model, psim = psim, data = data, ...)
    }
    
    resp.var <- all.vars(stats::formula(model))[1]
    bdat[[resp.var]] <- new.y
    
    assign(".bdat.nls", bdat, envir = nlraa.nls.env)
    umod <- tryCatch(update(model, data = get(".bdat.nls", envir = nlraa.nls.env)), error = function(e){invisible(e)})
    
    ## It should be rare for an lm model to fail to converge though
    ## Trying to catch any condition in which the model above does not converge
    if(inherits(umod, "error") || any(is.na(fn(umod))) || is.null(fn(umod)) || any(is.nan(fn(umod)))){
      out <- rep(NA, length(f0))
      NA.j <<- NA.j + 1
    }else{
      out <- fn(umod)
    }
    out
  }
  
  ans <- boot::boot(data = dat, 
                    sim = "parametric",
                    statistic = boot_fun_resid, 
                    R = R, fn = f,
                    model = object,
                    psim = psim, ...)
  
  ## There is no reason why a lm should not converge, but just in case
  if(sum(NA.j) > 0 && verbose) cat("Number of times model fit did not converge", NA.j, "out of", R, "\n")
  
  assign(".k.boot.nls", 0L, envir = nlraa.nls.env)
  assign(".bdat.nls", NULL, envir = nlraa.nls.env)
  return(ans)
}


### Just initializing variables in the environment
### Not sure if this is even necessary
nlraa.nls.env <- new.env(parent = emptyenv())
assign('.bdat.nls', NA, nlraa.nls.env)
assign('.k.boot.nls', 0L, nlraa.nls.env)