daris_label_y, info_threshold_label_y,
historical_label_y, endpoint_label_y and events_label_y
still override them. The vertical stacking order of the four DARIS
labels is mirrored, not changed.?plot.tsa_hr); a test covers both signs."*** NOTE (HKSJ INFERENCE): ..." caveat printed by tsa_hr(verbose = TRUE) under re_inference = "hksj"/"hksj_adhoc" is now separated
from the "=== Trial sequential monitoring boundaries ..." table above
it by a blank line (0.2.8.12/13 had deliberately printed it directly
after, with no blank line there -- see NEWS for those versions)..tsahr_color_enabled()/.tsahr_bg_note() (rtsa_engine.R). Colouring
is opt-in by detection: it only fires in a session that looks like an
interactive terminal or RStudio console with ANSI support, so captured
output (R CMD check, testthat, Rscript, sink()/
utils::capture.output()) stays plain text, matching prior behaviour.
options(tsahr.color = TRUE)/options(tsahr.color = FALSE) overrides
the auto-detection in either direction.test-re-inference-0.2.8.12.R updated for the new blank-line placement;
a new test-console-color-0.2.8.16.R covers the colour-detection
override and confirms the default (colour off) output has no ANSI codes."Add'l events to AR endpoint (hist. rate; median of the study-level information per event)" (and the analogous "design route"
DARIS row). Both now read "Add'l events to DARIS (hist. rate; <basis> IPE)" / "Add'l events to AR endpoint (hist. rate; <basis> IPE)",
where <basis> is "pooled" or the projection_stat statistic (e.g.
"median"), matching the style of the neighbouring
"...(<basis>-based proj.)" row. The new abbreviation "IPE"
("information per event") is added to the "Abbreviations" legend
returned as attr(summary_table, "abbreviations") and printed by
summary.tsa_hr(). The full-sentence description (e.g. "pooled ratio: total information / total events") still appears in the prose
projection_note; only the summary-table row label was shortened.test-re-inference-0.2.8.12.R (R CMD check --as-cran
failure on 0.2.8.13): the test checked for the substring
"not interpreted as directly comparable", but the actual, correct
printed HKSJ caveat says "not be interpreted as directly comparable"
(see ?tsa_hr and NEWS.md for 0.2.8.12) -- the test was missing "be".
No R/ code changed in this release.=== Trial sequential monitoring boundaries (alpha/beta spending) ===
table and the note itself (0.2.8.12 had a blank line there); a single
blank line still follows the note (or the table, when the note is not
shown), as before.NA first-look Z (under "hksj"/"Hksj_adhoc"; see 0.2.8.12)
removed explicitly, rather than relying only on the layers' own
na.rm = TRUE. This fixes an R CMD check --as-cran test failure where
the built plot still contained a point for the NA look on some
ggplot2/R combinations.Printed HKSJ caveat. With verbose = TRUE and a non-standard
re_inference ("hksj"/"knha" or "Hksj_adhoc"), a highlighted note is
now printed directly after the === Trial sequential monitoring boundaries (alpha/beta spending) === table:
HKSJ inference can be unstable at early cumulative looks because the degrees of freedom are k-1. In particular, the HKSJ statistic is undefined at k=1 ("NA" at first look) and is based on only one degree of freedom at k=2. Early cumulative HKSJ values should therefore not be interpreted as directly comparable in magnitude with conventional normal Z-statistics.
The same caveat is now in ?tsa_hr ("Random-effects inference" under
Details).
cumulative$Z is NA at the first look for "hksj"/"Hksj_adhoc"
(the HKSJ statistic is undefined at k = 1; a t distribution on 0 df has no
defined quantile). This shows as NA in both printed cumulative tables
and is simply not drawn on the plotted Z-curve (the point is skipped, the
line does not break elsewhere); cumulative$se/$pval/etc. at that look
are unaffected, only Z is withheld. results$crossed_conventional is
now computed with na.rm = TRUE so a first-look NA cannot propagate
into it.
Readable summary table. summary_table$Parameter now uses short
abbreviations (e.g. "DARIS", "AR endpoint" for analysis-route
endpoint, "RE" for random effects, "psi", "hist. rate", "Add'l")
instead of the long descriptive labels of earlier versions, so the table
no longer wraps illegibly at typical console widths. The expansions are
returned as a named character vector in
attr(summary_table, "abbreviations") and summary() now prints an
Abbreviations: ... line directly beneath the table. print() is
unaffected (it already used short, separate lines). See ?tsa_hr, "Value".
tsa_hr(re_inference = ) (new, last argument, so positional calls are
unaffected). "standard" (default) is the behaviour of all earlier
versions and gives identical results. "hksj" (alias "knha") applies the
Hartung-Knapp-Sidik-Jonkman adjustment (variance multiplied by q, t
distribution with k - 1 df); "Hksj_adhoc" (alias "knha_adhoc") uses
max(1, q) so the standard error is never smaller than the standard one.
Matching is case-insensitive; anything else is an error.Z is the normal-equivalent of the HKSJ t
statistic (same two-sided p-value), so the boundaries and the conventional
boundary keep their meaning; the raw t statistic is in cumulative$zval,
and non-standard options add cumulative$re_df and cumulative$re_scale.
The first look (one study) and looks with a zero/non-finite scale factor
keep the z-based values.method), pooled point
estimate, tau^2, I^2, Q, Diversity D^2, DARIS, information fractions and
therefore the alpha/beta boundaries. D^2 is deliberately computed from the
standard random-effects variance, not the HKSJ-scaled one.print(), summary(). With a non-standard option the
caption gets a line "Random-effects inference: ..." under the first
"Methods" line; print() prints the option and summary() gains a
"Random-effects inference" row. Nothing is added for "standard".parameters$re_inference (normalised) and parameters$re_inference_requested
are returned. If the HKSJ scale factor of the full data is zero or not
finite, tsa_hr() warns and falls back to "standard".?tsa_hr (argument, Details, Value) and ?plot.tsa_hr
updated; new tests in test-re-inference-0.2.8.11.R.useDynLib(tsahr, .registration = TRUE) and importFrom(Rcpp, sourceCpp) with no corresponding roxygen tags
anywhere in R/*.R, and an importFrom(utils, globalVariables) that
nothing requested and that isn't needed (the call in zzz.R is already
fully qualified as utils::globalVariables()). Added the missing
@useDynLib/@importFrom Rcpp sourceCpp tags (on a new "_PACKAGE"
sentinel in zzz.R) and regenerated NAMESPACE in standard roxygen2
output format; confirmed against devtools::document() with roxygen2
8.1.0 (matching Config/roxygen2/version in DESCRIPTION).xmax_mult's roxygen doc (?plot.tsa_hr) now lists the
historical-rate projection among the values that must fit in the plot's
x-axis; it was added to the plot in 0.2.8.7/0.2.8.8 but the doc line was
never updated.test-tsa_hr.R used to wrap
tsa_hr() in tryCatch(..., error = function(e) NULL) and skip_if on a
NULL result, so the suite could pass whether or not zero-event handling
was broken. Zero-event studies are a documented, deliberate feature (see
"Zero-event studies" in ?tsa_hr; they're excluded from the events
projection only, and the exclusion is reported), so both tests now call
tsa_hr() unconditionally and assert expect_s3_class(res, "tsa_hr") --
they fail instead of skipping if this behavior ever regresses.plot() draws as
the "DARIS (historical rate)" line (accrued events + estimated additional
events). summary() gets a matching row. This mirrors the analysis route's
"Historical information/event-rate projection" line.summary() gets a matching row.boundary_route = "analysis" and the analysis-route endpoint
(design_R x DARIS) not reached, plot() now also draws an
"Analysis-route endpoint (historical rate)" vertical line at accrued events
show_historical_daris, historical_label_size, historical_label_x and
historical_label_y now apply to both routes (design: "DARIS (historical
rate)").summary() (analysis route) gain the same
theoretical/historical-rate pair; the former "Estimated additional events
required" line and "Estimated additional events to reach the analysis-route
endpoint (projected, approximate)" row are renamed accordingly.projection$additional_events_theoretical for both routes:
max(ceiling(DARIS_events * route_endpoint - events_accrued), 0); for the
design route it equals additional_events_required_design.Design route output. When DARIS has not been reached, tsa_hr() now
reports three separate quantities (printed output, summary(), plot
caption):
Theoretical additional events to DARIS (Schoenfeld): xxx
Estimated additional events to DARIS (historical rate): ~XXX
Estimated additional studies required: N
The first is the deterministic DARIS_events - events_accrued (floored at
0); the second is the information-per-event projection introduced in
0.2.8.5 with I_required = DARIS. They rest on different assumptions and
can disagree (the second can be positive when the first is 0).
Plot. For the design route with DARIS not reached, plot() also draws a
"DARIS (historical rate)" vertical line at accrued events + estimated
additional events (projection$target_events_historical_rate), next to the
theoretical DARIS line. New arguments show_historical_daris,
historical_label_size, historical_label_x, historical_label_y. The
x-axis range accounts for it.
New tsa_hr() argument info_per_event_basis ("per_study", default,
or "pooled"). "per_study" uses projection_stat (median by default) of
each study's information per event; "pooled" uses total information /
total events. It only affects the events projection; the studies
projection always uses projection_stat. Both variants are always
returned (projection$additional_events_study_level,
projection$additional_events_pooled, and the corresponding
study_level_info_per_event / pooled_info_per_event).
Zero-event studies. They are still accepted, but are excluded from the
information-per-event (events) projection (they carry information but no
events). The number excluded is now stated explicitly in the printed
output, in summary() (new row) and in the plot caption, and returned as
projection$n_zero_event_studies / n_excluded_events_projection. They
stay in the studies projection and in all other results.
Caveat strengthened (no change to the mathematics). The projection note and the documentation now state that all projections are linear extrapolations on the fixed-effect, study-level inverse-variance information scale, and do not model how random-effects weights or tau^2 would change as studies are added; they are indicative only.
New projection fields: info_per_event_basis, events_accrued,
study_level_info_per_event, additional_events_study_level,
target_events_historical_rate, n_studies, n_zero_event_studies,
n_excluded_events_projection. The summary-table row "Additional events
required to reach DARIS (Schoenfeld-scale)" is renamed "Theoretical
additional events to DARIS (Schoenfeld)".
Unchanged: the required information, boundaries, decisions, and the analysis route's events and studies figures on the default basis.
ceiling(shortfall / typical study information) x typical events per study, i.e. it was a consequence of the rounded-up number of studies. That
answers "how many events would that many whole studies bring?", not "how
many events are needed?" (e.g. shortfall 150, typical study 100 information
units / 500 events: 2 studies -> 1,000 events, whereas the information
projection is 150 / (100/500) = 750 events).projection$additional_events_estimated = shortfall / typical
information per event, where the typical value is the projection_stat
(median by default, or mean) of each study's own
(1/Std_Error^2) / total_events. It is continuous and rounded up only for
display (print, summary(), plot caption). Studies with no events are
ignored for this figure only.projection$n_additional_studies is unchanged (shortfall / typical study
information, rounded up) and is now explicitly the secondary quantity.projection fields: central_info_per_event, pooled_info_per_event
(total information / total events), additional_events_pooled (the events
figure under that ratio-of-sums, as a sensitivity check) and
events_from_whole_studies (the old chained figure, kept and clearly
named).DARIS_events - events_accrued), the required
information, all boundaries and decisions.tsa_hr() now projects how many more studies, and
roughly how many more events, would be needed to reach it -- a
retrospective, HR-specific counterpart to RTSA's own prospective
minTrial() / ris(..., type = "retrospective") machinery, built
directly from tsahr's own observed study-level information increments
rather than reimplementing RTSA's.boundary_route = "design": I_required = DARIS.boundary_route = "analysis": I_required = design_R x DARIS (the
analysis-route endpoint).1/Std_Error^2) and event count. New argument
projection_stat = c("median", "mean") selects which; median is the
default because it is more robust to one unusually large/small study.
n_additional_studies is always rounded up to a natural number
(minimum 1 whenever a genuine shortfall exists).DARIS_events - events_accrued (no projection needed there); only the
additional-studies estimate uses the median/mean projection. For the
analysis route, both the additional-events and additional-studies
figures are projected (the events figure is marked with ~, since it is
an approximation, unlike the design route's figure).verbose = TRUE), printed right after the existing
"Required information size (DARIS) reached" block, and new rows in
summary()'s table, shown only when the relevant endpoint has not been
reached.plot() now appends a line to the methods caption, below everything
else, when the endpoint has not been reached: "Estimated additional events required: xxxx, Estimated additional studies required: xx".projection element of the
tsa_hr() result: method, I_required, info_accrued,
additional_info_required, central_info_increment,
central_event_increment, n_additional_studies,
additional_events_estimated, and (design route only)
additional_events_required_design. All fields are NA once the
route's own endpoint has already been reached.print()'s
existing lines, or to any existing summary()/plot() output when the
endpoint has already been reached -- this is purely additive.boundary_route = "analysis", the vertical line and label of
the analysis-route endpoint (design_R x DARIS) were drawn only when that
endpoint had been reached in the observed data (results$final_reached).
When design_R > 1 the endpoint lies beyond the theoretical DARIS
event-equivalent and is usually not yet reached, so neither the vertical nor
the label appeared, although the formal boundaries still ended at
DARIS_events * design_R. (When design_R < 1 the endpoint is reached as soon
as DARIS is, so the marker was always shown.)DARIS_events * design_R when the endpoint has not been reached, labelled
"Analysis-route endpoint (... x DARIS) not yet reached; theoretical ~ N
events". The reached case is unchanged. endpoint_label_x/_y/_size apply to
both wordings.plot() argument xmax_mult (default 1.15, the previous hard-coded
value): multiplier applied to the largest x element to set the upper x-axis
limit, e.g. plot(res, xmax_mult = 1.5) for more room on the right.print() or summary().The alpha- and beta-boundary engines were independently validated against the corresponding RTSA implementations.
For the reference validation design:
beta value.In the reference test, cumulative and incremental beta spending agreed with
RTSA to a tolerance of 1e-12, and the finite beta boundaries agreed
numerically at all monitored looks (n=40).
These tests establish numerical equivalence of the tsahr alpha- and
beta-engine calculations with the corresponding RTSA reference
implementation for the validated configuration. They should not be
interpreted as a formal proof of equivalence for every possible parameter
combination.
<https://ctu.dk/tsa/> returned a 404; replaced with the correct
<https://ctu.dk/tools> (also fixed in README.md).Author: field. It duplicated (and
had drifted from) the field R derives automatically from Authors@R,
triggering a checking DESCRIPTION meta-information ... NOTE. Authors@R
is now the sole source of truth for the Author field.plot(): pooled-effect line in the subtitle. Below the existing subtitle
line ("Random-effects model | Diversity D² | Anticipated HR = ... | ...") the
plot now shows a second line: Pooled HR = 0.51 [95% CI: 0.48, 0.54] | p < 0.001 | Tau² = 0.0136 | I² = 73.3%. The values are the ones print() and
summary() already report (res$res_re, res$heterogeneity); nothing is
recomputed. HR and CI use 2 decimals, tau² 4, I² 1; p-values below 0.001 are
shown as p < 0.001, others with 3 decimals. The "2" in Tau² and I² is the
Unicode superscript two, as for Diversity D². Internal helper
.tsahr_pooled_subtitle().plot(): endpoint_label_x, endpoint_label_y, endpoint_label_size.
Position (data coordinates: cumulative events, Z-score) and font size of the
"Analysis-route endpoint (Design_R x DARIS) reached" label, which is drawn
only for boundary_route = "analysis" when the endpoint was reached. The
defaults reproduce the previous plot exactly: the label sits just right of
its vertical line at 58% of the upper y limit, and endpoint_label_size = NULL still follows info_threshold_label_size (3.2 by default). Same
conventions as the other *_label_x/_y/_size arguments.inst/extdata/HR_meta_example.xlsx (10 studies) is replaced by HR_meta.xlsx
(20 studies) and HR_meta_2.xlsx (40 studies). tsahr_example_data() gains
a dataset argument, "HR_meta" (default) or "HR_meta_2"; calling it with
no argument now returns the 20-study file. The new sheets already use
underscores in their headers (Events_Treatment, ...). Examples in the
documentation and README use the new data.tests/testthat/testdata/HR_meta_legacy_10studies.xlsx and read
through legacy_example_data() (tests/testthat/helper-data.R); those tests
are otherwise untouched. New tests in test-plot-and-example-data-0.2.8.R
cover the two datasets (shape, columns, tsa_hr() and plot() run), the
subtitle line and the endpoint-label arguments. No numbers are pinned for the
new datasets.tsa_hr()'s
results or any decision field.spend == beta knife-edge fixed. When every interim look is suppressed
(rm_bs = nt - 1: small information fractions, "evidence still insufficient"),
the final look carries the whole beta spend, and RTSA's exact-float shortcut
else if (spend == beta) za = 0 fires or not depending on the last bit of
beta. RTSA is normally called with a literal beta = 0.2, which differs from
its own spend arithmetic in the last bit; tsahr computes beta = 1 - power
(0.19999999999999996), which equals it. When the shortcut fired, za was 0
whatever the information scale, so the calibration gap was constant and the
root search reported "no root bracket": a design that calibrates at
power 0.8 + 1e-12 failed at 0.8 (and at 0.95; 0.9 passed), the compiled engine
was lost and tsa_hr() fell back to the legacy engine (warning, banner, 12-23
s). The shortcut is no longer ported, at the first or later looks (a
deliberate, documented departure from RTSA; za = 0 has no statistical
meaning there, and it is unreachable in everything the frozen RTSA reference
covers, so parity is unchanged). Reproduced with the compiled core on a 30-look
schedule up to 26% of the required information: before, beta = 1 - 0.8 and
1 - 0.95 failed; after, powers 0.80-0.99 and their +/-1e-12 neighbours all
calibrate to the same root (1.0000404 at power 0.80). Tests:
test-robust-schedules.R (roots by power and last-bit perturbations; an
end-to-end low-information tsa_hr() run that must use the compiled engine at
powers 0.80/0.90/0.95); tools/cpp_vs_py.py section 11. The legacy R engine
keeps RTSA's shortcuts (it did not fail on the audited design; unaudited for
other powers) and is reached only if the compiled engine fails.
Found by an external audit of 0.2.7.21.?tsa_hr and pinned by a test. Reported values are unchanged.License: GPL (>= 2); the
MIT LICENSE file is removed; Authors@R credits Anne Lyngholm Soerensen,
Markus Harboe Olsen, Theis Lange and Christian Gluud (ctb, cph);
Copyright: inst/COPYRIGHTS documents provenance file by file, including RTSA's
origin in the Copenhagen Trial Unit's TSA software and its manual's authors;
the derived source files carry a copyright/GPL notice; README has an
attribution section.summary_table$Value is a character column. It was numeric only because
c() coerced logicals, so decision rows printed as 1/0. Numbers keep their
rounding; logicals print TRUE/FALSE/NA. (Anything that used
summary_table$Value numerically must convert.)%d overflow fixed. Event counts derived from the required information
can exceed the integer range (e.g. target_HR = 0.9999); the verbose output,
print() and the DARIS lines now use %.0f.DESCRIPTION: ggplot2 (>= 3.4.0) (the linewidth aesthetic needs it).legacy_fallback = TRUE). Reworded. The rtsa_core.h header no
longer claims "everything is bounds-checked" (array lengths are validated on
entry; inner loops use unchecked indexing).Documentation and test release for the LEGACY fallback engine. The default compiled paths are untouched, and no shipped code path changes.
.obf_alpha_boundary(), default n_grid = 16000) gives 4.332634 2.963388
2.358980 2.012955 against live RTSA 0.2.2's 4.332634 2.963131 2.359044
2.014090: errors 4e-7, 2.6e-4, -6.4e-5 and -1.14e-3, largest at the final
look (the value the futility root search targets). The legacy warp root is
1.132483 against RTSA's 1.133242 (error 7.6e-4), where before the 0.2.7.19
fix it was 1.0972 (error 3.6e-2) -- a ~47-fold reduction. Feeding the same
FFT alpha vector to an independent Simpson beta implementation gives root
1.132483339, i.e. exactly what the legacy R beta code reported: the legacy
beta engine is exact given its alpha input, and the whole residual error comes
from the FFT alpha. One schedule only -- not a general accuracy bound for the
fallback, which remains labelled approximate and "NOT comparable with RTSA".
This also answers the n_grid = 16000 question left open in the
.obf_alpha_boundary() notes for this schedule (the 5-look snippet there has
still not been run)..obf_alpha_boundary()
was called "the same validated side = 2 recursion"; it is the FFT
approximation of RTSA's Simpson recursion (error above), not a validated
reproduction. (2) The equality "final boundary = qnorm(1 - alpha/2)" was said to
hold "only in a continuous-monitoring limit"; it holds only for a single look.
With more looks the final wall is larger and keeps growing (2.014 for 4 even
looks, 2.185 for 100), which is exactly why assuming the constant was a bug..tsahr_legacy_boundaries() (it
still said the futility route substitutes qnorm(1 - alpha/2)), and the
headers of the two legacy test files (they said the qnorm(1 - alpha/2) final
boundary asserted in them is the legacy engine's convention; since 0.2.7.19 it
survives only for single-look schedules and the analysis wrapper's design_R
endpoint).test-rtsa-live-reference.R): the legacy engine on the frozen
reference schedule stays within 3e-3 of live RTSA for the alpha bounds and
futility bounds and within 2e-3 for the warp root, and its final wall is no
longer qnorm(1 - alpha/2). The tolerances are set from the measurements above
and separate the fixed engine (errors <= 1.14e-3) from the unfixed one (root
3.6e-2, futility bounds 3.4e-2 to 5.4e-2 off).0.2.7.19's production fix (.rtsa_beta_boundary() recomputing the true,
schedule-dependent final efficacy boundary via .obf_alpha_boundary()
instead of assuming qnorm(1-alpha/2)) was correct, but two tests that
still encoded the old assumption were not actually updated in the
uploaded package, and a real R CMD check run caught both:
test-boundaries-rtsa.R, "0.2.7.4: the information-scale root search
actually engages and hits its target" -- gap_ans was checked against
an alpha_ref vector whose last entry was still the stale
qnorm(1-alpha/2) constant, while ans$warp_root (from
.rtsa_beta_boundary()) was by then solving for the CORRECT,
recomputed final value -- two different equations, so their roots no
longer agreed (off by 0.072). Fixed by deriving alpha_ref from
.obf_alpha_boundary(), the same recursion .rtsa_beta_boundary()
itself now uses internally, so both sides check the same equation.test-tsa_hr.R, "LEGACY R-only boundary engine scales to many studies"
-- asserted the final futility boundary equals the fixed constant 1.96
for a 40-look schedule; the true value there is materially different
(~2.16). Fixed by asserting against the schedule's own true final
alpha boundary (c40's own last entry, already computed earlier in
the same test) instead of a hardcoded constant, so this test cannot
re-encode a wrong assumption about what that value is.Two more tests were fixed for correctness even though R CMD check
did not flag them as failing, because they were passing for the wrong
reason: .rtsa_beta_boundary()'s own defensive
pmin(boundary, c_vec_alpha) safeguard clips its (now correct, and
therefore LARGER) recomputed final boundary back down to whatever
c_vec_alpha's own last entry says, so a test that still supplied the
stale, smaller qnorm(1-alpha/2) as c_vec_alpha's last entry got that
same stale value back out of ans$boundary -- appearing to confirm the
old convention while actually only exercising the safety clip, not the
recomputation the test's own comment claimed to check:
test-boundaries-rtsa.R, "RTSA beta boundary returns NA only for
negligible early spend, not for every non-positive value"test-boundaries-rtsa.R, "0.2.7.5: the two-pass root search converges
on a realistic multi-look design"Both now derive their expected final value from .obf_alpha_boundary()
directly, matching the pattern the two already-fixed 0.2.7.19 tests
("RTSA retrospective inner-wedge engine has definitive alpha boundary"
and "RTSA over-powered analysis uses the true (recursion-based)
definitive boundary") already established. One further qnorm(1-alpha/2)
occurrence, in "0.2.7.5: a single-look schedule (DARIS already reached at
the first study) does not error", was checked and left as-is: with a
single look, RTSA's discretised recursion collapses to the ordinary
fixed-sample critical value, so qnorm(1-alpha/2) genuinely is the
correct value there (verified directly against the compiled engine),
not an unfixed instance of the bug. .rtsa_beta_boundary_analysis()'s
own, separately-disclosed use of the same constant (its harder
design_R-endpoint case) is unaffected, as previously noted.
A live RTSA cross-check caught a real bug in the LEGACY, R-only fallback
engine (.rtsa_beta_boundary() in R/obf_boundaries.R, used only when
legacy_fallback = TRUE and the compiled engine fails -- the default,
compiled boundary_route = "design"/"analysis" paths are NOT affected).
.rtsa_beta_boundary() assumed "by construction of any
properly normalised two-sided alpha-spending function, the efficacy
boundary at t = 1 is exactly qnorm(1-alpha/2)", and substituted that
constant for the final entry of the alpha wall it root-finds the
information-scale factor (warp_root) against. That assumption is wrong
for RTSA's actual discretised O'Brien-Fleming-type spending recursion --
the equality holds only for a SINGLE look (no earlier spending); with more looks the final wall is larger and grows with the number of looks (2.014 for 4 even looks, 2.185 for 100) -- it depends on the schedule, not on a continuous-monitoring limit. [Corrected in 0.2.7.21; this entry originally attributed the equality to a continuous-monitoring limit.] A live RTSA reconstruction
confirmed this concretely: for one real schedule, RTSA's own final alpha
boundary was 2.014090377368289, not qnorm(1-alpha/2) = 1.959963984540054 -- a difference large enough to materially shift the
solved root (1.133241903483384 vs. the wrong 1.097192717591548) and,
through it, every futility boundary derived from that root (e.g.
0.632531.../1.404162... vs. the wrong 0.598718.../1.362718...).
A second line (za_seq[length(za_seq)] <- final_alpha_bound) then
independently re-asserted the same wrong constant onto the reported
final futility boundary..rtsa_beta_boundary() no longer approximates the t = 1
efficacy boundary at all. It recomputes the TRUE, schedule-dependent
value via .obf_alpha_boundary() -- the same side = 2 FFT recursion used everywhere else in this legacy engine (an APPROXIMATION of RTSA's Simpson recursion, not a validated reproduction of it; see 0.2.7.21) -- run on timing_beta
itself (the observed pre-1 looks plus the definitive t = 1 point),
rather than reusing whichever c_vec_alpha the caller happened to
already have. This is correct and self-contained whether or not the
caller's t already contains an exact t = 1 entry, and the previously-
separate za_seq[...] <- final_alpha_bound override line is now
harmless (it re-asserts the same, now-correct value the root search
already converges to) rather than removed, since the bug was in what
final_alpha_bound held, not in tightening to it..rtsa_beta_boundary_analysis() (the legacy engine's
retrospective-analysis path) calls the now-fixed .rtsa_beta_boundary()
internally for its own design_R calibration, so that part is fixed
automatically. Its own, separate qnorm(1-alpha/2) approximation for
the design_R endpoint's alpha value (a materially harder problem,
since design_R essentially never coincides with an observed look, and
a correct fix requires porting RTSA's type = "analysis" alpha
recursion rather than reusing the type = "design" one) is UNCHANGED
and remains a disclosed approximation -- consistent with this engine's
existing "NOT comparable with RTSA" warning whenever it is actually
used.test-boundaries-rtsa.R that had encoded the old, wrong
qnorm(1-alpha/2) convention as correct behaviour are corrected to
check against a direct .obf_alpha_boundary() recomputation instead
(which is what the fixed function now does internally), and to assert
that the true value differs materially from the old approximation --
regression protection against reintroducing this exact bug.2.014090377368289/1.133241903483384
figures above are taken directly from the live RTSA reconstruction that
reported this bug. The corrected formula is a direct, minimal port of
the same .obf_alpha_boundary() call already used and validated
elsewhere in this file.Documentation/diagnostic-text release following the 0.2.7.17 assessment. No boundary value, root, test reference or code path changes.
qmax < target - tol (tol = 1e-15, plus floating-point summation error): the
classification is decided directly from the analytic limit qmax = sum(last),
up to the numerical tolerance, not mathematically exactly. Reworded in
src/rtsa_core.h, R/rtsa_engine.R, inst/REVERSE_ENGINEERING_RTSA.md and
the 0.2.7.15/16 NEWS entries.[lo, lo + 1e-8], and the
warning now says so. While checking this I found the collapse warning's "each
was widened automatically" was also inaccurate: only intervals with no positive
width are widened; a zero-width two-node grid is merely counted (arithmetic
unchanged) and a collapsed grid whose interval already has positive width is
used as is. The warning and the z_n_w() header comment now describe exactly
that.inst/extdata/rtsa_0.2.2_reference.R now records the
maintainer's R version and platform (R 4.6.1 (2026-06-24 ucrt),
x86_64-w64-mingw32), flagged for correction if the RTSA run was made in a
different session.inst/REVERSE_ENGINEERING_RTSA.md "Current implementation"
heading now spans 0.2.7.11 - 0.2.7.18.The open item from 0.2.7.12 onward is closed for the case supplied. Numbers
printed by RTSA 0.2.2 itself (RTSA::boundaries(), timing 0.25/0.50/0.75/1.00,
alpha 0.05, beta 0.20, side 2, non-binding, esOF) are frozen in
inst/extdata/rtsa_0.2.2_reference.R (provenance and call shapes documented in
the file) and tested in tests/testthat/test-rtsa-live-reference.R, which does
not need RTSA to be installed.
NA, 0.63253131244591, 1.40416178814087, 2.01409037736827
to 2e-16 (including the two-pass root search and the suppressed first look);
beta spend vectors (bs_cum, bs_incr) to 1e-9.boundaries(type = "analysis"), given a timing
that stops short of design_R, appends design_R to the information scale, so
the beta timeline has one more look than the alpha bounds it is handed (RTSA
silently recycles them; the harmless "zb * info$sd_proc" warning in the
reference run). The C++ beta_boundary() now accepts exactly that shape
(design-R mode, alpha bound one element short, padded with NaN at the appended
look, never read); every other length mismatch is still an error.
.rtsa_analysis_bounds() accepts the unextended timing, and its interior
"reaches the efficacy wall" check now counts looks from the beta timeline
rather than from the input timing.
RTSA()'s own internal call (timing already extended with design_R) is what
tsa_hr(boundary_route = "analysis") uses; it is covered by the design-R
parity tests and by the Python port, but is not frozen live.rm_bs > 1, or the extended-timing analysis call.
tools/cpp_vs_py.py section 10 repeats the comparison without R.inst/REVERSE_ENGINEERING_RTSA.md: "Validation status" rewritten around
this.Hardening release. The alpha/beta mathematics and every pinned reference value
are unchanged (the compiled core still agrees with the Python port to ~1e-14 and
all schedules in tools/cpp_vs_py.py give the same roots).
qout < as where an unconverged search happened to stop, which is not a
proof: a reachable target the search merely failed to find would also end
below the target. As the boundary goes to +Inf every Phi(...) in prob()
tends to 1, so the largest spend any boundary can produce at a look is
qmax = sum(last), the probability mass still alive. searchfunc() now
compares the target with qmax BEFORE iterating: target > qmax ->
SearchUnreachable (qmax < target - tol, tol = 1e-15; immediate instead of
400 rounds);
otherwise a search that fails to converge beyond the loose tolerance is a hard
error, never "beyond the wall". (An instrumented run over 303 schedules found
the old inference never misclassified -- 171 unreachable declarations, none
with qmax >= target -- so this makes the claim provable rather than
fixing an observed fault.)unreachable_look onward the raw za is the
deliberately non-physical zb + 1, an encoding of "beyond the efficacy wall"
for the sign of the root-search gap -- never a boundary. Stated in
BetaOut, .rtsa_beta_cpp() and on .rtsa_analysis_bounds()'s raw za
(only beta_ubound is meaningful), and asserted by a test..rtsa_alpha_cpp(), .rtsa_design_bounds(),
.rtsa_analysis_bounds() and .rtsa_beta_cpp() now reject Inf, -Inf
and NaN (!is.finite()), not only NA; anyNA() let Inf reach the
compiled recursion, which then returns garbage.RTSA()).tools/cpp_vs_py.py
section 9 (fixed seed, 3-45 looks) now reproduces "300 of 300 calibrate";
section 8 checks the reachability classification (converged / unreachable /
hard error). Both are standalone (no R) -- they are not part of the R test
suite, and the "186 of 300 under the 0.2.7.13/14 behaviour" figure in the
0.2.7.15 entry was measured on that older build, which is not shipped.packageVersion("RTSA")).searchfunc() throw when a search could not converge.
During the information-scale root search that routinely happens at the upper
edge of a bracketing window: the futility spend asked for at a look exceeds
the probability mass still alive, i.e. the futility bound would lie BEYOND the
efficacy wall. The candidate aborted, the bracket window slid on, no sign
change was ever found ("no root bracket"), and tsa_hr() fell back to the
legacy engine (or errored with legacy_fallback = FALSE). Measured with the
compiled core through the design-route orchestration: 114 of 300 random
3-45-look schedules failed, and so did a real 37-look (40-study, target HR
0.94) schedule and 100 evenly spaced looks; 0.2.7.12's silent cap-and-accept
succeeded on them only because the returned junk value happened to have the
right sign, and returns junk elsewhere.
Now an unreachable beta target is reported as such
(BetaOut::unreachable_look, exception SearchUnreachable): from that look on
the futility bounds are set to zb + 1 ("beyond the wall"), so the final gap
the root search works on is NEGATIVE -- the right sign -- and the search
continues. All 300 random schedules and the schedules above now calibrate, with
the same roots as before where 0.2.7.12 succeeded (e.g. 37-look schedule: root
1.2278, 14 suppressed early looks, matching the printed 0.2.7.12 results to
every digit). (The random-schedule figures come from a standalone harness, not
from the R tests; tools/cpp_vs_py.py section 9 reproduces the 300-schedule
calibration for the current build -- added in 0.2.7.16.)
Safeguards: a CONVERGED calibration pass must never be unreachable and its
residual final gap (futility vs efficacy at t = 1) is re-verified
(.rtsa_check_converged_pass()); in the analysis route an unreachable FINAL
look is handled by RTSA's own final clamp (beta_ubound[last] > alpha_ubound[last]) and flagged (final_beyond_wall), while an interior one is
an error. Genuine non-convergence (target reachable but not found) still throws.
Root-search failures now report the last real engine error instead of only
"no sign change".type = "design" timings that
add < 1% of the required information and drops such looks in RTSA(); tsahr
keeps every study. Measuring the default grid (r = 18 vs r = 72) shows the
recursion is accurate for even schedules down to increments of ~0.5% (error
<= 6e-4), for a single close pair down to ~0.25%, and numerically unreliable
below that (500 even looks: final wall 11.7 vs 2.216; one extra look 0.0005
after another: error 6e-3 -- measurements on even schedules and one close-pair
family, not a systematic benchmark of the failure region). tsa_hr() now WARNS (does not alter the schedule; the 0.25%
level is an empirical threshold for this implementation's grid, not an RTSA
rule) when the
increment between looks is < 0.25% of the required information
(.rtsa_check_look_spacing()). Merging near-simultaneous studies, as RTSA()
does below 1%, is left as a decision for the user.test-robust-schedules.R: real 37-look schedule (values checked
against a real 0.2.7.12 printout), 50/100 even looks, dense-at-start,
dense-at-end, tiny first fraction, many early looks (pinned root / rm_bs /
final wall); unreachable candidates (final and interior) give a negative gap
and are flagged; converged-pass checks; root-search error message; the
analysis-route final clamp vs interior error (mocked); the spacing diagnostic.
tools/cpp_vs_py.py (section 7) reproduces the design-route calibration on
these schedules without R.inst/REVERSE_ENGINEERING_RTSA.md is split into
"Current implementation (0.2.7.11-0.2.7.15)" and a clearly marked
"HISTORICAL" part; the three statements that had become false (final boundary
set directly to qnorm(1 - alpha/2), "tsahr has no compiled-code
dependency", "no compiled code was needed") are marked superseded. Test files
that exercise the legacy R-only engine (.obf_*, .rtsa_beta_boundary*) now say
so at the top and where they assert qnorm(1 - alpha/2).packageVersion("RTSA"));
schedules with looks closer than ~0.25% of the required information are only
warned about, not repaired.Implements four items from the 0.2.7.13 audit. The alpha/beta mathematics is
untouched: the compiled engine still agrees with the Python port
(tools/cpp_vs_py.py) to ~1e-14 and reproduces the same reference vectors.
final_non_efficacy fixed (was a semantic bug). 0.2.7.13 defined it as
!crossed_tsa, but crossed_tsa means "crossed at ANY formal look", so a
trial that crossed efficacy at an interim look and then fell back below the
boundary at the definitive look reported final_non_efficacy = FALSE.
New definitive-look fields refer to results$final_tsa_look only:
final_crossed_efficacy, final_non_efficacy (= !final_crossed_efficacy)
and final_entered_futility_region. All three are NA when the route
endpoint has not been reached (there is no definitive look then).
crossed_tsa and entered_futility_region are unchanged and remain the
"at any formal look" versions. Logic lives in the tested helper
.tsahr_definitive_look(); summary_table gained the matching rows and
relabelled the "at any formal look" ones; ?tsa_hr documents both families.boundary_route = "analysis" the formal endpoint is design_R * DARIS.
DARIS itself and the route endpoint are now kept apart:
results$daris_reached and information_size$DARIS_info_threshold_events
always refer to DARIS (t = 1); results$final_reached,
information_size$route_endpoint_info and route_endpoint_events refer to
the route endpoint (identical to DARIS for "design", so the default route
is unchanged). The verbose output, print(), summary_table and plot()
call it "analysis-route endpoint (x.xxx x DARIS)" and draw it as its own
marker; the not-yet-reached fallback endpoint is the event-equivalent of
design_R * DARIS, not of DARIS. The interpolation moved to
.tsahr_events_at_fraction() (identical arithmetic for t = 1).settings$fallback_used, fallback_route, fallback_reason,
route_used. fallback_route is "none", "design" (analysis route
failed; the RTSA-derived design-route result is returned) or "legacy"
(RTSA-derived engine failed; approximate legacy engine). route_used is
"design", "analysis" or "legacy"; boundary_route still records what
was requested. print()/summary() also flag the design-route fallback.
The fallback paths are now covered by mocked-failure tests.z_n_w(), a REVERSED interval (lower wall above the upper wall, za > zb)
is now counted in its own Diagnostics::grid_reversed -- no longer folded
into grid_collapses -- and surfaced by its own, louder warning ("NOT a
valid RTSA computation ... must not be trusted or reported as
RTSA-equivalent"). It deliberately does NOT throw: z_n_w() also runs at
every candidate information scale tried by the root searches, where an
interior look can be transiently reversed at a bad candidate without the
search (which targets the final look) being unable to continue and
converge; erroring there would turn recoverable searches into failures.
Also counted now: an exactly zero-width two-node grid (za == zb), as a
degenerate grid_collapses (arithmetic unchanged).
Because of the same transient-candidate argument, the root searches'
candidate evaluations no longer emit diagnostics warnings (warn = FALSE
in .rtsa_beta_cpp()); the converged passes always do, so a persistent
problem at the accepted root is still reported.DESCRIPTION/README.md now say "no fixed software-imposed limit
on the number of looks, subject to available computational resources"; the
README bullet describing the engine now describes the compiled RTSA-derived
engine (the R-only engine is the opt-in fallback).inst/extdata (needs one run of the real RTSA 0.2.2 calls -- design and
analysis -- with packageVersion("RTSA")), and the information-scale
convention (unchanged, as argued in 0.2.7.12/13).This release implements the five actioned points of the 0.2.7.12 audit, plus
an opt-in boundary_route argument the audit recommended over changing the
default. Two audit points were declined, with reasons (see below).
searchfunc() now errors on genuine non-convergence
(src/rtsa_core.h). Previously, hitting the (disclosed, RTSA-diverging)
400-round iteration cap silently returned whatever value the search had
reached. It now checks the residual against a loose tolerance
(kLooseSearchTol = 1e-6 -- far coarser than the 1e-9/1e-15 convergence
tolerances the recursion itself runs at, so it does not reject harmless
floating-point noise on an otherwise-converged search): within it, the
value is accepted and counted as having taken the "slow path"
(Diagnostics::slow_searches); beyond it, searchfunc() throws
std::runtime_error, which surfaces as a real R error. RTSA's own
while(cond) loop has no cap and would simply never terminate in this
situation, so erroring here is stricter than RTSA, not a departure from
it -- unlike returning a value RTSA's own algorithm would never have
produced.z_n_w()). The
degenerate-interval widening (a disclosed departure from RTSA, which
errors here) is now counted (Diagnostics::grid_collapses) rather than
passed through silently.
Both diagnostics are threaded through alpha_boundary()/beta_boundary()
into AlphaOut/BetaOut, returned from the compiled entry points
(rtsa_alpha_boundary_cpp()/rtsa_beta_boundary_cpp()), and turned into
immediate R warnings by every call site in R/rtsa_engine.R
(.rtsa_warn_diagnostics()), each naming which boundary/pass triggered
it (e.g. "the analysis-route beta (futility) boundary (pass 2)").results$final_non_efficacy. Because the futility and efficacy
boundaries are calibrated to meet exactly at the definitive final look,
entered_futility_region restricted to that look is, in practice, the
complement of crossed_tsa. final_non_efficacy <- !crossed_tsa names
that explicitly, and the "entered_futility_region is not a stopping rule"
documentation under ?tsa_hr is strengthened accordingly. Added to
results and to summary_table.DESCRIPTION,
README.md) and "no artificial limit on the number of looks" (a comment
in R/tsa_hr.R) are now uniformly "no fixed software limit on the number
of looks" -- the previously-stated claim was stronger than the package
actually guarantees.legacy_fallback argument, default TRUE. Governs what happens when
the compiled RTSA-derived engine fails to produce a result. TRUE
(unchanged default) falls back to the legacy, pre-0.2.7.11 R-only
approximate engine, with an immediate warning and a console/print/summary
banner -- a caller wrapping the call in suppressWarnings() will still
see the banner in print()/summary() output and the
beta_engine$engine == "legacy_r_fallback" flag, but not the warning
itself. legacy_fallback = FALSE instead makes tsa_hr() stop() with
an error, appropriate when a silently-substituted, non-RTSA-comparable
result would be worse than a hard failure. Applied both to the (always-run)
design pass and, separately, to the analysis pass when
boundary_route = "analysis" (whose fallback target is the already-
computed design-route result, not the legacy R engine, since the legacy
engine has no analysis-route equivalent).boundary_route = c("design", "analysis") argument, default
"design" (unchanged default, per the audit's recommendation not to
switch it). "analysis" runs RTSA's real
RTSA(type = "analysis", design = NULL) route
(.rtsa_analysis_bounds()/.rtsa_retrospective(), already implemented
in 0.2.7.11 but not previously wired into tsa_hr()): a design pass
first solves an inflation factor design_R, and the alpha/beta
boundaries actually reported are recomputed on the timeline scaled by
it. This is more than swapping the futility numbers -- the formal
endpoint moves from DARIS to design_R * DARIS, so final_reached,
the interpolated DARIS-threshold event count, boundary_timeline's
synthetic endpoint, and every decision-layer field in results are now
computed relative to a new route_endpoint variable (1 for
"design", design_R for "analysis") rather than the previously
hard-coded 1. settings$boundary_route/settings$route_endpoint in
the returned object record which was used. Both routes share the same
compiled recursion; only the orchestration differs, matching RTSA's own
two code paths.FALSE. The audit's
suppressWarnings() scenario is real, but the visible print()/
summary() banner and the beta_engine$engine flag remain even then.
Per direction, the compromise is the new legacy_fallback argument
above (default TRUE, i.e. the existing behaviour kept) rather than a
default change.sum(1/(SE^2 + tau^2))).
Not changed: DARIS already inflates the required information by the
diversity factor 1/(1-D^2), so accruing information on the same
(heterogeneity-inflated) scale it is compared against is the standard
TSA convention, not an inconsistency; option C would likely double-count
heterogeneity. The limitation remains documented, not fixed, per this
release's direction.beta_final is now the final efficacy
bound (RTSA's design pass makes the futility bound meet the efficacy bound at
t = 1). Versions 0.2.6.x-0.2.7.11 used min(qnorm(1 - alpha/2), final efficacy) (1.96 at alpha = 0.05). Consequence: at the definitive look,
results$entered_futility_region is TRUE exactly when the cumulative Z-curve
did not reach the final efficacy boundary. Pre-DARIS bounds are unchanged.
Docs (?tsa_hr) and test-daris-futility-stop.R updated.tsa_hr() still falls back to the pre-0.2.7.11 R-only engine,
but now emits an immediate warning and (if verbose) a console banner;
print()/summary() repeat it; the result is flagged
(res$beta_engine$engine == "legacy_r_fallback", engine_error).tsa_hr() runs RTSA's
type = "design" route on the observed information fractions (DARIS =
t = 1). RTSA's retrospective chain (design pass -> design_R -> analysis pass)
exists internally (.rtsa_retrospective()) but is not used by tsa_hr();
the analysis route also changes the efficacy bounds (alpha recomputed on
t / design_R)..Rbuildignore now excludes all of tools/; internal helper
.rtsa_na_to_null() renamed .rtsa_null_to_na() (it converts NULL to NA);
importFrom(Rcpp, sourceCpp) is deliberately kept (it is what keeps
R CMD check from flagging Rcpp in Imports as unused).The package now needs a C++ compiler (NeedsCompilation: yes, Rcpp).
Numerical comparison against a Python port of RTSA 0.2.2's R sources (which
reproduces RTSA's published futility-vignette output: SMA timing
0.541/0.812/1.083, futility 0.332/1.292/2.014) traced the reported "beta
bounds far from RTSA" symptom to orchestration, not to the C++ kernels
(init_int/recur_int/prob were already term-for-term equal; first,
other, fcab, qpos, trap are never called from any RTSA R file and are
not needed):
root) so the futility bound meets the efficacy bound at t = 1, where
that efficacy bound is the value of RTSA's alpha recursion at t = 1
(2.127 for the 9-look schedule 0.461 ... 0.961, 1). tsahr <= 0.2.7.10
used qnorm(1 - alpha/2) = 1.96 instead: root 1.1517 instead of
1.2107, first-look bound 0.479 instead of 0.531, and so on.type = "analysis" route
was run against alpha bounds computed for the design timeline; RTSA
recomputes them on t / design_R for that route. More importantly, the
RTSA futility numbers used as the reference for these fractions are the
type = "design" route's (reproduced to 4 decimals), not the analysis
route's (which gives ~0.13 ... 1.46 for the same fractions).src/rtsa_core.h, src/rtsa_engine.cpp, src/RcppExports.cpp,
R/RcppExports.R: end-to-end C++ port of RTSA's alpha_boundary(),
beta_boundary(), z_n_w(), searchfunc(), esOF(), sd_inf() and
the three used first.cpp functions. Calls R::dnorm/pnorm/qnorm.
Disclosed differences from RTSA: searchfunc() has an iteration cap
(RTSA converges in <= 9 rounds), and a collapsed integration grid is
widened instead of erroring.R/rtsa_engine.R: .rtsa_design_bounds() (RTSA
boundaries(side = 2, futility = "non-binding", type = "design")),
.rtsa_analysis_bounds() (type = "analysis") and .rtsa_retrospective()
(RTSA(type = "analysis", design = NULL) chain), all using stats::uniroot
as RTSA does.tsa_hr() now takes alpha bounds and pre-DARIS futility bounds from
.rtsa_design_bounds() on the (t < 1, 1) timeline. If that engine errors
it falls back, with a warning, to the pre-0.2.7.11 R-only engine
(.tsahr_legacy_boundaries()).res$beta_engine now carries engine = "rtsa_design_cpp", root,
warp_root, rm_bs, beta_ubound, boundary.test-rtsa-engine-parity.R (reference values from a Python port of
RTSA's R sources, tools/rtsa_port.py; optional live comparisons against
RTSA when it is installed); test-daris-futility-stop.R and the 0.2.7.7
beta_engine test updated.RTSA added to Suggests (used only by the skip_if_not_installed("RTSA")
parity tests).0.2.7.9 changed .rtsa_beta_boundary_analysis()'s design_R/delta/
rm_bs calibration from a side = 2, futility = "non-binding" basis
to a side = 1, futility = "none", right_power()-based one, on the
strength of an external "live RTSA 0.2.2 reconstruction" that reported
different numbers (delta = 2.486475, design_R = 1.057434,
rm_bs = 0) than 0.2.7.7/0.2.7.8 produced
(delta = 2.801585, design_R = 1.151571, rm_bs = 5).
That reconstruction queried the wrong branch of RTSA's source.
Re-reading RTSA's own top-level RTSA() wrapper (R/RTSA.R), not just
boundaries() in isolation, shows the internal call that manufactures
design_R when no design object is supplied uses side = side, futility = futility -- RTSA()'s own top-level arguments (side = 2, futility = "non-binding" for this package's design), not
hardcoded side = 1/futility = "none". The correct calibration is
therefore exactly .rtsa_beta_boundary()'s existing two-pass
warp_root search (side = 2, futility = "non-binding",
type = "design") -- what 0.2.7.7 already used.
design_R, delta, and rm_bs reverted to the 0.2.7.7/0.2.7.8
calculation, now with the source of that calculation directly
cited against R/RTSA.R's exact lines rather than an unverifiable
external reconstruction. .rtsa_beta_boundary_analysis() again
calls .rtsa_beta_boundary() for design_R, uses
delta = abs(qnorm(alpha/2) + qnorm(beta)) (side = 2, matching this
package's own alpha engine), and runs the 3-pass rm_bs
fixed-point iteration RTSA's own side == 2 analysis branch uses.side = 1/futility = "none"/right_power() functions 0.2.7.9
added (.rtsa_design_R(), .rtsa2_alpha_boundary_design_side1(),
.rtsa2_esOF(), .rtsa2_ma_power_upper(), .rtsa2_right_power())
are removed rather than left in as dead code.NA" symptom, without any separate display-layer change.
RTSA's own boundaries() converts every look pinned at the +/-20
sentinel (suppressed by rm_bs) to NA; with 0.2.7.9's rm_bs
fixed at 0, that sentinel was essentially never hit, so early
looks surfaced as small genuine finite negative numbers instead.
tsa_hr() has no NA-hiding logic beyond what
.rtsa_beta_boundary_analysis() returns, so restoring the correct
rm_bs iteration restores the correct NAs..rtsa2_init_int()/.rtsa2_recur_int()/.rtsa2_prob()) against
RTSA's actual src/first.cpp shows they already match exactly, term
for term. RTSA's own R code only calls three of first.cpp's eight
exported functions (init_int, recur_int, prob) -- the rest
(first, trap, fcab, other, qpos) are unused/orphaned. All
three functions that are used call R::dnorm()/R::pnorm()
internally, which is the identical Rmath C library stats::dnorm()/
stats::pnorm() call from R -- there is no floating-point gap
between the two, compiled or not. The "beta bounds too small"
symptom was fully explained by the design_R/delta/rm_bs
regression above.test-boundaries-analysis-rtsa.R updated to assert the
reverted (correct) behaviour instead of 0.2.7.9's; a new test
confirms sentinel-suppressed looks come back as NA.inst/REVERSE_ENGINEERING_RTSA.md rewritten to document this
regression-and-revert plainly, including the exact R/RTSA.R lines
that settle the side/futility question, so this mistake is not
repeated from an unverifiable external claim again without checking
the top-level wrapper first.A live R CMD check/devtools::test() run against this fix (on real
data, including the package's own bundled example dataset, which
reaches an info_fraction of 1.50 at its last study) surfaced two
further problems in the brand-new .rtsa_design_R(), both now fixed:
uniroot() bracket failure on any real, over-powered dataset
("f() values at end points not of opposite sign", raised from
every non-trivial call to tsa_hr(), including the package's own
example and every affected test). Root cause: .rtsa_design_R() was
extending the raw observed information fractions with a final 1
when needed, but never removed fractions that already exceed 1.
Retrospective meta-analyses routinely overshoot the required
information size at their last look or two -- but RTSA's own source
does not merely cap that overshoot; it drops every such look
outright before the design_R search (R/RTSA.R:
if (max(timing) > 1) { timing <- timing[timing <= 1] }). Passing
fractions greater than 1 into the side = 1 alpha-spending recursion
broke its monotonicity in the info-scale factor badly enough that
the search bracket [0.9, 1.2] no longer contained a sign change.
Fixed by filtering to t[t <= 1] before extending with the final 1,
exactly mirroring RTSA's own line.test-boundaries-analysis-rtsa.R "over-powered single early look"
test), the filtered set above is empty and there is no sub-design_R
data left to calibrate a power root from at all. RTSA's own source
has a sanctioned value for exactly this situation --
if (power_adj == FALSE) design_R <- 1 -- and .rtsa_design_R() now
falls back to it, letting the existing t > design_R "over-power"
routing in .rtsa_beta_boundary_analysis() correctly send every one
of those looks straight to the definitive final boundary (as it did
before this release), rather than raising an error.delta/design_R-derivation method itself, the
single-pass rm_bs = 0 beta calculation, or the alpha engine -- only
to what timing values are allowed into the design_R search.testthat cases (surfaced by the
devtools::test() run that caught this) exercise a path through
.rtsa_design_R(); this patch is expected to resolve them. No
further test files were changed for this patch beyond the fixes
already made for the original three bugs described below..rtsa_beta_boundary_analysis()
beta-futility engine, found by a direct numerical comparison against a
live RTSA 0.2.2 type = "analysis" reconstruction on real data, are
fixed:
delta used the wrong side. RTSA's type = "analysis"
branch is reached from inside boundaries()'s side == 1
code path, where the fixed theoretical drift is
abs(qnorm(alpha/side) + qnorm(beta)) evaluated at side = 1
(i.e. plain alpha) -- 0.2.486475 for alpha = 0.05, beta = 0.20. The previous code hard-coded abs(qnorm(alpha/2) + qnorm(beta)) (the design engine's side = 2 drift), giving
2.801585 instead -- about 13% too high, and the root cause of
every downstream discrepancy below.design_R was calibrated from the wrong branch. When no
design_R is supplied, RTSA's own RTSA() wrapper
manufactures one via a separate
boundaries(timing = <observed, extended to 1>, side = 1, futility = "none", type = "design") call -- a single
power-only uniroot(right_power, ...) search against a dummy
lower bound. The previous code instead called
.rtsa_beta_boundary(), which implements a different branch
(futility = "non-binding", a two-pass warp_root/inner-wedge
search) -- the right engine for prospective design work, but
the wrong one for calibrating design_R for a retrospective
analysis. This alone produced a design_R about 8-9% too high
in the case checked (1.151571 vs. the correct ~1.057434),
which shifted every futility boundary derived from it. Fixed by
adding .rtsa_design_R() (and its supporting port of RTSA's
alpha_boundary(..., side = 1, type = "design") and
right_power()/ma_power(), reusing the same
Simpson's-rule integration primitives already ported for the
futility engine) and calling it instead.rm_bs early-look-suppression loop. The
previous code ran RTSA's side = 2 non-binding-futility
type = "analysis" pattern (three passes of beta_boundary(),
re-deriving rm_bs -- the count of early looks whose beta
spend gets manually zeroed -- from the previous pass's
negative-boundary count each time). RTSA's actual side = 1
type = "analysis" branch (the one this package's side = 1
design reaches) calls beta_boundary() exactly once, with
rm_bs left at its default of 0 -- there is no early-look
suppression mechanism in that branch at all. The previous code
was accordingly producing rm_bs = 5 (five looks with their
beta spend forced to exactly zero) where the correct value is
rm_bs = 0 (a smoothly increasing beta-spend curve from the
first look onward).beta_engine$boundary, beta_engine$design_R,
beta_engine$delta, beta_engine$rm_bs, and
beta_engine$beta_spent/beta_spent_delta (and, downstream,
TSA_futility_upper/TSA_futility_lower on the returned
cumul_df) all change for any retrospective analysis. The alpha
(efficacy) engine (.obf_alpha_boundary()) is unchanged -- it
was independently validated against live RTSA output previously and
is not implicated in this bug.tests/testthat/test-boundaries-analysis-rtsa.R updated: the two
tests that previously asserted the old (buggy) behaviour
(design_R equal to .rtsa_beta_boundary()'s warp_root; rm_bs
re-derived by iteration) now assert the corrected behaviour
(design_R equal to .rtsa_design_R()'s root and measurably
different from the design-mode warp_root; rm_bs fixed at 0
with no exactly-zero early-look run). All other tests in that file,
and .rtsa_beta_boundary() (design mode) itself, are unchanged.delta formula is an exact,
closed-form match to the RTSA 0.2.2 reconstruction that prompted
this fix (2.486475). The recursive integration primitives
reused for the new design_R search
(.rtsa2_z_n_w/.rtsa2_init_int/.rtsa2_recur_int/.rtsa2_searchfunc/.rtsa2_prob)
are the same ones already used by, and validated for, the futility
engine, and were additionally spot-checked in an independent Python
re-implementation against RTSA's own published side = 2
esOF/alpha_boundary() reference values
(RTSA::boundaries(timing = c(0.2,0.4,0.6,0.8,1), alpha = 0.05, side = 2, es_alpha = "esOF") -> 4.877, 3.357, 2.680, 2.290, 2.031), which that re-implementation reproduces exactly, and,
separately, against the actual info_fraction values of this
package's own bundled example data (design_R root of ~1.045 with
power exactly matching the 0.80 target at that root, in the same
Python re-implementation). A live R CMD check/devtools::test()
run against it (see the post-check patch above) did catch two
further bugs in .rtsa_design_R()'s timing handling, now fixed.
Before relying on this for a real analysis, run the direct
comparison snippet in the block comment above
.rtsa_beta_boundary_analysis() in R/obf_boundaries.R against a
live RTSA installation on your own data, and please report anything
else R CMD check/testthat still catches..rtsa_beta_boundary_analysis() test suite
added that same release. The failure was in the test, not the
production code it was testing -- the .rtsa_beta_boundary_analysis()
code itself ran correctly end-to-end (including through a full
tsa_hr() call) in every other test, and all pre-existing tests
continued to pass unchanged.ans$rm_bs != 5L || length(t) == 5L -- i.e.
that rm_bs landing on exactly 5 should only be possible when the
input schedule also happened to have exactly 5 looks. That was never
a valid invariant: rm_bs is a count of early looks suppressed by
the fixed-point iteration, which can legitimately equal 5 for a
schedule of any length (the test's own 6-look example genuinely
suppressed 5 of them). The test conflated "happens to equal the old
hard-coded constant" with "is therefore still hard-coded", which
does not follow.rm_bs stays within [0, length(t_ext)], and two
differently-shaped schedules (tightly-spaced early looks vs.
well-separated looks) are confirmed to produce independently-derived
rm_bs values rather than both landing on a shared magic number.R/obf_boundaries.R or R/tsa_hr.R changed in
this release -- only the test file..rtsa_beta_boundary()) is a faithful port of RTSA's
boundaries(..., type = "design") branch: given only a timing
vector, it root-finds its own information-scale inflation
(warp_root) from scratch on every call. That is RTSA's
prospective-design math (solve everything given only alpha,
beta, and a planned schedule). tsa_hr() uses it retrospectively
-- on whatever information fractions the included studies actually
produced -- which is conceptually RTSA's type = "analysis" use
case, but type = "analysis" has a hard prerequisite (design_R, a
sample-size-inflation root from an earlier, separate design call)
that tsa_hr()'s architecture never had anywhere to get from, so it
ran the design-mode math on observed data instead. This was already
disclosed as an open item in the previous release's development
notes; it is a real, mechanical mismatch with RTSA's own
type = "analysis" code path, not just semantics..rtsa_beta_boundary() (design mode, warp_root) is unchanged
and still used internally, but only as step 1 of a new two-step
pipeline: a new .rtsa_beta_boundary_analysis() first calibrates a
design_R internally -- by running that same design-mode root
search once on the observed timing, exactly as RTSA's own RTSA()
wrapper does when no prior design object is supplied
(R/RTSA.R's type == "analysis"/design = NULL branch: it calls
boundaries(type = "design") once to get design_R <- bounds$root)
-- then runs RTSA's actual type = "analysis" beta_boundary() branch
against that fixed design_R:
warp_root/inf_warp() re-scaling at this stage);design_R) is appended as an
explicit endpoint whenever the observed data has not yet reached
it, exactly as RTSA's beta_boundary() does for its own
design_R branch;inf_frac / design_R (how much of the eventually-
planned total information has actually accrued), not the raw
observed fraction;rm_bs) is iterated to a
fixed point (three passes, matching RTSA's own three
beta_boundary() calls for this branch) rather than hard-coded
or re-derived via another root search -- design_R itself stays
fixed throughout, unlike the design-mode pipeline.tsa_hr()'s beta_engine now comes from
.rtsa_beta_boundary_analysis() instead of .rtsa_beta_boundary();
its $boundary output shape (aligned to the unique observed
information fractions) is unchanged, so no other code in the package
needed to change.beta_boundary()/boundaries() design_R branch
(RTSA 0.2.2 source), not a numerically confirmed match against a
live two-step RTSA::boundaries(type = "design") +
RTSA::boundaries(type = "analysis", design_R = ...) call -- no R
interpreter was available in the environment that wrote this port.
See the block comment on .rtsa_beta_boundary_analysis() in
R/obf_boundaries.R for a runnable comparison snippet.R CMD check on 0.2.7.5 came back with 242 passing tests and
exactly one failure, a substantial jump in confidence over the
three previous 0.2.7.x releases (each of which had shipped with a
real package bug). The failure itself was in the test, not the code
it was testing.ans$warp_root (the root found by
the FULL two-pass .rtsa_beta_boundary() pipeline, which may apply
rm_bs > 0 on its second pass) should nearly equal an isolated
root computed with rm_bs fixed at 0. For the specific design
this test uses, the pipeline's first pass suppressed one early look
(rm_bs = 1), which genuinely changes the root-search objective
function -- so the two roots are, correctly, solutions to two
different equations, and were never guaranteed to coincide. They
differed by about 4.6e-4, just past this test's (too tight,
wrongly-reasoned) 1e-4 tolerance.ans$warp_root should be a root of the pipeline's own
rm_bs-adjusted objective function (.rtsa2_inf_warp(ans$warp_root, ..., rm_bs = ans$rm_bs) should be numerically zero), which is what
the fixed test now checks, rather than requiring it to match a
different, unadjusted equation's root.R/obf_boundaries.R changed in this release --
only the test file.R CMD check on 0.2.7.4 showed PASS 2 of the two-pass root search (any
design where at least one look came back negative on PASS 1, i.e.
most realistic designs) failing with "Non-binding futility
boundaries could not be computed" -- not an edge case, essentially
every test that called tsa_hr() on real-shaped data.rm_bs suppression deliberately zeroes the first rm_bs entries of
its timing vector before calling .rtsa_beta_spend_OF() -- that is
the whole mechanism. .rtsa_beta_spend_OF() had an input-validation
guard, if (any(t <= 0)) stop(...), added well before the rm_bs
mechanism existed, that rejected exactly the t = 0 values rm_bs
produces. Every PASS-2 call with rm_bs > 0 therefore failed with
this function's own error on every evaluation, exhausted all 50
widening attempts, and surfaced one level up as a root-search
non-convergence -- a misleading symptom for what was actually an
unconditional, immediate failure. The underlying spending-function
formula was never the problem: qnorm(1-beta/2)/sqrt(0) is +Inf
in IEEE 754 (no error, no NaN), and pnorm(Inf, lower.tail = FALSE) = 0, so t = 0 already gave exactly the correct "no spend yet"
answer -- only the guard rejecting it was wrong. Fixed by relaxing
the guard to t < 0 (negative fractions are still, correctly,
rejected).timing_beta = 1. At that single
point the incremental beta-spend is exactly beta by construction
of any complete spending function, which hits the d1 == beta
special case (za = 0) regardless of the warp factor -- so the root
search's own objective function is constant and can never bracket a
sign change. This is a genuine mathematical degeneracy, not a
search-tuning problem. .rtsa_beta_boundary() now detects the
single-point case up front and skips the root search entirely,
going straight to the same definitive-look convention used in every
other case.RTSA::boundaries() call, which
remains outstanding -- see inst/REVERSE_ENGINEERING_RTSA.md..rtsa_beta_spend_OF() accepts t = 0
and still rejects t < 0; a full two-pass root search on a
design constructed to produce rm_bs > 0 converges rather than
erroring; a single-look (nn = 1) schedule resolves via the new
special case rather than exhausting the root search.beta_boundary(),
z_n_w(), searchfunc(), init_int(), recur_int(), prob())
faithfully, but silently assumed an implicit warp_root = 1: it used
the SAME information fractions both to look up how much beta had been
spent at each look and as the standard-deviation scale the recursion
runs on. RTSA's own boundaries() never does this. It always finds a
scalar inflation of the information SCALE alone (org_inf_frac = inf_frac * root; the beta-spending fractions themselves are never
warped) via root-finding (uniroot(inf_warp, ...)), chosen so the
resulting futility boundary lands exactly on the fixed efficacy
boundary at the definitive final look -- and then repeats that root
search a second time with any look that came back with a negative
futility boundary on the first pass suppressed (rm_bs = sum(lb$za < 0), its beta-spend fraction zeroed, which routes it through the
"negligible spend" sentinel instead).
rm_bs) root-search pass above. Without that pass,
values that real RTSA would have suppressed by construction were
left visible..rtsa2_inf_warp() and .rtsa2_find_warp_root() (new), and extends
.rtsa2_beta_boundary_core() with warp_root and rm_bs
parameters. .rtsa_beta_boundary() now runs the full two-pass
sequence RTSA's own boundaries() does: find a root with rm_bs = 0; run the recursion once to see which early looks came back
negative; re-find the root with those looks' beta-spend zeroed; run
the recursion once more for the final result. RTSA's own narrow,
incrementally-widening uniroot() bracket search ([start - step, start], start stepped by step on failure, up to 50 attempts) is
ported as-is rather than substituted with a generic wide-bracket
search, since a wide bracket risks landing on a spurious root far
from the physically sensible region RTSA's own narrow search is
specifically designed to stay within.za[last] to be able to reach -- and, for
uniroot() to bracket a sign change at all, briefly cross -- the
fixed efficacy boundary as the candidate warp factor is varied, a
fixed-margin ceiling on that one value would make the root search's
objective function unable to ever change sign. The clamp is kept, as
before, for every look before the final one, where the original
crash was actually triggered; that case is unaffected by this
change. .rtsa2_z_n_w()'s own degenerate-grid fallback (0.2.7.1/
0.2.7.2) remains as a second line of defense regardless..rtsa_beta_spend_OF() -- replacing
qnorm(1 - beta / 2) with qnorm(1 - beta), on the claim that
RTSA's esOF() itself omits the /2 -- was checked directly against
the actual RTSA 0.2.2 source and found to be incorrect.
R/RTSA_helperfunctions.R's esOF() is
as_cum[i] <- 2*(1 - pnorm(qnorm(1-alpha/2)/sqrt(timing[i]))) -- the
/2 is present. This formula was already correct in 0.2.7-0.2.7.3
and is unchanged here; applying the proposed fix would have halved
the effective beta-spend rate and reintroduced a real error,
understating early-look futility boundaries even further than the
bug this release actually fixes. A comment directly above
.rtsa_beta_spend_OF() now documents this check so the claim isn't
reintroduced without a fresh, direct read of the real source.stop("Non-binding futility boundaries could not be computed...") when no bracket is
found within 50 widening attempts) is itself ported from RTSA's own
code -- it is RTSA's
defined behaviour for a sufficiently unusual design, not necessarily
a tsahr defect, if it is ever hit. Please re-run the comparison in
inst/REVERSE_ENGINEERING_RTSA.md against a live RTSA/rpact call
before relying on this for anything beyond an approximate futility
band.warp_root = 1 case is measurably
different from it, i.e. that the fix is not a no-op.Re-running R CMD check after 0.2.7.2 reported 4 failures, all in the
one new test added in 0.2.7.2
("0.2.7.2: a 3-node (m = 3) Simpson grid produces correct weights, not an error") -- not in .rtsa2_z_n_w(), .rtsa2_seq_by2(), or any
other package code, and not a crash: expect_length()/expect_equal()
mismatches. The underlying 0.2.7.2 fix itself was already correct and is
unchanged in this release.
Cause: that test picked lo = 1.5 as a stand-in for a near-collapsed
[za[i], zb[i]] window. For the specific r = 18, delta = 0,
sd_incr = 1 inputs used in the test, 1.5 happens to land EXACTLY on
one of .rtsa2_z_n_w()'s own log-spaced grid nodes (node j = 72
evaluates to precisely -3 + (72 - 18)/12 = 1.5), so one extra original
node coincided with lo and survived trimming, correctly producing a
5-node grid (m = 5) instead of the 2-node/3-node (m = 3) case the
test's hard-coded expected values assumed. The code was doing the right
thing with the input it was given; the test's prediction of what that
input would produce was wrong.
Fix: rewrote the test (tests/testthat/test-boundaries-rtsa.R) to
(a) use lo = 1.234567, deliberately off the r = 18 grid's node
spacing, and, more importantly, (b) stop hard-coding an exact node count
or exact zj/wj values altogether, in favor of checking the
invariants that must hold for ANY valid Simpson grid built on
[lo, hi] regardless of how many interior nodes happen to survive
trimming: no error, an odd node count >= 3, endpoints exactly at
lo/hi, a non-decreasing sequence, and total Simpson weight equal to
the window width (the property Simpson's rule actually guarantees, and
what the downstream recursion actually relies on). This is more robust
against similar grid-alignment coincidences with any future choice of
r/lo, not just a fix for this one value.
No change to R/obf_boundaries.R in this release.
0.2.7.1 fixed the length(xi) <= 1 grid-collapse crash, but re-running
R CMD check surfaced a second, closely related crash in the exact same
function, .rtsa2_z_n_w(), now at seq(3, m - 2, 2): wrong sign in 'by' argument -- hit by the same target_HR = NA example-data scenario as
before (and by two of the new 0.2.7.1 regression tests themselves, which
is how it was caught).
Root cause: the 0.2.7.1 clamp keeps za[i] only 1e-6 below zb[i]
when the two would otherwise converge -- deliberately tiny, so as not to
perturb any well-separated boundary. But a 1e-6-wide interval is far
narrower than the spacing between .rtsa2_z_n_w()'s log-spaced grid
nodes, so trimming to [za[i], zb[i]] in that situation legitimately
leaves exactly 2 points (the two endpoints, no interior node survives) --
which is precisely what the 0.2.7.1 safety net widens degenerate cases
to as well. Both routes land on the same m = length(xi) * 2 - 1 = 3
grid. That case turns out to have been separately broken already, for
any m = 3 grid regardless of cause: the Simpson-weight loop's
k %in% seq(3, m - 2, 2) membership check evaluates seq(3, m - 2, 2)
eagerly, and for m = 3 that's seq(3, 1, 2) -- an empty index range,
but, as with 0.2.7.1's seq(1, length(xi) - 1, 1), R's seq() with an
explicit by throws rather than returning integer(0) when to < from.
So m = 3 (i.e. a 2-node grid) was never actually safe to reach, even
after 0.2.7.1 -- it was just newly reachable by 0.2.7.1's own fix,
where previously (pre-0.2.7.1) the crash at length(xi) <= 1 happened
first and masked it.
Fix: new internal .rtsa2_seq_by2(from, to) helper -- seq(from, to, 2), but returns integer(0) instead of erroring when to < from -- used
in place of the two raw seq(3, m - 2, 2) / seq(2, m - 1, 2) calls in
.rtsa2_z_n_w()'s Simpson-weight loop. For any m where those ranges
are non-empty (m >= 5), this is byte-for-byte the same seq() call as
before -- no change to any boundary value outside the previously-
crashing case. Hand-traced the m = 3 case through the full weight
loop to confirm it now produces the correct composite-Simpson weights
(1/6, 4/6, 1/6 of the interval width, summing to the full width) rather
than just failing to crash.
Also fixed the two 0.2.7.1 regression tests that this same bug broke
(test-boundaries-rtsa.R): both exercised exactly the m = 3 path (the
real tsa_hr(path, verbose = FALSE) call, and the direct
.rtsa2_z_n_w() unit test with equal/reversed za/zb) and so were
themselves failing under R CMD check, for the same underlying reason
they were written to catch the previous bug -- no test changes were
needed once the underlying .rtsa2_seq_by2() fix was in place; they
pass as originally written.
R CMD check on 0.2.7 failed 3 tests with Error in seq.default(1, length(xi) - 1, 1): wrong sign in 'by' argument,
raised from .rtsa2_z_n_w() via .rtsa2_beta_boundary_core() /
.rtsa_beta_boundary(), triggered by the package's own bundled example
data under target_HR = NA (the circular-target scenario exercised by
tsa_hr(path, verbose = FALSE) with no target_HR).
Root cause: .rtsa2_z_n_w() builds a Simpson's-rule integration grid
on [za[i], zb[i]] (the futility and efficacy boundaries at look i) by
trimming a fixed log-spaced node vector down to that interval. The
non-binding futility boundary za[i], returned unconstrained by
.rtsa2_searchfunc()'s root search, is expected to approach the fixed
efficacy wall zb[i] closely near the final look by design (the two are
constructed to meet at t = 1) -- but nothing kept it from landing at or
past zb[i], which collapses the trimmed grid to a single point (or
fewer). The two-argument form seq(1, length(xi) - 1, 1) then throws
"wrong sign in 'by' argument" once length(xi) <= 1, instead of quietly
returning an empty sequence -- this is a real edge case surfaced by this
specific port (RTSA's own equivalent R source has the identical
unguarded pattern), not a mistranslation.
Fix, two parts:
.rtsa2_beta_boundary_core()): new internal
.rtsa2_clamp_za_below_zb() helper clamps za[i] to stay at least a
small margin (1e-6 on the Z scale) below zb[i], applied right after
every branch that sets za[i] (the zninf sentinel and 0 branches
included, for uniformity, though the observed failure was in the
.rtsa2_searchfunc() branch). This keeps the grid-collapse condition
from arising during normal operation in the first place; the margin is
small enough that it does not perturb any well-separated boundary
value..rtsa2_z_n_w()): if the trimmed grid still collapses
to fewer than 2 points (e.g. from a future direct call with unclamped
za/zb), it is now widened back out to the two-point interval
[za[i], zb[i]] (with a minimal positive width inserted if the two are
equal or reversed) instead of proceeding to the vectorized Simpson's-
rule construction that assumed at least 2 points. The specific
seq(1, length(xi) - 1, 1) call is also replaced with
seq_len(length(xi) - 1), which returns an empty sequence rather than
erroring when length(xi) == 1 (note: 1:0 is not such an idiom -- it
evaluates to the length-2 sequence c(1, 0), not an empty one;
seq_len(0) is the correct empty-sequence idiom).No change to any boundary value for the normal (well-separated za/zb)
case tested elsewhere in the suite -- the clamp is a no-op unless za[i]
would otherwise land within 1e-6 of zb[i], and the safety net is only
reached if the clamp somehow didn't fire first.
getInnerWedge(),
betas_Obf, sdfunc, first_old, other_old, searchfunc_old,
qpos_old, fcab_old, gfunc, and tsa_beta_bound as its source.
Having now obtained and read the actual RTSA package (RTSA 0.2.2,
https://github.com/AnneLyng/RTSA) directly, none of those names
exist anywhere in it. 0.2.6.6-0.2.6.9 softened the surrounding
comment from "literal R implementation" to "adapted from" on an
external audit's advice, but the algorithm itself was never
re-examined or rebuilt until now.beta_boundary(), z_n_w(),
and searchfunc() (R/RTSA_helperfunctions.R, ported near-verbatim,
since they are already pure R), and init_int(), recur_int(), and
prob() (src/first.cpp, ported from C++ to pure R, since tsahr has
no compiled-code dependency) -- specialised to es_beta = "esOF",
the only spending family this package exposes. See
inst/REVERSE_ENGINEERING_RTSA.md for the full algorithmic write-up,
including a side-by-side of the old and new algorithms.delta = |qnorm(alpha/side) + qnorm(beta)|,
computed directly from alpha/beta/side before any recursion runs
-- independent of the data or of the futility construction's own
shape.+-za) null-referenced wedge and only
reconciled it with the real (generally asymmetric,
decreasing-then-flattening) efficacy boundary after the fact.
RTSA's real engine uses the actual efficacy boundary as the upper
wall of the recursive integration at every step from the start.b[b <= 0] <- NA_real_). RTSA returns NA only where the
incremental beta spend at a look is negligible enough to hit an
internal sentinel (zninf = -20); a finite, possibly negative,
non-binding futility boundary is a legitimate result (it says the
trial would only be flagged for advisory futility if the
cumulative Z-statistic had already crossed to the "wrong" side of
the null by that point) and is no longer hidden.qnorm(1-alpha/2), rather than porting RTSA's own
root-finding (uniroot(inf_warp,...)) machinery, which solves a
prospective design problem tsahr's retrospective architecture has
no counterpart for; (b) a defensive pmin(boundary, c_vec_alpha)
safeguard is still applied on top of the reconstructed recursion,
since RTSA's own algorithm does not hard-enforce that a futility
value can never exceed the efficacy wall at every intermediate step.RTSA::boundaries() call. Unlike the alpha
engine (0.2.6, checked against a live 5-look call and matching
essentially exactly), none of the new code in this port has actually
been executed. It is a careful, line-by-line reading of RTSA's
published source, not a numerically confirmed match.
inst/REVERSE_ENGINEERING_RTSA.md
gives the exact comparison to run before relying on this for
anything beyond an approximate, illustrative futility band (which is
how this package has always described its futility output)..rtsa_gfunc, .rtsa_trap, .rtsa_fcab,
.rtsa_qpos, .rtsa_first, .rtsa_other, .rtsa_searchfunc_old,
and .rtsa_get_inner_wedge are deleted; a new test
("legacy (non-RTSA-matching) beta engine internals have been
removed") checks they stay gone. .rtsa_beta_spend_OF() is
unchanged -- its formula already matched RTSA's real esOF(); only
the recursion built around it was wrong.delta to the fixed
theoretical formula; a new test confirms the dead functions are
gone. The pre-existing "definitive alpha boundary,"
"over-powered analysis," and "capped at the corresponding alpha
boundary" tests are unchanged and still pass under the reconstructed
engine. No downstream test (DARIS/futility-stop end-to-end tests,
the 40-study integration test) required changes: the only
beta-engine-derived quantity that reaches tsa_hr()'s formal output
for those scenarios is the single pre-DARIS interim boundary and the
independently-computed definitive final value, neither of which
depends on the specific numeric internals changed here.test-tsa_hr.R "non-numeric input columns are diagnosed as a type
problem" failed under R CMD check with
replacement has 0 rows, data has 10. The bug was in the test, not in
the package: it read the bundled example workbook directly with
readxl::read_excel() and then addressed d$Events_Treatment, but
several headers in that sheet are stored with spaces (Events Treatment, N treatment, Events controls, N controls; only
Study, log_HR, and Std_Error are already underscored). Space-to-
underscore normalisation happens inside tsa_hr(), so at that point
in the test d$Events_Treatment was NULL, and
as.character(NULL) is character(0) -- which cannot be assigned
into a 10-row data frame.gsub(" ", "_", names(d)) normalisation
that tsa_hr() performs before addressing columns, so it targets real
columns and stays correct if the example sheet's headers change.tsa_hr(). The other four tests that read the workbook
directly were checked for the same latent flaw and are unaffected:
they either add new columns or use Std_Error, which is genuinely
underscored in the source file.Prompted by an external (ChatGPT) audit of 0.2.6.7, which rated the
release otherwise ready and identified circularity_warning as the one
item to fix before finalising. Each item was checked against the code
before acting on it.
circularity_warning now means what its name and message say. It
was defined as
is.na(target_HR) && sum(total_events) / DARIS_events > 3 -- i.e. it
carried a severity condition on top of the circularity condition.
But summary.tsa_hr() gates a message worded purely for the general
case on that flag ("target_HR was not specified, so the observed
pooled HR was used... This is circular"). The consequence was a real
user-visible gap, not just a misnomer: a circular analysis sitting at,
say, 1.5x DARIS reported no circularity note at all, directly
contradicting ?tsa_hr, which correctly documents that any RIS
derived from the observed pooled effect is circular. Now split into
two fields:
circularity_warning <- is.na(target_HR) -- circular, full stop.circularity_severe -- circular and accrued events exceed three
times DARIS (the runaway case where the boundary collapses to the
conventional one almost immediately).
Both are returned in information_size and documented in ?tsa_hr.
The statistical calculation is unchanged; only the reporting is.circularity_severe, but a shorter note now
fires for any circular analysis, so the console and summary() agree
with each other and with the documentation. summary() prints the
general note whenever circular, and adds the collapse warning on top
when severe.method documentation no longer says "passed straight through"
to metafor::rma(), which stopped being exactly true in 0.2.6.7 when
"CO"/"VC" began being normalised to "HE" first. Now states that
it is passed after validation and, for those aliases, normalisation.is.finite() returns all-FALSE (rather than erroring) on a
character or factor column, so a column read in as text -- Excel cells
stored as strings, a stray footnote marker forcing the column to
character -- surfaced as "Column(s) must contain only finite values
(found NA/NaN/Inf)". True of the is.finite result, but a misleading
diagnosis. log_HR, Std_Error, and the four count columns are now
type-checked first, and the error names each offending column and its
actual class.circularity_severe never holds without circularity_warning), for
summary() reporting circularity in the non-severe case, and for the
numeric-type diagnosis.Not changed, deliberately: the beta engine remains frozen and
RTSA-derived/adapted, as in 0.2.6.7. The audit's suggestion to escalate
order_by NA handling from warning to error was not adopted -- the
current behaviour is transparent and reproducible, and the audit itself
scoped that as a possible future major release, not a fix for this one.
Prompted by an external (ChatGPT) review of 0.2.6.6, which flagged the
"CO" removal as a release blocker. The flag was correct; this release
reverses that removal and picks up several smaller items from the same
review. Each was checked against primary sources or the actual code
before acting on it.
method = "CO" (and "VC") are now accepted, and normalised to
"HE". 0.2.6.6 removed "CO" from the supported methods on the
basis of a direct test against an installed metafor, which threw
Unknown 'method' specified. That observation was real, but the
conclusion drawn from it -- that "CO" is "not a recognised method
string in current metafor at all" -- was wrong, and the
0.2.6.6 source comment and documentation asserting that have been
corrected. metafor's current documentation states that the Hedges
estimator is also called the variance-component or Cochran estimator
and that method = "VC" or method = "CO" may be used to select it.
The two observations are reconciled by version skew: the alias is
present in current metafor but absent from the older release the
0.2.6.6 test ran against (the CRAN 4.6-0 reference manual carries the
same sentence without the alias parenthetical, which is exactly the
boundary in question)."CO" to metafor::rma() would make
tsa_hr() work or fail depending on which metafor the user happens to
have installed. Mapping "CO"/"VC" to "HE" inside tsa_hr()
makes behaviour identical on every metafor version, and avoids
declaring a minimum metafor version in DESCRIPTION purely to pin
down an alias. All three strings denote the same estimator, so there
is no numerical consequence.parameters$method_requested (what the caller passed)
alongside parameters$method (the normalised string actually used).
method_requested is always populated, including when no aliasing
occurred, so downstream code can rely on the field existing."CO" and
"VC" must produce identical tau2, D2, and cumulative Z to
"HE". The 0.2.6.6 test asserting "CO" errors has been removed and
replaced with one asserting it does not."GENQ"/"GENQM" remain unsupported, unchanged and for the unchanged
reason: they require a user-supplied weights argument that
tsa_hr() does not collect. That part of the 0.2.6.6 finding was
correct and is not affected by the "CO" reversal.Std Error and Std_Error produced
two identically-named columns, after which data$Std_Error silently
resolved to whichever came first -- a wrong-column-used bug yielding a
plausible-looking but incorrect analysis with no error anywhere.
tsa_hr() now stops and names the offending columns instead of
guessing which was meant.order_by is normalised the same way as the column names.
Previously, passing a header exactly as it reads in the user's
spreadsheet (order_by = "Publication Year") failed with "not a
column in data" for a column that is plainly there, because matching
happened after normalisation. Both spaced and underscored forms now
work; documented in ?tsa_hr.order_by sorting is now tested, not just its tied-value warning.
The previous test only confirmed the warning fires; order_by could
have silently no-opped and only the warning path was covered. The new
test reverses a strictly-ordered dataset, sorts it back via
order_by, and requires identical Study order, Z, and
info_fraction to the already-sorted data. TSA is order-dependent, so
this is a reproducibility guarantee.pmin(boundary, c_vec_alpha), dynamic grid sizing, a convergence
fallback, defensive NA handling, post-DARIS mapping onto the
definitive t=1 boundary). It is now described as RTSA-derived and
adapted, with the open point stated plainly: no live multi-look RTSA
beta reference has yet been obtained, so numerical identity with RTSA
is not claimed. The engine itself is unchanged and remains frozen for
this release, as intended.D2 comment's claim that D2 is "mathematically bounded" in
[0,1) is softened. D2 is bounded by definition, and
var_random >= var_fixed holds for essentially all of metafor's
estimators, but the computed value is a ratio of two separately
estimated variances -- which is precisely why the max(0, .) and the
0.999 cap exist rather than being redundant.Prompted by a second external (ChatGPT) audit of 0.2.6.5, verified against the actual code and against current metafor documentation before acting on any of it (see below for what was independently corroborated and how).
method = "CO", "GENQ", "GENQM" never workedtsa_hr(method = ...) has advertised 13 metafor method strings as
supported since the method parameter was first added, but 3 of them
never actually worked -- this went unnoticed because only "DL",
"REML", and "ML" were ever exercised by the test suite. Found while
adding a test that runs every advertised method (itself prompted by
the audit's suggestion to test the full advertised set) -- independent
of anything the audit itself flagged.
"CO" is not a recognised method string in current metafor at all
(verified directly against an installed metafor: it throws
"Unknown 'method' specified", not a Hedges-estimator result).
It appears as a documented alias for "HE" in some older/secondary
sources, but is not accepted by rma() as currently shipped.
Independently corroborated: a search of metafor's own current
documentation and training materials turned up no mention of "CO"
in any enumerated method list."GENQ" / "GENQM" require the caller to also supply a weights
argument to metafor::rma() (per metafor's own documentation, e.g.
rma(yi, vi, weights = 1/vi, method = "GENQ")), which tsa_hr()
does not currently collect or pass through, so calling it with
these methods and no weights errors out inside rma().
Independently corroborated against metafor's own documentation and
Cochrane training materials.valid_methods corrected to the 10 that actually work standalone:
DL, HE, HS, HSk, SJ, ML, REML, EB, PM, PMM. Passing "GENQ" or
"GENQM" now raises a specific, explanatory error (these ARE real
metafor method names, just not usable here without a feature this
package doesn't yet have, so the generic "method must be one of"
list would be misleading); passing "CO" or any other unsupported
string still gets the generic list..tsahr_method_label)
to drop the now-invalid "CO" entry.@param method in R/tsa_hr.R and man/tsa_hr.Rd to list the
correct 10 methods and explain why the other 3 aren't supported."every advertised method value runs and returns a valid object") exercises all 10 currently-supported methods on the bundled
example data, plus explicit tests that "CO"/"GENQ"/"GENQM" each
fail with the right error message -- so a similar gap can't recur
silently.method = "CO", "GENQ", or "GENQM": those calls were already
failing before this release (with a less legible error from inside
metafor::rma()), so no one could have been relying on them
succeeding -- this release only makes the failure explicit and
immediate, with a clearer message.R/tsa_hr.R and man/tsa_hr.Rd incorrectly stated that the final
futility boundary "equals the same value as the final efficacy
boundary." It doesn't, and the actual code never implied it should:
beta_final <- min(qnorm(1 - alpha_two_sided / 2), tail(alpha_bounds_design, 1))
qnorm(1-alpha/2) value, min()
almost always selects the conventional value instead -- confirmed
numerically on the bundled example data (2.0786 for the final
efficacy boundary vs. 1.959964 for the final futility boundary).
Practical consequence: for 1.96 < |Z| < 2.08 (using that example),
conventional p<.05 is YES, the TSA efficacy boundary is NO, and the
futility region is also NO -- a legitimate intermediate state the old
wording obscured. Wording corrected in both files to describe the
final futility boundary as "the conventional two-sided alpha critical
value at the definitive analysis," not as equal to the efficacy
boundary. New regression test
("final futility boundary is the conventional critical value, not the final efficacy boundary", test-daris-futility-stop.R) pins
this numerically so it can't silently drift back.D2_raw and D2_was_capped added to the returned heterogeneity
list, alongside the existing (possibly-capped) D2 -- so a downstream
user can tell "D2 genuinely computed as 0.999" apart from "D2 was
capped here from something larger/degenerate" instead of the capping
being silent. Documented in @return/\value{}.Study identifiers are now validated as non-missing and
non-empty (previously only checked for duplicates).order_by now warns on tied values (previously only warned on
NAs and non-numeric/non-Date types) -- TSA is order-dependent, so
the relative order of tied studies can matter.target_HR very close to 1 (0.90-1.10, excluding exactly 1, which
is still a hard error) now triggers a warning, not silence --
required information size grows rapidly as log(target_HR) -> 0
(e.g. target_HR=0.95 needs roughly 4x the information of
target_HR=0.90 for an otherwise identical design), so this is a
nudge to confirm the target is intentional, not a defect.R/tsa_hr.R and man/tsa_hr.Rd
to the more precise framing: the canonical Lan-DeMets/O'Brien-Fleming
theory assumes a fixed, canonical information process with independent
Brownian-motion increments, and tsahr's random-effects (tau^2
re-estimated at every look) Z-process doesn't exactly satisfy that --
rather than the previous, slightly imprecise "holds exactly only for
a fixed-effect cumulative process."inst/REVERSE_ENGINEERING_RTSA.md), not a generic approximation.method section reworded to the more precise framing:
method does not change the mathematical alpha-spending function or
boundary-calculation algorithm, but it DOES change tau^2, the
cumulative information schedule, and therefore which study lands on
which information fraction -- so it can still change the practical
timing of a boundary crossing indirectly, which the previous wording
understated.Prompted by an external (ChatGPT) review; each item below was checked against the actual code before acting on it, not applied at face value -- see the per-item notes.
allocation_p validation hardened. allocation_source = "manual"
with allocation_p = NA, a length>1 vector, or a non-numeric value
previously could reach a generic/uninformative R error (e.g. "missing
value where TRUE/FALSE needed" for NA) instead of the package's own
informative message. Confirmed by direct testing before fixing, not
just taken on faith. Now validates is.numeric(), length() == 1L,
and is.finite() before the range check, with three new regression
tests (NA, c(0.5, 0.6), "0.5") added to test-tsa_hr.R.man/tsa_hr.Rd, out of sync with R/tsa_hr.R in three
compounding ways. The source's roxygen documentation had already
been updated in previous releases, but man/tsa_hr.Rd was never
regenerated to match, so it still said info_fraction "is capped at
1 for every subsequent look" (false -- only the boundary/decision
mapping is capped via pmin(info_fracs, 1); the stored
cumulative$info_fraction itself is not, and can exceed 1, e.g.
1.0126, after DARIS) and "...computed at the exact accrued
statistical information" (the source had already been corrected to
"observed"). The Rd was also entirely missing the source's current
"Retrospective boundary timeline" paragraph (about boundary_timeline
and the synthetic t=1 point), still showing an older, superseded
"Looks after DARIS is reached" paragraph instead. man/tsa_hr.Rd is
resynced to the current source. If you use roxygen2/devtools to
regenerate docs going forward, this class of drift won't recur; for
now the fix was done by hand.boundary_timeline and the caveated components of results are
now documented in @return/\value{} (both the roxygen source and
the Rd), including that results$entered_futility_region == TRUE at
the DARIS-reaching look reflects a comparison against the definitive
t=1 futility boundary, not an interim one, and is not itself a
formal stopping recommendation -- that caveat previously existed only
in the printed/verbose console output, not in the object a caller
gets back programmatically.crossed_conventional's full-cumulative-curve scope (as opposed to
crossed_tsa's decision-horizon-restricted scope) is now an explicit,
documented design choice, not left ambiguous. Checked directly:
yes, crossed_conventional uses the unrestricted cumul_df$Z
(including studies added after DARIS) while crossed_tsa is
restricted to decision_idx -- these do answer genuinely different
questions. Concluded this is intentional (the printed label already
says "Cumulative Z-curve crossed...", and contrasting it with the
properly-scoped crossed_tsa line illustrates exactly the
repeated-testing inflation risk TSA exists to guard against), so
behaviour is unchanged -- but this was a judgment call, not a
certainty, so it's now spelled out in both a code comment and
?tsa_hr rather than left for a future reader to guess at. If you
intended the decision-horizon-restricted reading instead, that's a
one-line change (see the code comment above crossed_conventional
for the exact replacement) -- flag it and it'll be made explicitly
rather than silently.README.md now states explicitly that method (the
heterogeneity-variance estimator) affects tau^2, the random-effects
cumulative Z-curve, and D-squared/DARIS-related quantities, but does
not change the alpha-spending function or TSA monitoring
boundaries themselves.rtsa_reported <- c(4.877, 3.357, 2.680, 2.290, 2.031) alpha test).
This needs an actual reference vector from a live
RTSA::boundaries() call with a beta-spending function, the same way
the alpha reference came from your own live RTSA session in 0.2.5.2.
Send the output of a call like
RTSA::boundaries(timing=c(...), alpha=0.05, beta=0.20, side=2, es_alpha="esOF", es_beta="esOF") and this test gets added properly.plot.tsa_hr()'s subtitle and methods caption now render "Diversity
D²" using the Unicode superscript-two character (U+00B2), matching
the existing convention elsewhere in the package (e.g. the tau²
console note added in 0.2.5.1), instead of adding a new dependency
(e.g. ggtext) just for this. Purely cosmetic -- no change to any
computed value, and x$heterogeneity$D2 (the underlying R variable
and list element name) is unchanged, since it must stay valid R
syntax. Not changed in this release: the "Diversity D2" wording in
summary.tsa_hr()'s printed table and in the near-boundary warning
message from tsa_hr() -- only the plot was in scope for this fix..beta_spend_OF() (R/obf_boundaries.R), the un-doubled
candidate beta-spending formula beta*(t) = 1 - Phi(z_beta/sqrt(t)).
It was flagged as dead code when the RTSA-matched beta/futility engine
was added in 0.2.4: the actual boundary-solving path has used
.rtsa_beta_spend_OF() (beta*(t) = 2*(1 - Phi(z_{beta/2}/sqrt(t))),
reverse-engineered directly from RTSA's getInnerWedge() -- see
inst/REVERSE_ENGINEERING_RTSA.md) ever since, and .beta_spend_OF()
was never called from anywhere in that path. Not a behaviour change --
tsa_hr()'s output is identical before and after this release..beta_spend_OF()'s (dead, un-doubled) formula as
though it were what the package actually computes. It now describes
.rtsa_beta_spend_OF()'s real (doubled) formula instead, with a
pointer to inst/REVERSE_ENGINEERING_RTSA.md.test-tsa_hr.R), which
exercised only the deleted dead function. The equivalent property for
the function actually in use is already covered by "RTSA beta
spending formula is used" in test-boundaries-rtsa.R.(Supersedes an interim 0.2.6.1 that introduced this fix but also introduced a second, separate documentation problem -- see below.)
VALIDATION comment block in R/obf_boundaries.R said two
incompatible things about the same question. One passage stated
the full 5-look boundary schedule had been checked end-to-end
against live RTSA::boundaries() output and matched "essentially
exactly"; a later passage in the same block said only the first look
had ever been checked (in closed form) and explicitly warned readers
that the rest was unverified. These cannot both be true, and
inst/REVERSE_ENGINEERING_RTSA.md disagreed with the second passage
too (it also claimed the full comparison was done) -- so the two
files disagreed with each other on top of the R file disagreeing
with itself.timing = c(0.2, 0.4, 0.6, 0.8, 1.0); it reports back a rounded
SMA_Timing column (0.205, 0.409, 0.614, 0.818, 1.023) reflecting its
internal event-count discretisation, and at one point during this
package's development that reported (rounded) column got used as the
comparison input instead of the nominal schedule that was actually
requested. Since the first-look boundary is an exact, closed-form
function of the timing value fed in, that mismatch alone produces a
spurious ~0.06 "discrepancy" that has nothing to do with the formula
or the recursive engine. Using the correct nominal timing, an
independent from-scratch Python port of .obf_alpha_boundary()
reproduces all 5 of RTSA's reported boundaries to within
~0.003-0.006 -- even at this engine's long-standing default of
n_grid=2000 -- consistent with "essentially exactly"..obf_alpha_boundary()'s default n_grid is still raised from
2000 to 16000 in this release, but the justification is narrower
than previously stated: not "fixes a confirmed ~0.03 error", but "a
conservative, essentially-free accuracy margin" -- FFT-based
convolution makes the extra grid resolution cheap for this package's
actual usage pattern (a handful of interim looks per tsa_hr() call),
so there's no real cost to preferring the larger grid even without a
confirmed problem at 2000. bmax is unchanged. Anyone who can run
the real R implementation is encouraged to confirm (or correct) the
Python-based numbers above -- see the VALIDATION note in
R/obf_boundaries.R for the exact snippet to run.test-alpha-spend-rtsa.R) pins the full
recursive-engine output against all 5 RTSA reference boundaries
(using the correct nominal timing) at a 0.01 tolerance. Based on the
Python check above this is expected to pass at both n_grid=2000 and
n_grid=16000, not only at the new default -- unlike what an earlier
draft of this test's own comment claimed..alpha_spend_OF() itself, the beta/futility engine, or
anything else; this release is scoped entirely to documentation
accuracy and the default grid resolution..alpha_spend_OF()'s two-sided alpha-spending formula was wrong for
the methodology this package documents itself as implementing, and
is corrected in this release. tsa_hr()'s O'Brien-Fleming-type
efficacy boundaries were computed with
alpha*(t) = 2*(1-Phi(z_{alpha/2}/sqrt(t))) in every prior tsahr
release (0.2.0-0.2.5.1). That form is a textbook Lan-DeMets
two-sided spending function — it's internally valid as a spending
function (it does reach exactly alpha at t=1) and is quoted as-is
in some general group-sequential-design references (e.g. gsDesign's
documentation) — but it does not match the specific TSA
methodology (Copenhagen Trial Unit / RTSA, Thorlund et al.) that this
package documents itself as following (Miladinovic et al. 2013,
Wetterslev et al. 2009). The corrected formula is
alpha*(t) = 4*(1-Phi(z_{alpha/4}/sqrt(t))) -- RTSA's own
side-parameterised spending function evaluated at side=2.RTSA::boundaries() call
(timing=c(0.2,0.4,0.6,0.8,1), alpha=0.05, side=2, es_alpha="esOF"),
which returns boundaries 4.877, 3.357, 2.680, 2.290, 2.031 -- matching
the corrected formula (run through tsahr's existing recursive
integration engine) essentially exactly, across all 5 boundaries of
that schedule, and clearly distinguishable from the pre-correction
formula's 4.383, 3.099, 2.554, 2.254, 2.063 for the same schedule. A
first-look regression test (test-alpha-spend-rtsa.R) locks in the
corrected value against this RTSA reference in closed form.tsa_hr() from
a version between 0.2.0 and 0.2.5.1 for TSA boundaries at an interim
look (not just a completed, fully-informed meta-analysis) should
re-run their analysis with this version.R/obf_boundaries.R confirm -- they check the total spend, not the
shape across interim looks. This bug is a genuinely separate
property from what those checks were ever positioned to catch, not a
contradiction of them.VALIDATION note at the top of R/obf_boundaries.R
for the full numeric account, and inst/REVERSE_ENGINEERING_RTSA.md
for the alpha-engine analogue of the beta-engine writeup added in
0.2.4..alpha_spend_OF() and its documentation changed. The recursive
integration engine that solves for the actual sequential boundaries
given a spending function (the FFT-based convolution machinery) is
untouched, as is the beta/futility engine (which is non-binding,
advisory-only, and doesn't affect type I error control), DARIS
calculation, and everything else in the package. No existing test
asserted a specific numeric alpha-boundary value tied to the old
formula (the one test that checks a specific constant, "internal OF
alpha boundary is consistent with the classical O'Brien-Fleming
boundary", uses a tolerance of 0.05 around the published ~2.040
reference and still passes under the corrected formula's ~2.031
final-look value); no existing test changes were needed.tsa_hr()'s verbose console output now ends with a final note on the
practical impact of the method (τ² estimator) choice: it may have
limited influence on the pooled effect-size when the evidence base is
substantial, but can materially influence heterogeneity-dependent
quantities, prediction intervals, DARIS, and the timing of TSA
conclusions -- particularly when cumulative information is near the
DARIS threshold.tsa_hr() gains a method argument to choose the heterogeneity-variance
(tau^2) estimator used for the random-effects meta-analysis and the
cumulative (sequential) TSA model, passed straight through to
\code{metafor::rma(method = ...)}. Accepts any of metafor's random-effects
estimators: \code{"DL"} (DerSimonian-Laird), \code{"HE"} (or its alias
\code{"CO"}), \code{"HS"}, \code{"HSk"}, \code{"SJ"}, \code{"ML"},
\code{"REML"}, \code{"EB"}, \code{"PM"}, \code{"GENQ"}, \code{"PMM"}, or
\code{"GENQM"}. Default remains \code{"DL"} for full backward
compatibility with all earlier tsahr versions (which always used
DerSimonian-Laird internally and offered no choice) -- this is a
deliberate difference from \code{metafor::rma()}'s own default of
\code{"REML"}. An invalid method value now fails fast with a clear
error listing the accepted values, rather than propagating into
metafor::rma()."=== Random-effects (...) meta-analysis ===") and
the plot caption ("Methods: Random-effects (...) model, ...") now
reflect whichever method was actually used, instead of always saying
"DerSimonian-Laird".method = "FE", regardless of the new method
argument, so D^2 retains its original definition. Alpha-/beta-spending
boundary calculations (the O'Brien-Fleming-type spending functions) are
also unaffected -- this release does not touch, and does not validate
against rpact, any of that machinery.tsa_hr()'s returned parameters list now includes the method that
was used, for downstream inspection/reporting.plot.tsa_hr() example crashed with conversion failure ... in 'mbcsToSbcs': for \u2248 (U+2248). The DARIS-information
marker label used the Unicode "almost equal to" character (\u2248) in
text handed to ggplot2::annotate(); on Windows, grid's text-drawing
routine tries to convert the string to the graphics device's native
single-byte codepage, and that conversion hard-fails for a codepage
that can't represent U+2248 (this is a real, environment-dependent
crash for affected users generally, not just a check artifact -- R CMD
check's "non-ASCII characters" step doesn't catch it because it checks
source-file encoding declarations, not runtime grid-rendering behavior
on a given codepage). Replaced with a plain ASCII ~, matching the
sibling "Theoretical DARIS event-equivalent ~ ..." label two lines
above, which already used ~ for the identical "approximately"
meaning. Confirmed no other non-ASCII characters (literal or \u
escapes) remain anywhere in R/.tail() was called bare (no namespace qualification) in
.rtsa_beta_boundary() (R/obf_boundaries.R, 2 call sites) and in
tsa_hr() (R/tsa_hr.R). Changed to utils::tail() at all three sites,
consistent with this package's existing convention of explicit
pkg::fn() calls rather than blanket importFrom — no NAMESPACE change
needed."caption size/face are customizable and actually applied" asserted expect_null(p_none$theme$plot.caption)
for plot(res, caption = FALSE). plot.tsa_hr()'s actual behavior here
was already correct (the whole labs(caption = ...) +
theme(plot.caption = ...) block is skipped when caption = FALSE) --
the test itself was the problem: newer, S7-based ggplot2 theme objects
resolve an unset plot.caption slot to an inherited element_text
object when accessed via $, rather than the literal NULL older
ggplot2's plain-list themes returned, so the assertion broke purely from
a ggplot2 version upgrade, not from any change in tsahr. Replaced the
check with expect_null(p_none$labels$caption), which tests the actual
intent (labs(caption = ...) was never called) and is stable across
ggplot2's internal theme representation.plot.tsa_hr() gains two new arguments for the methods caption at the
bottom of the plot (previously fixed at size 8, italic):
caption_size: font size (default 8, unchanged from before).caption_face: font face, e.g. "italic" (default, unchanged from
before) or "plain"; also accepts anything
ggplot2::element_text() understands for face (e.g. "bold",
"bold.italic").
Example: plot(res, caption_size = 10, caption_face = "plain").This release does not change any file under R/ -- tsa_hr.R,
methods.R, and obf_boundaries.R are byte-identical to 0.2.4.4 (this
was verified directly, not assumed). Both the alpha-boundary recursive
engine and the RTSA-matched beta-boundary engine are unchanged, per
explicit instruction that both are deliberate design choices, not bugs
to be second-guessed here.
A review of 0.2.4.4 claimed three things: (1) the alpha-spending formula
had "regressed" to a one-sided-parameterized bug, with a specific claim
that the correct first-look boundary was 4.88 rather than 4.38; (2) the
package fails R CMD check with real errors, not just warnings, due to
a test file that used readLines(file.path("R", "tsa_hr.R")) (a
relative path that does not exist when tests run against an installed
package) and a second test asserting stale (pre-RTSA-match) behavior at
the final look; (3) none of this had ever actually been run.
Each claim was checked independently before acting on it, via independent recomputation in Python rather than by running the package itself:
Claim 1 (alpha regression): not substantiated, and not changed.
.alpha_spend_OF() -- 2*(1-Phi(z_alpha/2/sqrt(t))) -- was
independently re-verified via a closed-form (non-simulation) 2D
integration at a schedule matching the disputed case (t1=0.0716): the
resulting first-look boundary was 7.3247, which (a) exactly matches
the elementary closed-form O'Brien-Fleming approximation
z_alpha/sqrt(t1) at that t1, (b) reproduces an overall two-sided alpha
of 0.0498 against a 0.05 target via the same independent integration,
and (c) matches the 7.32 tsahr originally reported for this same
scenario several revisions ago. No evidence for "4.38 vs the correct
4.88" could be reproduced. If a genuine discrepancy against a specific
rpact/RTSA run persists, the most likely explanation remains a
difference in the information-fraction SCHEDULE the two tools were
given (see "Information-fraction convention" in ?tsa_hr), not this
formula. See the VALIDATION note in R/obf_boundaries.R for the full
account, including this specific re-check.
Claim 2 (test failures): confirmed, and fixed. Both specific problems were independently reproduced by inspection/computation:
test-daris-futility-stop.R did use a working-directory-dependent
readLines() call that regex-matched literal source text rather
than testing behavior. This has been REPLACED (not merely patched)
with genuine behavioral tests that call tsa_hr() and check
properties of its actual return value. Note: an intermediate
"clean" revision had already replaced this file with a
behavioral-style test, but with a dataset that does not actually
reach DARIS under the stated design (3 studies of 30 information
units each, needing ~157.6 to reach DARIS at target_HR=0.80) --
this was caught by direct recomputation before being trusted, and
the dataset was corrected (100 information units/study) so the
test's own stated premise ("DARIS reached at look 2") is actually
true. That revision's alpha-boundary reference check was also
wrong in a different way: it compared against
.obf_alpha_boundary(1, alpha=0.05) -- a trivial SINGLE-look
calculation (numerically just qnorm(1-alpha/2)=1.960) -- against
the real MULTI-look recursive final boundary, which is a
different, larger quantity (verified numerically: ~1.992 for the
schedule in question, a difference far outside the test's stated
tolerance). Fixed to compute the reference via the actual
multi-look sequence the code path constructs.test-tsa_hr.R asserted is.na(b40[length(b40)]) ("no futility
decision at final look"), contradicting the RTSA-matched engine's
deliberate behavior (confirmed present in the actual code:
b[length(b)] <- stats::qnorm(1 - alpha / 2, ...) in the relevant
branch) of assigning a definitive final-look futility value rather
than leaving it undefined. Updated to assert the actual (correct,
deliberate) value.Corrected the retrospective boundary endpoint when DARIS is reached: the formal alpha/futility curves now terminate at the interpolated cumulative event coordinate corresponding to the observed DARIS information threshold, rather than at the theoretical Schoenfeld event-equivalent DARIS.
The theoretical DARIS event-equivalent remains displayed as a separate reference line.
The cumulative Z curve and observed-study data remain unchanged.
Preserved the 0.2.4.3 early futility-boundary values. The final-endpoint correction does not recompute the RTSA beta recursion with a synthetic t = 1 endpoint, because doing so would alter the earlier futility looks.
NA, matching RTSA's beta_ubound
post-processing rather than plotting them as negative boundaries.tsa_beta_bound / getInnerWedge()
algorithm in RTSA, which itself documents the old TSA functions as
translated from the original Copenhagen Trial Unit Java TSA software.2 * (1 - Phi(qnorm(1-beta/2) / sqrt(t))).t < 1 and the definitive
futility boundary is the conventional two-sided alpha quantile.TSA$beta_engine exposes the RTSA-derived diagnostics.Built from the 0.2.2 release (not from an interim 0.2.3 that had unrelated problems and was withdrawn), applying fixes from a further external statistical review of the 0.2.2 source.
info_fraction was already capped at 1 for every subsequent look --
but crossed_tsa and entered_futility_region were previously
computed with any() across ALL looks, meaning every look after DARIS
was reached got tested against the same final (t=1) boundary as if
each were an independent additional "final analysis". This is not
accounted for by the alpha-spending calculation, which treats t=1 as a
single final look. Fixed: both verdicts are now evaluated only through
the FIRST look at which DARIS was reached, exposed as the new
results$final_tsa_look (equal to the total number of looks if DARIS
was never reached). The full cumulative Z-curve, including any studies
added after DARIS, is still returned and plotted in full -- only the
formal sequential decision is restricted, not what is shown. verbose
output now prints a note identifying which look this was, when
relevant.inst/CITATION cited
Miladinovic et al. (2013) with DOI 10.1016/j.jclinepi.2013.01.007,
which belongs to an unrelated article ("How to write a research
paper"). Corrected to the actual DOI, 10.1016/j.jclinepi.2012.11.007
(verified against PubMed/ScienceDirect/Wikidata, not just taken on
trust).alpha_two_sided and power (must be
finite, length-1, strictly between 0 and 1) and for target_HR (must
be a single finite value, or NA) -- previously e.g.
alpha_two_sided = 2 would silently propagate into qnorm() rather
than being rejected with a clear message.tsa_hr() now requires at least two studies (stop), and warns (not
errors) below 10 studies, since heterogeneity/D2 -- and therefore
DARIS and the monitoring boundaries -- can be unstable with few
studies (per the Copenhagen TSA manual's own caution on this point).Events_Treatment, N_treatment,
Events_controls, N_controls) must now be whole numbers.final_tsa_look.seq(0.02, 1, length.out = 39) was commented
as "K=39, unequally spaced" but is in fact equally spaced (constant
step). Replaced with a genuinely unequal schedule,
(seq_len(39)/39)^1.4 with the last point fixed at 1.tests/testthat/Rplots.pdf build artifact from
the package source, and added a .Rbuildignore to keep future stray
plot output from being shipped.Per the same reviewer's own priority ranking, and consistent with prior NEWS entries: the two-sided beta/futility betaAdjustment geometry and independent numerical comparison against another group sequential package (e.g. rpact) both remain open. The next step recommended in the 0.2.2 review -- and still the recommendation here -- is external numerical validation against an established group-sequential implementation.
Beta/futility spending fixed. .beta_spend_OF() was reusing the
same "doubled" form as the alpha-spending function, which made its
cumulative target at t=1 equal 2*beta (e.g. 0.40 for the default
beta=0.20) rather than beta. The alpha function's doubling reflects
genuine two-tailed spending (symmetric outer rejection regions); the
futility construction instead targets a single central "inner wedge"
and should not be doubled. Fixed to beta*(t) = 1 - Phi(z_beta/sqrt(t)),
which reaches exactly beta at t=1. This only affects the non-binding,
advisory futility band -- it has no effect on type-I error control
(governed entirely by the independently-validated alpha engine).
Independent verification performed via a faithful line-by-line
Python/scipy reimplementation of R/obf_boundaries.R, cross-checked
three independent ways:
R/obf_boundaries.R rather than left as an implicit assumption.
See the VALIDATION comment block at the top of R/obf_boundaries.R
for the full write-up. No comparison against rpact specifically was
performed; this validation is independent of, and in addition to, the
package's own R test suite.Terminology tightened. info_accrued/info_fraction are now
documented explicitly as reported, study-level inverse-variance
information (sum(1/SE^2) from each study's reported standard error),
not the exact Fisher information of the cumulative random-effects
estimator -- since tau^2 is re-estimated at every cumulative look
under the random-effects model, these are not exactly the same
process. "Exact" language was removed from code comments, roxygen
docs, and the alpha-boundary test description; replaced with
"observed"/"reported" and a note that the Lan-DeMets alpha-spending
approximation to O'Brien-Fleming (used here) is closely related to,
but not numerically identical to, the original O'Brien-Fleming (1979)
construction.
D2 cap wording clarified: the 99.9% cap on Diversity D2 is documented as a purely numerical safeguard against division by (near-)zero, not a statistically justified correction.
Stronger input validation: tsa_hr() now explicitly checks
Events_Treatment, N_treatment, Events_controls, and N_controls
for finite values (previously NA/NaN/Inf could silently pass some of
the range checks, since comparisons against non-finite values return
NA rather than TRUE/FALSE).
New order_by argument for tsa_hr(): optionally names a column
to explicitly sort by (ascending, e.g. a publication-year column)
before the cumulative analysis, rather than relying solely on however
the input data happened to be ordered. Default NULL preserves the
previous behaviour (row order assumed chronological), with a printed
reminder that order_by is available.
Expanded the alpha-boundary Monte Carlo test to include a K=39
unequally-spaced configuration, and added a unit test asserting the
corrected beta-spending function reaches exactly beta (not 2*beta)
at t=1.
Two items from the external review are explicitly not resolved here, and are flagged for a future release rather than attempted without adequate tooling:
rpact (or another
independently validated group-sequential package) for a range of K
and information schedules.plot.tsa_hr() gains a show_theoretical_daris argument (default
TRUE). Set show_theoretical_daris = FALSE to hide the theoretical
DARIS event-equivalent reference line/label (e.g. when only the
observed-information "DARIS information reached" marker is of
interest, or to reduce clutter). This only affects what is drawn --
the underlying DARIS calculation and the printed "DARIS reached"
verdict are unchanged.Fixed a discrepancy where verbose/print() output could report
Required information size (DARIS) reached: YES while plot() still
showed the cumulative Z-curve short of the DARIS reference line (or vice
versa).
Cause: plot.tsa_hr() placed its DARIS vertical line at the
theoretical event-equivalent of the required information
(DARIS_events, derived from the Schoenfeld formula assuming a constant
psi*(1-psi) amount of information per event), while the "DARIS
information reached" verdict was based on the observed accrued
inverse-variance statistical information (sum(1/SE^2) across studies
vs. DARIS_info). These two quantities only agree when every included
study's actual information-per-event matches the pooled psi*(1-psi)
assumption; with heterogeneous studies (differing allocation ratios,
censoring patterns), they can diverge.
Fix: tsa_hr() now also computes DARIS_info_threshold_events, an
estimated cumulative-events position (via linear interpolation between
looks) at which the observed information first meets DARIS_info,
whenever DARIS has been reached. plot.tsa_hr() now shows both
reference lines, separately labelled, rather than substituting one for
the other:
"Theoretical DARIS event-equivalent ~ N", always shown), and"DARIS information reached \u2248 M events (est.)").This means the plot can no longer contradict the printed verdict, while
keeping the two quantities visually and terminologically distinct: the
interpolated crossing point is deliberately never called "DARIS events"
(no study was actually observed at that exact event count), and the
observed accrued information is no longer referred to as "exact" in
code comments/docs, since tau^2 is re-estimated at every cumulative
look under the random-effects model.
verbose = TRUE output and the summary table now report both the
theoretical DARIS event-equivalent and (when applicable) the estimated
event count at which the DARIS information threshold was reached, so
both figures are visible together rather than one replacing the other.
New plot() arguments info_threshold_label_size,
info_threshold_label_x, info_threshold_label_y control the position
of the new label.