Skip to content

Commit

Permalink
Finish documentation
Browse files Browse the repository at this point in the history
  • Loading branch information
@jan-imbi
jan-imbi committed Sep 21, 2023
1 parent 47e6749 commit 2d1ab0f
Show file tree
Hide file tree
Showing 28 changed files with 900 additions and 183 deletions.
4 changes: 4 additions & 0 deletions NAMESPACE
View file
Expand Up @@ -2,6 +2,7 @@

S3method("[",EstimatorScoreResultList)
S3method(format,EstimatorScoreResultList)
export(AdaptivelyWeightedSampleMean)
export(Bias)
export(BiasReduced)
export(Centrality)
Expand Down Expand Up @@ -76,8 +77,11 @@ importFrom(ggplot2,ggplot)
importFrom(ggplot2,scale_color_manual)
importFrom(ggplot2,scale_fill_gradient)
importFrom(ggplot2,scale_x_continuous)
importFrom(ggpubr,ggarrange)
importFrom(ggpubr,theme_pubclean)
importFrom(ggpubr,theme_pubr)
importFrom(grDevices,xy.coords)
importFrom(graphics,plot.default)
importFrom(latex2exp,TeX)
importFrom(progressr,progressor)
importFrom(scales,percent)
Expand Down
31 changes: 18 additions & 13 deletions R/adestr_package.R
View file
@@ -1,20 +1,25 @@
#' adestr
#' @description Point estimates, confidence intervals, and p-values for optimal adaptive two-stage designs
#' @description Point estimates, confidence intervals, and p-values for optimal adaptive two-stage designs.
#'
#' @details This package implements methods to evaluate the performance characteristics of
#' various point and interval estimators for optimal adaptive two-stage designs.
#' @details This package implements methods to \link[adestr:evaluate_estimator]{evaluate the performance characteristics} of
#' various \link[adestr:PointEstimator]{point} and \link[adestr:IntervalEstimator]{interval} estimators for optimal adaptive two-stage designs.
#' Specifically, this package is written to interface with trial designs created by the \code{\link{adoptr}} package
#' \insertCite{kunzmann2021adoptr,pilz2021optimal}{adestr}.
#' Apart from the a priori evaluation of performance characteristics, this package also allows for the
#' evaluation of the implemented estimators on real datasets, and it implements methods
#' to calculate p-values.
#' \link[adestr:analyze]{calculation of the values of the estimators} given real datasets, and it implements methods
#' to calculate \link[adestr:PValue]{p-values}.
#'
#' @docType package
#' @name adestr
#' @import adoptr methods
#' @importFrom stats dnorm pnorm qnorm dt pt qt dchisq pchisq qchisq integrate uniroot var
#' @importFrom cubature hcubature
#' @importFrom Rdpack reprompt
#' @seealso \code{\link[adestr:evaluate_estimator]{evaluate_estimator}}
#' @seealso \code{\link[adestr:analyze]{analyze}}
#' @seealso \code{\link[adestr:Statistic]{Statistic}} \code{\link[adestr:PointEstimator]{PointEstimator}} \code{\link[adestr:IntervalEstimator]{IntervalEstimator}} \code{\link[adestr:PValue]{PValue}}
#' @seealso \code{\link[adestr:plot,EstimatorScoreResultList-method]{plot}} \code{\link[adestr:plot_p]{plot_p}}
#' @seealso \url{https://jan-imbi.github.io/adestr/}
#'
#' @references
#' \insertAllCited{}
Expand All @@ -25,16 +30,16 @@

.adestr_options <- list(
# Root finding inside estimators
adestr_tol_roots = 1e-3,
adestr_maxiter_roots = 1e3,
adestr_tol_roots = 2e-3,
adestr_maxiter_roots = 5e2,
# Integrals used inside estimators
adestr_tol_inner = 5e-3,
adestr_maxEval_inner = 1e4,
adestr_absError_inner = 1e-5,
adestr_tol_inner = 2e-3,
adestr_maxEval_inner = 5e2,
adestr_absError_inner = 5e-5,
# Integrals to evaluate estimators
adestr_tol_outer = 5e-6,
adestr_maxEval_outer = 3e5,
adestr_absError_outer = 1e-8
adestr_tol_outer = 1e-3,
adestr_maxEval_outer = 1e3,
adestr_absError_outer = 1e-5
)

#### Universal argument order for all functions: ####
Expand Down
53 changes: 35 additions & 18 deletions R/analyze.R
View file
Expand Up @@ -11,8 +11,18 @@ Results <- setClass("Results", slots = c(data ="data.frame",
#' \link[adestr:ConfidenceInterval]{confidence intervals},
#' and \link[adestr:PValue]{p-values} for a given dataset.
#'
#' Note that in \code{\link{adestr}}, statistics are codes as functions of the
#' stage-wise sample means (and stage-wise sample variances if data_distribution is
#' \code{\link{Student}}). In a first-step, the data is summarized to produce these
#' parameters. Then, the list of statistics are evaluated at the values of these parameters.
#'
#' The output of the \code{analyze} function also displays information on the hypothesis
#' test and the interim decision. If the \code{\link[adestr:Statistic-class]{statistics}} list is empty, this will
#' be the only information displayed.
#'
#' @param data a data.frame containing the data to be analyzed.
#' @param statistics a list of objects of class \code{\link{PointEstimator}}, \code{\link{ConfidenceInterval}} or
#' @param statistics a list of objects of class \code{\link{PointEstimator}},
#' \code{\link{ConfidenceInterval}} or
#' \code{\link{PValue}}.
#' @inheritParams evaluate_estimator
#'
Expand All @@ -21,21 +31,28 @@ Results <- setClass("Results", slots = c(data ="data.frame",
#' @export
#'
#' @examples
#' set.seed(321)
#' set.seed(123)
#' dat <- data.frame(
#' endpoint = c(rnorm(28, .2, 1), rnorm(28, 0, 1),
#' rnorm(23, .2, 1), rnorm(23, 0, 1)),
#' group = factor(rep(c("ctl", "trt", "ctl", "trt"),
#' c(28,28,23,23))),
#' stage = rep(c(1L, 2L), c(56, 46))
#' )
#' analyze(
#' data = dat,
#' statistics = c(get_example_statistics()),
#' data_distribution = Normal(),
#' sigma = 1,
#' design = get_example_design()
#' endpoint = c(rnorm(28, 0.3)),
#' stage = rep(1, 28)
#' )
#' analyze(data = dat,
#' statistics = list(),
#' data_distribution = Normal(FALSE),
#' design = get_example_design(),
#' sigma = 1)
#'
#' # The results suggest recruiting 32 patients for the second stage
#' dat <- rbind(
#' dat,
#' data.frame(
#' endpoint = rnorm(32, mean = 0.3),
#' stage = rep(2, 32)))
#' analyze(data = dat,
#' statistics = get_example_statistics(),
#' data_distribution = Normal(FALSE),
#' design = get_example_design(),
#' sigma = 1)
setGeneric("analyze", function(data,
statistics = list(),
data_distribution,
Expand Down Expand Up @@ -78,18 +95,18 @@ setMethod("analyze", signature("data.frame"),
if(data_distribution@two_armed) c(sdata$n_s1_g1, sdata$n_s1_g2) else sdata$n1,
data_distribution@two_armed)
if (length(statistics)>1) {
if (abs(sdata@n_s1_g1 - design@n1)/ (design@n1) > 0.1)
if (abs(sdata$n_s1_g1 - design@n1)/ (design@n1) > 0.1)
warning("Planned first-stage sample size in group 1 differs from actually observed sample size by more than 10%. Results may be unreliable.")
if (data_distribution@two_armed){
if (abs(sdata@n_s1_g2 - design@n1)/ (design@n1) > 0.1)
if (abs(sdata$n_s1_g2 - design@n1)/ (design@n1) > 0.1)
warning("Planned first-stage sample size in group 2 differs from actually observed sample size by more than 10%. Results may be unreliable.")
}
calc_n2 <- n2(design, test_val, round=FALSE)
if (sdata$n_stages==2L){
if (abs(sdata@n_s2_g1 - calc_n2 )/ (calc_n2) > 0.1)
if (abs(sdata$n_s2_g1 - calc_n2 )/ (calc_n2) > 0.1)
warning("Planned second-stage sample size in group 1 differs from actually observed sample size by more than 10%. Results may be unreliable.")
if (data_distribution@two_armed) {
if (abs(sdata@n_s2_g2 - calc_n2 )/ (calc_n2) > 0.1)
if (abs(sdata$n_s2_g2 - calc_n2 )/ (calc_n2) > 0.1)
warning("Planned second-stage sample size in group 2 differs from actually observed sample size by more than 10%. Results may be unreliable.") }
if (test_val > design@c1e | test_val < design@c1f)
warning("Second-stage data was recorded but trial should have been stopped at interim. Results may be unreliable.")
Expand Down
120 changes: 109 additions & 11 deletions R/estimators.R
View file
Expand Up @@ -26,11 +26,58 @@ setClass("Estimator", contains = "Statistic")

#' Point estimators
#'
#' This is the parent class for all point estimators implemented in this package.
#' Currently, only estimators for the parameter \eqn{\mu} of a normal distribution
#' are implemented.
#'
#' @details
#' Details about the point estimators can be found in (our upcoming paper).
#' ## Sample Mean (\code{SampleMean()})
#' The sample mean is the maximum likelihood estimator for the mean and probably the
#' 'most straightforward' of the implemented estimators.
#' ## Fixed weighted sample means (\code{WeightedSampleMean()})
#' The first- and second-stage (if available) sample means are combined via fixed, predefined
#' weights. See \insertCite{brannath2006estimation}{adestr} and \insertCite{@Section 8.3.2 in @wassmer2016group}{adestr}.
#' ## Adaptively weighted sample means (\code{AdaptivelyWeightedSampleMean()})
#' The first- and second-stage (if available) sample means are combined via a combination of
#' fixed and adaptively modified
#' weights that depend on the standard error.
#' See \insertCite{@Section 8.3.4 in @wassmer2016group}{adestr}.
#' ## Minimizing peak variance in adaptively weighted sample means (\code{MinimizePeakVariance()})
#' For this estimator, the weights of the adaptively weighted sample mean are chosen to
#' minimize the variance of the estimator for the value of \eqn{\mu} which maximizes
#' the expected sample size.
#' ## (Pseudo) Rao-Blackwell estimators (\code{RaoBlackwell} and \code{PseudoRaoBlackwell})
#' The conditional expectation of the first-stage sample mean given the overall sample
#' mean and the second-stage sample size. See \insertCite{emerson1997computationally}{adestr}.
#' ## A bias-reduced estimator (\code{BiasReduced()})
#' This estimator is calculated by subtracting an estimate of the bias from the MLE.
#' See \insertCite{whitehead1986bias}{adestr}.
#' ## Median-unbiased estimators
#' The implemented median-unbiased estimators are:
#' * \code{MedianUnbiasedMLEOrdering()}
#' * \code{MedianUnbiasedLikelihoodRatioOrdering()}
#' * \code{MedianUnbiasedScoreTestOrdering()}
#' * \code{MedianUnbiasedStagewiseCombinationFunctionOrdering()}
#'
#' These estimators are constructed by specifying an ordering of the sample space
#' and finding the value of \eqn{\mu}, such that the observed sample is the
#' median of the sample space according to the chosen ordering.
#' Some of the implemented orderings are based on the work presented in
#' \insertCite{emerson1990parameter}{adestr},
#' \insertCite{@Sections 8.4 in @jennison1999group}{adestr},
#' and \insertCite{@Sections 4.1.1 and 8.2.1 in @wassmer2016group}{adestr}.
#' @md
#'
#' @param g1 functional representation of the estimator in the early futility and efficacy regions.
#' @param g2 functional representation of the estimator in the continuation region.
#' @param label name of the estimator. Used in printing methods.
#' @seealso \code{\link{evaluate_estimator}}
#'
#' @return An object of class \code{PointEstimator}.
#' @return an object of class \code{PointEstimator}.
#'
#' @references
#' \insertAllCited{}
#'
#' @export
#'
Expand All @@ -46,18 +93,42 @@ VirtualPointEstimator <- function() stop("Cannot create instance of class Virtua

#' P-values
#'
#' This is the parent class for all p-values implemented in this package.
#' Details about the methods for calculating p-values can be found in
#' (our upcoming paper).
#' @param g1 functional representation of the p-value in the early futility and efficacy regions.
#' @param g2 functional representation of the p-value in the continuation region.
#' @param label name of the p-value. Used in printing methods.
#' @seealso [plot_p]
#'
#' @return An object of class \code{PValue}.
#' @details
#' The implemented p-values are:
#' * \code{MLEOrderingPValue()}
#' * \code{LikelihoodRatioOrderingPValue()}
#' * \code{ScoreTestOrderingPValue()}
#' * \code{StagewiseCombinationFunctionOrderingPValue()}
#'
#' The p-values are calculated by specifying an ordering of the sample space
#' calculating the probability that a random sample under the null hypothesis is
#' larger than the observed sample.
#' Some of the implemented orderings are based on the work presented in
#' \insertCite{emerson1990parameter}{adestr},
#' \insertCite{@Sections 8.4 in @jennison1999group}{adestr},
#' and \insertCite{@Sections 4.1.1 and 8.2.1 in @wassmer2016group}{adestr}.
#' @md
#'
#' @references
#' \insertAllCited{}
#'
#' @export
#'
#' @examples
#' # This is the definition of a 'naive' p-value based on a Z-test for a one-armed trial
#' PValue(
#' g1 = \(smean1, ...) runif(length(smean1)),
#' g2 = \(smean2, ...) runif(length(smean2)),
#' g1 = \(smean1, n1, sigma, ...) pnorm(smean1*sqrt(n1)/sigma, lower.tail=FALSE),
#' g2 = \(smean1, smean2, n1, n2, ...) pnorm((n1 * smean1 + n2 * smean2)/(n1 + n2) *
#' sqrt(n1+n2)/sigma, lower.tail=FALSE),
#' label="My custom p-value")
setClass("PValue", slots = c(g1 = "function", g2 = "function"), contains = "Statistic")
#' @rdname PValue-class
Expand All @@ -69,26 +140,52 @@ VirtualPValue <- function() stop("Cannot create instance of class VirtualPValue.

#' Interval estimators
#'
#' This is the parent class for all confidence intervals implemented in this package.
#' Currently, only confidence intervals for the parameter \eqn{\mu} of a normal distribution
#' are implemented. Details about the methods for calculating confidence intervals can be found in
#' (our upcoming paper).
#' @param l1 functional representation of the lower boundary of the interval in the early futility and efficacy regions.
#' @param u1 functional representation of the upper boundary of the interval in the early futility and efficacy regions.
#' @param l2 functional representation of the lower boundary of the interval in the continuation region.
#' @param u2 functional representation of the upper boundary of the interval in the continuation region.
#' @param two_sided logical indicating whether the confidence interval is two-sided.
#' @param label name of the estimator. Used in printing methods.
#' @seealso \code{\link{evaluate_estimator}}
#'
#' @return An object of class \code{IntervalEstimator}.
#' @return an object of class \code{IntervalEstimator}.
#'
#' @export
#' @aliases ConfidenceInterval ConfidenceInterval-class
#'
#' @details
#' The implemented confidence intervals are:
#' * \code{MLEOrderingCI()}
#' * \code{LikelihoodRatioOrderingCI()}
#' * \code{ScoreTestOrderingCI()}
#' * \code{StagewiseCombinationFunctionOrderingCI()}
#'
#' These confidence intervals are constructed by specifying an ordering of the sample space
#' and finding the value of \eqn{\mu}, such that the observed sample is the
#' \eqn{\alpha/2} (or (\eqn{1-\alpha/2})) quantile of the sample space according to the
#' chosen ordering.
#' Some of the implemented orderings are based on the work presented in
#' \insertCite{emerson1990parameter}{adestr},
#' \insertCite{@Sections 8.4 in @jennison1999group}{adestr},
#' and \insertCite{@Sections 4.1.1 and 8.2.1 in @wassmer2016group}{adestr}.
#' @md
#'
#' @references
#' \insertAllCited{}
#'
#' @examples
#' # This is the definition of the 'naive' confidence interval for one-armed trials
#' IntervalEstimator(
#' two_sided = FALSE,
#' l1 = \(smean1, ...) smean1 - 1,
#' u1 = \(smean1, ...) smean1 + 1,
#' l2 = \(smean2, ...) smean2 - 1,
#' u2 = \(smean2, ...) smean2 + 1,
#' label="My custom p-value")
#' two_sided = TRUE,
#' l1 = \(smean1, n1, sigma, ...) smean1 - qnorm(.95, sd = sigma/sqrt(n1)),
#' u1 = \(smean1, n1, sigma, ...) smean1 + qnorm(.95, sd = sigma/sqrt(n1)),
#' l2 = \(smean1, smean2, n1, n2, sigma, ...) smean2 - qnorm(.95, sd = sigma/sqrt(n1 + n2)),
#' u2 = \(smean1, smean2, n1, n2, sigma, ...) smean2 + qnorm(.95, sd = sigma/sqrt(n1 + n2)),
#' label="My custom CI")
setClass("IntervalEstimator", slots = c(two_sided = "logical", l1 = "function", u1 = "function", l2 = "function", u2 = "function"), contains = "Estimator")
#' @rdname IntervalEstimator-class
#' @export
Expand Down Expand Up @@ -415,7 +512,7 @@ FirstStageSampleMean <- function() new("PointEstimator", g1 = \(smean1, ...) sme
setClass("WeightedSampleMean", contains = "PointEstimator")

#' @rdname PointEstimator-class
#' @param w1 weight of the first-stage sample mean.
#' @param w1 weight of the first-stage data.
#' @importFrom scales percent
#' @export
WeightedSampleMean <- function(w1=.5) new("PointEstimator", g1 = \(smean1, ...) smean1, g2 = \(smean1, smean2, n1, n2, ...) w1 * smean1 + (1-w1) * smean2,
Expand All @@ -424,6 +521,7 @@ WeightedSampleMean <- function(w1=.5) new("PointEstimator", g1 = \(smean1, ...)
setClass("AdaptivelyWeightedSampleMean", contains = "VirtualPointEstimator", slots = c(w1 = "numeric"))
#' @rdname PointEstimator-class
#' @importFrom scales percent
#' @export
AdaptivelyWeightedSampleMean <- function(w1 = 1/sqrt(2)) {
new(
"AdaptivelyWeightedSampleMean",
Expand Down

0 comments on commit 2d1ab0f

Please sign in to comment.