R/f.quad.plateau.R

In FertBoot: Fertilizer Response Curve Analysis by Bootstrapping Residuals

#'
#' Fitting quadratic plateau model using multiple initial values
#'
#' \code{f.quad.plateau} fits quadratic plateau model using multiple initial values. The multiple initial values are randomly sampled in  a "cube" of parameter space. More precisely, quadratic plateau model assumes
#'     y ~ (a + b * x + c *x^2) * (x <= -0.5*b/c) + (a + -b^2/(4 * c)) * (x > -0.5 * b/c).
#'
#' @param d data frame with two columns (\code{x} and \code{y})
#' @param start initial estimate for non-linear least square (default value: \code{ list(a = 1, b = 1, c = 1)})
#' @param plus_minus radius of random initial values (default: \code{100})
#' @param n.start total number of initial points considered (default: \code{1000})
#' @param msg logical flag whether printing progress
#'
#' @return \code{f.quad.plateau} returns a list of two components (if converged): \code{nls.summary}: summary of the fitted model; \code{nls.model}: nls object
#' @import stats nls.multstart simpleboot
#' @examples
#'
#' set.seed(3)
#' x <- rep(1:300, each=4)
#' a <- 8; b <- 0.05; c <- -1e-4
#' y <- (a + b * x + c *x^2) * (x <= -0.5*b/c) + (a + -b^2/(4 * c)) * (x > -0.5 * b/c) +
#'     rnorm(length(x), sd=0.1)
#' d <- cbind(x,y)
#'
#' # a converged example:
#' ans <- f.quad.plateau(d, start=list(a = 7, b = 0.02, c = 1e-5),
#'     plus_minus=10, n.start=10, msg=FALSE)
#'
#' summary(ans$nls.model)
#'
#'
#'
#' @export
#'
f.quad.plateau <- function(d, start = list(a = 1, b = 1, c = 1),
                          plus_minus = 1e2, n.start=1000, msg=FALSE) {

  if ("matrix" %in% class(d)) {d <- data.frame(d)}

  names(d) = c("x", "y")
  mq = lm(y ~ x + I(x^2), data=d)
  c3 = coef(mq)[[1]]
  c4 = coef(mq)[[2]]
  c5 = coef(mq)[[3]]
  a = start[[1]]
  b = start[[2]]
  c = start[[3]]
  a11 = ifelse(a == 1, c3, a)
  b11 = ifelse(b == 1, c4, b)
  c11 = ifelse(c == 1, c5, c)
  m <-  nls.multstart::nls_multstart(y ~ (a + b * x + c * I(x^2)) * (x <= -0.5 *
                                                        b/c) + (a + I(-b^2/(4 * c))) * (x > -0.5 * b/c),
                      start_lower=list(a = a11 - plus_minus, b=b11 - plus_minus, c=c11 - plus_minus),
                      start_upper=list(a = a11 + plus_minus, b=b11 + plus_minus, c=c11 + plus_minus),
                      iter=n.start, data = d, supp_errors ="Y")
  if(!is.null(m)) {
    c = coef(m)
    a = c[1]
    b = c[2]
    c = c[3]
    pmm = (a + I(-b^2/(4 * c)))
    pcc = -0.5 * b/c
    nls.summary = c(a, b, c, summary(m)[11][[1]][10], summary(m)[11][[1]][11],
                    summary(m)[11][[1]][12], AIC(m), BIC(m),
                    pmm, pcc)
    # nls.summary = round(nls.summary, 8)
    nls.summary = as.data.frame(nls.summary)
    rownames(nls.summary) = c("coefficient a", "coefficient b", "coefficient c",
                              "p-value t.test for a", "p-value t.test for b", "p-value t.test for c",
                              "AIC", "BIC",
                              "maximum or minimum value for y", "critical point in x")
    res <- list(nls.summary=nls.summary, nls.model=m)
  } else{
    res <- NULL
    if (msg) cat("Not converged! may try larger n.start!\n")
  }
  res
}