Skip to content

Commit

Permalink
Make plot function and fix analayze output
Browse files Browse the repository at this point in the history
  • Loading branch information
@jan-imbi
jan-imbi committed Sep 19, 2023
1 parent abf9933 commit 30d8ea4
Show file tree
Hide file tree
Showing 4 changed files with 174 additions and 14 deletions.
27 changes: 26 additions & 1 deletion R/analyze.R
View file
Expand Up @@ -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,
Expand All @@ -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,
Expand Down
103 changes: 103 additions & 0 deletions R/evaluate_estimator.R
View file
Expand Up @@ -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
#'
Expand Down Expand Up @@ -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)
Expand Down
40 changes: 33 additions & 7 deletions R/print.R
View file
Expand Up @@ -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,
Expand All @@ -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)
Expand All @@ -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 {
Expand All @@ -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")
}
Expand Down
18 changes: 12 additions & 6 deletions vignettes/Introduction.Rmd
View file
Expand Up @@ -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)
```



Expand Down

0 comments on commit 30d8ea4

Please sign in to comment.