#' Confidence interval for the distance between two response surface optima (2 regressors)
#'
#' Computes bootstrapped confidence intervals for the mean and median
#' distance between the optima of two response surface models. Models can be thin
#' plate splines or quadratic polynomials
#' \insertCite{DelCastilloCR}{OptimaRegion}.
#'
#' Computes distribution-free bootstrapped confidence intervals on the mean and
#' median distance between the optima of two different responses. The responses can
#' be Thin Plate Spline models or Quadratic polynomial models. Program calls
#' each response, next computes all pairwise distances between points in each CR,
#' and finally bootstraps the distances to compute bca bootstrapped confidence
#' intervals for the mean and median distance.
#'
#' @param X1 nx2 matrix with the values of the 2 regressors (experimental factors)
#' corresponding to the first response. Note: can have replicates. They
#' will be eliminated by the program and the corresponding y-values averaged
#' @param y1 nx1 vector of values for the first response corresponding to X1
#' @param X2 nx2 matrix with the values of the 2 regressors (experimental factors)
#' corresponding to the second response. Note: can have replicates. They
#' will be eliminated by the program and the corresponding y-values averaged
#' @param y2 nx1 vector of values for the second response corresponding to X2
#' @param lambda psmoothing penalty if a TPS model is selected (default=0.04)
#' @param responseType use 'TPS' if fitting thin plate spline responses, 'Quad'
#' if fitting quadratic polynomials
#' @param nosim1and2 number of simulations(default = 200) used to find each of the two
#' confidence regions of optima
#' @param alpha confidence level (0 < alpha < 1; default = 0.05)
#' @param LB1 vector of lower bounds for x (2x1 vector) above which the optimum
#' is sought for the first response
#' @param LB2 vector of lower bounds for x (2x1 vector) above which the optimum
#' is sought for the second response
#' @param UB1 vector of upper bounds for x (2x1 vector) below which the optimum
#' is sought for the first response
#' @param UB2 vector of upper bounds for x (2x1 vector) below which the optimum
#' is sought for the second response
#' @param triangularRegion1 logical: if TRUE it will constrain the maximum points of
#' response 1 to lie inside a triangle defined by the coordinates (0,0), and those in
#' "vertex11", and "vertex21", see below (in addition to being constrained to
#' lie inside the region defined by LB1 and UB1). NOTE: use TRUE when the treatments
#' form a triangular experimental region in shape. If FALSE, optima will only be
#' constrained to lie inside the rectangular region defined by LB1 and UB1.
#' Default is FALSE.
#' @param vertex11 2 times 1 vector with coordinates defining one of the
#' 3 vertices of the triangular region where the first response is being
#' optimized. Must be provided if triangularRegion1 is TRUE
#' (NOTE: vertices numbered clockwise, with vertex0 fixed to (0,0))
#' @param vertex21 2 times 1 vector with coordinates defining a second
#' vertex of a triangular region where the first response is being
#' optimized. Must be provided if triangularRegion1 is TRUE
#' @param triangularRegion2 logical: if TRUE it will constrain the maximum points of
#' response 2 to lie inside a triangle defined by the coordinates (0,0), and those in
#' "vertex12", and "vertex22", see below (in addition to being constrained to
#' lie inside the region defined by LB2 and UB2).NOTE: use TRUE when the treatments
#' form a triangular experimental region in shape. If FALSE, optima will only be
#' constrained to lie inside the rectangular region defined by LB2 and UB2.
#' Default is FALSE.
#' @param vertex12 2 times 1 vector with coordinates defining one of the
#' 3 vertices of the triangular region where the second response is being
#' optimized. Must be provided if triangularRegion2 is TRUE
#' (NOTE: vertices numbered clockwise, with vertex0 fixed to (0,0))
#' @param vertex22 2 times 1 vector with coordinates defining a second
#' vertex of a triangular region where the second response is being
#' optimized. Must be provided if triangularRegion2 is TRUE
#' @param maximization1 logical: if TRUE (default) it maximizes response 1,
#' if FALSE it minimizes it
#' @param maximization2 logical: if TRUE (default) it maximizes response 2,
#' if FALSE it minimizes it
#' @param xlab1and2 text label for x axis in both confidence region plots
#' (default: "Protein eaten (mg)")
#' @param ylab1and2 text label for y axis in both confidence region plots
#' (default: "Carbohydrates eaten (mg)")
#' @param outputPDFFile1 name of the PDF file where the CR plot of the
#' first response is saved (default: "CR_plot.pdf")
#' @param outputPDFFile2 name of the PDF file where the CR plot of the
#' second response is saved (default: "CR_plot.pdf")
#' @param outputOptimaFile1 name of the text file containing the coordinates of
#' all the simulated optima of the first response
#' @param outputOptimaFile2 name of the text file containing the coordinates of
#' all the simulated optima of the second response
#' @return Upon completion, two PDF files with the CR plots and two text files
#' with the coordinates of each set of optima are created, and the
#' function also returns a list consisting of the following 5 components:
#' \describe{
#' \item{dist}{vector of distances between pairs of points taken
#' from each set of optima}
#' \item{mean}{mean of dist}
#' \item{median}{median of dist}
#' \item{ciMean}{95% confidence interval for the mean of dist using
#' bca bootstrapping; it is a vector with 5 columns,
#' containing the signicance level, the next two
#' containing the indices of the order statistics used
#' in the calculations and the final two the calculated
#' endpoints of the CI's.}
#' \item{ciMEdian}{95% confidence interval for the median of dist
#' using bca bootstrapping; it is a vector with 5 columns,
#' containing the signicance level, the next two
#' containing the indices of the order statistics used
#' in the calculations and the final two the calculated
#' endpoints of the CI's.}
#' }
#' @inheritSection OptRegionQuad Author(s)
#' @references{
#' \insertAllCited{}
#' }
#' @examples
#' \dontrun{
#' # Example: two randomly generated data sets, quadratic polynomial responses.
#' X1 <- cbind(runif(100, -2, 2), runif(100, -2, 2))
#' y1 <- as.matrix(72 - 11.78 * X1[, 1] + 0.74 * X1[, 2] - 7.25 * X1[, 1]^2 -
#' 7.55 * X1[, 2]^2 - 4.85 * X1[, 1] * X1[, 2] + rnorm(100, 0, 8))
#' X2 <- cbind(runif(100, -2, 2), runif(100, -2, 2))
#' y2 <- as.matrix(72 - 11.78 * X2[, 1] + 0.74 * X2[, 2] - 7.25 * X2[, 1]^2 -
#' 7.55 * X2[, 2]^2 - 4.85 * X2[, 1] * X2[, 2] + rnorm(100, 0, 8))
#' out <- CRcompare(
#' X1 = X1, y1 = y1, X2 = X2, y2 = y2, responseType = "Quad", nosim1and2 = 200,
#' alpha = 0.05, LB1 = c(-2, -2), UB1 = c(2, 2), LB2 = c(-2, -2), UB2 = c(2, 2)
#' )
#' }
#' @importFrom stats median
#' @export
CRcompare <- function(X1, y1, X2, y2, responseType = "TPS",
lambda = 0.04, nosim1and2 = 200, alpha = 0.05,
LB1, LB2, UB1, UB2,
triangularRegion1 = FALSE, vertex11 = NULL, vertex21 = NULL,
triangularRegion2 = FALSE, vertex12 = NULL, vertex22 = NULL,
maximization1 = TRUE, maximization2 = TRUE,
xlab1and2 = "Protein eaten (mg)",
ylab1and2 = "Carbohydrates eaten (mg)",
outputPDFFile1 = "CR_plot1.pdf", outputOptimaFile1 = "Optima1.txt",
outputPDFFile2 = "CR_plot2.pdf", outputOptimaFile2 = "Optima2.txt") {
# Load required libraries
# t<-.libPaths()
# library("spam", lib.loc=t)
# library("boot", lib.loc=t)
# library("OptimaRegion", lib.loc=t)
if (responseType == "TPS") { # fit TPS models
# Run OptRegionTps (from package OptimaRegion) twice
out1 <- OptRegionTps(
X = X1, y = y1, nosim = nosim1and2, lambda = lambda, alpha = alpha, LB = LB1, UB = UB1, triangularRegion = triangularRegion1, vertex1 = vertex11, vertex2 = vertex21, maximization = maximization1, xlab = xlab1and2, ylab = ylab1and2,
outputPDFFile = outputPDFFile1, outputOptimaFile = outputOptimaFile1
)
out2 <- OptRegionTps(
X = X2, y = y2, nosim = nosim1and2, lambda = lambda, alpha = alpha, LB = LB2, UB = UB2, triangularRegion = triangularRegion2, vertex1 = vertex12, vertex2 = vertex22, maximization = maximization2, xlab = xlab1and2, ylab = ylab1and2,
outputPDFFile = outputPDFFile2, outputOptimaFile = outputOptimaFile2
)
} else if (responseType == "Quad") # fit quadratic polynomails instead
{
# Run OptRegionQuad (from package OptimaRegion) twice
out1 <- OptRegionQuad(
X = X1, y = y1, nosim = nosim1and2, alpha = alpha, LB = LB1, UB = UB1, triangularRegion = triangularRegion1, vertex1 = vertex11, vertex2 = vertex21, maximization = maximization1, xlab = xlab1and2, ylab = ylab1and2,
outputPDFFile = outputPDFFile1
)
out2 <- OptRegionQuad(
X = X2, y = y2, nosim = nosim1and2, alpha = alpha, LB = LB2, UB = UB2, triangularRegion = triangularRegion2, vertex1 = vertex12, vertex2 = vertex22, maximization = maximization2, xlab = xlab1and2, ylab = ylab1and2,
outputPDFFile = outputPDFFile2
)
}
# Form as many pairs of points as in the smallest set and compute their euclidean distances
Dmatrix <- spam::nearest.dist(x = out1$xin, y = out2$xin, upper = TRUE, delta = 1e20)
dist <- spam::diag(Dmatrix)
mean <- mean(dist)
median <- median(dist)
# Bootstrap the distances between points in each CR
bmean <- boot::boot(dist, sampleMeans, R = 1000)
bmedian <- boot::boot(dist, sampleMedians, R = 1000)
ciMean <- boot::boot.ci(bmean, type = c("bca"))
ciMedian <- boot::boot.ci(bmedian, type = c("bca"))
return(list(dist = dist, mean = mean, median = median, ciMean = ciMean$bca, ciMedian = ciMedian$bca))
}
sampleMeans <- function(W, i) {
mean(W[i])
}
sampleMedians <- function(W, i) {
median(W[i])
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.