More documentation work

NAMESPACE

@@ -44,7 +44,6 @@ export(StagewiseCombinationFunctionOrderingPValue)
 export(TestAgreement)
 export(UniformPrior)
 export(Variance)
-export(VirtualIntervalEstimator)
 export(WeightedSampleMean)
 export(Width)
 export(analyze)
@@ -58,6 +57,7 @@ exportClasses(EstimatorScore)
 exportClasses(IntervalEstimator)
 exportClasses(PValue)
 exportClasses(PointEstimator)
+exportClasses(Statistic)
 import(adoptr)
 import(ggplot2)
 import(ggpubr)

R/analyze.R

@@ -13,7 +13,7 @@ Results <- setClass("Results", slots = c(data ="data.frame",
 #'
 #' @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
-#' \code{\link{PValues}}.
+#' \code{\link{PValue}}.
 #' @inheritParams evaluate_estimator
 #'
 #' @return \code{Results} object containing the values of the statistics

R/estimators.R

@@ -1,3 +1,26 @@
+#' Statistics and Estimators of the adestr package
+#'
+#' The \code{\link{Statistic}} class is a parent class for the classes
+#' \code{\link{Estimator}} and \code{\link{PValue}}. The \code{\link{Estimator}} class is a parent
+#' for the classes \code{\link{PointEstimator}} and \code{\link{ConfidenceInterval}}.
+#'
+#' The function \code{\link{analyze}} can be used to calculate the value
+#' of a \code{\link{Statistic}} for a given dataset.
+#'
+#' The function \code{\link{evaluate_estimator}} can be used to evaluate
+#' \link[EstimatorScore]{distributional quantities} of an \code{\link{Estimator}}
+#' like the \code{\link{MSE}} for a \code{\link{PointEstimator}} or the
+#' \code{\link{Coverage}} for a \code{\link{ConfidenceInterval}}.
+#'
+#'
+#' @param label name of the statistic. Used in printing methods.
+#'
+#' @export
+#' @rdname Statistic-class
+#' @aliases Statistic Statistics Estimator
+#' @seealso \code{\link{PointEstimator}} \code{\link{ConfidenceInterval}} \code{\link{PValue}}
+#' @seealso \code{\link{analyze}} \code{\link{evaluate_estimator}}
+#' @seealso \code{\link{EstimatorScore}}
 setClass("Statistic", slots = c(label = "character"))
 setClass("Estimator", contains = "Statistic")
 
@@ -18,7 +41,6 @@ setClass("PointEstimator", slots = c(g1 = "function", g2 = "function"), contains
 #' @export
 PointEstimator <- function(g1, g2, label) new("PointEstimator", g1 = g1, g2 = g2, label = label)
 setClass("VirtualPointEstimator", contains = "PointEstimator")
-#' @rdname PointEstimator-class
 VirtualPointEstimator <- function() stop("Cannot create instance of class VirtualPointEstimator.")
 
 
@@ -42,7 +64,6 @@ setClass("PValue", slots = c(g1 = "function", g2 = "function"),  contains = "Sta
 #' @export
 PValue <- function(g1, g2, label) new("PValue", g1 = g1, g2 = g2, label = label)
 setClass("VirtualPValue", contains = "PValue")
-#' @rdname PValue-class
 VirtualPValue <- function() stop("Cannot create instance of class VirtualPValue.")
 
 
@@ -58,6 +79,7 @@ VirtualPValue <- function() stop("Cannot create instance of class VirtualPValue.
 #' @return An object of class \code{IntervalEstimator}.
 #'
 #' @export
+#' @aliases ConfidenceInterval ConfidenceInterval-class
 #'
 #' @examples
 #' IntervalEstimator(
@@ -72,8 +94,6 @@ setClass("IntervalEstimator", slots = c(two_sided = "logical", l1 = "function",
 #' @export
 IntervalEstimator <- function(two_sided, l1, u1, l2, u2, label) new("IntervalEstimator", two_sided = two_sided, l1 = l1, u1 = u1, l2 = l2, u2 = u2, label = label)
 setClass("VirtualIntervalEstimator", contains = "IntervalEstimator")
-#' @rdname IntervalEstimator-class
-#' @export
 VirtualIntervalEstimator <- function() stop("Cannot create instance of class VirtualIntervalEstimator")
 
 

R/evaluate_estimator.R

@@ -1,14 +1,14 @@
-
 #' Performance scores for point and interval estimators
 #'
-#' These functions encode various metrics which can be used to evaluate
+#' These classes encode various metrics which can be used to evaluate
 #' the performance characteristics of point and interval estimators.
 #'
 #'
 #' @slot label name of the performance score. Used in printing methods.
 #'
 #' @return an \code{EstimatorScore} object.
 #' @export
+#' @aliases EstimatorScore
 #' @seealso \code{\link{evaluate_estimator}}
 #'
 #' @inherit evaluate_estimator examples

R/helper_functions.R

@@ -102,6 +102,8 @@ get_overall_svar_twoarm <- function(smean1, smean1T, svar1, smean2, smean2T, sva
 #' for a normally distributed test statistic (i.e. known variance).
 #' For an alternative hypothesis of mu=0.4, the overall power is 80\%.
 #'
+#' @param label (optional) label to be assigned to the design.
+#'
 #' @return an exmplary design of class \code{TwoStageDesign}.
 #' @export
 #'

R/reference_implementation.R

@@ -19,26 +19,26 @@
   x2 <- .z_to_x(z = z2, n = n2, mu0 = mu0, sigma = sigma)
   (n1 * x1 + n2 * x2) / (n1 + n2)
 }
-# n2_extrapol <- function(design, x1) {
-#   if (length(design@n2_pivots)>1){
-#     h <- (design@c1e - design@c1f) / 2
-#     return(stats::splinefun(
-#       h * design@x1_norm_pivots + (h + design@c1f),
-#       design@n2_pivots,
-#       method = "monoH.FC"
-#     )(x1))
-#   } else{
-#     return(design@n2_pivots)
-#   }
-# }
-# c2_extrapol <- function(design, x1) {
-#   h <- (design@c1e - design@c1f) / 2
-#   return(stats::splinefun(
-#     h * design@x1_norm_pivots + (h + design@c1f),
-#     design@c2_pivots,
-#     method = "monoH.FC"
-#   )(x1))
-# }
+.n2_extrapol <- function(design, x1) {
+  if (length(design@n2_pivots)>1){
+    h <- (design@c1e - design@c1f) / 2
+    return(stats::splinefun(
+      h * design@x1_norm_pivots + (h + design@c1f),
+      design@n2_pivots,
+      method = "monoH.FC"
+    )(x1))
+  } else{
+    return(design@n2_pivots)
+  }
+}
+.c2_extrapol <- function(design, x1) {
+  h <- (design@c1e - design@c1f) / 2
+  return(stats::splinefun(
+    h * design@x1_norm_pivots + (h + design@c1f),
+    design@c2_pivots,
+    method = "monoH.FC"
+  )(x1))
+}
 
 ## The densities for integration.
 .f1 <- function(z1, n1, mu, mu0, sigma) {
@@ -61,7 +61,7 @@
     if (z1 < c1f || z1 > c1e) {
       .f1(z1 = z1, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       .f2(z1 = z1, z2 = z2, n1 = n1, n2 = n2, mu = mu, mu0 = mu0, sigma = sigma)
     }
   ret
@@ -115,7 +115,7 @@
   # Cannot use 2d integral because boundaries of integral are function of z1
   continuation_int <- hcubature(
     \(z1) {
-       n2 <- n2_extrapol(design, z1)
+       n2 <- .n2_extrapol(design, z1)
        x1 <- .z_to_x(z = z1, n = n1, mu0 = mu0, sigma = sigma)
        .f1(z1 = z1, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma) *
        hcubature(
@@ -136,7 +136,7 @@
   n1 <- design@n1
   int <- hcubature(
     \(x, n1, mu, mu0, sigma, design) {
-      n2 <- n2_extrapol(design, x[1,,drop=FALSE])
+      n2 <- .n2_extrapol(design, x[1,,drop=FALSE])
       z2 <- ( ((n1 + n2) * y - n1 * .z_to_x(z = x[1,,drop=FALSE], n = n1, mu0 = mu0, sigma = sigma)) /n2 - mu0) * sqrt(n2)/sigma
       (n1 + n2)/n2 * sqrt(n2)/sigma * .f2(z1 = x[1,,drop=FALSE], z2 = z2, n1 = n1, n2 = n2, mu = mu, mu0 = mu0, sigma = sigma)
     },
@@ -177,7 +177,7 @@
   )$integral
   continuation <- hcubature(
     f = \(x, n1, mu, mu0, sigma, design){
-      n2 <- n2_extrapol(design, x[1,,drop=FALSE])
+      n2 <- .n2_extrapol(design, x[1,,drop=FALSE])
       .z1_z2_to_x(z1 = x[1, , drop = FALSE], z2 = x[2, , drop = FALSE], n1 = n1, n2 = n2, mu0 = mu0, sigma = sigma ) *
         .f2(z1 = x[1, , drop = FALSE], z2 = x[2, , drop = FALSE], n1 = n1, n2 = n2, mu = mu, mu0 = mu0, sigma = sigma)
     },
@@ -221,7 +221,7 @@
   )$integral
   continuation <- hcubature(
     f = \(x, n1, mu, mu0, sigma, design){
-      n2 <- n2_extrapol(design, x[1,,drop=FALSE])
+      n2 <- .n2_extrapol(design, x[1,,drop=FALSE])
       x1 <- .z_to_x(z = x[1,,drop=FALSE], n = n1, mu0 = mu0, sigma = sigma)
       x2 <- .z_to_x(z = x[2,,drop=FALSE], n = n2, mu0 = mu0, sigma = sigma)
       .adaptively_weighted_mle(x1 = x1, x2 = x2, n1 = n1, n2 = n2, sigma = sigma, w = w) *
@@ -259,7 +259,7 @@
   )$integral
   continuation <- hcubature(
     f = \(x, n1, mu, mu0, sigma, design){
-      n2 <- n2_extrapol(design, x[1,,drop=FALSE])
+      n2 <- .n2_extrapol(design, x[1,,drop=FALSE])
       x1 <- .z_to_x(x[1,,drop=FALSE], n = n1, mu0 = mu0, sigma = sigma)
       x2 <- .z_to_x(x[2,,drop=FALSE], n = n2, mu0 = mu0, sigma = sigma)
       (.adaptively_weighted_mle(x1 = x1, x2 = x2, n1 = n1, n2 = n2, sigma = sigma, w = w) - E)^2 *
@@ -283,8 +283,8 @@
     sigma = sigma,
     design = design
   )$minimum
-  max_n2 <- n2_extrapol(design, design@c1f)
-  min_n2 <- n2_extrapol(design, design@c1e)
+  max_n2 <- .n2_extrapol(design, design@c1f)
+  min_n2 <- .n2_extrapol(design, design@c1e)
   max_tau <- (minw * sqrt(design@n1) / sigma) / (minw * sqrt(design@n1) / sigma + sqrt(1 - minw^2) * sqrt(min_n2) / sigma)
   min_tau <- (minw * sqrt(design@n1) / sigma) / (minw * sqrt(design@n1) / sigma + sqrt(1 - minw^2) * sqrt(max_n2) / sigma)
   message(sprintf("Minimum weight w = %.3f, which means tau ranges from %.3f to %.3f.", minw, min_tau, max_tau))
@@ -307,7 +307,7 @@
     if (z1 < c1f || z1 > c1e) {
       x1
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       x <- .x1_x2_to_x(x1 = x1, x2 = x2, n1 = n1, n2 = n2)
       numerator <-
         hcubature(
@@ -353,9 +353,9 @@
       # The difference between the pseudo Rao-Blackwell and the Rao-Blackwell
       # estimator is that the (normal) Rao-Blackwell estimator uses the actual
       # n2 preimage.
-      n2 <- ceiling(n2_extrapol(design, z1))
-      lower_z1 <- uniroot(\(x) n2_extrapol(design, x) - n2 - -2, c(c1f, c1e))$root
-      upper_z1 <- uniroot(\(x) n2_extrapol(design, x) - n2 - -1, c(c1f, c1e))$root
+      n2 <- ceiling(.n2_extrapol(design, z1))
+      lower_z1 <- uniroot(\(x) .n2_extrapol(design, x) - n2 - -2, c(c1f, c1e))$root
+      upper_z1 <- uniroot(\(x) .n2_extrapol(design, x) - n2 - -1, c(c1f, c1e))$root
       x <- .x1_x2_to_x(x1 = x1, x2 = x2, n1 = n1, n2 = n2)
       numerator <-
         hcubature(
@@ -412,7 +412,7 @@
   )$integral
   continuation <- hcubature(
     f = \(x, n1, mu, mu0, sigma, design){
-      n2 <- n2_extrapol(design, x[1,,drop=FALSE])
+      n2 <- .n2_extrapol(design, x[1,,drop=FALSE])
       .z1_z2_to_x(z1 = x[1, , drop = FALSE], z2 = x[2, , drop = FALSE], n1 = n1, n2 = n2, mu0 = mu0, sigma = sigma) *
         ((x[1,,drop=FALSE] + (mu0 - mu) * sqrt(n1)/sigma) * sqrt(n1)/sigma +
            (x[2,,drop=FALSE] + (mu0 - mu) * sqrt(n2)/sigma) * sqrt(n2)/sigma) *
@@ -495,10 +495,10 @@
   2 + pnorm(.x_to_z(x = x1, n = n1, mu0 = mu0, sigma = sigma))
 }
 .swcf <- function(x1, x2, n1, n2, mu, sigma, design){
-  .x_to_z(x = x2, n = n2, mu0 = mu, sigma = sigma) - c2_extrapol(design, .x_to_z(x = x1, n = n1, mu0 = mu, sigma = sigma))
+  .x_to_z(x = x2, n = n2, mu0 = mu, sigma = sigma) - .c2_extrapol(design, .x_to_z(x = x1, n = n1, mu0 = mu, sigma = sigma))
 }
 .rho_swcf_B <- function(cf_value, x1, n1, n2, mu, sigma, design){
-  sigma / sqrt(n2) * (cf_value + c2_extrapol(design, (x1 - mu)*sqrt(n1)/sigma)) + mu
+  sigma / sqrt(n2) * (cf_value + .c2_extrapol(design, (x1 - mu)*sqrt(n1)/sigma)) + mu
 }
 .calculate_A <- function(y, n1, mu, sigma, design){
   c1f <- design@c1f
@@ -533,14 +533,14 @@
     if (z1 < c1f || z1 > c1e) {
       .rho_ml1(x1 = x1)
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       .rho_ml2(x1 = x1, x2 = x2, n1 = n1, n2 = n2)
     }
   y <- .rho_ml_y(rho_ml = rho)
   A <- .calculate_A(y = y, n1 = n1, mu = mu, sigma = sigma, design = design)
   int <- hcubature(
     f = \(z_) {
-      n2_ <- n2_extrapol(design, z_)
+      n2_ <- .n2_extrapol(design, z_)
       x1_ <- .z_to_x(z_, n = n1, mu0 = mu0, sigma = sigma)
       (1 - pnorm((.rho_ml_B(rho_ml = rho, x1 = x1_, n1 = n1, n2 = n2_, mu = mu) - mu) * sqrt(n2_) / sigma)) *
         .f1(z1 = z_, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
@@ -585,14 +585,14 @@
     if (z1 < c1f || z1 > c1e) {
       .rho_lr1(x1 = x1, n1 = n1, mu = mu)
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       .rho_lr2(x1 = x1, x2 = x2, n1 = n1, n2 = n2, mu = mu)
     }
   y <- .rho_lr_y(rho_lr = rho, n1 = n1, mu = mu)
   A <- .calculate_A(y = y, n1 = n1, mu = mu, sigma = sigma, design = design)
   int <- hcubature(
     f = \(z_) {
-      n2_ <- n2_extrapol(design, z_)
+      n2_ <- .n2_extrapol(design, z_)
       x1_ <- .z_to_x(z_, n = n1, mu0 = mu0, sigma = sigma)
       (1 - pnorm((.rho_lr_B(rho_lr = rho, x1 = x1_, n1 = n1, n2 = n2_, mu = mu) - mu) * sqrt(n2_) / sigma)) *
         .f1(z1 = z_, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
@@ -637,14 +637,14 @@
     if (z1 < c1f || z1 > c1e) {
       .rho_st1(x1 = x1, n1 = n1, mu = mu)
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       .rho_st2(x1 = x1, x2 = x2, n1 = n1, n2 = n2, mu = mu)
     }
   y <- .rho_st_y(rho_st = rho, n1 = n1, mu = mu)
   A <- .calculate_A(y = y, n1 = n1, mu = mu, sigma = sigma, design = design)
   int <- hcubature(
     f = \(z_) {
-      n2_ <- n2_extrapol(design, z_)
+      n2_ <- .n2_extrapol(design, z_)
       x1_ <- .z_to_x(z_, n = n1, mu0 = mu0, sigma = sigma)
       (1 - pnorm((.rho_st_B(rho_st = rho, x1 = x1_, n1 = n1, n2 = n2_, mu = mu) - mu) * sqrt(n2_) / sigma)) *
         .f1(z1 = z_, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
@@ -691,15 +691,15 @@
     if (z1 < c1f || z1 > c1e) {
       .rho_np1(x1 = x1, n1 = n1, mu0 = mu0, mu1 = mu1)
     } else {
-      n2 <- n2_extrapol(design, z1)
+      n2 <- .n2_extrapol(design, z1)
       .rho_np2(x1 = x1, x2 = x2, n1 = n1, n2 = n2, mu0 = mu0, mu1 = mu1)
     }
   y <- .rho_np_y(rho_np = rho, n1 = n1, mu0 = mu0, mu1 = mu1)
 
   A <- .calculate_A(y = y, n1 = n1, mu = mu, sigma = sigma, design = design)
   int <- hcubature(
     f = \(z_) {
-      n2_ <- n2_extrapol(design, z_)
+      n2_ <- .n2_extrapol(design, z_)
       x1_ <- .z_to_x(z_, n = n1, mu0 = mu0, sigma = sigma)
       (1 - pnorm((.rho_np_B(rho_np = rho, x1 = x1_, n1 = n1, n2 = n2_, mu0 = mu0, mu1 = mu1) - mu) * sqrt(n2_) / sigma)) *
         .f1(z1 = z_, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
@@ -746,12 +746,12 @@
     if (z1 < c1f || z1 > c1e) {
       1 - pnorm( (x1 - mu)*sqrt(n1)/sigma )
     } else {
-      n2 <- n2_extrapol(design, z1)
-      cf_value <- (x2 - mu) * sqrt(n2) / sigma - c2_extrapol(design, (x1 - mu) * sqrt(n1) / sigma)
+      n2 <- .n2_extrapol(design, z1)
+      cf_value <- (x2 - mu) * sqrt(n2) / sigma - .c2_extrapol(design, (x1 - mu) * sqrt(n1) / sigma)
       A <- 1 - pnorm(c1e -mu*sqrt(n1)/sigma)
       int <- hcubature(
         f = \(z_) {
-          n2_ <- n2_extrapol(design, z_)
+          n2_ <- .n2_extrapol(design, z_)
           x1_ <- .z_to_x(z_, n = n1, mu0 = mu0, sigma = sigma)
           (1 - pnorm((.rho_swcf_B(cf_value = cf_value, x1 = x1_, n1 = n1, n2 = n2_, mu = mu, sigma = sigma, design = design) - mu) * sqrt(n2_) / sigma)) *
             .f1(z1 = z_, n1 = n1, mu = mu, mu0 = mu0, sigma = sigma)
@@ -802,20 +802,20 @@
   c1f <- design@c1f
   c1e <- design@c1e
   z1 <- .x_to_z(x = x1, n = n1, mu0 = mu0, sigma = sigma)
-  n2 <- n2_extrapol(design, z1)
+  n2 <- .n2_extrapol(design, z1)
   lower_l <- x1 - c1e * sigma / sqrt(n1)
   upper_l <- x1 - c1f * sigma / sqrt(n1)
   lower_z2 <- (x2 - lower_l) * sqrt(n2) / sigma
   upper_z2 <- (x2 - upper_l) * sqrt(n2) / sigma
-  minc2 <- c2_extrapol(design, c1e)
-  maxc2 <- c2_extrapol(design, c1f)
+  minc2 <- .c2_extrapol(design, c1e)
+  maxc2 <- .c2_extrapol(design, c1f)
   # Assumes c2 is monotonically decreasing
   if (maxc2 < lower_z2) {
     ret <- upper_l
   } else if (minc2 >upper_z2) {
     ret <- lower_l
   } else{
-    ret <- uniroot(\(x) c2_extrapol(design, (x1 - x)*sqrt(n1)/sigma) - (x2 - x)*sqrt(n2)/sigma,
+    ret <- uniroot(\(x) .c2_extrapol(design, (x1 - x)*sqrt(n1)/sigma) - (x2 - x)*sqrt(n2)/sigma,
                    interval = c(lower_l, upper_l))$root
   }
   ret