Nothing
#' Confidence or prediction interval limit
#'
#' The function \code{get_intvl_limit()} calculates the upper or lower
#' confidence or prediction interval limit(s) at a given value of \eqn{x}.
#'
#' @param x_new A numeric value of \eqn{x} at which the upper or lower
#' confidence or prediction interval limit(s) should be calculated.
#' @param model A linear model object of type \sQuote{\code{lm}}.
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_intvl_limit()} calculates the upper or lower
#' confidence or prediction interval limit(s) for the linear regression model
#' provided by \code{model}.
#'
#' @return A numeric value or, if \code{model} contains a categorical variable,
#' a numeric vector of the predicted upper or lower confidence or prediction
#' interval limit(s) at a given value of \eqn{x}.
#'
#' @seealso \code{\link[stats]{predict.lm}}.
#'
#' @importFrom stats predict
#'
#' @keywords internal
#' @noRd
get_intvl_limit <- function(x_new, model, alpha = 0.05, ivl = "confidence",
ivl_type = "one.sided", ivl_side = "lower") {
if (!is.numeric(x_new) && !is.na(x_new)) {
stop("x_new must be a numeric value.")
}
if (!inherits(model, "lm")) {
stop("Please provide a model of type \"lm\".")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (!(ivl %in% c("confidence", "prediction", "fitted.line"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\" or ",
"\"fitted.line\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper"))) {
stop("Please specify ivl_side either as \"lower\" or \"upper\".")
}
# If x_new is NA, terminate with NA
if (is.na(x_new)) {
return(NA)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Note: since the predict.lm() function from the 'stats' package always
# calculates two-sided limits the value of alpha must be doubled in case that
# the ivl_type is "one-sided".
if (ivl_type == "one.sided") {
alpha <- alpha * 2
}
variable_names <- names(attr(model$terms, which = "dataClasses"))
is_factor <- sum(attr(model$terms, which = "dataClasses") %in% "factor")
if (is_factor == 0 && length(variable_names) == 2) {
l_newdata <- list(x_new)
names(l_newdata) <- variable_names[2]
if (ivl == "fitted.line") {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha, interval = "confidence"))[, "fit"]
} else {
switch(ivl_side,
"lower" = {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha,
interval = ivl))[, "lwr"]
},
"upper" = {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha,
interval = ivl))[, "upr"]
})
}
} else {
grouping_variables <-
variable_names[attr(model$terms, which = "dataClasses") %in% "factor"]
n_levels <- vapply(grouping_variables,
function(x) {
nlevels(model$model[, grouping_variables])
},
numeric(1))
n_levels <- prod(n_levels)
l_newdata <-
vapply(grouping_variables,
function(x) {
list(rep(levels(model$model[, grouping_variables]),
each = n_levels /
nlevels(model$model[, grouping_variables]) *
length(x_new)))
},
list(1))
l_newdata[[length(l_newdata) + 1]] <- rep(x_new, n_levels)
names(l_newdata) <- c(grouping_variables,
variable_names[!(variable_names %in%
grouping_variables)][-1])
if (ivl == "fitted.line") {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha, interval = "confidence"))[, "fit"]
} else {
switch(ivl_side,
"lower" = {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha,
interval = ivl))[, "lwr"]
},
"upper" = {
res <- do.call("predict",
list(object = model, newdata = l_newdata,
level = 1 - alpha,
interval = ivl))[, "upr"]
})
}
names(res) <- l_newdata[[grouping_variables]]
}
return(res)
}
#' Determine distance of lines
#'
#' The function \code{get_distance()} calculates the distance between two
#' lines at a given value of \eqn{x}.
#'
#' @param x_new A numeric value of \eqn{x} for which the distance between two
#' lines is sought, e.g. the distance of upper or lower confidence/prediction
#' interval limits from the upper or lower specification or expiry limits,
#' respectively.
#' @param sl A numeric variable that specifies the \dQuote{specification limit}
#' (SL) or, for determinations according to ARGPM guidance \dQuote{Stability
#' testing for prescription medicines}, the \dQuote{expiry limit} (EL). The
#' EL is defined as the intercept \eqn{\pm} the difference between the
#' specification limit and the release limit (RL). If it is the upper limit
#' which is the relevant limit, it is added (\code{+}) to the intercept,
#' otherwise it is subtracted (\code{-}) from the intercept.
#' @param mode A character string of either \code{"minimal"} or \code{"all"},
#' that specifies if only the minimal distance of a factor regression model
#' is returned or if the distances of all lines belonging to the different
#' factor levels is returned. The default is \code{"minimal"}.
#' @inheritParams get_intvl_limit
#' @inheritParams expirest_osle
#'
#' @details The function \code{find_get_distance()} estimates the distance
#' between the upper or lower confidence or prediction interval and the upper
#' or lower acceptance criterion (i.e. the specification or the expiry limit)
#' at a certain value of \code{x_new}. The confidence or prediction interval
#' is calculated for the linear regression model provided by \code{model}.
#' Recommendations on how to estimate shelf life or expiry can be found in the
#' corresponding section below.
#'
#' @section How to estimate shelf life or expiry:
#' ICH Q1E recommends that \dQuote{\emph{For an attribute known to decrease with
#' time, the lower one-sided 95 percent confidence limit should be compared
#' to the acceptance criterion. For an attribute known to increase with time,
#' the upper one-sided 95 percent confidence limit should be compared to the
#' acceptance criterion. For an attribute that can either increase or decrease,
#' or whose direction of change is not known, two-sided 95 percent confidence
#' limits should be calculated and compared to the upper and lower acceptance
#' criteria.}} Since attributes often either decrease or increase, the default
#' for \code{ivl_type} is \code{one.sided}.
#'
#' According to the ARGPM guidance \dQuote{Stability testing for prescription
#' medicines}, the shelf life or expiry limit is estimated as the point where
#' the upper or lower limit of the 95\% confidence interval of the linear model
#' fitted to the data intersects the worst case scenario limit. The worst case
#' scenario limit is obtained by adding/subtracting the absolute difference of
#' specification limit and release limit to/from the common intercept of the
#' test batches or the intercept of the worst performing batch.
#'
#' @return A numeric value representing the distance of the respective lines
#' is returned.
#'
#' @inherit expirest_wisle references
#'
#' @seealso \code{\link{get_intvl_limit}}, \code{\link{find_poi}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}},
#' \code{\link[stats]{lm}}.
#'
#' @keywords internal
#' @noRd
get_distance <- function(x_new, model, sl, mode = "minimal", alpha = 0.05,
ivl = "confidence", ivl_type = "one.sided",
ivl_side = "lower") {
if (!is.numeric(x_new)) {
stop("x_new must be a numeric value.")
}
if (!inherits(model, "lm")) {
stop("Please provide a model of type \"lm\".")
}
if (!is.numeric(sl) || length(sl) > 1) {
stop("sl must be a numeric value of length 1.")
}
if (!(mode %in% c("minimal", "all"))) {
stop("Please specify mode either as \"minimal\" or \"all\".")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (!(ivl %in% c("confidence", "prediction", "fitted.line"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\" or ",
"\"fitted.line\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper"))) {
stop("Please specify ivl_side either as \"lower\" or \"upper\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
pred_lim <-
get_intvl_limit(x_new = x_new, model = model, alpha = alpha,
ivl = ivl, ivl_type = ivl_type, ivl_side = ivl_side)
if (mode == "minimal") {
switch(ivl_side,
"lower" = {
res <- sl - min(pred_lim, na.rm = TRUE)
},
"upper" = {
res <- sl - max(pred_lim, na.rm = TRUE)
})
} else {
res <- sl - pred_lim
}
return(res)
}
#' Find the point of intersection (POI)
#'
#' The function \code{find_poi()} determines the point where the distance
#' between two lines is minimal, e.g. the distance between a specification or
#' expiry limit and a confidence or prediction interval, or in other words
#' the point of intersection (POI). The estimation is done by aid of
#' \code{\link[stats]{uniroot}()} from the \sQuote{\code{stats}} package.
#'
#' @param srch_range A vector of length \code{2} that specifies the endpoints
#' of the (time) range within which the minimum distance is expected.
#' @param sl A numeric variable that specifies the \dQuote{specification limit}
#' (SL). Another kind of acceptance criterion may be regarded as SL.
#' @param ... Additional named or unnamed arguments passed on to
#' \code{\link[stats]{uniroot}()}.
#' @inheritParams expirest_osle
#' @inheritParams get_distance
#'
#' @details The function \code{find_poi()} (find the \dQuote{point of
#' intersection}) estimates the value of \eqn{x} (e.g. the time) where the
#' difference between the upper or lower confidence or prediction interval and
#' the upper or lower acceptance criterion (e.g. the specification or the
#' expiry limit) is minimal, or in other words where the corresponding lines
#' intersect each other. Confidence or prediction intervals are calculated
#' for the \code{model} provided. The POI is determined by aid of the
#' \code{\link[stats]{uniroot}()} function from the \sQuote{\code{stats}}
#' package. The distance between the two lines of interest is calculated using
#' the function \code{\link{get_distance}()}, and it is this distance which
#' \code{uniroot()} tries to minimise. Recommendations on how to estimate shelf
#' life or expiry can be found in the corresponding section below.
#'
#' @inheritSection get_distance How to estimate shelf life or expiry
#'
#' @return A numeric value representing the value of \eqn{x} where the distance
#' between the two lines of interest is minimal is returned.
#'
#' @inherit expirest_wisle references
#'
#' @seealso \code{\link{get_distance}}, \code{\link{get_poi_list}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}},
#' \code{\link[stats]{uniroot}}.
#'
#' @importFrom stats uniroot
#' @importFrom stats setNames
#'
#' @keywords internal
#' @noRd
find_poi <- function(srch_range, model, sl, mode = "minimal", alpha = 0.05,
ivl = "confidence", ivl_type = "one.sided",
ivl_side = "lower", ...) {
if (!is.numeric(srch_range) || length(srch_range) != 2) {
stop("The parameter srch_range must be a vector of length 2.")
}
if (!inherits(model, "lm")) {
stop("Please provide a model of type \"lm\".")
}
if (!is.numeric(sl) || length(sl) > 1) {
stop("sl must be a numeric value of length 1.")
}
if (!(mode %in% c("minimal", "all"))) {
stop("Please specify mode either as \"minimal\" or \"all\".")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (!(ivl %in% c("confidence", "prediction", "fitted.line"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\" or ",
"\"fitted.line\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper"))) {
stop("Please specify ivl_side either as \"lower\" or \"upper\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (mode == "minimal") {
tmp <- try_get_model(
uniroot(f = get_distance, interval = srch_range, model = model,
sl = sl, mode = mode, alpha = alpha, ivl = ivl,
ivl_type = ivl_type, ivl_side = ivl_side, ...)[["root"]]
)
if (is.null(tmp[["Error"]])) {
res <- tmp[["Model"]]
} else {
stop("Error in uniroot (find_poi)!: ", tmp[["Error"]])
}
} else {
t_dist <-
get_distance(x_new = 0, model = model, sl = sl, mode = mode,
alpha = alpha, ivl = ivl, ivl_type = ivl_type,
ivl_side = ivl_side)
res <- setNames(rep(NA, length(t_dist)), names(t_dist))
for (i in seq_along(t_dist)) {
tmp <- try_get_model(
uniroot(f = function(x, ii, ...) get_distance(x_new = x, ...)[ii],
interval = srch_range, ii = i, model = model, sl = sl,
mode = "all", alpha = alpha, ivl = ivl, ivl_type = ivl_type,
ivl_side = ivl_side, ...)[["root"]]
)
if (is.null(tmp[["Error"]])) {
res[i] <- tmp[["Model"]]
}
}
}
return(res)
}
#' List of points of intersection
#'
#' The function \code{get_poi_list()} prepares a list of points of intersection
#' (POI) for multiple regression models using the \code{find_poi()} function.
#'
#' @param data The data frame that was used for fitting the models of parameter
#' \code{model_list}.
#' @param batch_vbl A character string that specifies the column in \code{data}
#' with the grouping information (i.e. a categorical variable) for the
#' differentiation of the observations of the different batches.
#' @param model_list A list of regression models of different type. Usually,
#' it is a list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids}, where the first three elements contain
#' \sQuote{\code{lm}} objects of the \dQuote{common intercept / common slope}
#' (\code{cics}), \dQuote{different intercept / common slope} (\code{dics})
#' and \dQuote{different intercept / different slope} (\code{dids.pmse}) type.
#' The fourth element with the label \code{dids} is usually a list of the
#' \sQuote{\code{lm}} objects that is obtained from fitting a regression
#' model to the data of each level of the categorical variable separately.
#' The \code{dids.pmse} model differs from the \code{dids} model in that it
#' is a model with the categorical variable as a fixed main effect and with
#' an interaction term of the categorical variable with the time variable,
#' i.e. a model where the mean square error is pooled across batches (thus
#' the \dQuote{pmse} suffix meaning \dQuote{pooled mean square error}). The
#' \code{cics}, \code{dics} and \code{dids.pmse} elements are \code{NA} if
#' data of only a single batch is available.
#' @param sl A numeric variable that specifies the \dQuote{specification limit}
#' (SL). Another kind of acceptance criterion may be regarded as SL.
#' @inheritParams find_poi
#'
#' @details The function \code{get_poi_list()} applies the \code{find_poi()}
#' function (find the \dQuote{point of intersection}) on all the models that
#' are provided.
#'
#' @return A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} is returned. Each of them contains a named
#' vector of the POI values estimated for each batch and named accordingly.
#'
#' @seealso \code{\link{get_model_list}}, \code{\link{find_poi}},
#' \code{\link{get_distance}}, \code{\link{get_osle_poi_list}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}},
#' \code{\link[stats]{uniroot}}.
#'
#' @keywords internal
#' @noRd
get_poi_list <- function(data, batch_vbl, model_list, sl, srch_range,
mode = "minimal", alpha = 0.05, ivl = "confidence",
ivl_type = "one.sided", ivl_side = "lower", ...) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
if (!is.list(model_list)) {
stop("The parameter model_list must be a list.")
}
if (sum(names(model_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The parameter model_list must have four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.numeric(sl) || length(sl) > 1) {
stop("The parameter sl must be a numeric value of length 1.")
}
if (!is.numeric(srch_range) || length(srch_range) != 2) {
stop("The parameter srch_range must be a vector of length 2.")
}
if (!(mode %in% c("minimal", "all"))) {
stop("Please specify mode either as \"minimal\" or \"all\".")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (!(ivl %in% c("confidence", "prediction"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper", "both"))) {
stop("Please specify ivl_side either as \"lower\", \"upper\" or \"both\".")
}
if (ivl_side == "both" && length(sl) == 1) {
stop("Since ivl_side = \"both\", a specification with two sides is ",
"expected. Only one side has been specified, though, i.e. ",
"sl = ", sl, ".\nPlease provide a specification with two sides.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
d_dat <- droplevels(data)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of points of intersection (POIs)
l_poi <- lapply(model_list, function(x) NA)
l_poi[["dids"]] <-
vapply(model_list[["dids"]],
function(x) {
tmp <- try_get_model(
find_poi(srch_range = srch_range, model = x, sl = sl,
mode = "minimal", alpha = alpha, ivl_type = ivl_type,
ivl_side = ivl_side, ivl = ivl)
)
ifelse(is.null(tmp[["Error"]]), tmp[["Model"]], NA)
},
numeric(1))
if (nlevels(d_dat[, batch_vbl]) > 1) {
for (variety in names(l_poi)[names(l_poi) != "dids"]) {
if (variety == "cics") {
tmp <- try_get_model(
find_poi(srch_range = srch_range, model = model_list[[variety]],
sl = sl, mode = "minimal", alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl))
}
if (variety %in% c("dics", "dids.pmse")) {
tmp <- try_get_model(
find_poi(srch_range = srch_range, model = model_list[[variety]],
sl = sl, mode = "all", alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl))
}
if (is.null(tmp[["Error"]])) {
l_poi[[variety]] <- tmp[["Model"]]
}
}
}
return(l_poi)
}
#' Getting the intercept(s) of a linear model
#'
#' The function \code{get_icpt()} determines the intercept(s) of the provided
#' model.
#'
#' @param response_vbl A character string that specifies the response variable
#' name that must be a column of the data frame that was used for model
#' fitting.
#' @param time_vbl A character string that specifies the time variable name
#' that must be a column of data frame that was used for model fitting.
#' @param batch_vbl A character string that specifies the column of the data
#' frame that was used for model fitting with the grouping information (i.e.
#' a categorical variable) for the differentiation of the observations from
#' the different batches.
#' @inheritParams get_intvl_limit
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_icpt()} determines the intercept(s) of the
#' model that has been handed over via the \code{model} parameter.
#'
#' @return A list with a single element containing the numeric value or a
#' numeric vector of the intercept(s) or, if the data have been transformed,
#' a list with an additional element that contains the numeric value or
#' numeric vector on the original scale is returned.
#'
#' @seealso \code{\link{get_icpt_list}}, \code{\link{expirest_osle}},
#' \code{\link{expirest_wisle}}, \code{\link[stats]{lm}}.
#'
#' @importFrom stats coef
#' @importFrom stats formula
#'
#' @keywords internal
#' @noRd
get_icpt <- function(model, response_vbl, time_vbl, batch_vbl,
xform = c("no", "no"), shift = c(0, 0)) {
if (!inherits(model, "lm")) {
stop("Please provide a model of type \"lm\".")
}
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (length(grep(response_vbl, as.character(formula(model$terms)))) == 0) {
stop("The parameter response_vbl was not found in the provided model.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (length(grep(time_vbl, as.character(formula(model$terms)))) == 0) {
stop("The parameter time_vbl was not found in the provided model.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (length(attr(model$terms, "order")) > 1) {
if (length(grep(batch_vbl, as.character(formula(model$terms)))) == 0) {
stop("The parameter batch_vbl was not found in the provided model.")
}
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of intercept(s)
if (length(coef(model)) == 2) {
intercept <- coef(model)["(Intercept)"]
} else {
t_coef <- coef(model)
t_batch_levels <- levels(model[["model"]][, batch_vbl])
t_batch_id <- grep(batch_vbl, names(t_coef))
t_time_id <- grep(time_vbl, names(t_coef))
# Calculate the intercepts
intercept <- c(t_coef["(Intercept)"], t_coef["(Intercept)"] +
t_coef[t_batch_id[t_batch_id < t_time_id[1]]])
# Rename the vector
t_names <- names(intercept)
t_names <- sub(batch_vbl, "", t_names)
t_names[which(!(t_names %in% t_batch_levels))] <-
t_batch_levels[which(!(t_names %in% t_batch_levels))]
names(intercept) <- t_names
}
# ---------
# Depending on the transformation of the response variable the intercept(s)
# has (have) to be back-transformed to the original scale.
if (xform[2] != "no") {
switch(xform[2],
"log" = {
icpt_orig <- exp(intercept) - shift[2]
},
"sqrt" = {
icpt_orig <- intercept^2 - shift[2]
},
"sq" = {
icpt_orig <- rep(NA, length(intercept))
names(icpt_orig) <- names(intercept)
ok <- unname(which(intercept >= 0))
if (length(ok) > 0) {
icpt_orig[ok] <- sqrt(intercept[ok]) - shift[2]
}
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Return of the intercept(s)
if (xform[2] == "no") {
l_res <- list(icpt = intercept)
} else {
l_res <- list(icpt = intercept,
icpt.orig = icpt_orig)
}
return(l_res)
}
#' List of intercepts
#'
#' The function \code{get_icpt_list()} prepares a list of the intercepts
#' of the regression models fitted to the data.
#'
#' @param data The data frame that was used for fitting the models of parameter
#' \code{model_list}.
#' @inheritParams expirest_osle
#' @inheritParams get_poi_list
#'
#' @details The function \code{get_icpt_list()} extracts the intercepts of
#' the various regression models (fitted by aid of the \code{lm()} function)
#' that are passed in via the \code{model_list} parameter.
#'
#' @return A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} is returned. Each of them contains a list
#' element named \code{icpt} with a vector of the intercepts. If the data
#' have been transformed, each of the primary list elements contains a further
#' list element called \code{icpt.orig} with a numeric vector of the intercepts
#' on the original scale.
#'
#' @seealso \code{\link{get_model_list}}, \code{\link{get_icpt}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_icpt_list <- function(data, response_vbl, time_vbl, batch_vbl, model_list,
xform = c("no", "no"), shift = c(0, 0)) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (!(response_vbl %in% colnames(data))) {
stop("The response_vbl was not found in the provided data frame.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (!(time_vbl %in% colnames(data))) {
stop("The time_vbl was not found in the provided data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
if (!is.list(model_list)) {
stop("The parameter model_list must be a list.")
}
if (sum(names(model_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The parameter model_list must have four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
d_dat <- droplevels(data)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of the intercepts
if (nlevels(d_dat[, batch_vbl]) > 1) {
l_icpt <- vapply(model_list[c("cics", "dics", "dids.pmse")], function(x) {
list(get_icpt(model = x, response_vbl = response_vbl,
time_vbl = time_vbl, batch_vbl = batch_vbl,
xform = xform, shift = shift))
},
list(1))
} else {
l_icpt <- list(cics = NA, dics = NA, dids.pmse = NA)
}
tmp <- vapply(model_list[["dids"]], function(x) {
list(get_icpt(model = x, response_vbl = response_vbl,
time_vbl = time_vbl, batch_vbl = batch_vbl,
xform = xform, shift = shift))
},
list(1))
tmp <- unlist(tmp)
names(tmp) <- sub("\\.\\(Intercept\\)", "", names(tmp))
if (xform[2] == "no") {
names(tmp) <- sub("\\.icpt", "", names(tmp))
l_icpt <- c(l_icpt, list(dids = list(icpt = tmp)))
} else {
t_i_orig <- grep("orig", names(tmp))
names(tmp) <- sub("\\.icpt", "", names(tmp))
names(tmp) <- sub("\\.orig", "", names(tmp))
l_icpt <- c(l_icpt, list(dids = list(icpt = tmp[-t_i_orig],
icpt.orig = tmp[t_i_orig])))
}
return(l_icpt)
}
#' Determination of the \dQuote{worst case scenario} (wcs) limit
#'
#' The function \code{get_wcs_limit()} calculates \dQuote{worst case scenario}
#' (wcs) limit following the ARGPM Guidance \dQuote{Stability testing for
#' prescription medicines}.
#'
#' @param rl A numeric value that specifies the release specification limit(s),
#' on the same scale as \code{sl} and \code{intercept}.
#' @param sl A numeric value that specifies the specification limit, on the same
#' scale as \code{rl} and \code{intercept}.
#' @param intercept A numeric value representing the intercept of a linear
#' regression model fitted to sample data, on the same scale as \code{rl}
#' and \code{sl}.
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_wcs_limit()} determines the \dQuote{worst
#' case scenario} (wcs) limit as is proposed by the Australian Regulatory
#' Guidelines for Prescription Medicines (ARGPM) guidance \dQuote{Stability
#' testing for prescription medicines}. According to this guideline, the shelf
#' life or expiry limit is estimated as the point where the upper or lower
#' limit of the 95\% confidence interval of the linear model fitted to the
#' data intersects the wcs limit. The wcs limit is obtained by
#' adding/subtracting the absolute difference of specification limit and
#' release limit to/from the common intercept of the test batches or the
#' intercept of the worst performing batch.
#'
#' If data have been linearised by transformation, all elements, i.e. \code{rl},
#' \code{sl} and \code{intercept} must be on the same, i.e. transformed, scale.
#' The results are returned on the transformed scale and on the original scale.
#'
#' @return A list with the following elements is returned:
#' \item{xform}{A vector of two character strings that specifies the
#' transformation of the response and the time variable.}
#' \item{shift}{A vector of two values which has been added to the values of
#' the transformed \eqn{x} and/or \eqn{y} variables (specified via the
#' \code{xform} parameter).}
#' \item{delta.lim}{A numeric value or a numeric vector of the absolute
#' difference(s) between \code{rl} and \code{sl}, if \code{xform[2] != "no"}
#' on the transformed scale.}
#' \item{delta.lim.orig}{A numeric value or a numeric vector of the absolute
#' difference(s) between \code{rl} and \code{sl} on the original scale.}
#' \item{wcs.lim}{A numeric value or a numeric vector of the worst case
#' scenario (wcs) limit(s), if \code{xform[2] != "no"} on the transformed
#' scale.}
#' \item{wcs.lim.orig}{A numeric value or a numeric vector of the worst case
#' scenario (wcs) limit(s) on the original scale.}
#'
#' @references
#' Therapeutic Goods Administration (TGA) of the Department of Health of the
#' Australian Government, Australian Regulatory Guidelines for Prescription
#' Medicines (ARGPM), Stability testing for prescription medicines,
#' Version 1.1, March 2017
#'
#' @seealso \code{\link{extract_from_ll_wcsl}},
#' \code{\link{get_wisle_poi_list}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_wcs_limit <- function(rl, sl, intercept, xform = c("no", "no"),
shift = c(0, 0), ivl_side = "lower") {
if (!is.numeric(rl) || length(rl) > 1) {
stop("The parameter rl must be a numeric value of length 1.")
}
if (!is.numeric(sl) || length(sl) > 1) {
stop("The parameter sl must be a numeric value of length 1.")
}
if (!is.numeric(intercept)) {
stop("The parameter intercept must be a numeric value of length 1.")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
if (!(ivl_side %in% c("lower", "upper"))) {
stop("Please specify ivl_side either as \"lower\" or \"upper\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of worst case scenario (wcs) limit(s)
switch(ivl_side,
"lower" = {
delta_lim <- rl - sl
wcs_lim <- intercept - delta_lim
},
"upper" = {
delta_lim <- sl - rl
wcs_lim <- intercept + delta_lim
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Based on the transformation of the response variable delta_lim and
# wcs_lim have to be back-transformed to the original scale.
if (xform[2] != "no") {
switch(xform[2],
"log" = {
delta_lim_orig <- exp(delta_lim) - shift[2]
wcs_lim_orig <- exp(wcs_lim) - shift[2]
},
"sqrt" = {
delta_lim_orig <- delta_lim^2 - shift[2]
wcs_lim_orig <- wcs_lim^2 - shift[2]
},
"sq" = {
delta_lim_orig <- sqrt(delta_lim) - shift[2]
if (wcs_lim < 0) {
wcs_lim_orig <- NA
} else {
wcs_lim_orig <- sqrt(wcs_lim) - shift[2]
}
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Collect and return the data
if (xform[2] == "no") {
l_res <- list(xform = xform,
shift = shift,
delta.lim = delta_lim,
delta.lim.orig = delta_lim,
wcs.lim = wcs_lim,
wcs.lim.orig = wcs_lim)
} else {
l_res <- list(xform = xform,
shift = shift,
delta.lim = delta_lim,
delta.lim.orig = delta_lim_orig,
wcs.lim = wcs_lim,
wcs.lim.orig = wcs_lim_orig)
}
return(l_res)
}
#' Get intercepts of the worst case batches
#'
#' The function \code{get_wc_icpt()} prepares a vector of the intercepts
#' of the worst case batches of all the regression models fitted to the data.
#'
#' @param icpt_list A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} with the intercepts of each linear
#' regression model and batch. The \code{cics}, \code{dics} and
#' \code{dids.pmse} elements are \code{NA} if data of only a single batch
#' is available.
#' @param poi_list A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} with the points of intersection (POI)
#' of each linear regression model and batch. The \code{cics}, \code{dics}
#' and \code{dids.pmse} elements are \code{NA} if data of only a single
#' batch is available.
#' @param wc_batch A numeric vector of the indices of the worst case batches
#' of each model type, i.e. a vector of four elements named \code{cics},
#' \code{dics}, \code{dids.pmse} and \code{dids}. The \code{cics} element
#' is \code{NA} because in the \dQuote{common intercept / common slope} model
#' the data from different batches is pooled.
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_wc_icpt()} extracts the worst case batches
#' from the list of intercepts (\code{icpt_list}), given that the estimation
#' of the corresponding POIs (\code{poi_list}) was successful.
#'
#' @return A named vector of the intercepts of the worst case batches is
#' returned.
#'
#' @seealso \code{\link{get_icpt_list}}, \code{\link{get_poi_list}},
#' \code{\link{get_osle_poi_list}}, \code{\link{get_wisle_poi_list}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @importFrom stats setNames
#'
#' @keywords internal
#' @noRd
get_wc_icpt <- function(data, batch_vbl, icpt_list, poi_list, wc_batch, xform) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
if (!is.list(icpt_list)) {
stop("The parameter icpt_list must be a list.")
}
if (sum(names(icpt_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The parameter icpt_list must have four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.list(poi_list)) {
stop("The parameter poi_list must be a list.")
}
if (sum(names(poi_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The parameter poi_list must have four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.numeric(wc_batch)) {
stop("The parameter wc_batch must be a numeric vector.")
}
if (sum(names(wc_batch) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The parameter wc_batch must have four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
d_dat <- droplevels(data)
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (nlevels(d_dat[, batch_vbl]) > 1) {
wc_icpt <-
vapply(names(icpt_list), function(nn) {
if (nn == "cics") {
if (xform[2] != "no") {
unname(icpt_list[[nn]][["icpt.orig"]])
} else {
unname(icpt_list[[nn]][["icpt"]])
}
} else {
if (xform[2] != "no") {
icpt_list[[nn]][["icpt.orig"]][wc_batch[nn]]
} else {
icpt_list[[nn]][["icpt"]][wc_batch[nn]]
}
}
},
numeric(1))
} else {
wc_icpt <-
setNames(rep(NA, length(wc_batch)), names(wc_batch))
if (!any(is.na(poi_list[["dids"]]))) {
if (xform[2] != "no") {
wc_icpt["dids"] <-
icpt_list[["dids"]][["icpt.orig"]][wc_batch["dids"]]
} else {
wc_icpt["dids"] <-
icpt_list[["dids"]][["icpt"]][wc_batch["dids"]]
}
}
}
return(wc_icpt)
}
#' Transformation of variables
#'
#' The function \code{get_xformed_variables()} transforms the variables as
#' needed.
#'
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_xformed_variables()} transforms the variables
#' (\code{response_vbl} and/or \code{time_vbl}) as specified by the parameters
#' \code{xform} and \code{shift}.
#'
#' @return The provided data frame with (a) new column(s) of the transformed
#' variable(s).
#'
#' @seealso \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_xformed_variables <- function(data, response_vbl, time_vbl,
xform = c("no", "no"), shift = c(0, 0)) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (!(response_vbl %in% colnames(data))) {
stop("The response_vbl was not found in the provided data frame.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (!(time_vbl %in% colnames(data))) {
stop("The time_vbl was not found in the provided data frame.")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Linearisation of data by variable transformation
# Transformations:
# log: natural logarithm of the variable
# sqrt: square root of the variable variable
# sq: square of the variable
#
# Note: The log and sqrt transformations include adding the value defined by
# the shift parameter before performing the transformation.
d_dat <- droplevels(data)
if (xform[1] != "no") {
switch(xform[1],
"log" = {
trfmd_time_vbl <- paste("log", time_vbl, sep = ".")
d_dat <- cbind(d_dat, log(d_dat[, time_vbl] + shift[1]))
colnames(d_dat)[ncol(d_dat)] <- trfmd_time_vbl
},
"sqrt" = {
trfmd_time_vbl <- paste("sqrt", time_vbl, sep = ".")
d_dat <- cbind(d_dat, sqrt(d_dat[, time_vbl] + shift[1]))
colnames(d_dat)[ncol(d_dat)] <- trfmd_time_vbl
},
"sq" = {
trfmd_time_vbl <- paste("sq", time_vbl, sep = ".")
d_dat <- cbind(d_dat, (d_dat[, time_vbl] + shift[1])^2)
colnames(d_dat)[ncol(d_dat)] <- trfmd_time_vbl
})
}
if (xform[2] != "no") {
switch(xform[2],
"log" = {
trfmd_response_vbl <- paste("log", response_vbl, sep = ".")
d_dat <- cbind(d_dat, log(d_dat[, response_vbl] + shift[2]))
colnames(d_dat)[ncol(d_dat)] <- trfmd_response_vbl
},
"sqrt" = {
trfmd_response_vbl <- paste("sqrt", response_vbl, sep = ".")
d_dat <- cbind(d_dat, sqrt(d_dat[, response_vbl] + shift[2]))
colnames(d_dat)[ncol(d_dat)] <- trfmd_response_vbl
},
"sq" = {
trfmd_response_vbl <- paste("sq", response_vbl, sep = ".")
d_dat <- cbind(d_dat, (d_dat[, response_vbl] + shift[2])^2)
colnames(d_dat)[ncol(d_dat)] <- trfmd_response_vbl
})
}
return(d_dat)
}
#' Listing of variable names
#'
#' The function \code{get_variable_list()} makes a list of the (original) and,
#' if applicable, the transformed variable names.
#'
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_variable_list()} makes a list of the
#' variable names. If data have been transformed, the list comprises the
#' original variable name(s) (with suffix .orig in the corresponding list
#' element names) and the transformed variable name(s).
#'
#' @return A list with the variable names. If the data have been transformed,
#' the list element names of the original variables have the suffix
#' \code{".orig"}.
#'
#' @seealso \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_variable_list <- function(response_vbl, time_vbl, batch_vbl,
xform = c("no", "no")) {
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (xform[1] != "no") {
switch(xform[1],
"log" = {
new_time_vbl <- paste("log", time_vbl, sep = ".")
},
"sqrt" = {
new_time_vbl <- paste("sqrt", time_vbl, sep = ".")
},
"sq" = {
new_time_vbl <- paste("sq", time_vbl, sep = ".")
})
}
if (xform[2] != "no") {
switch(xform[2],
"log" = {
new_response_vbl <- paste("log", response_vbl, sep = ".")
},
"sqrt" = {
new_response_vbl <- paste("sqrt", response_vbl, sep = ".")
},
"sq" = {
new_response_vbl <- paste("sq", response_vbl, sep = ".")
})
}
if (xform[1] == "no" && xform[2] == "no") {
l_variables <- list(batch = batch_vbl,
response = response_vbl,
time = time_vbl)
}
if (xform[1] != "no" && xform[2] != "no") {
l_variables <- list(batch = batch_vbl,
response.orig = response_vbl,
response = new_response_vbl,
time.orig = time_vbl,
time = new_time_vbl)
} else {
if (xform[1] != "no") {
l_variables <- list(batch = batch_vbl,
response = response_vbl,
time.orig = time_vbl,
time = new_time_vbl)
}
if (xform[2] != "no") {
l_variables <- list(batch = batch_vbl,
response.orig = response_vbl,
response = new_response_vbl,
time = time_vbl)
}
}
return(l_variables)
}
#' Adjustment of limits
#'
#' The function \code{set_limits()} adjusts the limits according to the number
#' of relevant decimal places and according to the transformation requirement.
#'
#' @inheritParams expirest_osle
#'
#' @details The function \code{set_limits()} adjusts the limits according to
#' \code{rl_sf} and \code{sl_sf} and, if necessary, transforms the limits
#' (\code{rl} and \code{sl}) as specified by the parameters \code{xform} and
#' \code{shift}.
#'
#' @return A list with the following elements is returned:
#' \item{sf.option}{A character string that specifies the option concerning the
#' significant figures.}
#' \item{xform}{A vector of two character strings that specifies the
#' transformation of the response and the time variable.}
#' \item{shift}{A vector of two values to be added to the values of the
#' transformed \eqn{x} and/or \eqn{y} variables (specified via the
#' \code{xform} parameter).}
#' \item{rl.orig}{A numeric value or a numeric vector of the release limit(s)
#' on the original scale. If \code{NA} was passed in, \code{NA} is returned.}
#' \item{rl.sf}{A numeric value or a numeric vector that specifies the
#' significant figures of the release limit(s). If \code{NA} was passed in,
#' \code{NA} is returned.}
#' \item{rl}{A numeric value or a numeric vector of the adjusted (as specified
#' by the \code{sf_option} parameter) release limit(s). If \code{NA} was
#' passed in, \code{NA} is returned.}
#' \item{rl.trfmd}{A numeric value or a numeric vector of the adjusted and
#' transformed, if applicable (as specified by the the \code{sf_option}
#' parameter and the second element of the \code{xform} parameter,
#' respectively), release limit(s), otherwise the same as \code{rl}.}
#' \item{sl.orig}{A numeric value or a numeric vector of length \code{2} of the
#' specification limit(s) on the original scale.}
#' \item{sl.sf}{A numeric value or a numeric vector of length \code{2}
#' that specifies the significant figures of the specification limit(s).}
#' \item{sl}{A numeric value or a numeric vector of length \code{2} of the
#' adjusted (as specified by the \code{sf_option} parameter) specification
#' limit(s).}
#' \item{sl.trfmd}{A numeric value or a numeric vector of length \code{2} of
#' the adjusted and transformed, if applicable (as specified by the the
#' \code{sf_option} parameter and the second element of the \code{xform}
#' parameter, respectively) specification limit(s), otherwise the same as
#' \code{sl}.}
#'
#' @seealso \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
set_limits <- function(rl, rl_sf, sl, sl_sf, sf_option = "tight",
xform = c("no", "no"), shift = c(0, 0),
ivl_side = "lower") {
if (!is.numeric(rl) && all(!is.na(rl))) {
stop("The parameter rl must be a numeric value or NA.")
}
if (!is.numeric(rl_sf) && all(!is.na(rl_sf))) {
stop("The parameter rl_sf must be a positive integer of length rl, or NA.")
}
if (sum(rl_sf < 0) > 0 && all(!is.na(rl_sf))) {
stop("The parameter rl_sf must be a positive integer of length rl, or NA.")
}
if (length(rl_sf) != length(rl)) {
stop("The parameter rl_sf must be a positive integer of length rl, or NA.")
}
if (!isTRUE(all.equal(rl_sf, as.integer(rl_sf))) && all(!is.na(rl_sf))) {
stop("The parameter rl_sf must be a positive integer of length rl, or NA.")
}
if (!is.numeric(sl) || length(sl) > 2) {
stop("The parameter sl must be a numeric value or vector of length 1 or 2.")
}
if (length(sl) == 2) {
if (sl[2] < sl[1]) {
stop("The parameter sl must be of the form c(lower, upper).")
}
}
if (!is.numeric(sl_sf) && all(!is.na(sl_sf))) {
stop("The parameter sl_sf must be a positive integer of length sl.")
}
if (sum(sl_sf < 0) > 0) {
stop("The parameter sl_sf must be a positive integer of length sl.")
}
if (length(sl_sf) != length(sl)) {
stop("The parameter sl_sf must be a positive integer of length sl.")
}
if (!isTRUE(all.equal(sl_sf, as.integer(sl_sf)))) {
stop("The parameter sl_sf must be a positive integer of length sl.")
}
if (!(sf_option %in% c("tight", "loose"))) {
stop("Please specify sf_option either as \"tight\" or \"loose\".")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
if (!(ivl_side %in% c("lower", "upper", "both"))) {
stop("Please specify ivl_side either as \"lower\", \"upper\" or \"both\".")
}
if (ivl_side == "both" && length(sl) == 1) {
stop("Since ivl_side = \"both\", a specification with two sides is ",
"expected. Only one side has been specified, though, i.e. ",
"sl = ", sl, ".\nPlease provide a specification with two sides.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Save the original values of rl and sl.
rl_orig <- rl
sl_orig <- sl
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# If necessary, determine the adjusted (based on rl_sf and sl_sf) values of
# rl and sl.
if (sf_option == "loose") {
# Specification limit(s)
sl_factor <- 10^(get_n_whole_part(sl) - 1)
sl_std <- signif(sl / sl_factor, sl_sf)
sl_size <- as.numeric(sl < 1)
n_zero_sl <- attr(regexpr("(?<=\\.)0+|$",
as.character(format(sl, scientific = FALSE,
drop0trailing = TRUE)),
perl = TRUE),
"match.length")
if (length(sl) == 2) {
sl[1] <- (sl_std[1] - 5 / 10^(sl_sf[1] + n_zero_sl[1] + sl_size[1])) *
sl_factor[1]
sl[2] <- (sl_std[2] + 4 / 10^(sl_sf[2] + n_zero_sl[2] + sl_size[2])) *
sl_factor[2]
} else {
switch(ivl_side,
"lower" = {
sl <- (sl_std - 5 / 10^(sl_sf + n_zero_sl + sl_size)) *
sl_factor
},
"upper" = {
sl <- (sl_std + 4 / 10^(sl_sf + n_zero_sl + sl_size)) *
sl_factor
})
}
# Release limit(s)
rl_factor <- 10^(get_n_whole_part(rl) - 1)
rl_std <- signif(rl / rl_factor, rl_sf)
rl_size <- as.numeric(rl < 1)
n_zero_rl <- attr(regexpr("(?<=\\.)0+|$",
as.character(format(rl, scientific = FALSE,
drop0trailing = TRUE)),
perl = TRUE),
"match.length")
if (ivl_side != "both") {
switch(ivl_side,
"lower" = {
rl <- (rl_std - 5 / 10^(rl_sf + n_zero_rl + rl_size)) *
rl_factor
},
"upper" = {
rl <- (rl_std + 4 / 10^(rl_sf + n_zero_rl + rl_size)) *
rl_factor
})
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Based on the transformation of the response variable the rl and sl limits
# also have to be transformed.
if (xform[2] != "no") {
switch(xform[2],
"log" = {
rl_trfmd <- log(rl + shift[2])
sl_trfmd <- log(sl + shift[2])
},
"sqrt" = {
rl_trfmd <- sqrt(rl + shift[2])
sl_trfmd <- sqrt(sl + shift[2])
},
"sq" = {
rl_trfmd <- (rl + shift[2])^2
sl_trfmd <- (sl + shift[2])^2
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Collect and return the data
if (xform[2] != "no") {
l_res <- list(sf.option = sf_option,
xform = xform,
shift = shift,
rl.orig = rl_orig,
rl.sf = rl_sf,
rl = rl,
rl.trfmd = rl_trfmd,
sl.orig = sl_orig,
sl.sf = sl_sf,
sl = sl,
sl.trfmd = sl_trfmd)
} else {
l_res <- list(sf.option = sf_option,
xform = xform,
shift = shift,
rl.orig = rl_orig,
rl.sf = rl_sf,
rl = rl,
rl.trfmd = rl,
sl.orig = sl_orig,
sl.sf = sl_sf,
sl = sl,
sl.trfmd = sl)
}
return(l_res)
}
#' Get relevant limits
#'
#' The function \code{get_relevant_limits()} expects a list returned by the
#' function \code{set_limits} and returns a list of only the relevant limits,
#' i.e. those that are relevant with respect to transformation.
#'
#' @param limits_list A list returned by the \code{set_limits()} function.
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_relevant_limits()} makes a subset of the
#' list returned by the \code{set_limits()} function.
#'
#' @return A list with the following elements is returned:
#' \item{sl.orig}{A numeric value or a numeric vector of length \code{2} of the
#' specification limit(s) on the original scale.}
#' \item{sl}{A numeric value or a numeric vector of length \code{2} of adjusted
#' specification limit(s).}
#' \item{rl.orig}{A numeric value or a numeric vector that specifies the release
#' limit(s) on the original scale. If \code{NA} was passed in, \code{NA} is
#' returned.}
#' \item{rl}{A numeric value or a numeric vector of adjusted release limit(s).
#' If \code{NA} was passed in, \code{NA} is returned.}
#' \item{sl.bt}{A numeric value or a numeric vector of length \code{2} of
#' adjusted specification limit(s) before transformation. If no
#' transformation has been performed it is \code{NA}.}
#' \item{rl.bt}{A numeric value or a numeric vector of adjusted release
#' limit(s) before transformation. If no transformation has been performed
#' it is \code{NA}.}
#'
#' @seealso \code{\link{set_limits}}
#'
#' @keywords internal
#' @noRd
get_relevant_limits <- function(limits_list, xform = c("no", "no"),
ivl_side = "lower") {
if (!is.list(limits_list)) {
stop("The parameter limits_list must be a list.")
}
if (sum(names(limits_list) %in%
c("sf.option", "xform", "shift", "sl.orig", "sl.sf", "sl",
"sl.trfmd")) != 7) {
stop("The limits_list must have at least the elements \"sf.option\",
\"xform\", \"shift\", \"sl.orig\", \"sl.sf\", \"sl\" and ",
"\"sl.trfmd\".")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!(ivl_side %in% c("lower", "upper", "both"))) {
stop("Please specify ivl_side either as \"lower\", \"upper\" or \"both\".")
}
if (ivl_side == "both" && length(limits_list[["sl"]]) == 1) {
stop("Since ivl_side = \"both\", a specification with two sides is ",
"expected. Only one side has been specified, though, i.e. ",
"sl = ", limits_list[["sl"]], ".\nPlease provide a specification ",
"with two sides.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extract the relevant limits
# ---------
# Specification limits
if (length(limits_list[["sl"]]) == 2 && ivl_side != "both") {
switch(ivl_side,
"lower" = {
sl_orig <- limits_list[["sl.orig"]][1]
if (xform[2] == "no") {
sl <- limits_list[["sl"]][1]
sl_bt <- NA
} else {
sl <- limits_list[["sl.trfmd"]][1]
sl_bt <- limits_list[["sl"]][1]
}
},
"upper" = {
sl_orig <- limits_list[["sl.orig"]][2]
if (xform[2] == "no") {
sl <- limits_list[["sl"]][2]
sl_bt <- NA
} else {
sl <- limits_list[["sl.trfmd"]][2]
sl_bt <- limits_list[["sl"]][2]
}
})
} else {
sl_orig <- limits_list[["sl.orig"]]
if (xform[2] == "no") {
sl <- limits_list[["sl"]]
sl_bt <- NA
} else {
sl <- limits_list[["sl.trfmd"]]
sl_bt <- limits_list[["sl"]]
}
}
# ---------
# Release limits
rl_orig <- limits_list[["rl.orig"]]
if (xform[2] == "no") {
rl <- limits_list[["rl"]]
rl_bt <- NA
} else {
rl <- limits_list[["rl.trfmd"]]
rl_bt <- limits_list[["rl"]]
}
# ---------
# Compile and return results
return(list(sl.orig = sl_orig,
sl = sl,
rl.orig = rl_orig,
rl = rl,
sl.bt = sl_bt,
rl.bt = rl_bt))
}
#' Result of ANCOVA model check
#'
#' The function \code{check_ancova()} fits an ANalysis of COVAriance (ANCOVA)
#' model to figure out which kind of linear regression model suits the
#' (historical) data best.
#'
#' @param alpha A numeric value that specifies the significance level for the
#' decision which model is appropriate, i.e. if the assumption of
#' \emph{common slope} or \emph{common intercept} is appropriate or not.
#' The default is \code{0.05}.
#' @inheritParams expirest_osle
#'
#' @details The function \code{check_ancova()} fits an ANCOVA (ANalyis of
#' COVAriance) model to the data contained in the provided data frame. Based
#' on \code{alpha}, it checks if the intercepts and/or slopes between the
#' groups differ significantly or not. \cr
#' Three possible models are regarded as appropriate,
#' i.e.
#' \itemize{
#' \item a \emph{common intercept / common slope} model (cics),
#' \item a \emph{different intercept / common slope} model (dics) or
#' \item a \emph{different intercept / different slope} model (dids).
#' }
#' The \emph{common intercept / different slope} model (cids) is of limited
#' practical relevance because the corresponding model is missing an effect.
#' When slopes exhibit significant differences, comparing intercepts becomes
#' inconsequential. Moreover, while initial levels of different batches in a
#' manufacturing process might be relatively well-controlled, it is improbable
#' that these levels are identical. Consequently, if the model probabilities
#' associated with the intercepts and slopes suggest the appropriateness of
#' the cids model, the decision is taken in favour of a dids model.
#'
#' @return A list of two elements is returned that specifies which model, based
#' on the ANCOVA analysis, suits best. The first element (\code{type.spec})
#' is a numeric vector of length 2 that specifies the best model accepted at
#' the significance level specified by \code{alpha}. It has the form
#' \code{c(ci, cs)}, where \code{ci} specifies if a common intercept is
#' appropriate (\code{ci = 1}) or not (\code{ci = 0}) and \code{cs} specifies
#' if a common slope is appropriate (\code{cs = 1}) or not (\code{cs = 0}).
#' The second element (\code{type.acronym}) is an acronym representing the
#' first item. In case of a linear model including only a single batch,
#' all elements are \code{NA}.
#'
#' @seealso \code{\link{expirest_osle}}, \code{\link{expirest_wisle}},
#' \code{\link[stats]{aov}}.
#'
#' @importFrom stats aov
#' @importFrom stats as.formula
#' @importFrom stats summary.aov
#' @importFrom stats setNames
#'
#' @keywords internal
#' @noRd
check_ancova <- function(data, response_vbl, time_vbl, batch_vbl,
alpha = 0.05) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (!(response_vbl %in% colnames(data))) {
stop("The response_vbl was not found in the provided data frame.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (!(time_vbl %in% colnames(data))) {
stop("The time_vbl was not found in the provided data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Remove unused factor levels
data <- droplevels(data)
if (nlevels(data[, batch_vbl]) > 1) {
t_formula <-
paste(response_vbl, "~", paste(batch_vbl, time_vbl, sep = " * "))
lm_ancova <-
do.call("aov", list(as.formula(t_formula), data = as.name("data")))
slm_ancova <- summary(lm_ancova)[[1]]
p_batch <-
slm_ancova[grepl(batch_vbl, rownames(slm_ancova)) &
!grepl(time_vbl, rownames(slm_ancova)), "Pr(>F)"]
p_interaction <-
slm_ancova[grepl(batch_vbl, rownames(slm_ancova)) &
grepl(time_vbl, rownames(slm_ancova)), "Pr(>F)"]
# Store the outcome of the test in two logical parameters, i.e.
# common_icpt: yes or no (common intercept model)
# common_slp: yes or no (common slope model)
ifelse(p_batch > alpha, common_icpt <- 1L, common_icpt <- 0L)
ifelse(p_interaction > alpha, common_slp <- 1L, common_slp <- 0L)
# In case of a cids model, set the common_icpt parameter to 0 to obtain
# a dids model because the cids model is regarded as not relevant
# (See # LeBlond, D., Griffith, D. and Aubuchon, K. Linear Regression 102:
# Stability Shelf Life Estimation Using Analysis of Covariance.
# J Valid Technol (2011) 17(3): 47-68.)
if (common_icpt == 1 && common_slp == 0) {
common_icpt <- 0
}
# ---------
# Compile and return results
l_model_type <-
list(type.spec = setNames(c(common_icpt, common_slp),
c("common.icpt", "common.slp")),
type.acronym =
paste0(c("di", "ci")[common_icpt + 1],
c("ds", "cs")[common_slp + 1]))
} else {
l_model_type <-
list(type.spec = setNames(c(NA, NA), c("common.icpt", "common.slp")),
type.acronym = "n.a.")
}
return(l_model_type)
}
#' Linear model fitting
#'
#' The function \code{get_model_list()} fits four types of linear regression
#' models that are often used for the assessment of stability data, e.g. for
#' the estimation of the shelf life or retest period.
#'
#' @inheritParams expirest_osle
#'
#' @details The function \code{get_model_list()} expects a data frame with
#' a response variable, a time variable and a categorical variable which
#' usually has factor levels of multiple batches of a drug product that was
#' assessed over a certain period of time with respect to the time-dependent
#' behaviour of characteristic parameters. Using these results, the function
#' fits
#' \itemize{
#' \item a \emph{common intercept / common slope} model (cics),
#' \item a \emph{different intercept / common slope} model (dics) or
#' \item a \emph{different intercept / different slope} model with pooled
#' mean square error (dids.pmse) and
#' \item a \emph{different intercept / different slope} model (dids) in which
#' individual models are fitted to each level of the categorical variable.
#' }
#'
#' If the categorical variable has only a single factor level, then the first
#' three models are \code{NA} and only a single regression model is fitted.
#'
#' @return A list of three elements is returned,
#' containing the following elements:
#' \item{Models}{A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids}. The first three elements contain the
#' \sQuote{\code{lm}} objects of the \dQuote{common intercept / common slope}
#' (\code{cics}), \dQuote{different intercept / common slope} (\code{dics})
#' and \dQuote{different intercept / different slope} (\code{dids}) models.
#' The fourth element is a list of the \sQuote{\code{lm}} objects that is
#' obtained from fitting a regression model to the data of each level of the
#' categorical variable separately. The \code{cics}, \code{dics} and
#' \code{dids.pmse} elements are \code{NA} if data of only a single batch
#' is available.}
#' \item{AIC}{A numeric named vector of the Akaike Information Criterion (AIC)
#' values of the \code{cics}, \code{dics} and \code{dids.pmse} models.}
#' \item{BIC}{A numeric named vector of the Bayesian Information Criterion (BIC)
#' values of each of the \code{cics}, \code{dics} and \code{dids.pmse}
#' models.}
#'
#' @seealso \code{\link{expirest_osle}}, \code{\link{expirest_wisle}},
#' \code{\link[stats]{lm}}, \code{\link[stats]{AIC}}, \code{\link[stats]{BIC}}.
#'
#' @importFrom stats lm
#' @importFrom stats as.formula
#' @importFrom stats coef
#' @importFrom stats AIC
#' @importFrom stats BIC
#' @importFrom stats setNames
#'
#' @keywords internal
#' @noRd
get_model_list <- function(data, response_vbl, time_vbl, batch_vbl) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(response_vbl)) {
stop("The parameter response_vbl must be a string.")
}
if (!(response_vbl %in% colnames(data))) {
stop("The response_vbl was not found in the provided data frame.")
}
if (!is.character(time_vbl)) {
stop("The parameter time_vbl must be a string.")
}
if (!(time_vbl %in% colnames(data))) {
stop("The time_vbl was not found in the provided data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Remove unused factor levels
d_dat <- droplevels(data)
l_models <- setNames(vector(mode = "list", length = 4),
c("cics", "dics", "dids.pmse", "dids"))
if (nlevels(d_dat[, batch_vbl]) > 1) {
# ---------
# Common Intercept / Common Slope
t_formula <- paste(response_vbl, "~", time_vbl)
l_models[["cics"]] <-
do.call("lm", list(as.formula(t_formula), data = as.name("d_dat")))
# ---------
# Different Intercept / Common Slope
t_formula <-
paste(response_vbl, "~", paste(batch_vbl, time_vbl, sep = " + "))
l_models[["dics"]] <-
do.call("lm", list(as.formula(t_formula), data = as.name("d_dat")))
# ---------
# Different Intercept / Different Slope (pooled mean square error)
t_formula <-
paste(response_vbl, "~", paste(batch_vbl, time_vbl, sep = " * "))
l_models[["dids.pmse"]] <-
do.call("lm", list(as.formula(t_formula), data = as.name("d_dat")))
# ---------
# Different Intercept / Different Slope (individual models)
t_formula <- paste(response_vbl, "~", time_vbl)
tmp <- lapply(levels(d_dat[, batch_vbl]),
function(batch) {
t_dat <- d_dat[d_dat[, batch_vbl] == batch, ]
do.call("lm", list(as.formula(t_formula),
data = as.name("t_dat")))
})
names(tmp) <- levels(d_dat[, batch_vbl])
l_models[["dids"]] <- tmp
# ---------
# Determination of the Akaike Information Criterion (AIC) and Bayesian
# Information Criterion (BIC) of each of the relevant models
t_AIC <- vapply(l_models[c("cics", "dics", "dids.pmse")], AIC, numeric(1))
t_BIC <- vapply(l_models[c("cics", "dics", "dids.pmse")], BIC, numeric(1))
} else {
t_formula <- paste(response_vbl, "~", time_vbl)
l_models[names(l_models) != "dids"] <- NA
tmp <- lapply(levels(d_dat[, batch_vbl]),
function(batch) {
t_dat <- d_dat[d_dat[, batch_vbl] == batch, ]
do.call("lm", list(as.formula(t_formula),
data = as.name("t_dat")))
})
names(tmp) <- levels(d_dat[, batch_vbl])
l_models[["dids"]] <- tmp
t_AIC <- t_BIC <- setNames(rep(NA, 3), c("cics", "dics", "dids.pmse"))
}
return(list(Models = l_models,
AIC = t_AIC,
BIC = t_BIC))
}
#' Extract information from \dQuote{worst case scenario} (wcs) limit lists list
#'
#' The function \code{extract_from_ll_wcsl()} extracts specific elements from
#' a list of lists returned by the \code{\link{get_wcs_limit}()} function.
#'
#' @param ll A list of lists returned by the \code{\link{get_wcs_limit}()}
#' function. The list must have four elements named \code{"cics"},
#' \code{"dics"}, \code{dids.pmse} and \code{"dids"}. Each of these elements
#' has a sub-list of the same length as the number of intercepts. And each
#' of these elements has a sub-sub-list of the same length as the number of
#' release limits (\code{rl}).
#' @param element A character string that specifies the element to be extracted,
#' i.e. either one of \code{"delta.lim"}, \code{"delta.lim.orig"},
#' \code{"wcs.lim"} or \code{"wcs.lim.orig"}.
#'
#' @details Information in a bulk list of lists that has been obtained by
#' aid of the function \code{\link{get_wcs_limit}()} for a set of release
#' limit values (\code{rl}) and intercepts.
#'
#' @return A list of the same length as \code{ll_wcsl} is returned. The
#' individual elements of the list are matrices of the values specified by
#' \code{element} that have been extracted from \code{ll}.
#'
#' @seealso \code{\link{get_wcs_limit}}.
#'
#' @keywords internal
#' @noRd
extract_from_ll_wcsl <- function(ll, element) {
if (sum(names(ll) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The list ll must have four elements named \"cics\", \"dics\" ",
"\"dids.pmse\" and \"dids\".")
}
if (!(element %in%
c("delta.lim", "delta.lim.orig", "wcs.lim", "wcs.lim.orig"))) {
stop("Please specify element either as \"delta.lim\", \"delta.lim.orig\", ",
"\"wcs.lim\" or \"wcs.lim.orig\".")
}
if (get_n_list_levels(ll) != 4) {
stop("The parameter ll must be a list of lists returned by ",
"get_wcs_limit() at level three.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of worst case scenario (wcs) limit(s)
t_expected_names <-
c("delta.lim", "delta.lim.orig", "wcs.lim", "wcs.lim.orig")
t_model_names <- names(ll)
l_res <- lapply(t_model_names, function(mm) {
if (mm == "cics") {
if (get_n_list_levels(ll[[mm]]) == 0) {
NA
} else {
if (any(vapply(ll[[mm]], function(x) {
sum(t_expected_names %in% names(x[[1]])) != 4
},
logical(1)))) {
stop("The element was not found in one of the sub-elements of ",
"model ", mm, ". Please provide a list returned ",
"from get_wcs_limit().")
} else {
matrix(vapply(seq_along(ll[[mm]][[1]]), function(j) {
ll[[mm]][[1]][[j]][[element]]
},
numeric(1)),
nrow = length(ll[[mm]][[1]]), ncol = length(ll[[mm]]),
dimnames = list(NULL, names(ll[[mm]])))
}
}
} else {
if (get_n_list_levels(ll[[mm]]) == 0) {
NA
} else {
if (any(vapply(ll[[mm]], function(x) {
sum(t_expected_names %in% names(x[[1]])) != 4
},
logical(1)))) {
stop("The element was not found in one of the sub-elements of ",
"model ", mm, ". Please provide a list returned ",
"from get_wcs_limit().")
} else {
matrix(vapply(seq_along(ll[[mm]]), function(bb) {
vapply(seq_along(ll[[mm]][[1]]), function(j) {
ll[[mm]][[bb]][[j]][[element]]
},
numeric(1))
}, numeric(length(ll[[mm]][[1]]))),
nrow = length(ll[[mm]][[1]]), ncol = length(ll[[mm]]),
dimnames = list(NULL, names(ll[[mm]])))
}
}
}
})
names(l_res) <- t_model_names
return(l_res)
}
#' Extract worst case x value
#'
#' The function \code{extract_wc_x()} extracts the worst case \eqn{x} value
#' from a list of matrices of possible \eqn{x} values based on a list of
#' vectors of indices which specify the worst case elements.
#'
#' @param l1 A list of matrices of \eqn{x} values or a list of lists of one
#' vectors of \eqn{x} values. The list must have four elements named
#' \code{"cics"}, \code{"dics"}, \code{"dids.pmse"} and \code{"dids"}.
#' @param l2 A list of vectors of indices. As \code{l1}, the list must have
#' four elements named \code{"cics"}, \code{"dics"}, \code{"dids.pmse"} and
#' \code{"dids"}. The length of the vectors of \code{l2} must be equal to
#' the number of rows of the matrices in \code{l1} or the length of the
#' vectors in \code{l1}.
#'
#' @details Information from a list of matrices or a list of lists of vectors
#' is extracted by aid of a list of vectors of indices which specify which
#' elements per row of each matrix or which elements of the vectors have to
#' be returned.
#'
#' @return A matrix of the worst case values with the number of rows
#' corresponding to the length of the vectors in \code{l2} and the number of
#' columns corresponding to the length of \code{l1} or \code{l2} is returned.
#'
#' @keywords internal
#' @noRd
extract_wc_x <- function(l1, l2) {
if (!is.list(l1)) {
stop("The parameter l1 must be a list.")
}
if (!is.list(l2)) {
stop("The parameter l2 must be a list.")
}
if (sum(names(l1) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The list l1 must have four elements named \"cics\", \"dics\" ",
"\"dids.pmse\" and \"dids\".")
}
if (sum(names(l2) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The list l2 must have four elements named \"cics\", \"dics\" ",
"\"dids.pmse\" and \"dids\".")
}
if (get_n_list_levels(l1) == 1) {
if (sum(vapply(l1, function(x) {
!is.matrix(x)
},
logical(1))) > 0) {
stop("The elements of l1 must be matrices or lists of vectors.")
}
} else {
if (sum(vapply(l1, function(x) {
!is.numeric(x[[1]]) & !is.logical(x[[1]])
},
logical(1))) > 0) {
stop("The elements of l1 must be matrices or lists of vectors.")
}
}
if (sum(vapply(l2, function(x) {
!is.numeric(x) & !is.logical(x)
},
logical(1))) > 0) {
stop("The elements of l2 must be numeric vectors or vectors of NA.")
}
if (sum(vapply(l1, function(x) is.matrix(x), logical(1))) == 4) {
if (!isTRUE(all.equal(vapply(l1, function(x) nrow(x), numeric(1)),
vapply(l2, function(x) length(x), numeric(1))))) {
stop("The number of rows of the matrices in l1 must be equal ",
"to the length of the vectors in l2.")
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
m_res <- matrix(NA, nrow = length(l2[[1]]), ncol = length(l1))
colnames(m_res) <- names(l1)
for (i in seq_along(l1)) {
if (names(l1)[i] == "cics") {
if (is.matrix(l1[["cics"]])) {
m_res[, "cics"] <- l1[["cics"]]
}
if (is.list(l1[["cics"]])) {
if (is.vector(l1[["cics"]][[length(l1[["cics"]])]]) &&
length(l1[["cics"]][[length(l1[["cics"]])]]) == 1) {
m_res[, "cics"] <- rep(unname(l1[[i]][[1]]), length(l2[[1]]))
}
}
} else {
if (is.matrix(l1[[i]])) {
m_res[, i] <- vapply(seq_along(l2[[1]]), function(j) {
ifelse(!is.na(l2[[i]][j]),
l1[[i]][j, l2[[i]][j]],
NA)
}, numeric(1))
}
if (is.list(l1[[i]])) {
m_res[, i] <- vapply(seq_along(l2[[1]]), function(j) {
ifelse(!is.na(l2[[i]][j]),
l1[[i]][[1]][l2[[i]][j]],
NA)
},
numeric(1))
}
}
}
return(m_res)
}
#' Compile information on worst case batches for ordinary shelf life estimation
#'
#' The function \code{get_osle_poi_list()} prepares a list of points of
#' intersection (POI) for multiple regression models using the
#' \code{find_poi()} function.
#'
#' @inheritParams expirest_osle
#' @inheritParams get_wc_icpt
#' @inheritParams get_poi_list
#'
#' @details The function \code{get_osle_poi_list()} applies the
#' \code{find_poi()} function (find the \dQuote{point of intersection}) on
#' all the models and for each release limit (\code{rl}) provided. With respect
#' to the latter it differs from the \code{\link{get_poi_list}()} function.
#'
#' @return A list with the following elements is returned:
#' \item{all.poi}{A list of the POI values, i.e. a list with one or two list
#' elements for the side (i.e. \code{lower} or \code{upper}) of the
#' corresponding specification limit, each containing a list returned by the
#' \code{\link{get_poi_list}()} function, containing of four elements named
#' \code{cics}, \code{dics}, \code{dids.pmse} and \code{dids}. Each of them
#' contains a named vector of the POI values estimated for each batch and
#' named accordingly.}
#' \item{poi}{A named vector of the worst case POI values of each model, i.e.
#' named \code{cics}, \code{dics}, \code{dids.pmse} and \code{dids}. In
#' addition, the vector has an attribute called \code{side} that specifies
#' the side of the specification limit which is crossed by the confidence or
#' prediction interval at the corresponding POI value.}
#' \item{wc.icpt}{A named vector of the intercepts of the worst case batches of
#' each model, i.e. named \code{cics}, \code{dics}, \code{dids.pmse} and
#' \code{dids}. In addition, the vector has an attribute called \code{side}
#' that specifies the side of the specification limit which is crossed by the
#' confidence or prediction interval at the POI value of the corresponding
#' worst case batch.}
#' \item{which.wc.batch}{A named vector of the indices of the worst case
#' batches of each model, i.e. named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids}. In addition, the vector has an attribute
#' called \code{side} that specifies the side of the specification limit
#' which is crossed by the confidence or prediction interval at the POI value
#' of the corresponding worst case batch.}
#'
#' @seealso \code{\link{get_icpt_list}}, \code{\link{get_model_list}},
#' \code{\link{get_poi_list}}, \code{\link{get_wc_icpt}},
#' \code{\link{expirest_osle}}, \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_osle_poi_list <- function(data, batch_vbl, icpt_list, model_list, sl,
srch_range, alpha = 0.05, xform = c("no", "no"),
shift = c(0, 0), ivl = "confidence",
ivl_type = "one.sided", ivl_side = "lower", ...) {
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.list(icpt_list) ||
sum(names(icpt_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The icpt_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.list(model_list) ||
sum(names(model_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The model_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.numeric(sl) || length(sl) > 2) {
stop("The parameter sl must be a numeric or vector of length 1 or 2.")
}
if (length(sl) == 2) {
if (sl[2] < sl[1]) {
stop("The parameter sl must be of the form c(lower, upper).")
}
}
if (!is.numeric(srch_range) || length(srch_range) != 2) {
stop("The parameter srch_range must be a vector of length 2.")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
if (!(ivl %in% c("confidence", "prediction"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper", "both"))) {
stop("Please specify ivl_side either as \"lower\", \"upper\" or \"both\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
d_dat <- droplevels(data)
l_icpt <- icpt_list
l_models <- model_list
t_sides <- c("lower", "upper")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of POI values of all models
if (ivl_side == "both" && length(sl) == 2) {
l_poi <-
lapply(seq_along(t_sides), function(i) {
get_poi_list(data = d_dat, batch_vbl = batch_vbl, model_list = l_models,
srch_range = srch_range, sl = sl[i], mode = "minimal",
alpha = alpha, ivl = ivl, ivl_type = ivl_type,
ivl_side = t_sides[i])
})
names(l_poi) <- t_sides
# Extraction of all worst case POI values
d_poi <-
rbind(lower = vapply(l_poi$lower, function(x) {
ifelse(all(is.na(x)),
NA,
min(x, na.rm = TRUE))
},
numeric(1)),
upper = vapply(l_poi$upper, function(x) {
ifelse(all(is.na(x)),
NA,
min(x, na.rm = TRUE))
},
numeric(1)))
d_poi[is.infinite(d_poi)] <- NA
# Determination of the side of the worst case POI value
t_poi_side <- vapply(colnames(d_poi), function(cn) {
ifelse(all(is.na(d_poi[, cn])),
"NA",
t_sides[which.min(d_poi[, cn])])
},
character(1))
# Summary vector of the worst case POI values
t_poi <- vapply(colnames(d_poi), function(cn) {
ifelse(t_poi_side[cn] == "NA",
NA,
d_poi[t_poi_side[cn], cn])
},
numeric(1))
attr(t_poi, "side") <- t_poi_side
} else {
l_poi <- setNames(list(NA), ivl_side)
l_poi[[ivl_side]] <-
get_poi_list(data = d_dat, batch_vbl = batch_vbl, model_list = l_models,
srch_range = srch_range, sl = sl, mode = "minimal",
alpha = alpha, ivl = ivl, ivl_type = ivl_type,
ivl_side = ivl_side)
# "Determination" of the side of the worst case POI value
t_poi_side <- rep(ivl_side, length(l_models))
names(t_poi_side) <- names(l_models)
# Summary vector of the worst case POI values
t_poi <- vapply(l_poi[[ivl_side]], function(x) {
ifelse(all(is.na(x)),
NA,
min(x, na.rm = TRUE))
},
numeric(1))
t_poi[is.infinite(t_poi)] <- NA
attr(t_poi, "side") <- t_poi_side
}
if (nlevels(d_dat[, batch_vbl]) > 1) {
if (sum(is.na(t_poi)) != 0) {
if (min(srch_range) == 0) {
warning("Not for all model types POI values obtained. ",
"Possibly, changing srch_range could solve the issue ",
"(a lower limit > 0 might be a solution).")
} else {
warning("Not for all model types POI values obtained. ",
"Possibly, changing srch_range could solve the issue. ")
}
}
} else {
if (is.na(t_poi["dids"])) {
if (min(srch_range) == 0) {
warning("No POI value was obtained. ",
"Possibly, changing srch_range could solve the issue ",
"(a lower limit > 0 might be a solution).")
} else {
warning("No POI value was obtained. ",
"Possibly, changing srch_range could solve the issue. ")
}
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of worst case batch (wc_batch) and its intercept (wc_icpt)
# In case of cics model: wc_icpt is the common intercept of all batches
# and none of the batches is the worst case batch and thus NA.
# Worst case batch (of each model)
if (ivl_side == "both" && length(sl) == 2) {
l_wc_batch <- lapply(l_poi, function(ll) {
vapply(ll, function(x) {
ifelse(!length(which.min(x)),
NA,
which.min(x))
},
numeric(1))
})
l_wc_batch$lower["cics"] <- NA
l_wc_batch$upper["cics"] <- NA
# Summary vector of the worst case batches
wc_batch <- vapply(names(t_poi_side), function(cn) {
ifelse(t_poi_side[cn] == "NA",
NA,
l_wc_batch[[t_poi_side[cn]]][cn])
},
numeric(1))
attr(wc_batch, "side") <- t_poi_side
} else {
wc_batch <- vapply(l_poi[[ivl_side]], function(x) {
ifelse(!length(which.min(x)),
NA,
which.min(x))
},
numeric(1))
wc_batch["cics"] <- NA
attr(wc_batch, "side") <- t_poi_side
}
# Intercept of the worst case batch (of each model)
if (ivl_side == "both" && length(sl) == 2) {
l_wc_icpt <-
lapply(t_sides, function(side) {
get_wc_icpt(data = d_dat, batch_vbl = batch_vbl,
icpt_list = l_icpt, poi_list = l_poi[[side]],
wc_batch = l_wc_batch[[side]], xform = xform)
})
names(l_wc_icpt) <- t_sides
# Summary vector of the intercepts of the worst case batches
wc_icpt <- vapply(names(t_poi_side), function(cn) {
ifelse(t_poi_side[cn] == "NA",
NA,
l_wc_icpt[[t_poi_side[cn]]][cn])
},
numeric(1))
attr(wc_icpt, "side") <- t_poi_side
} else {
wc_icpt <- get_wc_icpt(data = d_dat, batch_vbl = batch_vbl,
icpt_list = l_icpt, poi_list = l_poi[[ivl_side]],
wc_batch = wc_batch, xform = xform)
attr(wc_icpt, "side") <- t_poi_side
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compile and return results
return(list(all.poi = l_poi,
poi = t_poi,
wc.icpt = wc_icpt,
which.wc.batch = wc_batch))
}
#' Compile information on worst case batches for what-if shelf life estimation
#'
#' The function \code{get_wisle_poi_list()} prepares a list of points of
#' intersection (POI) for multiple regression models and release limits
#' using the \code{find_poi()} function.
#'
#' @inheritParams expirest_wisle
#' @inheritParams get_wc_icpt
#' @inheritParams get_poi_list
#'
#' @details The function \code{get_wisle_poi_list()} applies the
#' \code{find_poi()} function (find the \dQuote{point of intersection}) on
#' all the models and for each release limit (\code{rl}) provided. With respect
#' to the latter it differs from the \code{\link{get_poi_list}()} function.
#'
#' @return A list with the following elements is returned:
#' \item{all.wcsl}{A list of the worst case scenario (wcs) limits with a list
#' of four elements for each linear model named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids}. Each of these elements contains a list
#' element for each batch (intercept) containing itself a list element for
#' each release limit. The wcs limit is obtained by adding/subtracting the
#' absolute difference of specification limit and release limit to/from the
#' common intercept of the test batches or the intercept of the worst
#' performing batch.}
#' \item{all.poi}{A list of the POI values, i.e. a list with a list of four
#' elements for each linear model named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids}. Each of these elements contains a
#' matrix with columns for each intercept and rows for each release limit.}
#' \item{all.prl}{A list of the confidence or prediction interval limits that
#' are associated with a POI value, i.e. a list with a list of four elements
#' for each linear model named \code{cics}, \code{dics}, \code{dids.pmse}
#' and \code{dids}. Each of these elements contains an array with a level for
#' each batch, containing matrices with columns for each batch and rows for
#' each release limit, where the matrices contain the estimated interval
#' limits at each POI value per batch.}
#' \item{which.min.dist}{A list of four elements for each linear model named
#' \code{cics}, \code{dics}, \code{dids.pmse} and \code{dids}. Each of these
#' list elements contains a matrix of the indices of the batches with the
#' minimal intercept in the \code{all.prl} list. The matrices have a column
#' for each batch and a row for each release limit.}
#' \item{which.min.poi}{A list of four elements for each linear model named
#' \code{cics}, \code{dics}, \code{dids.pmse} and \code{dids}. Each of these
#' list elements contains a numeric vector with the minimal POI value of
#' associated with each release limit.}
#' \item{which.wc.batch}{A list of four elements for each linear model named
#' \code{cics}, \code{dics}, \code{dids.pmse} and \code{dids}. Each of these
#' list elements contains a numeric vector with the indices of the worst
#' case batches associated with each release limit.}
#'
#' @seealso \code{\link{get_icpt_list}}, \code{\link{get_model_list}},
#' \code{\link{get_wcs_limit}}, \code{\link{find_poi}},
#' \code{\link{get_intvl_limit}}, \code{\link{expirest_osle}},
#' \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
get_wisle_poi_list <- function(icpt_list, model_list, rl, sl, srch_range,
alpha = 0.05, xform = c("no", "no"),
shift = c(0, 0), ivl = "confidence",
ivl_type = "one.sided",
ivl_side = "lower", ...) {
if (!is.list(icpt_list) ||
sum(names(icpt_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The icpt_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids.pmse\" and \"dids\".")
}
if (!is.list(model_list) ||
sum(names(model_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The model_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids\" and \"dids.pmse\".")
}
if (!is.numeric(rl)) {
stop("The parameter rl must be a numeric.")
}
if (!is.numeric(sl) || length(sl) > 1) {
stop("The parameter sl must be a numeric value of length 1.")
}
if (!is.numeric(srch_range) || length(srch_range) != 2) {
stop("The parameter srch_range must be a vector of length 2.")
}
if (alpha <= 0 || alpha > 1) {
stop("Please specify alpha as (0, 1].")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
if (!(ivl %in% c("confidence", "prediction"))) {
stop("Please specify ivl either as \"confidence\" or \"prediction\".")
}
if (!(ivl_type %in% c("one.sided", "two.sided"))) {
stop("Please specify ivl_type either as \"one.sided\" or \"two.sided\".")
}
if (!(ivl_side %in% c("lower", "upper"))) {
stop("Please specify ivl_side either as \"lower\" or \"upper\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
l_icpt <- icpt_list
l_models <- model_list
# Preliminary definition of the lists that will be required below
l_poi <- l_prl <- l_wc_batch <- setNames(rep(list(NA), 4), names(l_models))
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of worst case scenario (wcs) limits for all intercepts of
# all models (on the transformed scale, if data have been transformed)
# List of all wcs_limit lists
ll_wcsl <- lapply(seq_along(l_icpt),
function(i) {
if (get_n_list_levels(l_icpt[[i]]) != 0) {
lapply(l_icpt[[i]]$icpt, function(xx) {
lapply(rl, function(j) {
get_wcs_limit(rl = j, sl = sl, intercept = xx,
xform = xform, shift = shift,
ivl_side = ivl_side)
})
})
} else {
NA
}
})
names(ll_wcsl) <- names(l_icpt)
l_wcsl <- extract_from_ll_wcsl(ll_wcsl, "wcs.lim")
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Calculation of POI values for all wcs limits
# Determination of worst case POI values
# The worst case POI value is the POI value obtained with the batch whose
# lower or upper confidence or prediction interval limit is closest to the
# corresponding specification limit, i.e. the worst case batch.
# Example: the response is the assay, and the lower specification limit is
# the relevant limit. A batch may have a shorter POI than the other batches,
# but because it has a higher intercept or/and smaller variability than one
# or more of the other batches, the lower confidence or prediction interval
# limit of one of the other batches still may be closer to the lower
# specification limit so that their POI values are the POI values of
# relevance.
for (variety in names(l_wcsl)) {
if (get_n_list_levels(l_icpt[[variety]]) != 0) {
# Initialise empty arrays
m_poi <- matrix(NA,
nrow = length(rl),
ncol = length(l_icpt[[variety]][["icpt"]]))
colnames(m_poi) <- names(l_icpt[[variety]][["icpt"]])
a_prl <- array(NA,
dim = c(length(rl), length(l_icpt[[variety]][["icpt"]]),
length(l_icpt[[variety]][["icpt"]])),
dimnames = list(as.character(seq_along(rl)),
names(l_icpt[[variety]][["icpt"]]),
names(l_icpt[[variety]][["icpt"]])))
# Fill arrays
for (j in seq_along(rl)) {
for (k in seq_len(ncol(l_wcsl[[variety]]))) {
if (variety != "dids") {
tmp_poi <- try_get_model(
find_poi(srch_range = srch_range,
model = l_models[[variety]],
sl = l_wcsl[[variety]][j, k], alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl))
} else {
tmp_poi <- try_get_model(
find_poi(srch_range = srch_range,
model = l_models[["dids"]][[k]],
sl = l_wcsl[[variety]][j, k], alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl))
}
if (is.null(tmp_poi[["Error"]])) {
m_poi[j, k] <- tmp_poi[["Model"]]
if (variety != "dids") {
tmp_prl <- try_get_model(
get_intvl_limit(
x_new = tmp_poi[["Model"]],
model = l_models[[variety]], alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl)
)
if (is.null(tmp_prl[["Error"]])) {
a_prl[j, k, ] <- tmp_prl[["Model"]]
}
} else {
t_prl <- rep(NA, ncol(l_wcsl[[variety]]))
for (kk in seq_len(ncol(l_wcsl[["dids"]]))) {
tmp_prl <- try_get_model(
get_intvl_limit(
x_new = tmp_poi[["Model"]],
model = l_models[["dids"]][[kk]], alpha = alpha,
ivl_type = ivl_type, ivl_side = ivl_side, ivl = ivl)
)
if (is.null(tmp_prl[["Error"]])) {
t_prl[kk] <- tmp_prl[["Model"]]
}
}
a_prl[j, k, ] <- t_prl
}
}
}
}
# Put the resulting arrays into the corresponding list entries
l_poi[[variety]] <- m_poi
l_prl[[variety]] <- a_prl
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determination of the batches with the confidence or prediction interval
# limits that are closest to the respective specification limit for each
# model and each POI
switch(ivl_side,
"lower" = {
l_min_dist <- lapply(l_prl, FUN = function(x) {
if (is.logical(x)) {
NA
} else {
apply(x, c(1, 2), FUN = function(y) {
ifelse(length(which.min(y)) != 0,
which.min(abs(y)),
NA)
})
}
})
},
"upper" = {
l_min_dist <- lapply(l_prl, FUN = function(x) {
if (is.logical(x)) {
NA
} else {
apply(x, c(1, 2), FUN = function(y) {
ifelse(length(which.max(y)) != 0,
which.max(abs(y)),
NA)
})
}
})
})
# Determination of the smallest POI value for each model and each rl value
l_min_poi <- lapply(l_poi, FUN = function(x) {
if (is.logical(x)) {
NA
} else {
apply(x, 1, function(y) {
ifelse(length(which.min(y)) != 0,
which.min(y),
NA)
})
}
})
# Determination of the worst case batches for each model and each rl value:
# The worst case batches are the ones with the confidence or prediction
# interval limits that are closest to the respective specification limit
# where the POI values are smallest.
# In case of cics model: wc_icpt is the common intercept of all batches
# and none of the batches is the worst case batch.
for (i in seq_along(l_min_dist)) {
if (!is.logical(l_min_dist[[i]])) {
if (names(l_min_dist)[i] == "cics") {
l_wc_batch[[i]] <- rep(NA, length(rl))
} else {
l_wc_batch[[i]] <-
vapply(seq_along(rl), function(j) {
ifelse(!is.na(l_min_poi[[i]][j]),
l_min_dist[[i]][j, l_min_poi[[i]][j]],
NA)
},
numeric(1))
}
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Compile and return results
return(list(all.wcsl = ll_wcsl,
all.poi = l_poi,
all.prl = l_prl,
which.min.dist = l_min_dist,
which.min.poi = l_min_poi,
which.wc.batch = l_wc_batch))
}
#' Compile \dQuote{what-if shelf life estimation} (wisle) assessment results
#'
#' The function \code{compile_wisle_summary()} extracts results from various
#' lists that are generated during the wisle estimation and compiles a
#' summary data frame.
#'
#' @param wcsl_list A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} with the worst case scenario limits
#' of each batch and linear regression model. The \code{cics}, \code{dics}
#' and \code{dids.pmse} elements are \code{NA} if data of only a single
#' batch is available.
#' @param wcb_list A list of four elements named \code{cics}, \code{dics},
#' \code{dids.pmse} and \code{dids} with the indices of the worst case
#' batches. The \code{cics}, \code{dics} and \code{dids.pmse} elements are
#' \code{NA} if data of only a single batch is available.
#' @param poi_ich A numeric named vector of the POI values of the worst case
#' batches of each model.
#' @inheritParams expirest_osle
#' @inheritParams get_wc_icpt
#' @inheritParams get_relevant_limits
#'
#' @details Information stored in multiple lists that are generated during the
#' \dQuote{what-if shelf life estimation} is extracted and compiled in a single
#' data frame.
#'
#' @return A list with two element is returned, containing the following
#' elements:
#' \item{wc.icpt}{A data frame of the worst case intercepts of each of the
#' four fitted models.}
#' \item{POI}{A data frame of the intercepts, the differences between release
#' and shelf life limits, the WCSLs, the expiry and release specification
#' limits, the shelf lives and POI values.}
#'
#' Structure of the \code{POI} data frame:
#' \item{Intercept.cics}{The intercept of the worst case batch of the cics
#' model.}
#' \item{Intercept.dics}{The intercept of the worst case batch of the dics
#' model.}
#' \item{Intercept.dids.pmse}{The intercept of the worst case batch of the dids
#' model with pooled mean square error (pmse).}
#' \item{Intercept.dids}{The intercept of the worst case batch of the dids
#' model obtained by fitting individual models to the data of each batch.}
#' \item{Delta.cics}{Absolute difference between the release and and the shelf
#' life specification of the cics model.}
#' \item{Delta.dics}{Absolute difference between the release and and the shelf
#' life specification of the dics model.}
#' \item{Delta.dids.pmse}{Absolute difference between the release and and the
#' shelf life specification of the dids model with pooled mean square error
#' (pmse).}
#' \item{Delta.dids}{Absolute difference between the release and and the shelf
#' life specification of the dids model obtained by fitting individual
#' models to the data of each batch.}
#' \item{WCSL.cics}{WCSL of the cics model.}
#' \item{WCSL.dics}{WCSL of the dics model.}
#' \item{WCSL.dids.pmse}{WCSL of the dids model with pooled mean square error
#' (pmse).}
#' \item{WCSL.dids}{WCSL of the dids model obtained by fitting individual
#' models to the data of each batch.}
#' \item{Exp.Spec}{The (expiry) specification, i.e. the specification which is
#' relevant for the determination of the expiry.}
#' \item{Rel.Spec}{The calculated release specification.}
#' \item{Shelf.Life.cics}{The estimated shelf life of the cics model.}
#' \item{Shelf.Life.dics}{The estimated shelf life of the dics model.}
#' \item{Shelf.Life.dids.pmse}{The estimated shelf life of the dids model with
#' pooled mean square error (pmse).}
#' \item{Shelf.Life.dids}{The estimated shelf life of the dids model obtained
#' by fitting individual models to the data of each batch.}
#' \item{POI.Model.cics}{The POI of the cics model.}
#' \item{POI.Model.dics}{The POI of the dics model.}
#' \item{POI.Model.dids.pmse}{The POI of the dids model with pooled mean
#' square error (pmse).}
#' \item{POI.Model.dids}{The POI of the dids model obtained by fitting
#' individual models to the data of each batch.}
#'
#' @seealso \code{\link{extract_wc_x}}, \code{\link{extract_from_ll_wcsl}},
#' \code{\link{expirest_wisle}}.
#'
#' @keywords internal
#' @noRd
compile_wisle_summary <- function(data, batch_vbl, rl, poi_list, icpt_list,
wcsl_list, wcb_list, limits_list, poi_ich,
xform = c("no", "no"), shift = c(0, 0)) {
# This function replaces the following section in expirest_wisle()
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Collection of data and compilation of summary data frame
if (!is.data.frame(data)) {
stop("The data must be provided as data frame.")
}
if (!is.character(batch_vbl)) {
stop("The parameter batch_vbl must be a string.")
}
if (!(batch_vbl %in% colnames(data))) {
stop("The batch_vbl was not found in the provided data frame.")
}
if (!is.factor(data[, batch_vbl])) {
stop("The column in data specified by batch_vbl must be a factor.")
}
if (!is.numeric(rl)) {
stop("The parameter rl must be a numeric.")
}
if (!is.list(poi_list) ||
sum(names(poi_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The poi_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids.pmse\" and \"dids\".")
}
if (!is.list(icpt_list) ||
sum(names(icpt_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The icpt_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids.pmse\" and \"dids\".")
}
if (!is.list(wcsl_list) ||
sum(names(wcsl_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The wcsl_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids.pmse\" and \"dids\".")
}
if (!is.list(wcb_list) ||
sum(names(wcb_list) %in% c("cics", "dics", "dids.pmse", "dids")) != 4) {
stop("The wcb_list must be a list with four elements named \"cics\", ",
"\"dics\", \"dids.pmse\" and \"dids\".")
}
if (!is.list(limits_list) ||
sum(names(limits_list) %in% c("sl.orig", "sl", "rl.orig", "rl")) != 4) {
stop("The limits_list must be a list with at least the four elements ",
"named \"sl.orig\", \"sl\", \"rl.orig\" and \"rl\".")
}
if (!is.numeric(poi_ich) || length(poi_ich) != 4) {
stop("The parameter poi_ich must be a vector of length 4.")
}
if (!all((names(poi_ich) %in% c("cics", "dics", "dids", "dids.pmse")))) {
stop("The parameter poi_ich must be a vector with the names \"cics\", ",
"\"dics\", \"dids\" or \"dids.pmse\".")
}
if (length(xform) != 2) {
stop("Please specify xform appropriately.")
}
if (!(xform[1] %in% c("no", "log", "sqrt", "sq")) ||
!(xform[2] %in% c("no", "log", "sqrt", "sq"))) {
stop("Please specify xform appropriately.")
}
if (!is.numeric(shift) || length(shift) != 2) {
stop("The parameter shift must be a numeric vector of length 2.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of data
d_dat <- droplevels(data)
l_poi <- poi_list
l_icpt <- icpt_list
ll_wcsl <- wcsl_list
l_wc_batch <- wcb_list
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Matrix of the worst case POI values for each model and each rl value
if (nlevels(d_dat[, batch_vbl]) > 1) {
m_poi <- extract_wc_x(l1 = l_poi, l2 = l_wc_batch)
} else {
m_poi <- matrix(NA, nrow = length(rl), ncol = length(l_wc_batch))
colnames(m_poi) <- names(l_wc_batch)
m_poi[, "dids"] <- as.numeric(l_poi[["dids"]])
}
# Depending on the transformation of the time variable the POI values have to
# be back-transformed.
if (xform[1] != "no") {
switch(xform[1],
"log" = {
m_poi <- exp(m_poi) - shift[1]
},
"sqrt" = {
m_poi <- m_poi^2 - shift[1]
},
"sq" = {
m_poi <- sqrt(m_poi) - shift[1]
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Worst case intercepts (wc_icpt_argpm) (on the original scale)
if (nlevels(d_dat[, batch_vbl]) > 1) {
if (xform[2] == "no") {
wc_icpt_argpm <- extract_wc_x(l1 = l_icpt, l2 = l_wc_batch)
} else {
l_icpt_sub <- lapply(l_icpt, function(x) list(x$icpt.orig))
wc_icpt_argpm <- extract_wc_x(l1 = l_icpt_sub, l2 = l_wc_batch)
}
} else {
wc_icpt_argpm <- matrix(NA, nrow = length(rl), ncol = length(l_wc_batch))
colnames(wc_icpt_argpm) <- names(l_wc_batch)
if (xform[2] == "no") {
wc_icpt_argpm[, "dids"] <- as.numeric(l_icpt[["dids"]][["icpt"]])
} else {
wc_icpt_argpm[, "dids"] <- as.numeric(l_icpt[["dids"]][["icpt.orig"]])
}
}
# ---------
# Delta and WCSL
if (xform[2] == "no") {
l_delta <- extract_from_ll_wcsl(ll_wcsl, "delta.lim")
l_wcsl <- extract_from_ll_wcsl(ll_wcsl, "wcs.lim")
if (nlevels(d_dat[, batch_vbl]) > 1) {
m_delta <- extract_wc_x(l1 = l_delta, l2 = l_wc_batch)
m_wcsl <- extract_wc_x(l1 = l_wcsl, l2 = l_wc_batch)
} else {
m_delta <- matrix(NA, nrow = length(rl), ncol = length(l_wc_batch))
colnames(m_delta) <- names(l_wc_batch)
m_wcsl <- m_delta
m_delta[, "dids"] <- as.numeric(l_delta[["dids"]])
m_wcsl[, "dids"] <- as.numeric(l_wcsl[["dids"]])
}
} else {
l_delta_orig <- extract_from_ll_wcsl(ll_wcsl, "delta.lim.orig")
l_wcsl_orig <- extract_from_ll_wcsl(ll_wcsl, "wcs.lim.orig")
if (nlevels(d_dat[, batch_vbl]) > 1) {
m_delta <- extract_wc_x(l1 = l_delta_orig, l2 = l_wc_batch)
m_wcsl <- extract_wc_x(l1 = l_wcsl_orig, l2 = l_wc_batch)
} else {
m_delta <- matrix(NA, nrow = length(rl), ncol = length(l_wc_batch))
colnames(m_delta) <- names(l_wc_batch)
m_wcsl <- m_delta
m_delta[, "dids"] <- as.numeric(l_delta_orig[["dids"]])
m_wcsl[, "dids"] <- as.numeric(l_wcsl_orig[["dids"]])
}
}
# ---------
# Summary data frame compilation
d_poi <- data.frame(
Intercept.cics = wc_icpt_argpm[, "cics"],
Intercept.dics = wc_icpt_argpm[, "dics"],
Intercept.dids = wc_icpt_argpm[, "dids"],
Intercept.dids.pmse = wc_icpt_argpm[, "dids.pmse"],
Delta.cics = m_delta[, "cics"],
Delta.dics = m_delta[, "dics"],
Delta.dids = m_delta[, "dids"],
Delta.dids.pmse = m_delta[, "dids.pmse"],
WCSL.cics = m_wcsl[, "cics"],
WCSL.dics = m_wcsl[, "dics"],
WCSL.dids = m_wcsl[, "dids"],
WCSL.dids.pmse = m_wcsl[, "dids.pmse"],
Exp.Spec.Report = rep(limits_list[["sl.orig"]], nrow(m_poi)),
Exp.Spec = rep(limits_list[["sl"]], nrow(m_poi)),
Rel.Spec.Report = limits_list[["rl.orig"]],
Rel.Spec = limits_list[["rl"]],
Shelf.Life.cics = m_poi[, "cics"],
Shelf.Life.dics = m_poi[, "dics"],
Shelf.Life.dids = m_poi[, "dids"],
Shelf.Life.dids.pmse = m_poi[, "dids.pmse"],
POI.Model.cics = rep(poi_ich["cics"], nrow(m_poi)),
POI.Model.dics = rep(poi_ich["dics"], nrow(m_poi)),
POI.Model.dids = rep(poi_ich["dids"], nrow(m_poi)),
POI.Model.dids.pmse = rep(poi_ich["dids.pmse"], nrow(m_poi)))
if (xform[2] != "no") {
d_poi[, "Exp.Spec"] <- rep(limits_list[["sl.bt"]], nrow(m_poi))
d_poi[, "Rel.Spec"] <- limits_list[["rl.bt"]]
}
rownames(d_poi) <- NULL
# ---------
# Compile and return results
return(list(wc.icpt = wc_icpt_argpm,
POI = d_poi))
}
#' Print value(s)
#'
#' The function \code{print_val()} generates a character string for the purpose
#' to print a value on a plot (together with associated information).
#'
#' @param val_name A character string that specifies the text preceding the
#' value of the parameter to be displayed.
#' @param val_value A numeric value that specifies the value of the parameter
#' to be displayed.
#' @param val_unit A character string that specifies the text following the
#' value of the parameter to be displayed.
#' @param val_sf A positive integer that specifies the number of significant
#' figures for the display of the limit.
#' @param prefix A character string at the beginning of the whole text. The
#' default is an empty string, i.e. \code{""}.
#' @param suffix A character string at the end of the whole text. The default
#' is an empty string, i.e. \code{""}.
#'
#' @details The function \code{print_val()} generates a character string that
#' is based on the provided information. The string is used as label of a
#' corresponding graph element. For the number formatting, the
#' \code{\link[base]{sprintf}()} function from the \sQuote{\code{base}} package
#' is used. For concatenation of the various elements, the
#' \code{\link[base]{paste}()} function from the \sQuote{\code{base}} package
#' is used.
#'
#' @return A single character string of the form \dQuote{val_name: val_value
#' (with the number of specified decimal places) val_unit}.
#'
#' @seealso \code{\link{get_text_annotation}}, \code{\link{get_n_whole_part}},
#' \code{\link[base]{formatC}}, \code{\link[base]{paste}}.
#'
#' @keywords internal
#' @noRd
print_val <- function(val_name, val_value, val_unit, val_sf,
prefix = "", suffix = "") {
if (!is.character(val_name)) {
stop("The parameter val_name must be a string.")
}
if (length(val_value) > 1) {
stop("The parameter val_value must be a numeric value of length 1.")
}
if (!is.numeric(val_value) && !is.na(val_value)) {
stop("The parameter val_value must be a numeric value of length 1.")
}
if (!is.character(val_unit)) {
stop("The parameter val_unit must be a string.")
}
if (!is.numeric(val_sf) || length(val_sf) > 1) {
stop("The parameter val_sf must be a positive integer of length 1.")
}
if (val_sf != as.integer(val_sf)) {
stop("The parameter val_sf must be a positive integer of length 1.")
}
if (val_sf < 0) {
stop("The parameter val_sf must be a positive integer of length 1.")
}
if (!is.character(prefix)) {
stop("The parameter prefix must be a string.")
}
if (!is.character(suffix)) {
stop("The parameter suffix must be a string.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Formatting of information
if (is.na(val_value)) {
paste(prefix, val_name,
formatC(signif(val_value, val_sf), digits = val_sf, format = "fg",
flag = "#"),
val_unit, suffix, sep = "")
} else {
if (val_sf <= get_n_whole_part(val_value)) {
paste(prefix, val_name,
formatC(signif(val_value, val_sf), digits = val_sf, format = "fg"),
val_unit, suffix, sep = "")
} else {
paste(prefix, val_name,
formatC(signif(val_value, val_sf), digits = val_sf, format = "fg",
flag = "#"),
val_unit, suffix, sep = "")
}
}
}
#' Get number of digits of whole part (of a decimal number)
#'
#' The function \code{get_n_whole_part()} counts the number of digits of the
#' whole number portion of a decimal number.
#'
#' @param x A decimal number (or an integer) or a vector of decimal numbers (or
#' of integers).
#'
#' @details The function \code{get_n_whole_part()} counts the number of digits
#' of the whole number portion of a decimal number.
#'
#' @return An integer representing the number of digits of the whole number
#' portion of the decimal number that was handed over.
#'
#' @seealso \code{\link{print_val}}, \code{\link{set_limits}},
#' \code{\link{get_text_annotation}}.
#'
#' @keywords internal
#' @noRd
get_n_whole_part <- function(x) {
if (!is.numeric(x) && all(!is.na(x))) {
stop("The parameter x must be a numeric value or NA.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Determine n
check <- function(xx) {
if (is.na(xx)) {
NA
} else {
if (xx <= 1 && xx >= -1) {
1
} else {
floor(log10(abs(xx))) + 1
}
}
}
vapply(x, function(xx) check(xx), numeric(1))
}
#' Determine the level of nesting of a list
#'
#' The function \code{get_n_list_levels()} determines the number of levels of
#' a nested list.
#'
#' @param x A list.
#'
#' @details The function \code{get_n_list_levels()} determines the number of
#' levels of a (nested) list.
#'
#' @return An integer representing the number of levels of the list. If an
#' object is passed on to \code{x} that is not a list \code{0} is returned.
#'
#' @seealso \code{\link{extract_from_ll_wcsl}}, \code{\link{extract_wc_x}},
#' \code{\link{get_wisle_poi_list}}.
#'
#' @keywords internal
#' @noRd
get_n_list_levels <- function(x) {
if (is.list(x)) {
1L + max(vapply(x, get_n_list_levels, numeric(1)))
} else {
0L
}
}
#' Prepare text annotation
#'
#' The function \code{get_text_annotation()} prepares a data frame for putting
#' text on a plot prepared by the \code{ggplot()} function from the
#' \sQuote{\code{ggplot2}} package.
#'
#' @param model An \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_Wisle}} object, i.e. a list returned
#' by the \code{\link{expirest_osle}()} or by the
#' \code{\link{expirest_wisle}()} function.
#' @param rvu A character string that specifies the unit associated with the
#' response variable.
#' @param x_range A numeric vector of the form \code{c(min, max)} that
#' specifies the range of the time variable to be plotted.
#' @param y_range A numeric vector of the form \code{c(min, max)} that
#' specifies the range of the response variable to be plotted.
#' @param rl_index A positive integer that specifies which of the release limit
#' values that have been handed over to \code{\link{expirest_wisle}()} should
#' be displayed. The default is \code{NULL}.
#' @param plot_option A character string of either \code{"full"}, \code{"lean"},
#' \code{"lean1"}, \code{"lean2"}, \code{"basic1"} and \code{"basic2"},
#' that specifies if additional information should be shown in the plot
#' (option \code{"full"}) or only basic information (options \code{"lean"}
#' and \code{"basic"}). Full means the data points, the fitted regression
#' line with the confidence or prediction interval, the specification
#' limit(s) and the estimated shelf life. For \sQuote{expirest_osle} objects,
#' only the options \code{"full"} and \code{"lean"} are relevant. The default
#' is \code{"full"}.
#' @inheritParams plot_expirest_osle
#'
#' @details The function \code{get_text_annotation()} expects various pieces
#' of information characterising an \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_wisle}} model. With this information, the function
#' prepares a data frame that that is used by the functions
#' \code{\link{plot_expirest_osle}()} or \code{\link{plot_expirest_wisle}(}))
#' to put text annotations on the graph that is prepared by these functions.
#'
#' @return A data frame with the columns \sQuote{Time}, \sQuote{Response},
#' Label and Colour is returned, where the column names \sQuote{Time} and
#' \sQuote{Response} are placeholders for the corresponding variable names.
#' If \code{model} is an \sQuote{expirest_osle} object, the data frame has up
#' to three rows representing the relevant specification limit(s) (row 1 or
#' rows 1 and 2) and the POI obtained from ordinary shelf life estimation (row
#' 2 or 3). If \code{model} is an \sQuote{expirest_wisle} object, the data
#' frame has up to seven rows representing the relevant specification limit(s),
#' the worst case scenario limit (row 2 or 3), the intercept (row 3 or 4), the
#' POI of the worst case scenario model (row 4 or 5), the POI obtained from
#' ordinary shelf life estimation (row 5 or 6) and the release limit (row 6
#' or 7).
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}},
#' \code{\link{print_val}}, \code{\link[ggplot2]{ggplot}},
#' \code{\link[ggplot2]{geom_text}}.
#'
#' @keywords internal
#' @noRd
get_text_annotation <- function(model, rvu, x_range, y_range, rl_index = NULL,
plot_option = "full", mtbs = "verified") {
if (!inherits(model, "expirest_osle") && !inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_osle or ",
"expirest_wisle.")
}
if (!is.character(rvu)) {
stop("The parameter rvu must be a string.")
}
if (!is.null(x_range)) {
if (!is.numeric(x_range) || length(x_range) != 2) {
stop("The parameter x_range must be a numeric vector of length 2.")
}
}
if (!is.null(y_range)) {
if (!is.numeric(y_range) || length(y_range) != 2) {
stop("The parameter y_range must be a numeric vector of length 2.")
}
}
if (!is.null(rl_index)) {
if (!is.numeric(rl_index) || length(rl_index) > 1) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
if (rl_index != as.integer(rl_index)) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
}
if (!(mtbs %in% c("verified", "cics", "dics", "dids", "dids.pmse"))) {
stop("Please specify mtbs either as \"verified\", \"cics\", \"dics\", ",
"\"dids\" or \"dids.pmse\".")
}
if (!(plot_option %in%
c("full", "lean", "lean1", "lean2", "basic1", "basic2"))) {
stop("Please specify plot_option either as \"full\", \"lean\", ",
"\"lean1\", \"lean2\", \"basic1\" or \"basic2\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
# General parameters
expob <- model
t_exp <- expob[["POI"]]
sl <- expob[["Limits"]][["sl"]]
if (expob[["Limits"]][["sf.option"]] == "tight") {
sl_sf <- expob[["Limits"]][["sl.sf"]]
} else {
sl_sf <- expob[["Limits"]][["sl.sf"]] + 1
}
xform <- expob[["Limits"]][["xform"]]
ivl_side <- expob[["Parameters"]][["ivl.side"]]
# Set model_name to dids if it is n.a.
model_name <- ifelse(mtbs == "verified",
expob[["Model.Type"]][["type.acronym"]],
mtbs)
model_name <- ifelse(model_name == "n.a.", "dids", model_name)
# Setting of y breaks
y_breaks <- pretty(y_range, 5)
switch(class(model),
"expirest_osle" = {
poi_model <- t_exp[model_name]
},
"expirest_wisle" = {
if (expob[["Limits"]][["sf.option"]] == "tight") {
rl_sf <- expob[["Limits"]][["rl.sf"]]
} else {
rl_sf <- expob[["Limits"]][["rl.sf"]] + 1
}
poi_model_name <- paste0("POI.Model.", model_name)
sl_model_name <- paste0("Shelf.Life.", model_name)
wcsl_model_name <- paste0("WCSL.", model_name)
poi_model <- t_exp[rl_index, poi_model_name]
poi_woca <- t_exp[rl_index, sl_model_name]
wc_icpt <- expob[["wc.icpt"]][, model_name]
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
switch(
class(model),
"expirest_osle" = {
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Display of text elements for plot_expirest_osle()
# The rows in data frame d_text have the following meaning and position
# (position in brackets):
# LSL (lower right), USL (upper right), POI model (low at poi.model)
if (length(sl) == 2) {
d_text <- data.frame(
Time = c(rep(x_range[2], 2), poi_model),
Response = c(sl, sl[1]),
Label = c(print_val("LSL: ", sl[1], rvu, sl_sf[1]),
print_val("USL: ", sl[2], rvu, sl_sf[2]),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1)),
Colour = c("black", "black", "forestgreen"),
stringsAsFactors = FALSE)
d_text$Response <- d_text$Response +
rep(diff(y_breaks[1:2]), 3) * 1 / c(-10, 10, -2)
} else {
switch(ivl_side,
"lower" = {
d_text <- data.frame(
Time = c(x_range[2], poi_model),
Response = rep(sl, 2),
Label =
c(print_val("LSL: ", sl, rvu, sl_sf),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1)),
Colour = c("black", "forestgreen"),
stringsAsFactors = FALSE)
d_text$Response <- d_text$Response +
rep(diff(y_breaks[1:2]), 2) * 1 / c(-10, -2)
},
"upper" = {
d_text <- data.frame(
Time = c(x_range[2], poi_model),
Response = rep(sl, 2),
Label =
c(print_val("USL: ", sl, rvu, sl_sf),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1)),
Colour = c("black", "forestgreen"),
stringsAsFactors = FALSE)
d_text$Response <- d_text$Response +
rep(diff(y_breaks[1:2]), 2) * 1 / c(10, 2)
})
}
}, "expirest_wisle" = {
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Display of text elements for plot_expirest_wisle()
# The rows in data frame d_text have the following meaning and position
# (position in brackets):
# LSL (lower right), USL (upper right),
# WCSL (left, at RL), Intercept (left, at intercept),
# POI worst case (low at poi.woca), POI model (low at poi.model)
# RL (lower right or upper right)
if (length(sl) == 2) {
d_text <- data.frame(
Time = c(rep(x_range[2], 2), 0, 0, poi_woca, poi_model),
Response = c(sl, t_exp[rl_index, wcsl_model_name], wc_icpt[rl_index],
rep(sl[1], 2)),
Label = c(print_val("LSL: ", sl[1], rvu, sl_sf[1]),
print_val("USL: ", sl[2], rvu, sl_sf[2]),
print_val("", t_exp[rl_index, wcsl_model_name], rvu,
rl_sf[rl_index], suffix = " "),
print_val("", wc_icpt[rl_index], rvu, rl_sf[rl_index]),
ifelse(plot_option %in% "lean2",
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1),
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1,
suffix = "\n(worst case\nscenario)")),
ifelse(plot_option %in% "lean2",
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1,
suffix = "\n(standard\nscenario)"))),
Colour = c("black", "black", "red", "royalblue", "forestgreen",
"grey50"),
stringsAsFactors = FALSE)
switch(ivl_side,
"lower" = {
d_text <- rbind(d_text, d_text[nrow(d_text), ])
d_text[nrow(d_text), "Time"] <- x_range[2]
d_text[nrow(d_text), "Response"] <- t_exp[rl_index, "Rel.Spec"]
d_text[nrow(d_text), "Label"] <-
print_val("LRL: ", t_exp[rl_index, "Rel.Spec"], rvu,
rl_sf[rl_index])
d_text[nrow(d_text), "Colour"] <- "grey0"
d_text$Response <- d_text$Response +
c(rep(diff(y_breaks[1:2]), 2), 0, 0,
rep(diff(y_breaks[1:2]), 2),
diff(y_breaks[1:2])) * 1 / c(-10, 10, 1, 1, -2, -2, -10)
},
"upper" = {
d_text <- rbind(d_text, d_text[nrow(d_text), ])
d_text[nrow(d_text), "Time"] <- x_range[2]
d_text[nrow(d_text), "Response"] <- t_exp[rl_index, "Rel.Spec"]
d_text[nrow(d_text), "Label"] <-
print_val("URL: ", t_exp[rl_index, "Rel.Spec"], rvu,
rl_sf[rl_index])
d_text[nrow(d_text), "Colour"] <- "grey0"
d_text[5:6, "Response"] <- rep(sl[2], 2)
d_text$Response <- d_text$Response +
c(rep(diff(y_breaks[1:2]), 2), 0, 0,
rep(diff(y_breaks[1:2]), 2),
diff(y_breaks[1:2])) * 1 / c(10, 10, 1, 1, 2, 2, 1)
})
} else {
switch(ivl_side,
"lower" = {
d_text <- data.frame(
Time = c(x_range[2], 0, 0, poi_woca, poi_model, x_range[2]),
Response = c(sl, t_exp[rl_index, wcsl_model_name],
wc_icpt[rl_index], rep(sl, 2),
t_exp[rl_index, "Rel.Spec"]),
Label =
c(print_val("LSL: ", sl, rvu, sl_sf),
print_val("", t_exp[rl_index, wcsl_model_name], rvu,
rl_sf[rl_index], suffix = " "),
print_val("", wc_icpt[rl_index], rvu, rl_sf[rl_index],
suffix = " "),
ifelse(plot_option %in% "lean2",
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1),
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1,
suffix = "\n(worst case\nscenario)")),
ifelse(plot_option %in% "lean2",
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1,
suffix = "\n(standard\nscenario)")),
print_val("LRL: ", t_exp[rl_index, "Rel.Spec"], rvu,
rl_sf[rl_index])),
Colour = c("black", "red", "royalblue", "forestgreen",
"grey50", "grey0"),
stringsAsFactors = FALSE)
d_text$Response <- d_text$Response +
c(diff(y_breaks[1:2]), 0, 0,
rep(diff(y_breaks[1:2]), 2),
diff(y_breaks[1:2])) * 1 / c(-10, 1, 1, -2, -2, -10)
},
"upper" = {
d_text <- data.frame(
Time = c(x_range[2], 0, 0, poi_woca, poi_model, x_range[2]),
Response = c(sl, t_exp[rl_index, wcsl_model_name],
wc_icpt[rl_index], rep(sl, 2),
t_exp[rl_index, "Rel.Spec"]),
Label =
c(print_val("USL: ", sl, rvu, sl_sf),
print_val("", t_exp[rl_index, wcsl_model_name], rvu,
rl_sf[rl_index], suffix = " "),
print_val("", wc_icpt[rl_index], rvu, rl_sf[rl_index],
suffix = " "),
ifelse(plot_option %in% "lean2",
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1),
print_val("", poi_woca, "",
get_n_whole_part(poi_woca) + 1,
suffix = "\n(worst case\nscenario)")),
ifelse(plot_option %in% "lean2",
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1),
print_val("", poi_model, "",
get_n_whole_part(poi_model) + 1,
suffix = "\n(standard\nscenario)")),
print_val("URL: ", t_exp[rl_index, "Rel.Spec"], rvu,
rl_sf[rl_index])),
Colour = c("black", "red", "royalblue", "forestgreen",
"grey50", "grey0"),
stringsAsFactors = FALSE)
d_text$Response <- d_text$Response +
c(diff(y_breaks[1:2]), 0, 0,
rep(diff(y_breaks[1:2]), 2),
diff(y_breaks[1:2])) * 1 / c(10, 1, 1, 2, 2, 10)
})
}
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Renaming of columns for plotting on the original scale
if (sum(xform %in% "no") == 2) {
colnames(d_text) <- c(expob[["Variables"]][["time"]],
expob[["Variables"]][["response"]],
"Label", "Colour")
}
if (sum(xform %in% "no") == 0) {
colnames(d_text) <- c(expob[["Variables"]][["time.orig"]],
expob[["Variables"]][["response.orig"]],
"Label", "Colour")
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
colnames(d_text) <- c(expob[["Variables"]][["time.orig"]],
expob[["Variables"]][["response"]],
"Label", "Colour")
}
if (xform[2] != "no") {
colnames(d_text) <- c(expob[["Variables"]][["time"]],
expob[["Variables"]][["response.orig"]],
"Label", "Colour")
}
}
return(d_text)
}
#' Prepare horizontal lines
#'
#' The function \code{get_hlines()} prepares a data frame for putting
#' horizontal lines on a plot prepared by the \code{ggplot()} function from
#' the \sQuote{\code{ggplot2}} package.
#'
#' @inheritParams get_text_annotation
#'
#' @details The function \code{get_hlines()} expects various pieces
#' of information characterising an \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_wisle}} model. Based on the information provided,
#' the function prepares a data frame that that is used by the functions
#' \code{\link{plot_expirest_osle}()} or \code{\link{plot_expirest_wisle}()}
#' to put lines on the graph that is prepared by these functions.
#'
#' @return A data frame with the columns \sQuote{Response}, Item, Colour and
#' Type is returned, where the column name \sQuote{Response} is a placeholder
#' for the corresponding variable name. The data frame has up to two rows
#' representing the specification limit(s).
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}},
#' \code{\link[ggplot2]{ggplot}}, \code{\link[ggplot2]{geom_text}}.
#'
#' @keywords internal
#' @noRd
get_hlines <- function(model) {
if (!inherits(model, "expirest_osle") && !inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_osle or ",
"expirest_wisle.")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
expob <- model
sl <- expob[["Limits"]][["sl"]]
xform <- expob[["Limits"]][["xform"]]
ivl_side <- expob[["Parameters"]][["ivl.side"]]
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (length(sl) == 2) {
d_hlines <- data.frame(Response = sl,
Item = c("LSL", "USL"),
Colour = as.character(c("black", "black")),
Type = as.character(c("dotted", "dotted")),
stringsAsFactors = FALSE)
} else {
switch(ivl_side,
"lower" = {
d_hlines <- data.frame(Response = sl,
Item = c("LSL"),
Colour = as.character(c("black")),
Type = as.character(c("dotted")),
stringsAsFactors = FALSE)
},
"upper" = {
d_hlines <- data.frame(Response = sl,
Item = c("USL"),
Colour = as.character(c("black")),
Type = as.character(c("dotted")),
stringsAsFactors = FALSE)
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Renaming of columns for plotting on the original scale
if (sum(xform %in% "no") == 2) {
colnames(d_hlines) <- c(expob[["Variables"]][["response"]],
"Item", "Colour", "Type")
}
if (sum(xform %in% "no") == 0) {
colnames(d_hlines) <- c(expob[["Variables"]][["response.orig"]],
"Item", "Colour", "Type")
}
if (sum(xform %in% "no") == 1) {
if (xform[2] != "no") {
colnames(d_hlines) <- c(expob[["Variables"]][["response.orig"]],
"Item", "Colour", "Type")
} else {
colnames(d_hlines) <- c(expob[["Variables"]][["response"]],
"Item", "Colour", "Type")
}
}
return(d_hlines)
}
#' Prepare vertical lines
#'
#' The function \code{get_vlines()} prepares a data frame for putting
#' vertical lines on a plot prepared by the \code{ggplot()} function from
#' the \sQuote{\code{ggplot2}} package.
#'
#' @inheritParams get_text_annotation
#'
#' @details The function \code{get_vlines()} expects various pieces
#' of information characterising an \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_wisle}} model. Based on the information provided,
#' the function prepares a data frame that that is used by the functions
#' \code{\link{plot_expirest_osle}()} or \code{\link{plot_expirest_wisle}(}))
#' to put text annotations on the graph that is prepared by these functions.
#'
#' @return A data frame with the columns \sQuote{Time}, Item, Colour and Type
#' is returned, where the column name \sQuote{Time} is a placeholder for the
#' corresponding variable name. If \code{model} is an \sQuote{expirest_osle}
#' object, the data frame has one row representing the POI obtained from
#' ordinary shelf life estimation. If \code{model} is an \sQuote{expirest_wisle}
#' object, the data frame has two rows representing the POI obtained from
#' ordinary shelf life estimation (row 1) and the POI of the What-if approach
#' for shelf life estimation (row 2).
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}},
#' \code{\link[ggplot2]{ggplot}}, \code{\link[ggplot2]{geom_text}}.
#'
#' @keywords internal
#' @noRd
get_vlines <- function(model, rl_index = NULL, mtbs = "verified") {
if (!inherits(model, "expirest_osle") && !inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_osle or ",
"expirest_wisle.")
}
if (!is.null(rl_index)) {
if (!is.numeric(rl_index) || length(rl_index) > 1) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
if (rl_index != as.integer(rl_index)) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
}
if (!(mtbs %in% c("verified", "cics", "dics", "dids", "dids.pmse"))) {
stop("Please specify mtbs either as \"verified\", \"cics\", \"dics\", ",
"\"dids\" or \"dids.pmse\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
expob <- model
t_exp <- expob[["POI"]]
xform <- expob[["Limits"]][["xform"]]
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of models and of the model type
# Set model_name to dids if it is n.a.
model_name <- ifelse(mtbs == "verified",
expob[["Model.Type"]][["type.acronym"]],
mtbs)
model_name <- ifelse(model_name == "n.a.", "dids", model_name)
switch(class(model),
"expirest_osle" = {
poi_model <- t_exp[model_name]
},
"expirest_wisle" = {
# Most appropriate model based on the ANCOVA analysis, or as
# specified via the mtbs parameter
poi_model_name <- paste0("POI.Model.", model_name)
sl_model_name <- paste0("Shelf.Life.", model_name)
# POI with the upper or lower confidence or prediction interval of
# the linear regression model representing the worst case scenario
# (woca) case
poi_model <- t_exp[rl_index, poi_model_name]
poi_woca <- t_exp[rl_index, sl_model_name]
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Setting variable names
if (sum(xform %in% "no") == 2) {
time_vbl <- expob[["Variables"]][["time"]]
}
if (sum(xform %in% "no") == 0) {
old_time_vbl <- expob[["Variables"]][["time.orig"]]
time_vbl <- expob[["Variables"]][["time"]]
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
old_time_vbl <- expob[["Variables"]][["time.orig"]]
time_vbl <- expob[["Variables"]][["time"]]
}
if (xform[2] != "no") {
time_vbl <- expob[["Variables"]][["time"]]
}
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
switch(class(model),
"expirest_osle" = {
d_vlines <- data.frame(Time = c(poi_model),
Item = c("poi.model"),
Colour = c("forestgreen"),
Type = c("dotdash"),
stringsAsFactors = FALSE)
},
"expirest_wisle" = {
d_vlines <- data.frame(Time = c(poi_woca, poi_model),
Item = c("poi.woca", "poi.model"),
Colour = c("forestgreen", "grey50"),
Type = c("dashed", "dotdash"),
stringsAsFactors = FALSE)
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Renaming of columns for plotting on the original scale
if (sum(xform %in% "no") == 2) {
colnames(d_vlines) <- c(time_vbl, "Item", "Colour", "Type")
}
if (sum(xform %in% "no") == 0) {
colnames(d_vlines) <- c(old_time_vbl, "Item", "Colour", "Type")
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
colnames(d_vlines) <- c(old_time_vbl, "Item", "Colour", "Type")
} else {
colnames(d_vlines) <- c(time_vbl, "Item", "Colour", "Type")
}
}
return(d_vlines)
}
#' Prepare segments explaining graphical elements
#'
#' The function \code{get_segments()} prepares a data frame for putting
#' segments on a plot prepared by the \code{ggplot()} function from the
#' \sQuote{\code{ggplot2}} package.
#'
#' @inheritParams get_text_annotation
#'
#' @details The function \code{get_segments()} expects various pieces
#' of information characterising an \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_wisle}} model. Based on the information provided,
#' the function prepares a data frame that that is used by the function
#' \code{\link{plot_expirest_wisle}(})) to put line segments on the graph that
#' is prepared by this function.
#'
#' @return A data frame with the columns \sQuote{Time.1}, \sQuote{Time.2},
#' \sQuote{Response.1}, \sQuote{Response.2}, Item, Colour, Type and Size is
#' returned, where the column names \sQuote{Time.1}, \sQuote{Time.2},
#' \sQuote{Response.1} and \sQuote{Response.2} are placeholders for the
#' corresponding variable names. The items in the four rows represent the
#' maximal allowed difference over time from the intercept (shown as red
#' dashed horizontal line), the release limit (shown as black dotted horizontal
#' line), the maximal allowed difference over time from the specification limit
#' (shown as grey solid vertical line) and the maximal allowed difference over
#' time from the intercept (shown as grey solid vertical line).
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}},
#' \code{\link[ggplot2]{ggplot}}, \code{\link[ggplot2]{geom_text}}.
#'
#' @keywords internal
#' @noRd
get_segments <- function(model, x_range, rl_index, mtbs = "verified") {
if (!inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_wisle.")
}
if (!is.null(x_range)) {
if (!is.numeric(x_range) || length(x_range) != 2) {
stop("The parameter x_range must be a numeric vector of length 2.")
}
}
if (!is.null(rl_index)) {
if (!is.numeric(rl_index) || length(rl_index) > 1) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
if (rl_index != as.integer(rl_index)) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
}
if (!(mtbs %in% c("verified", "cics", "dics", "dids", "dids.pmse"))) {
stop("Please specify mtbs either as \"verified\", \"cics\", \"dics\", ",
"\"dids\" or \"dids.pmse\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
# General parameters
expob <- model
t_exp <- expob[["POI"]]
sl <- expob[["Limits"]][["sl"]]
rl <- expob[["Limits"]][["rl"]]
xform <- expob[["Limits"]][["xform"]]
ivl_side <- expob[["Parameters"]][["ivl.side"]]
# Set model_name to dids if it is n.a.
model_name <- ifelse(mtbs == "verified",
expob[["Model.Type"]][["type.acronym"]],
mtbs)
model_name <- ifelse(model_name == "n.a.", "dids", model_name)
sl_model_name <- paste0("Shelf.Life.", model_name)
wcsl_model_name <- paste0("WCSL.", model_name)
poi_woca <- t_exp[rl_index, sl_model_name]
wc_icpt <- expob[["wc.icpt"]][, model_name]
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (length(sl) == 2) {
switch(ivl_side,
"lower" = {
d_seg <- data.frame(
Time.1 =
c(0, 0,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Time.2 =
c(poi_woca, x_range[2],
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Response.1 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], sl[1],
t_exp[rl_index, wcsl_model_name]),
Response.2 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], rl[rl_index],
wc_icpt[rl_index]),
Item = c("x.delta", "x.delta.shifted", "y.delta",
"y.delta.shifted"),
Colour = c("red", "grey0", "grey50", "grey50"),
Type = c("dashed", "dotted", rep("solid", 2)),
Size = c(rep(0.5, 2), rep(1, 2)),
stringsAsFactors = FALSE)
},
"upper" = {
d_seg <- data.frame(
Time.1 =
c(0, 0,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Time.2 =
c(poi_woca, x_range[2],
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Response.1 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], sl[2],
t_exp[rl_index, wcsl_model_name]),
Response.2 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], rl[rl_index],
wc_icpt[rl_index]),
Item = c("x.delta", "x.delta.shifted", "y.delta",
"y.delta.shifted"),
Colour = c("red", "grey0", "grey50", "grey50"),
Type = c("dashed", "dotted", rep("solid", 2)),
Size = c(rep(0.5, 2), rep(1, 2)),
stringsAsFactors = FALSE)
})
} else {
d_seg <- data.frame(
Time.1 = c(0, 0,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Time.2 = c(poi_woca, x_range[2],
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 3,
-round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9),
Response.1 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], sl,
t_exp[rl_index, wcsl_model_name]),
Response.2 =
c(t_exp[rl_index, wcsl_model_name], rl[rl_index], rl[rl_index],
wc_icpt[rl_index]),
Item = c("x.delta", "x.delta.shifted", "y.delta", "y.delta.shifted"),
Colour = c("red", "grey0", "grey50", "grey50"),
Type = c("dashed", "dotted", rep("solid", 2)),
Size = c(rep(0.5, 2), rep(1, 2)),
stringsAsFactors = FALSE)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Renaming of columns for plotting on the original scale
if (sum(xform %in% "no") == 2) {
colnames(d_seg) <-
c(paste(expob[["Variables"]][["time"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response"]], 1:2, sep = "."),
"Item", "Colour", "Type", "Size")
}
if (sum(xform %in% "no") == 0) {
colnames(d_seg) <-
c(paste(expob[["Variables"]][["time.orig"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response.orig"]], 1:2, sep = "."),
"Item", "Colour", "Type", "Size")
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
colnames(d_seg) <-
c(paste(expob[["Variables"]][["time.orig"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response"]], 1:2, sep = "."),
"Item", "Colour", "Type", "Size")
}
if (xform[2] != "no") {
colnames(d_seg) <-
c(paste(expob[["Variables"]][["time"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response.orig"]], 1:2, sep = "."),
"Item", "Colour", "Type", "Size")
}
}
return(d_seg)
}
#' Prepare arrow supporting the explanation of graphical elements
#'
#' The function \code{get_arrows()} prepares a data frame for putting an
#' arrow on a plot prepared by the \code{ggplot()} function from the
#' \sQuote{\code{ggplot2}} package.
#'
#' @inheritParams get_segments
#'
#' @details The function \code{get_arrows()} expects various pieces
#' of information that characterises an \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_wisle}} model. Based on the information provided,
#' the function prepares a data frame that that is used by the function
#' \code{\link{plot_expirest_wisle}(})) to put an arrow on the graph that
#' is prepared by this function.
#'
#' @return A data frame with the columns \sQuote{Time.1}, \sQuote{Time.2},
#' \sQuote{Response.1}, \sQuote{Response.2}, Item, Colour, Line.Type,
#' Arrow.Type, Size, Curvature, Angle and Length is returned, where the column
#' names \sQuote{Time.1}, \sQuote{Time.2}, \sQuote{Response.1} and
#' \sQuote{Response.2} are placeholders for the corresponding variable names.
#' The data frame has a single row representing the arrow that is put on the
#' graphical illustration.
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}},
#' \code{\link[ggplot2]{ggplot}}, \code{\link[ggplot2]{geom_text}}.
#'
#' @keywords internal
#' @noRd
get_arrow <- function(model, x_range, rl_index, mtbs = "verified") {
if (!inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_wisle.")
}
if (!is.null(x_range)) {
if (!is.numeric(x_range) || length(x_range) != 2) {
stop("The parameter x_range must be a numeric vector of length 2.")
}
}
if (!is.null(rl_index)) {
if (!is.numeric(rl_index) || length(rl_index) > 1) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
if (rl_index != as.integer(rl_index)) {
stop("The parameter rl_index must be a positive integer of length 1.")
}
}
if (!(mtbs %in% c("verified", "cics", "dics", "dids", "dids.pmse"))) {
stop("Please specify mtbs either as \"verified\", \"cics\", \"dics\", ",
"\"dids\" or \"dids.pmse\".")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
# General parameters
expob <- model
t_exp <- expob[["POI"]]
sl <- expob[["Limits"]][["sl"]]
rl <- expob[["Limits"]][["rl"]]
xform <- expob[["Limits"]][["xform"]]
ivl_side <- expob[["Parameters"]][["ivl.side"]]
# Set model_name to dids if it is n.a.
model_name <- ifelse(mtbs == "verified",
expob[["Model.Type"]][["type.acronym"]],
mtbs)
model_name <- ifelse(model_name == "n.a.", "dids", model_name)
sl_model_name <- paste0("Shelf.Life.", model_name)
wcsl_model_name <- paste0("WCSL.", model_name)
wc_icpt <- expob[["wc.icpt"]][, model_name]
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
if (length(sl) == 2) {
switch(ivl_side,
"lower" = {
d_arr <- data.frame(
Time.1 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9,
Time.2 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 2,
Response.1 = c((t_exp[rl_index, wcsl_model_name] +
wc_icpt[rl_index]) / 2),
Response.2 = c((sl[1] + rl[rl_index]) / 2),
Item = c("arrow"),
Colour = c("grey50"),
Line.Type = c("solid"),
Arrow.Type = c("closed"),
Size = 0.5,
Curvature = 0.5,
Angle = 90,
Length = ceiling(log2(x_range[2])),
stringsAsFactors = FALSE)
},
"upper" = {
d_arr <- data.frame(
Time.1 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9,
Time.2 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 2,
Response.1 = c((t_exp[rl_index, wcsl_model_name] +
wc_icpt[rl_index]) / 2),
Response.2 = c((sl[2] + rl[rl_index]) / 2),
Item = c("arrow"),
Colour = c("grey50"),
Line.Type = c("solid"),
Arrow.Type = c("closed"),
Size = 0.5,
Curvature = -0.5,
Angle = 90,
Length = ceiling(log2(x_range[2])),
stringsAsFactors = FALSE)
})
} else {
switch(ivl_side,
"lower" = {
d_arr <- data.frame(
Time.1 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9,
Time.2 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 2,
Response.1 = c((t_exp[rl_index, wcsl_model_name] +
wc_icpt[rl_index]) / 2),
Response.2 = c((sl + rl[rl_index]) / 2),
Item = c("arrow"),
Colour = c("grey50"),
Line.Type = c("solid"),
Arrow.Type = c("closed"),
Size = 0.5,
Curvature = 0.5,
Angle = 90,
Length = ceiling(log2(x_range[2])),
stringsAsFactors = FALSE)
},
"upper" = {
d_arr <- data.frame(
Time.1 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 9,
Time.2 = -round(sqrt(t_exp[rl_index, sl_model_name]) / 3, 0) / 2,
Response.1 = c((t_exp[rl_index, wcsl_model_name] +
wc_icpt[rl_index]) / 2),
Response.2 = c((sl + rl[rl_index]) / 2),
Item = c("arrow"),
Colour = c("grey50"),
Line.Type = c("solid"),
Arrow.Type = c("closed"),
Size = 0.5,
Curvature = -0.5,
Angle = 90,
Length = ceiling(log2(x_range[2])),
stringsAsFactors = FALSE)
})
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Renaming of columns for plotting on the original scale
if (sum(xform %in% "no") == 2) {
colnames(d_arr) <-
c(paste(expob[["Variables"]][["time"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response"]], 1:2, sep = "."),
"Item", "Colour", "Line.Type", "Arrow.Type",
"Size", "Curvature", "Angle", "Length")
}
if (sum(xform %in% "no") == 0) {
colnames(d_arr) <-
c(paste(expob[["Variables"]][["time.orig"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response.orig"]], 1:2, sep = "."),
"Item", "Colour", "Line.Type", "Arrow.Type",
"Size", "Curvature", "Angle", "Length")
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
colnames(d_arr) <-
c(paste(expob[["Variables"]][["time.orig"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response"]], 1:2, sep = "."),
"Item", "Colour", "Line.Type", "Arrow.Type",
"Size", "Curvature", "Angle", "Length")
}
if (xform[2] != "no") {
colnames(d_arr) <-
c(paste(expob[["Variables"]][["time"]], 1:2, sep = "."),
paste(expob[["Variables"]][["response.orig"]], 1:2, sep = "."),
"Item", "Colour", "Line.Type", "Arrow.Type",
"Size", "Curvature", "Angle", "Length")
}
}
return(d_arr)
}
#' Get the predicted values
#'
#' The function \code{get_predictions()} makes the predictions to be shown
#' on the graphs prepared by the function \code{\link{plot_expirest_osle}()}
#' or \code{\link{plot_expirest_wisle}()}.
#'
#' @param model An \sQuote{\code{expirest_osle}} or an
#' \sQuote{\code{expirest_Wisle}} object, i.e. a list returned
#' by the \code{\link{expirest_osle}()} or by the
#' \code{\link{expirest_wisle}()} function.
#' @param model_name A character string representing the acronym that specifies
#' which model, based on the ANCOVA analysis, suits best.
#' @param x_range A numeric vector of the form \code{c(min, max)} that
#' specifies the range of the time variable to be plotted.
#'
#' @details The function \code{get_predictions()} prepares a data frame that
#' contains the predicted values .
#'
#' @return A data frame with the columns named after the \code{batch_vbl}, the
#' \code{time_vbl} and the \code{response_vbl} (see function descriptions of
#' \code{\link{plot_expirest_osle}()} or \code{\link{plot_expirest_wisle}()})
#' and further two columns named \code{LL} and \code{UL} which represent the
#' lower or the upper confidence or prediction interval limits, respectively.
#'
#' @seealso \code{\link{plot_expirest_osle}}, \code{\link{plot_expirest_wisle}}.
#'
#' @importFrom stats predict
#'
#' @keywords internal
#' @noRd
get_predictions <- function(model, model_name, x_range) {
if (!inherits(model, "expirest_osle") && !inherits(model, "expirest_wisle")) {
stop("The model must be an object of class expirest_osle or ",
"expirest_wisle.")
}
if (!(model_name %in% c("cics", "dics", "dids", "dids.pmse"))) {
stop("Please specify model_name either as \"cics\", \"dics\", ",
"\"dids\" or \"dids.pmse\".")
}
if (!is.numeric(x_range) || length(x_range) != 2) {
stop("The parameter x_range must be a vector of length 2.")
}
if (x_range[2] < x_range[1]) {
stop("The parameter x_range must be of the form c(lower, upper).")
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Extraction of data and parameters
expob <- model
l_models <- expob[["Models"]]
model <- l_models[[model_name]]
d_dat <- expob[["Data"]]
alpha <- expob[["Parameters"]][["alpha"]]
ivl <- expob[["Parameters"]][["ivl"]]
ivl_type <- expob[["Parameters"]][["ivl.type"]]
# Note: since the predict.lm() function from the 'stats' package always
# calculates two-sided limits the value of alpha must be doubled in case that
# the ivl_type is "one-sided".
if (ivl_type == "one.sided") {
alpha <- alpha * 2
}
xform <- expob[["Limits"]][["xform"]]
shift <- expob[["Limits"]][["shift"]]
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation the prediction
# Generate the new x values (on the original scale) for prediction
x_new <- seq(min(x_range), max(x_range), length.out = 100)
# Transformation of new x values, if necessary
switch(xform[1],
"no" = {
},
"log" = {
x_new_trfmd <- log(x_new + shift[1])
},
"sqrt" = {
x_new_trfmd <- sqrt(x_new + shift[1])
},
"sq" = {
x_new_trfmd <- (x_new + shift[1])^2
})
# Creation of data frame for the predict.lm() function
t_batches <- levels(d_dat[, expob[["Variables"]][["batch"]]])
if (xform[1] == "no") {
d_new <- data.frame(rep(t_batches, each = length(x_new)),
rep(x_new, times = length(t_batches)))
} else {
d_new <- data.frame(rep(t_batches, each = length(x_new_trfmd)),
rep(x_new_trfmd, times = length(t_batches)))
}
colnames(d_new) <- c(expob[["Variables"]][["batch"]],
expob[["Variables"]][["time"]])
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Performing the prediction
if (model_name != "dids") {
m_pred <- predict(model, newdata = d_new, interval = ivl, level = 1 - alpha)
} else {
l_pred <- lapply(t_batches, function(x) {
predict(l_models$dids[[x]],
newdata = d_new[d_new[, expob[["Variables"]][["batch"]]] == x, ],
interval = ivl, level = 1 - alpha)
})
m_pred <- do.call(rbind, l_pred)
}
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Preparation of the results
# Back-transformation of predicted (response) values, if necessary
switch(xform[2],
"no" = {
},
"log" = {
m_pred <- exp(m_pred) - shift[2]
},
"sqrt" = {
m_pred <- m_pred^2 - shift[2]
},
"sq" = {
m_pred[m_pred < 0] <- NA
m_pred <- sqrt(m_pred) - shift[2]
})
# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
# Generation of data frame for plotting (with values in original scale)
if (sum(xform %in% "no") == 2) {
d_pred <- as.data.frame(cbind(d_new, m_pred))
colnames(d_pred) <- c(expob[["Variables"]][["batch"]],
expob[["Variables"]][["time"]],
expob[["Variables"]][["response"]],
"LL", "UL")
}
if (sum(xform %in% "no") == 0) {
d_pred <- data.frame(rep(t_batches, each = length(x_new)),
rep(x_new, times = length(t_batches)),
m_pred)
colnames(d_pred) <- c(expob[["Variables"]][["batch"]],
expob[["Variables"]][["time.orig"]],
expob[["Variables"]][["response.orig"]],
"LL", "UL")
}
if (sum(xform %in% "no") == 1) {
if (xform[1] != "no") {
d_pred <- data.frame(rep(t_batches, each = length(x_new)),
rep(x_new, times = length(t_batches)),
m_pred)
colnames(d_pred) <- c(expob[["Variables"]][["batch"]],
expob[["Variables"]][["time.orig"]],
expob[["Variables"]][["response"]],
"LL", "UL")
}
if (xform[2] != "no") {
d_pred <- as.data.frame(cbind(d_new, m_pred))
colnames(d_pred) <- c(expob[["Variables"]][["batch"]],
expob[["Variables"]][["time"]],
expob[["Variables"]][["response.orig"]],
"LL", "UL")
}
}
return(d_pred)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.