version 0.4.0

Author: Ting Fung (Ralph) Ma

DESCRIPTION

@@ -0,0 +1,21 @@
+Package: FertBoot
+Type: Package
+Title: Fertilizer Response Curve Analysis by Bootstrapping Residuals
+Version: 0.4.0
+Authors@R: c(person("Ting Fung (Ralph)", "Ma", email = "[email protected]", role = c("cre", "aut")),
+    person("Hannah", "Francis", email = "[email protected]", role = "aut"),
+    person("Matt", "Ruark", role = "ctb"))
+Maintainer: Ting Fung (Ralph) Ma <[email protected]>
+Description: Quantify variability (such as confidence interval) of fertilizer response curves and optimum fertilizer rates using bootstrapping residuals with several popular non-linear and linear models.
+Imports: stats, nls.multstart, simpleboot
+License: GPL-2
+Encoding: UTF-8
+LazyData: true
+RoxygenNote: 7.0.2
+NeedsCompilation: no
+Packaged: 2020-07-14 17:17:17 UTC; Ralph
+Author: Ting Fung (Ralph) Ma [cre, aut],
+  Hannah Francis [aut],
+  Matt Ruark [ctb]
+Repository: CRAN
+Date/Publication: 2020-07-18 09:32:17 UTC

MD5

@@ -0,0 +1,18 @@
+cf46442d80c8f64d5c98d450cabd65e6 *DESCRIPTION
+c90f950aef0038c9f28f95b0278153a3 *NAMESPACE
+e51d94089628154d5b91a780fb362c08 *R/boot.CI.R
+a790ca6b477fa3eac45312dc1f910acc *R/boot.resid.linear.plateau.R
+ad128e3018f4614f772f1ab873582fb3 *R/boot.resid.quad..R
+27a49d0e3aea09e667e574d8354cbcae *R/boot.resid.quad.plateau.R
+afc7548a77a84882e90c54b9da65300c *R/compare.two.sample.R
+78200e1f92ae26d08ad46c594e994c2f *R/f.linear.plateau.R
+f812fc47d9a078d5931d3705b546f82f *R/f.quad.R
+f29ad8c1d0c22cdc0f802389f6a5547d *R/f.quad.plateau.R
+c8b1617d446d92b00279b22588aa4ec9 *man/boot.CI.Rd
+2c303a1c98a115a3dbc6dd6ea76464d0 *man/boot.resid.linear.plateau.Rd
+bf0f9c3c8e610bbc5ab8cc9f114ed0b2 *man/boot.resid.quad.Rd
+0b6da7c23d44131a62073514cb0d5ecf *man/boot.resid.quad.plateau.Rd
+4674a98de8082bf9e84e8bcaea26e91a *man/compare.two.sample.Rd
+c289928245408a5784a745e3527c7607 *man/f.linear.plateau.Rd
+894f7c8610174e917ca7cf9f209be6f6 *man/f.quad.Rd
+25e15a97f7e2873c5acc1115a4b50607 *man/f.quad.plateau.Rd

NAMESPACE

@@ -0,0 +1,13 @@
+# Generated by roxygen2: do not edit by hand
+
+export(boot.CI)
+export(boot.resid.linear.plateau)
+export(boot.resid.quad)
+export(boot.resid.quad.plateau)
+export(compare.two.sample)
+export(f.linear.plateau)
+export(f.quad)
+export(f.quad.plateau)
+import(nls.multstart)
+import(simpleboot)
+import(stats)

R/boot.CI.R

@@ -0,0 +1,58 @@
+#' Bootstrap confidence intervals of mean
+#' @param x a vector of observation
+#' @param alpha significance level (default: 0.05)
+#' @param CI.type type of CI required (default: "all")
+#'
+#' @return \code{boot.CI} return list of confidence intervals of mean (\code{CI.percent}: percentile, \code{CI.BC}:bias-corrected and \code{CI.BCa}: bias-corrected and accelerated)
+#' @import stats nls.multstart simpleboot
+#' @examples
+#'
+#' set.seed(12)
+#' boot.CI(rnorm(1000, mean=0, sd=1), alpha=0.05, CI.type="per") # example of wrong input for type
+#' boot.CI(rnorm(1000, mean=0, sd=1), alpha=0.05, CI.type="all") # require all type
+#'
+#' @export
+#'
+boot.CI <- function(x, alpha=0.05, CI.type="all") {
+  B <- length(x)
+
+  result <- NULL
+
+  if ("percentile" %in% CI.type || "all" %in% CI.type) {
+  # Percentile CI
+  CI.percent <- quantile(x, c(alpha/2, 1 - alpha/2))
+
+  result <- c(result, list(CI.percent=CI.percent))
+  }
+
+
+  if ("BC" %in% CI.type || "all" %in% CI.type) {
+  # Bias-corrected CI
+  z0 <- qnorm(sum(x < mean(x))/B)
+  CI.BC <- quantile(x, c(pnorm(qnorm(alpha/2) + 2*z0), pnorm(qnorm(1-alpha/2) + 2*z0)))
+
+  result <- c(result, list(CI.BC=CI.BC))
+  }
+
+  if ("BCa" %in% CI.type || "all" %in% CI.type) {
+  # Accelerated bias-corrected CI
+  theta <- rep(NA, B)
+  for (i in 1:B) {theta[i] <- mean(x[-i])}
+  theta.mean <- mean(theta)
+
+  a <- sum((theta.mean - theta)^3)/(6*sum((theta.mean - theta)^2)^1.5)
+
+  alpha1 <- pnorm(z0 + (z0 + qnorm(alpha/2))/(1-a*(z0 + qnorm(alpha/2))))
+  alpha2 <- pnorm(z0 + (z0 + qnorm(1 - alpha/2))/(1-a*(z0 + qnorm(1 - alpha/2))))
+  CI.BCa <- quantile(x, c(alpha1, alpha2))
+
+  result <- c(result, list(CI.BCa=CI.BCa))
+  }
+
+  if (is.null(result)) {
+    cat("Auto-corrected with percentile CI!\n")
+    result <- boot.CI(x, alpha=0.05, CI.type="percentile")
+  }
+
+  result
+}

R/boot.resid.linear.plateau.R

@@ -0,0 +1,81 @@
+#' Linear plateau model estimation by bootstrapping residuals
+#'
+#' \code{boot.resid.linear.plateau} is the core function that implement bootstrapping residuals on quadratic plateau model, which assumes
+#'     y ~ a + b * (x - c) * (x <= c). Note that this functions may take minutes up to days. Parallel computing could be considered when necessary. We suggest users start with a smaller \code{B} and moderate \code{n.start} to see if the bootstrap models can convergence.  In general, increasing \code{n.start} and \code{plus_minus} may help convergence. For rigorous statistical inference, \code{B} should be reach order of thousand.
+#'
+#'
+#' @param mod a full model list, probably from \code{f.linear.plateau()}
+#' @param data data drame with two columns (\code{x} and \code{y})
+#' @param x.range vector of data.frame with one column for range of N rate of interested for prediction interval
+#' @param B bootstrap sample size
+#' @param plus_minus radius of random initial values (default: \code{100})
+#' @param n.start total number of initial points considered (deafult: \code{1000})
+#' @param print.progress logical flag whehter printing progress
+#'
+#' @return \code{boot.resid.linear.plateau} returns a list of two elements:
+#' \code{result}: matrix with B rows and columns containing bootstrap sample for parameter (\code{a,b,c}), optimal N and yield (\code{max_x, max_y}), log-likelihood (\code{logLik}) and N values of interest;
+#' \code{x.range}: range of x considered for prediction interval (same as \code{x.range} in vector form)
+#'
+#' @import stats nls.multstart simpleboot
+#'
+#' @examples
+#'
+#'\donttest{
+#' set.seed(1)
+#' x <- rep(1:300, each=4)
+#' a <- 8; b <- 0.05; c <- 100
+#' y <- a + b * (x - c) * (x <= c) +
+#'     rnorm(length(x), sd=1)
+#' d <- cbind(x,y)
+#'
+#' # a converged example:
+#' ans <- f.linear.plateau(d, start=list(a = 7, b = 0.1, c = 150),
+#'     plus_minus=10, n.start=10, msg=FALSE)
+#'
+#'
+#' ans.boot <- boot.resid.linear.plateau(ans, d, x.range=seq(0,280,by=40),
+#'     B=99, plus_minus = 1e2, n.start=1000, print.progress=TRUE)
+#'
+#' }
+#'
+#'
+#'
+#' @export
+#'
+boot.resid.linear.plateau <- function(mod, data, x.range=data.frame(x=seq(0,280,by=40)), B=1e2-1, plus_minus = 1e2, n.start=5000, print.progress=TRUE) {
+
+  if (class(x.range) == "numeric")  x.range <- data.frame(x=x.range)
+
+  # a, b, c, max x value, max y value, logLik
+  result <- data.frame(matrix(NA, nrow= B + 1, ncol=6 + NROW(x.range)))
+  names(result) <- c("a", "b", "c", "max_x", "max_y", "logLik", paste0("x_", x.range[,1]))
+
+  data <- data.frame(data)
+
+  data.tmp <- data ; names(data.tmp) <- c("x", "y")
+
+
+  fit.value <- fitted(mod$nls.model)
+  res.value <- residuals(mod$nls.model)
+  start.value <- list(a=coef(mod$nls.model)[1], b=coef(mod$nls.model)[2], c=coef(mod$nls.model)[3])
+
+  result[1, ] <- c(mod$nls.summary[c(1:3,10,9),], logLik(mod$nls.model),
+                   predict(mod$nls.model, newdata=x.range))
+  i <- 1
+  while (i <= B) {
+    if (print.progress) cat("Bootstrap residuals:", i,"out of",B,"\n")
+
+    # Bootstrap residual
+    data.tmp$y <- as.numeric(fit.value + sample(res.value, replace=TRUE))
+
+    m.tmp <- f.linear.plateau(d=data.tmp, start=start.value, plus_minus=plus_minus, n.start=n.start)
+
+    if(!is.null(m.tmp$nls.model)) {
+      result[i+1, ] <- c(m.tmp$nls.summary[c(1:3,10,9),], logLik(m.tmp$nls.model),
+                         predict(m.tmp$nls.model, newdata=x.range))
+      i <- i + 1
+    }else{ if (print.progress) cat("Not converged! Retry...\n")}
+  }
+  list(result=result, x.range = x.range[,1])
+}
+

R/boot.resid.quad..R

@@ -0,0 +1,82 @@
+#' Fitting quadratic plateau model using mutiple initial vaues
+#'
+#' \code{boot.resid.quad.plateau} is the core function that implement bootstrapping residuals on quadratic plateau model, which assumes
+#'     y = (a + b * x + c *x^2) * (x <= -0.5*b/c) + (a + -b^2/(4 * c)) * (x > -0.5 * b/c). Note that this functions may take minutes up to days. Parallel computing could be considered when necessary. We suggest users start with a smaller \code{B} and moderate \code{n.start} to see if the bootstrap models can convergence.  In general, increasing \code{n.start} and \code{plus_minus} may help convergence. For rigorous statistical inference, \code{B} should be reach order of thousand.
+#'
+#'
+#' @param mod a full model list, probably from \code{f.quad.plateau()}
+#' @param data data drame with two columns (\code{x} and \code{y})
+#' @param x.range vector of data.frame with one column for range of N rate of interested for prediction interval
+#' @param B bootstrap sample size
+#' @param plus_minus radius of random initial values (default: \code{100})
+#' @param n.start total number of initial points considered (deafult: \code{1000})
+#' @param print.progress logical flag whehter printing progress
+#'
+#' @return \code{boot.resid.quad.plateau} returns a list of two elements:
+#' \code{result}: matrix with B rows and columns containing bootstrap sample for parameter (\code{a,b,c}), optimal N and yield (\code{max_x, max_y}), log-likelihood (\code{logLik}) and N values of interest;
+#' \code{x.range}: range of x considered for prediction interval (same as \code{x.range} in vector form)
+#'
+#' @import stats nls.multstart simpleboot
+#'
+#' @examples
+#'
+#' \donttest{
+#'
+#' set.seed(1)
+#' x <- rep(1:300, each=5)
+#' a <- 8; b <- 0.05; c <- -1e-3
+#' 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)
+#'
+#' ans <- f.quad.plateau(d, start=list(a = 7, b = 0.02, c = 1e-5),
+#'     plus_minus=10, n.start=10, msg=FALSE)
+#'
+#'
+#' ans.boot <- boot.resid.quad.plateau(ans, d, x.range=seq(0,280,by=40),
+#'     B=1e1, plus_minus = 1e2, n.start=1e3, print.progress=TRUE)
+#'
+#' }
+#'
+#'
+#'
+#' @export
+#'
+boot.resid.quad <- function(mod, data, x.range=data.frame(x=seq(0,280,by=40)),
+                                  B=1e2-1, plus_minus = 1e1, n.start=20, print.progress=TRUE) {
+
+  if (class(x.range) == "numeric")  x.range <- data.frame(x=x.range)
+
+  # a, b, c, max x value, max y value, logLik
+  result <- data.frame(matrix(NA, nrow= B + 1, ncol=6 + NROW(x.range)))
+  names(result) <- c("a", "b", "c", "max_x", "max_y", "logLik", paste0("x_", x.range[,1]))
+
+  data.tmp <- data ; names(data.tmp) <- c("x", "y")
+
+
+
+  fit.value <- fitted(mod$nls.model)
+  res.value <- residuals(mod$nls.model)
+  start.value <- list(a=coef(mod$nls.model)[1], b=coef(mod$nls.model)[2], c=coef(mod$nls.model)[3])
+
+  result[1, ] <- c(mod$nls.summary[c(1:3,10,9),], logLik(mod$nls.model),
+                   predict(mod$nls.model, newdata=x.range))
+  i <- 1
+  while (i <= B) {
+    if (print.progress) cat("Bootstrap residuals:", i,"out of",B,"\n")
+
+    # Bootstrap residual
+    data.tmp$y <- as.numeric(fit.value + sample(res.value, replace=TRUE))
+
+    m.tmp <- f.quad(d=data.tmp, start=start.value, plus_minus=plus_minus, n.start=n.start)
+
+    if(!is.null(m.tmp$nls.model)) {
+
+      result[i+1, ] <- c(m.tmp$nls.summary[c(1:3,10,9),], logLik(m.tmp$nls.model),
+                         predict(m.tmp$nls.model, newdata=x.range))
+      i <- i + 1
+    }else{ if (print.progress) cat("Not converged! Retry...\n")}
+  }
+  list(result=result, x.range=x.range[,1])
+}
+

R/boot.resid.quad.plateau.R

@@ -0,0 +1,80 @@
+#' Quadratic plateau model estimation by bootstrapping residuals
+#'
+#' \code{boot.resid.quad.plateau} is the core function that implement bootstrapping residuals on quadratic plateau model, which assumes
+#'     y = (a + b * x + c *x^2) * (x <= -0.5*b/c) + (a + -b^2/(4 * c)) * (x > -0.5 * b/c). Note that this functions may take minutes up to days. Parallel computing could be considered when necessary. We suggest users start with a smaller \code{B} and moderate \code{n.start} to see if the bootstrap models can convergence.  In general, increasing \code{n.start} and \code{plus_minus} may help convergence. For rigorous statistical inference, \code{B} should be reach order of thousand.
+#'
+#'
+#' @param mod a full model list, probably from \code{f.quad.plateau()}
+#' @param data data drame with two columns (\code{x} and \code{y})
+#' @param x.range vector of data.frame with one column for range of N rate of interested for prediction interval
+#' @param B bootstrap sample size
+#' @param plus_minus radius of random initial values (default: \code{100})
+#' @param n.start total number of initial points considered (deafult: \code{1000})
+#' @param print.progress logical flag whehter printing progress
+#'
+#' @return \code{boot.resid.quad.plateau} returns a list of two elements:
+#' \code{result}: matrix with B rows and columns containing bootstrap sample for parameter (\code{a,b,c}), optimal N and yield (\code{max_x, max_y}), log-likelihood (\code{logLik}) and N values of interest;
+#' \code{x.range}: range of x considered for prediction interval (same as \code{x.range} in vector form)
+#'
+#' @import stats nls.multstart simpleboot
+#'
+#' @examples
+#'
+#'\donttest{
+#' set.seed(1)
+#' x <- rep(1:300, each=5)
+#' 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)
+#'
+#' ans <- f.quad.plateau(d, start=list(a = 7, b = 0.02, c = 1e-5),
+#'     plus_minus=10, n.start=10, msg=FALSE)
+#'
+#'
+#' boot.resid.quad.plateau(ans, d, x.range=seq(0,280,by=40),
+#'     B=1e2-1, plus_minus = 1e2, n.start=1000, print.progress=TRUE)
+#'
+#' }
+#'
+#'
+#'
+#' @export
+#'
+boot.resid.quad.plateau <- function(mod, data, x.range=data.frame(x=seq(0,280,by=40)), B=1e2-1, plus_minus = 1e2, n.start=5000, print.progress=TRUE) {
+
+  if (class(x.range) == "numeric")  x.range <- data.frame(x=x.range)
+
+  # a, b, c, max x value, max y value, logLik
+  result <- data.frame(matrix(NA, nrow= B + 1, ncol=6 + NROW(x.range)))
+  names(result) <- c("a", "b", "c", "max_x", "max_y", "logLik", paste0("x_", x.range[,1]))
+
+  data <- data.frame(data)
+
+  data.tmp <- data ; names(data.tmp) <- c("x", "y")
+
+
+  fit.value <- fitted(mod$nls.model)
+  res.value <- residuals(mod$nls.model)
+  start.value <- list(a=coef(mod$nls.model)[1], b=coef(mod$nls.model)[2], c=coef(mod$nls.model)[3])
+
+  result[1, ] <- c(mod$nls.summary[c(1:3,10,9),], logLik(mod$nls.model),
+                   predict(mod$nls.model, newdata=x.range))
+  i <- 1
+  while (i <= B) {
+    if (print.progress) cat("Bootstrap residuals:", i,"out of",B,"\n")
+
+    # Bootstrap residual
+    data.tmp$y <- as.numeric(fit.value + sample(res.value, replace=TRUE))
+
+    m.tmp <- f.quad.plateau(d=data.tmp, start=start.value, plus_minus=plus_minus, n.start=n.start)
+
+    if(!is.null(m.tmp$nls.model)) {
+      result[i+1, ] <- c(m.tmp$nls.summary[c(1:3,10,9),], logLik(m.tmp$nls.model),
+                         predict(m.tmp$nls.model, newdata=x.range))
+      i <- i + 1
+    }else{ if (print.progress) cat("Not converged! Retry...\n")}
+  }
+  list(result=result, x.range = x.range[,1])
+}
+

R/compare.two.sample.R

@@ -0,0 +1,29 @@
+#' Two sample bootstrap test for comparing different in \code{sample1} and \code{sample2}, not necessary with same sample size
+#' @param sample1 first sample
+#' @param sample2 second sample
+#' @param fun statistic (univariate) to be compared (default: \code{mean})
+#' @param R number of resample (deafult: \code{1000})
+#' @return \code{compare.two.sample} return a list with two components, namely,
+#' \code{p.value}: two tailed p-value for the bootstrap test
+#' \code{object}: a "\code{simpleboot}" object allowing further analysis using other R packages, such as \code{boot})
+#'
+#' @import stats nls.multstart simpleboot
+#' @examples
+#'
+#' set.seed(1203)
+#' # compare median of two expontential r.v.
+#' compare.two.sample(rexp(100, rate=1), rexp(100, rate=2), fun=median, R=1e3)$p.value
+#'
+#' f.Q1 <- function(x) quantile(x, probs=0.25)
+#' compare.two.sample(rnorm(100, mean=0), rnorm(200, mean=0.5), fun=f.Q1, R=1e3)$p.value
+#'
+#' @export
+#'
+compare.two.sample <- function(sample1, sample2, fun=mean, R=1000) {
+
+  b <- simpleboot::two.boot(sample1, sample2, fun, R = R)
+
+  p.value <- max(min(mean(b$t >=0), mean(b$t <= 0)),1/R)*2
+
+  list(object=b, p.value=p.value)
+}