Update for CRAN upload

jan-imbi · Sep 22, 2023 · 807c52a · 807c52a
1 parent 2d1ab0f
commit 807c52a
Show file tree

Hide file tree

Showing 31 changed files with 447 additions and 293 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,7 +1,7 @@
 Package: adestr
 Type: Package
 Title: Adaptive Design Estimation in R
-Version: 0.0.1
+Version: 0.5.0
 Authors@R: c(person("Jan", "Meis", role = c("aut", "cre"), email="meis@imbi.uni-heidelberg.de", comment = c(ORCID = "0000-0001-5407-7220")))
 Description:
   This package implements methods to evaluate the performance characteristics of
@@ -17,8 +17,7 @@ LazyData: true
 VignetteBuilder: knitr
 RoxygenNote: 7.2.3
 Depends:
-    R (>= 4.0.0),
-    adoptr
+    R (>= 4.0.0)
 Imports: 
     methods,
     stats,
@@ -41,6 +40,7 @@ Suggests:
 Config/testthat/edition: 3
 Collate: 
     'adestr_package.R'
+    'twostagedesign_with_cache.R'
     'analyze.R'
     'estimators.R'
     'densities.R'
@@ -54,12 +54,10 @@ Collate:
     'mle_distribution.R'
     'mlmse_score.R'
     'n2c2_helpers.R'
-    'twostagedesign_with_cache.R'
     'plot.R'
     'priors.R'
     'print.R'
 URL: https://jan-imbi.github.io/adestr/
 RdMacros: Rdpack
-Remotes: github::jan-imbi/adoptr@dissertation
 
 
diff --git a/NAMESPACE b/NAMESPACE
@@ -30,6 +30,7 @@ export(MinimizePeakVariance)
 export(NaiveCI)
 export(NeymanPearsonOrderingCI)
 export(NeymanPearsonOrderingPValue)
+export(Normal)
 export(NormalPrior)
 export(OverestimationProbability)
 export(PValue)
@@ -43,6 +44,7 @@ export(ScoreTestOrderingPValue)
 export(SoftCoverage)
 export(StagewiseCombinationFunctionOrderingCI)
 export(StagewiseCombinationFunctionOrderingPValue)
+export(Student)
 export(TestAgreement)
 export(UniformPrior)
 export(Variance)
@@ -61,8 +63,8 @@ exportClasses(IntervalEstimator)
 exportClasses(PValue)
 exportClasses(PointEstimator)
 exportClasses(Statistic)
+exportClasses(TwoStageDesign)
 exportMethods(plot)
-import(adoptr)
 import(ggplot2)
 import(methods)
 importFrom(Rdpack,reprompt)

diff --git a/NEWS.md b/NEWS.md
@@ -1,3 +1,7 @@
+# adestr 0.5.0
+
+* First CRAN submission.
+
 # adestr 0.0.1
 
 * Added a `NEWS.md` file to track changes to the package.
diff --git a/R/adestr_package.R b/R/adestr_package.R
@@ -3,15 +3,15 @@
 #'
 #' @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
+#' Specifically, this package is written to interface with trial designs created by the \code{adoptr} package
 #' \insertCite{kunzmann2021adoptr,pilz2021optimal}{adestr}.
 #' Apart from the a priori evaluation of performance characteristics, this package also allows for the
 #' \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
+#' @import methods
 #' @importFrom stats dnorm pnorm qnorm dt pt qt dchisq pchisq qchisq integrate uniroot var
 #' @importFrom cubature hcubature
 #' @importFrom Rdpack reprompt

diff --git a/R/analyze.R b/R/analyze.R
@@ -101,7 +101,7 @@ setMethod("analyze", signature("data.frame"),
                 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)
+              calc_n2 <- .n2_extrapol(design, test_val)
               if (sdata$n_stages==2L){
                 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.")

diff --git a/R/estimators.R b/R/estimators.R
@@ -566,12 +566,13 @@ setMethod("get_stagewise_estimators", signature("MinimizePeakVariance", "Normal"
                    design,
                    sigma,
                    exact) {
-            get_ess <- function(mu) {
-              H <- PointMassPrior(mu, 1)
-              ess <- ExpectedSampleSize(data_distribution, H)
-              evaluate(ess, design, optimization=TRUE)
-            }
-            max_ess_mu <- optimize(get_ess, z_to_smean(c(design@c1f, design@c1e), n1(design, round=FALSE), sigma, data_distribution@two_armed), maximum = TRUE)$maximum
+            # get_ess <- function(mu) {
+            #   H <- PointMassPrior(mu, 1)
+            #   ess <- ExpectedSampleSize(data_distribution, H)
+            #   evaluate(ess, design, optimization=TRUE)
+            # }
+            # max_ess_mu <- optimize(get_ess, z_to_smean(c(design@c1f, design@c1e), n1(design, round=FALSE), sigma, data_distribution@two_armed), maximum = TRUE)$maximum
+            max_ess_mu <- mean(z_to_smean(c(design@c1f, design@c1e), n1(design, round=FALSE), sigma, data_distribution@two_armed))
             get_var <- function(w1) {
               est <- AdaptivelyWeightedSampleMean(w1)
               evaluate_estimator(score = Variance(),

diff --git a/R/evaluate_estimator.R b/R/evaluate_estimator.R
@@ -81,8 +81,8 @@ setMethod("c", signature("EstimatorScoreResultList"), definition =
 #' First, a functional representation of the integrand is created by combining information
 #' from the \code{\link{EstimatorScore}} object (\code{score}) and the \code{\link{PointEstimator}} or
 #' \code{\link{IntervalEstimator}} object (\code{estimator}).
-#' The sampling distribution of a design is determined by the \code{\link{TwoStageDesign}} object
-#' (\code{design}) and the \code{\link{DataDistribution}} object (\code{data_distribution}),
+#' The sampling distribution of a design is determined by the \code{TwoStageDesign} object
+#' (\code{design}) and the \code{DataDistribution} object (\code{data_distribution}),
 #' as well as the assumed parameters \eqn{\mu} (mu) and \eqn{\sigma} (sigma).
 #' The other parameters control various details of the integration problem.
 #'
@@ -98,7 +98,7 @@ setMethod("c", signature("EstimatorScoreResultList"), definition =
 #' estimating the standard deviation.
 #'
 #' If the parameter \code{exact} is set to \code{FALSE}
-#' (the default), the continuous version of the second-stage sample-size function \code{\link{n2}}
+#' (the default), the continuous version of the second-stage sample-size function \code{n2}
 #' is used. Otherwise, an integer valued version of that function will be used,
 #' though this is considerably slower.
 #'

diff --git a/R/helper_functions.R b/R/helper_functions.R
@@ -102,46 +102,68 @@ 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 two_armed (logical) determins whether the design is for one- or
+#' two-armed trials.
 #' @param label (optional) label to be assigned to the design.
 #'
 #' @return an exmplary design of class \code{TwoStageDesign}.
 #' @export
 #'
 #' @examples
-#' design <- get_example_design()
-#' # Type I error
-#' evaluate(Power(Normal(FALSE), PointMassPrior(0, 1)), design)
-#' # Power
-#' evaluate(Power(Normal(FALSE), PointMassPrior(.4, 1)), design)
-#' # Expected sample size under the alternative
-#' evaluate(ExpectedSampleSize(Normal(FALSE), PointMassPrior(.4, 1)), design)
-#' # Expected sample size under the null
-#' evaluate(ExpectedSampleSize(Normal(FALSE), PointMassPrior(0, 1)), design)
-get_example_design <- function(label = NULL) {
-  d <- TwoStageDesign(
-    n1 = 28.16834031633078083701,
-    c1f = 0.7907304707554818623549,
-    c1e = 2.291260947864900643367,
-    n2_pivots = c(
-      39.39294353955478555918,
-      37.23397813905835818105,
-      33.27173714438612961430,
-      27.77227568901122012335,
-      21.41776450755991234587,
-      15.17163280081247300757,
-      10.25508398663193787570
-    ),
-    c2_pivots = c(
-      2.16914648055318837194250,
-      2.02493357331804890719695,
-      1.77299079624771049878973,
-      1.42524439642541422834654,
-      0.99916431580845133098023,
-      0.52325801518650127963639,
-      0.07133753446126563091401
-    ),
-    7
-  )
+#' get_example_design()
+#'
+get_example_design <- function(two_armed = FALSE, label = NULL) {
+  if (two_armed) {
+    d <- TwoStageDesign(
+      n1 = 56.33739084822602904978,
+      c1f = 0.7907356135206976555097,
+      c1e = 2.291313615804561720779,
+      n2_pivots = c(
+        78.74984462914770233510,
+        74.46976145811771630179,
+        66.54436357139142899086,
+        55.54462857388867291775,
+        42.83473938241603207189,
+        30.34059894227031151104,
+        20.48959489543554823854
+      ),
+      c2_pivots = c(
+        2.16905734410577055726321,
+        2.02492349183474207308109,
+        1.77298374394898416994693,
+        1.42521477360223225439029,
+        0.99909735679793421070372,
+        0.52314051699418129270924,
+        0.07058637352917693230658
+      ),
+      7
+    )
+  } else {
+    d <- TwoStageDesign(
+      n1 = 28.16834031633078083701,
+      c1f = 0.7907304707554818623549,
+      c1e = 2.291260947864900643367,
+      n2_pivots = c(
+        39.39294353955478555918,
+        37.23397813905835818105,
+        33.27173714438612961430,
+        27.77227568901122012335,
+        21.41776450755991234587,
+        15.17163280081247300757,
+        10.25508398663193787570
+      ),
+      c2_pivots = c(
+        2.16914648055318837194250,
+        2.02493357331804890719695,
+        1.77299079624771049878973,
+        1.42524439642541422834654,
+        0.99916431580845133098023,
+        0.52325801518650127963639,
+        0.07133753446126563091401
+      ),
+      7
+    )
+  }
   d <- TwoStageDesignWithCache(d)
   attr(d, "label") <- label
   d
@@ -219,13 +241,3 @@ get_statistics_from_paper <- function(point_estimators = TRUE,
     ret <- c(ret, p)
   return(ret)
 }
-
-
-
-
-
-
-
-
-
-
diff --git a/R/mlmse_score.R b/R/mlmse_score.R
@@ -1,38 +1,40 @@
-setClass("MLMSE",
-         slots = c(
-           mu = "numeric",
-           sigma = "numeric",
-           two_armed = "logical",
-           tol = "numeric",
-           maxEval = "numeric",
-           absError = "numeric"
-         ),
-         contains = "UnconditionalScore")
-MLMSE <- function(mu = 0, sigma = 1, two_armed = FALSE, tol = 1e-5, maxEval = 1e7, absError = 1e-7) {
-  new("MLMSE",
-      mu = mu,
-      sigma = sigma,
-      two_armed = two_armed,
-      tol = tol,
-      maxEval = maxEval,
-      absError = absError,
-      label = "E[|mu_ML - mu|^2]")
-}
-setMethod("evaluate",
-          signature("MLMSE", "TwoStageDesign"),
-          function(s, design, ...) {
-            mean(
-              evaluate_estimator(
-                score = MSE(),
-                estimator = SampleMean(),
-                data_distribution = Normal(two_armed = s@two_armed),
-                design = design,
-                mu = s@mu,
-                sigma = s@sigma,
-                tol = s@tol,
-                maxEval = s@maxEval
-              )@results$MSE
-            )
-          })
-
+### Uncomment this once adoptr is back on CRAN ###
 
+# setClass("MLMSE",
+#          slots = c(
+#            mu = "numeric",
+#            sigma = "numeric",
+#            two_armed = "logical",
+#            tol = "numeric",
+#            maxEval = "numeric",
+#            absError = "numeric"
+#          ),
+#          contains = "UnconditionalScore")
+# MLMSE <- function(mu = 0, sigma = 1, two_armed = FALSE, tol = 1e-5, maxEval = 1e7, absError = 1e-7) {
+#   new("MLMSE",
+#       mu = mu,
+#       sigma = sigma,
+#       two_armed = two_armed,
+#       tol = tol,
+#       maxEval = maxEval,
+#       absError = absError,
+#       label = "E[|mu_ML - mu|^2]")
+# }
+# setMethod("evaluate",
+#           signature("MLMSE", "TwoStageDesign"),
+#           function(s, design, ...) {
+#             mean(
+#               evaluate_estimator(
+#                 score = MSE(),
+#                 estimator = SampleMean(),
+#                 data_distribution = Normal(two_armed = s@two_armed),
+#                 design = design,
+#                 mu = s@mu,
+#                 sigma = s@sigma,
+#                 tol = s@tol,
+#                 maxEval = s@maxEval
+#               )@results$MSE
+#             )
+#           })
+#
+#
diff --git a/R/n2c2_helpers.R b/R/n2c2_helpers.R
@@ -74,7 +74,6 @@ get_c2_coefficients <- function(design){
 #' @param design an object of class \code{\link{TwoStageDesignWithCache}}.
 #' @param x1 first-stage test statistic
 #'
-#' @seealso \link[adoptr]{n2}
 n2_extrapol <- function(design, x1) {
   if (length(design@n2_pivots) > 1L){
     fastmonoH.FC_evaluate(x1, design@n2_coefficients)
@@ -89,7 +88,6 @@ n2_extrapol <- function(design, x1) {
 #' @param design an object of class \code{\link{TwoStageDesignWithCache}}.
 #' @param x1 first-stage test statistic
 #'
-#' @seealso \link[adoptr]{c2}
 c2_extrapol <- function(design, x1) {
   fastmonoH.FC_evaluate(x1, design@c2_coefficients)
 }

diff --git a/R/plot.R b/R/plot.R
@@ -76,7 +76,7 @@ setMethod("plot", signature = "list", definition =
 #' of values for the first and second-stage test statistics.
 #'
 #' When the first-stage test statistic lies below the futility threshold (c1f) or
-#' above the early efficacy threshold (c1e) of the \code{\link{TwoStageDesign}},
+#' above the early efficacy threshold (c1e) of the \code{TwoStageDesign},
 #' there is no second-stage test statistics. The p-values in these regions are only
 #' based on the first-stage values.
 #' For first-stage test statistic values between c1f and c1e, the first and second-stage
@@ -379,7 +379,7 @@ setMethod("plot_sample_mean", signature("DataDistribution", "TwoStageDesign"),
               plt <- ggplot(data = smean_dat, aes(x = .data$`smean`, y = .data$`Density`)) +
                 geom_line(size = 1) +
                 scale_x_continuous(name = "Sample mean") +
-                facet_wrap(vars(n))
+                facet_wrap(vars(.data$n))
               return(plt)
             }
           })