diff --git a/R/analyze.R b/R/analyze.R index 3fb78ee..37c7e6c 100644 --- a/R/analyze.R +++ b/R/analyze.R @@ -37,7 +37,7 @@ Results <- setClass("Results", slots = c(data ="data.frame", #' design = get_example_design() #' ) setGeneric("analyze", function(data, - statistics, + statistics = list(), data_distribution, use_full_twoarm_sampling_distribution = FALSE, design, @@ -63,6 +63,31 @@ setMethod("analyze", signature("data.frame"), if (!is.list(statistics)){ statistics <- list(statistics) } + test_val <- + if (is(data_distribution, "Normal")) + z_test(sdata$smean1, + sdata$n1, + sigma, + data_distribution@two_armed) + else + t_test(sdata$smean1, + sdata$svar1, + sdata$n1, + data_distribution@two_armed) + if (length(statistics)>1) { + if (abs(sdata@n1 - design@n1 * (1L + data_distribution@two_armed))/ (design@n1 * (1L + data_distribution@two_armed)) > 0.1) + warning("Planned first-stage sample size 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){ + obs_n2 <- sdata@n2 + if (abs(obs_n2 - calc_n2 * (1L + data_distribution@two_armed))/ (calc_n2 * (1L + data_distribution@two_armed)) > 0.1) + warning("Planned second-stage sample size 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.") + } + if (test_val <= design@c1e & test_val >= design@c1f & sdata$n_stages==1L) + warning("Calculating statics for interim data, although the trial should continue to the second stage. Results may be unreliable.") + } results <- lapply(statistics, .analzye, data_distribution = data_distribution, use_full_twoarm_sampling_distribution = use_full_twoarm_sampling_distribution, diff --git a/R/evaluate_estimator.R b/R/evaluate_estimator.R index fcf2afd..9c76ba3 100644 --- a/R/evaluate_estimator.R +++ b/R/evaluate_estimator.R @@ -19,6 +19,19 @@ EstimatorScoreResult <- setClass("EstimatorScoreResult", slots = c(score = "list design = "TwoStageDesign", mu = "ANY", sigma = "numeric", results = "list", integrals = "list")) setClass("EstimatorScoreResultList", contains = "list") +EstimatorScoreResultList <- function(...) { + r <- list(...) + class(r) <- c("EstimatorScoreResultList", "list") + r +} +setMethod("c", signature("EstimatorScoreResult"), definition = + function(x, ...) { + EstimatorScoreResultList(x, ...) + }) +setMethod("c", signature("EstimatorScoreResultList"), definition = + function(x, ...) { + EstimatorScoreResultList(x, ...) + }) #' Evaluate performance characteristics of an estimator #' @@ -819,6 +832,96 @@ setMethod("evaluate_estimator", signature("TestAgreement", "IntervalEstimator"), conditional_integral) }) + +setClass("Centrality", contains = "PointEstimatorScore", + slots = list(interval = "IntervalEstimator")) +#' @rdname EstimatorScore-class +#' @export +Centrality <- function(interval = NULL) new("Centrality", label = sprintf("Centrality with respect to %s", toString(interval)), interval = interval) +#' @rdname evaluate_estimator-methods +setMethod("evaluate_estimator", signature("Centrality", "PointEstimator"), + function(score, + estimator, + data_distribution, + use_full_twoarm_sampling_distribution, + design, + true_parameter, + mu, + sigma, + tol, + maxEval, + absError, + exact, + early_futility_part, + continuation_part, + early_efficacy_part, + conditional_integral) { + design <- TwoStageDesignWithCache(design) + stagewise_estimators <- get_stagewise_estimators(estimator = estimator, + data_distribution = data_distribution, + use_full_twoarm_sampling_distribution = use_full_twoarm_sampling_distribution, + design = design, sigma = sigma, exact = exact) + stagewise_intervals <- get_stagewise_estimators(estimator = score@interval, + data_distribution = data_distribution, + use_full_twoarm_sampling_distribution = use_full_twoarm_sampling_distribution, + design = design, sigma = sigma, exact = exact) + g1 <- stagewise_estimators[[1L]] + g2 <- stagewise_estimators[[2L]] + l1 <- stagewise_intervals[[1L]] + u1 <- stagewise_intervals[[2L]] + l2 <- stagewise_intervals[[3L]] + u2 <- stagewise_intervals[[4L]] + generate_g1 = \(tp) \(...) ((g1(...) - l1(...)) + (g1(...) - u1(...))) + generate_g2 = \(tp) \(...) ((g2(...) - l2(...)) + (g2(...) - u2(...))) + .evaluate_estimator( + score, + estimator, + data_distribution, + use_full_twoarm_sampling_distribution, + design, + generate_g1, + generate_g2, + true_parameter, + mu, + sigma, + tol, + maxEval, + absError, + exact, + early_futility_part, + continuation_part, + early_efficacy_part, + conditional_integral) + }) + +#' @importFrom ggplot2 ggplot scale_x_continuous geom_line facet_wrap +#' @importFrom latex2exp TeX +setMethod("plot", signature = "EstimatorScoreResultList", definition = + function(x, ...) { + dat <- data.frame() + for (estimatorscoreresult in x) { + plot_list <- names(estimatorscoreresult@results) + for (score in plot_list){ + dat <- rbind(dat, + data.frame(mu = estimatorscoreresult@mu, + Score = estimatorscoreresult@results[[score]], + Estimator = toString(estimatorscoreresult@estimator), + score_name = score + ) + ) + } + } + ggplot(data = dat, mapping = aes(x = mu, y = Score, col = Estimator)) + + scale_x_continuous(name = TeX("$\\mu$")) + + geom_line() + + facet_wrap(vars(score_name), scales = "free_y") + }) +setMethod("plot", signature = "EstimatorScoreResult", definition = + function(x, ...) { + l <- EstimatorScoreResultList(x) + plot(l, ...) + }) + #' @importFrom latex2exp TeX plot_rejection_regions <- function(estimators, data_distribution, design, mu, sigma, subdivisions = 100, ...){ design <- TwoStageDesignWithCache(design) diff --git a/R/print.R b/R/print.R index fbf6b3a..cad7a5d 100644 --- a/R/print.R +++ b/R/print.R @@ -65,11 +65,17 @@ setMethod("toString", signature("Results"), right <- x@summary_data$n_stages lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + left <- "Observed n1 (in total)" + nval1 <- x@summary_data$n1 + right <- format(nval1) + lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + test_str <- if (is(x@data_distribution, "Normal")) "Z1" else "T1" if (x@summary_data$n_groups == 2L) n1 <- c(x@summary_data$n_s1_g1, x@summary_data$n_s1_g2) else n1 <- x@summary_data$n1 + test_val <- if (is(x@data_distribution, "Normal")) z_test(x@summary_data$smean1, @@ -85,7 +91,32 @@ setMethod("toString", signature("Results"), right <- format(test_val, digits=3) lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + left <- "Interim decision:" + if (test_val > x@design@c1e) { + right <- "reject null (early efficacy stop)" + } else if (test_val< x@design@c1f) { + right <- "accept null (early futility stop)" + } else { + right <- "continue to second stage" + } + lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + + left <- "Calculated n2(Z1) (per group)" + nval2 <- n2(x@design, test_val) + right <- format(nval2) + lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + + left <- "Calculated c2(Z1)" + cval2 <- c2(x@design, test_val) + right <- format(cval2, digits=3) + lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + if (x@summary_data$n_stages==2L){ + left <- "Observed n2 (in total)" + nval2 <- x@summary_data$n2 + right <- format(nval2) + lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") + test_str <- if (is(x@data_distribution, "Normal")) "Z2" else "T2" if (x@summary_data$n_groups == 2L) n2 <- c(x@summary_data$n_s2_g1, x@summary_data$n_s2_g2) @@ -106,12 +137,7 @@ setMethod("toString", signature("Results"), right <- format(test_val2, digits=3) lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") - left <- "c2(Z1)" - cval2 <- c2(x@design, test_val) - right <- format(cval2, digits=3) - lines[[length(lines)+1L]] <- paste0(pad_middle(left, right), "\n") - - left <- "Test decision:" + left <- "Final test decision:" if (test_val2 > cval2) { right <- "reject null" } else { @@ -128,7 +154,7 @@ setMethod("toString", signature("Results"), lines[[length(lines)+1L]] <- paste0("Stage 1 results:\n") print_header <- FALSE } - left <- paste0(" ", res$estimator@label, collapse="") + left <- paste0(" ", res$statistic@label, collapse="") right <- format(res$stage1, ...) lines[[length(lines)+1L]] <- paste0(pad_middle(paste0(left, ":"), right), "\n") } diff --git a/vignettes/Introduction.Rmd b/vignettes/Introduction.Rmd index 5195355..79fcb85 100644 --- a/vignettes/Introduction.Rmd +++ b/vignettes/Introduction.Rmd @@ -205,18 +205,24 @@ Keep in mind that a difference of $\mu=0.3$ was used in the simulation. Note that while the median unbiased estimator performs well in this particular example, this is not universally true. - # More details on the evaluation of performance characteristics +```{r} +a <- evaluate_estimator(MSE(), estimator = SampleMean(), data_distribution = Normal(), design = design, mu=mu, sigma=1) +b <- evaluate_estimator(MSE(), estimator = AdaptivelyWeightedSampleMean(), data_distribution = Normal(), design = design, mu=mu, sigma=1) +plot(c(a,b)) +``` - - - - - +```{r} +plot_p(estimator = LikelihoodRatioOrderingPValue(), + data_distribution = Normal(), + design = design, + mu = 0, + sigma = 1) +```