Skip to content

Commit

Permalink
TestAgreement for p-values
Browse files Browse the repository at this point in the history
  • Loading branch information
@jan-imbi
jan-imbi committed Mar 25, 2024
1 parent 6d4fd65 commit 13823c8
Show file tree
Hide file tree
Showing 2 changed files with 104 additions and 11 deletions.
39 changes: 29 additions & 10 deletions R/estimators.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1265,8 +1265,8 @@ setMethod("get_stagewise_estimators", signature("MLEOrderingPValue", "Normal"),
design,
sigma,
exact) {
g1 <- Vectorize(\(design, smean1, n1, sigma, two_armed, ...) p1_ml(design, smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, sigma, two_armed, ...) p2_ml(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
g1 <- Vectorize(\(design, smean1, n1, mu, sigma, two_armed, ...) p1_ml(design, smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) p2_ml(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
Expand All @@ -1282,8 +1282,8 @@ setMethod("get_stagewise_estimators", signature("LikelihoodRatioOrderingPValue",
design,
sigma,
exact) {
g1 <- Vectorize(\(design, smean1, n1, sigma, two_armed, ...) p1_lr(design, smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, sigma, two_armed, ...) p2_lr(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
g1 <- Vectorize(\(design, smean1, n1, mu, sigma, two_armed, ...) p1_lr(design, smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) p2_lr(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
Expand All @@ -1299,8 +1299,8 @@ setMethod("get_stagewise_estimators", signature("ScoreTestOrderingPValue", "Norm
design,
sigma,
exact) {
g1 <- Vectorize(\(smean1, n1, sigma, two_armed, ...) p1_st(smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, sigma, two_armed, ...) p2_st(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
g1 <- Vectorize(\(smean1, n1, mu, sigma, two_armed, ...) p1_st(smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) p2_st(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
Expand All @@ -1316,8 +1316,8 @@ setMethod("get_stagewise_estimators", signature("StagewiseCombinationFunctionOrd
design,
sigma,
exact) {
g1 <- Vectorize(\(smean1, n1, sigma, two_armed, ...) p1_sw(smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, sigma, two_armed, ...) p2_sw(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
g1 <- Vectorize(\(smean1, n1, mu, sigma, two_armed, ...) p1_sw(smean1, n1, mu = 0, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) p2_sw(design, smean1, smean2, n1, n2, mu = 0, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
Expand All @@ -1335,8 +1335,27 @@ setMethod("get_stagewise_estimators", signature("NeymanPearsonOrderingPValue", "
design,
sigma,
exact) {
g1 <- Vectorize(\(design, smean1, n1, sigma, two_armed, ...) p1_np(design, smean1, n1, mu = 0, mu0=estimator@mu0, mu1=estimator@mu1, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, sigma, two_armed, ...) p2_np(design, smean1, smean2, n1, n2, mu = 0, mu0=estimator@mu0, mu1=estimator@mu1, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
g1 <- Vectorize(\(design, smean1, n1, mu, sigma, two_armed, ...) p1_np(design, smean1, n1, mu = 0, mu0=estimator@mu0, mu1=estimator@mu1, sigma, two_armed, ...), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) p2_np(design, smean1, smean2, n1, n2, mu = 0, mu0=estimator@mu0, mu1=estimator@mu1, sigma, two_armed, ...), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
setClass("NaivePValue", contains = "VirtualPValue")
#' @rdname PValue-class
#' @export
NaivePValue <- function() new("NaivePValue", label = "Naive p-value")
#' @rdname get_stagewise_estimators
setMethod("get_stagewise_estimators", signature("NaivePValue", "Normal"),
function(estimator,
data_distribution,
use_full_twoarm_sampling_distribution,
design,
sigma,
exact) {
g1 <- Vectorize(\(smean1, n1, mu, sigma, two_armed, ...) pnorm(smean_to_z(smean1, n1, sigma, two_armed), lower.tail=FALSE), c("smean1"))
g2 <- Vectorize(\(design, smean1, smean2, n1, n2, mu, sigma, two_armed, ...) pnorm(smean_to_z(smeans_to_smean(smean1 = smean1, smean2 = smean2, n1 = n1, n2 = n2),
n = n1 + n2, sigma = sigma, two_armed = two_armed),
lower.tail=FALSE), c("smean1", "smean2", "n2"))
list(g1 = g1,
g2 = g2)
})
Expand Down
76 changes: 75 additions & 1 deletion R/evaluate_estimator.R
View file
Original file line number Diff line number Diff line change
Expand Up @@ -42,7 +42,7 @@
setClass("EstimatorScore", slots = c(label = "character"))
setClass("PointEstimatorScore", contains = "EstimatorScore")
setClass("IntervalEstimatorScore", contains = "EstimatorScore")
EstimatorScoreResult <- setClass("EstimatorScoreResult", slots = c(score = "list", estimator = "Estimator", data_distribution = "DataDistribution",
EstimatorScoreResult <- setClass("EstimatorScoreResult", slots = c(score = "list", estimator = "Statistic", data_distribution = "DataDistribution",
design = "TwoStageDesign", mu = "ANY", sigma = "numeric",
results = "list", integrals = "list"))
setClass("EstimatorScoreResultList", contains = "list")
Expand Down Expand Up @@ -930,6 +930,80 @@ setMethod("evaluate_estimator", signature("TestAgreement", "IntervalEstimator"),
conditional_integral)
})

#' @rdname evaluate_estimator-methods
setMethod("evaluate_estimator", signature("TestAgreement", "PValue"),
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) {
if (missing(true_parameter))
true_parameter <- 0.05
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)
g1 <- stagewise_estimators[[1L]]
g2 <- stagewise_estimators[[2L]]

if (is(data_distribution, "Student")) {
generate_g1 = \(tp) \(design, smean1, svar1, n1, two_armed, ...) {
t1 <- smean_to_t(smean1, svar1, n1, two_armed)
!xor(design@c1e < t1, g1(design = design, smean1 = smean1, svar1 = svar1, n1 = n1, two_armed = two_armed, ...) < true_parameter)
}
generate_g2 <- \(tp) \(design, smean1, svar1, smean2, svar2, n1, n2, two_armed, ...) {
t1 <- smean_to_t(smean1, svar1, n1, two_armed)
t2 <- smean_to_t(smean2, svar2, n2, two_armed)
c2 <- c2_extrapol(design, t1)
!xor(c2 < t2, g2(design = design, smean1 = smean1, svar1 = svar1, smean2 = smean2, svar2 = svar2, n1 = n1, n2 = n2, two_armed = two_armed, ...) < true_parameter )
}
} else if (is(data_distribution, "Normal")) {
generate_g1 = \(tp) \(design, smean1, n1, sigma, two_armed, ...) {
z1 <- smean_to_z(smean1, n1, sigma, two_armed)
!xor(design@c1e < z1, g1(design = design, smean1 = smean1, n1 = n1, sigma = sigma, two_armed = two_armed, ...) < true_parameter)
}
generate_g2 <- \(tp) \(design, smean1, smean2, n1, n2, sigma, two_armed, ...) {
z1 <- smean_to_z(smean1, n1, sigma, two_armed)
z2 <- smean_to_z(smean2, n2, sigma, two_armed)
c2 <- c2_extrapol(design, z1)
!xor(c2 < z2, g2(design = design, smean1 = smean1, smean2 = smean2, n1 = n1, n2 = n2, sigma = sigma, two_armed = two_armed, ...) < true_parameter)
}
} else {
stop("This data_distribution class is not supported.")
}
.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)
})


setClass("Centrality", contains = "PointEstimatorScore",
slots = list(interval = "IntervalEstimator"))
Expand Down

0 comments on commit 13823c8

Please sign in to comment.