Getting started with infercalcR
Raban Heller
Source:vignettes/getting-started.Rmd
getting-started.Rmd
infercalcR carries out statistical procedures that the analyst names in advance. It never chooses a test from the data, for example after a normality pretest. Every procedure returns the same kind of object: a table of finalized rows (estimates, intervals, p-values, counts and a status), a list of metadata (method identity, settings, population and provenance) and a method-specific payload.
This tutorial follows that workflow on a small synthetic data set: describe the data, specify a method, compute it and read the result. It then converts results produced by base R and adjusts a declared family of p-values. The code runs offline and, apart from one optional chunk, needs only the packages that ship with R.
Synthetic data
The data imitate a two-arm study with three continuous outcomes and one binary characteristic. Two blood-pressure values are missing on purpose.
set.seed(2026)
trial <- data.frame(
id = 1:24,
arm = factor(rep(c("control", "treatment"), each = 12)),
sbp = round(c(rnorm(12, mean = 142, sd = 11), rnorm(12, mean = 132, sd = 11)), 1),
dbp = round(c(rnorm(12, mean = 88, sd = 8), rnorm(12, mean = 84, sd = 8)), 1),
heart_rate = round(c(rnorm(12, mean = 74, sd = 9), rnorm(12, mean = 73, sd = 9))),
smoker = factor(sample(c("no", "yes"), 24, replace = TRUE, prob = c(0.7, 0.3)))
)
trial$sbp[c(5, 17)] <- NA
head(trial)
#> id arm sbp dbp heart_rate smoker
#> 1 1 control 147.7 88.4 81 yes
#> 2 2 control 130.1 103.3 78 no
#> 3 3 control 143.5 101.8 63 no
#> 4 4 control 141.1 88.5 72 no
#> 5 5 control NA 93.2 66 no
#> 6 6 control 114.3 101.8 78 noDescribe the data
describe_data() summarizes variables, optionally within
groups. It computes no hypothesis tests. Continuous variables get the
mean, standard deviation, median, quartiles and range; categorical
variables get counts and proportions.
desc <- describe_data(
trial,
variables = c("sbp", "smoker"),
group = "arm",
labels = c(sbp = "Systolic blood pressure", smoker = "Current smoker"),
units = c(sbp = "mmHg")
)
class(desc)
#> [1] "descriptive_result" "inference_result"
desc_rows <- result_data(desc)
desc_rows[, c("group", "variable_id", "estimand", "level", "estimate",
"proportion", "n", "n.missing", "n.total")]
#> group variable_id estimand level estimate proportion n n.missing
#> 1 control sbp mean <NA> 135.718182 NA 11 1
#> 2 control sbp sd <NA> 9.023061 NA 11 1
#> 3 control sbp median <NA> 136.800000 NA 11 1
#> 4 control sbp q1 <NA> 132.350000 NA 11 1
#> 5 control sbp q3 <NA> 142.150000 NA 11 1
#> 6 control sbp min <NA> 114.300000 NA 11 1
#> 7 control sbp max <NA> 147.700000 NA 11 1
#> 8 control smoker count no 8.000000 0.6666667 12 0
#> 9 control smoker count yes 4.000000 0.3333333 12 0
#> 10 treatment sbp mean <NA> 130.036364 NA 11 1
#> 11 treatment sbp sd <NA> 12.747021 NA 11 1
#> 12 treatment sbp median <NA> 129.600000 NA 11 1
#> 13 treatment sbp q1 <NA> 125.750000 NA 11 1
#> 14 treatment sbp q3 <NA> 137.000000 NA 11 1
#> 15 treatment sbp min <NA> 104.000000 NA 11 1
#> 16 treatment sbp max <NA> 148.100000 NA 11 1
#> 17 treatment smoker count no 10.000000 0.8333333 12 0
#> 18 treatment smoker count yes 2.000000 0.1666667 12 0
#> n.total
#> 1 12
#> 2 12
#> 3 12
#> 4 12
#> 5 12
#> 6 12
#> 7 12
#> 8 12
#> 9 12
#> 10 12
#> 11 12
#> 12 12
#> 13 12
#> 14 12
#> 15 12
#> 16 12
#> 17 12
#> 18 12Missing values are counted in n.missing and kept apart
from the n used for each summary. Category proportions use
the non-missing observations within a group as their denominator. The
result has the same structure as an inferential result, so the same
accessors work on it.
Choose a method explicitly
method_catalogue() lists the method identities that
run_inference() accepts, the package that provides the
computation and the option names each method understands.
catalogue <- method_catalogue()
dim(catalogue)
#> [1] 69 6
table(catalogue$implemented)
#>
#> FALSE TRUE
#> 1 68
catalogue[c("t_means_two_independent_welch", "rank_sum_mann_whitney",
"regression_negative_binomial"), c("backend", "implemented")]
#> backend implemented
#> t_means_two_independent_welch stats TRUE
#> rank_sum_mann_whitney stats TRUE
#> regression_negative_binomial MASS,stats TRUEMost methods use only stats and utils. A
few rely on an optional package, such as MASS for negative binomial
regression. Such a package is used only when it is installed; if it is
missing, the method stops with an explicit error and no other model is
fitted in its place. ?observed-methods lists the arguments
that each method takes.
A specification fixes the method, the alternative hypothesis, the confidence level, the missing-data policy and any method options before the data are seen.
spec <- method_spec("t_means_two_independent_welch",
alternative = "two.sided", conf.level = 0.95)
str(unclass(spec))
#> List of 8
#> $ method_id : chr "t_means_two_independent_welch"
#> $ variant_id : chr "t_means_two_independent_welch"
#> $ alternative : chr "two.sided"
#> $ conf.level : num 0.95
#> $ missing : chr "error"
#> $ options : list()
#> $ engine : chr "base"
#> $ contract_version: chr "1.0.0"Unknown options and impossible settings are rejected when the specification is created, not silently ignored later.
tryCatch(
method_spec("t_means_two_independent_welch", options = list(var.equal = TRUE)),
error = conditionMessage
)
#> [1] "unknown options: var.equal"
tryCatch(
method_spec("t_means_two_independent_welch", conf.level = 1),
error = conditionMessage
)
#> [1] "conf.level must be in (0, 1)"Compute the result
The default missing-data policy is "error". Because the
blood-pressure values contain missing entries, the Welch test refuses to
run on them.
treated <- trial$sbp[trial$arm == "treatment"]
control <- trial$sbp[trial$arm == "control"]
tryCatch(run_inference(spec, x = treated, y = control), error = conditionMessage)
#> [1] "all(is.finite(x)) is not TRUE"missing = "omit" does not change this for the Welch
test: the package does not pick a missing-data assumption on the
analyst’s behalf. Here we decide to analyze the complete observations
and do so visibly in the code. A result_id gives the result
a stable name for later joins.
x_complete <- treated[!is.na(treated)]
y_complete <- control[!is.na(control)]
welch <- run_inference(spec, x = x_complete, y = y_complete,
result_id = "sbp-welch")
welch
#> <inference> t_means_two_independent_welch
#> row_id estimand
#> t_means_two_independent_welch:1 mean difference (x minus reference or y)
#> estimate p.reported status
#> -5.681818 0.2431838 okRead the result contract
result_data() returns the finalized rows as a data
frame. Every result has the same mandatory columns; unavailable
quantities are NA.
rows <- result_data(welch)
rows[, c("estimand", "estimate", "conf.low", "conf.high", "conf.level",
"interval.method")]
#> estimand estimate conf.low conf.high
#> 1 mean difference (x minus reference or y) -5.681818 -15.57427 4.210632
#> conf.level interval.method
#> 1 0.95 Welch t
rows[, c("statistic", "df", "p.value", "p.adjusted", "p.reported",
"alternative", "n", "n.excluded", "status", "reason")]
#> statistic df p.value p.adjusted p.reported alternative n n.excluded
#> 1 -1.206635 18.01017 0.2431838 NA 0.2431838 two.sided 22 0
#> status reason
#> 1 ok <NA>The estimand names the quantity, including its
direction: here the mean of x (treatment) minus the mean of
y (control). The interval columns record the level, the
method and the scale. p.value is the raw p-value,
p.adjusted stays empty until a multiplicity adjustment is
requested, and p.reported is the value meant for display.
Check status before using an estimate: besides
"ok", a row can be "undefined",
"not_applicable", "unsupported",
"invalid_input", "numerical_failure",
"infeasible", "reference_pending" or
"unresolved", with an explanation in
reason.
result_metadata() records how the result came about.
meta <- result_metadata(welch)
meta[c("method_id", "variant_id", "result_id", "kind", "contract_version")]
#> $method_id
#> [1] "t_means_two_independent_welch"
#>
#> $variant_id
#> [1] "t_means_two_independent_welch"
#>
#> $result_id
#> [1] "sbp-welch"
#>
#> $kind
#> [1] "inference"
#>
#> $contract_version
#> [1] "1.0.0"
str(meta$settings[c("alternative", "conf.level", "missing", "options")])
#> List of 4
#> $ alternative: chr "two.sided"
#> $ conf.level : num 0.95
#> $ missing : chr "error"
#> $ options : list()
str(meta$population)
#> List of 3
#> $ n_input : int 22
#> $ n_used : num 22
#> $ n_excluded: num 0
str(meta$provenance)
#> List of 3
#> $ package: chr "infercalcR"
#> $ version: chr "0.1.0"
#> $ R : chr "4.6.2"The payload keeps the complete output of the underlying computation, including quantities that do not fit in a single row.
str(welch$payload$native, give.attr = FALSE)
#> List of 8
#> $ method : chr "Welch Two Sample t-test"
#> $ n : Named int [1:2] 11 11
#> $ mean_difference : num -5.68
#> $ confidence_interval: num [1:2] -15.57 4.21
#> $ statistic : num -1.21
#> $ df : num 18
#> $ p_value : num 0.243
#> $ alternative : chr "two.sided"The numbers agree with a direct call to stats::t.test()
on the same vectors.
Prespecified alternatives use the same pattern
A different procedure is another specification, not a fallback chosen
after looking at the data. The Mann-Whitney rank-sum test, for example,
reports a probability of superiority as its descriptive estimate. The
exact option is declared in the specification.
rank_spec <- method_spec("rank_sum_mann_whitney", options = list(exact = FALSE))
rank_sum <- run_inference(rank_spec, x = x_complete, y = y_complete,
result_id = "sbp-rank-sum")
result_data(rank_sum)[, c("estimand", "estimate", "p.value", "n", "status")]
#> estimand
#> 1 identical group distributions; descriptive probability of superiority
#> estimate p.value n status
#> 1 0.3305785 0.1890818 22 okChoose the engine
Thirteen methods, the Brunner-Munzel test among them, have two
engines. engine = "base", the default, uses infercalcR’s
own implementation; engine = "reference" computes the same
procedure with the validated original package, here brunnermunzel. The
original is optional: this chunk runs only when it is installed. Without
it, the reference engine stops with an error instead of using the base
engine.
bm_base <- run_inference(method_spec("brunner_munzel"), x = x_complete, y = y_complete)
bm_reference <- run_inference(method_spec("brunner_munzel", engine = "reference"),
x = x_complete, y = y_complete)
rbind(base = result_data(bm_base)[1, c("estimate", "p.value")],
reference = result_data(bm_reference)[1, c("estimate", "p.value")])
#> estimate p.value
#> base 0.3305785 0.2027142
#> reference 0.3305785 0.2027142
result_metadata(bm_reference)$engine
#> $engine
#> [1] "reference"
#>
#> $package
#> [1] "brunnermunzel"
#>
#> $version
#> [1] "2.0"
#>
#> $call
#> [1] "brunnermunzel::brunnermunzel.test"Convert results computed elsewhere
as_inference() turns a supported htest,
lm or glm object into the same contract. It
does not refit anything. Sample sizes are not recovered from degrees of
freedom, so they are supplied together with a description of the
analysis population.
imported <- as_inference(
reference,
n = c(length(x_complete), length(y_complete)),
population = "Synthetic trial, complete observations",
groups = c("treatment", "control")
)
result_data(imported)[, c("estimand", "estimate", "conf.low", "conf.high",
"p.value", "n")]
#> estimand estimate conf.low conf.high p.value n
#> 1 population mean difference -5.681818 -15.57427 4.210632 0.2431838 22A fitted linear model is summarized once, with classical t intervals for each coefficient.
model <- lm(sbp ~ arm, data = trial)
coefficients <- result_data(
as_inference(model, population = "Synthetic trial, complete observations")
)
print(coefficients[, c("contrast", "estimate", "conf.low", "conf.high", "p.value")],
row.names = FALSE)
#> contrast estimate conf.low conf.high p.value
#> (Intercept) 135.718182 128.77269 142.663677 1.004296e-20
#> armtreatment -5.681818 -15.50423 4.140595 2.416520e-01lm() drops the two incomplete rows through its own
na.action; the population label states this explicitly.
Adjust a declared family of p-values
Suppose the three continuous outcomes form a prespecified family of
secondary comparisons. Each is analyzed with the same specification, and
the rows are collected into one result with
new_inference_result().
secondary <- lapply(c("sbp", "dbp", "heart_rate"), function(variable) {
values <- trial[[variable]]
keep <- !is.na(values)
fit <- run_inference(spec,
x = values[keep & trial$arm == "treatment"],
y = values[keep & trial$arm == "control"])
out <- result_data(fit)[, c("estimand", "estimate", "conf.low", "conf.high",
"conf.level", "p.value", "alternative", "n")]
out$row_id <- variable
out
})
family <- new_inference_result(
do.call(rbind, secondary),
metadata = list(method_id = "t_means_two_independent_welch",
result_id = "secondary-outcomes")
)
class(validate_inference(family))
#> [1] "inference_result"adjust_inference() adjusts the complete family once and
records the family identifier and size. The raw p-values and the
confidence intervals stay unchanged.
adjusted <- adjust_inference(family, method = "holm", family_id = "secondary")
result_data(adjusted)[, c("row_id", "estimate", "p.value", "p.adjusted",
"p.reported", "p.method", "family_id", "family_size")]
#> row_id estimate p.value p.adjusted p.reported p.method family_id
#> 1 sbp -5.68181818 0.24318377 0.486367549 0.486367549 holm secondary
#> 2 dbp -9.42500000 0.00304647 0.009139409 0.009139409 holm secondary
#> 3 heart_rate -0.08333333 0.97732161 0.977321607 0.977321607 holm secondary
#> family_size
#> 1 3
#> 2 3
#> 3 3A second adjustment of the same rows is refused, so a reported p-value cannot be adjusted twice by accident.
tryCatch(adjust_inference(adjusted, family_id = "secondary-again"),
error = conditionMessage)
#> [1] "family already has adjustment metadata; do not adjust twice"Keep a record
A result is an ordinary S3 list and can be saved as RDS. Load RDS files only from trusted sources.
path <- tempfile(fileext = ".rds")
saveRDS(welch, path)
restored <- readRDS(path)
stopifnot(identical(result_data(restored), result_data(welch)))
unlink(path)Archive the specification and the result together with the package
version (recorded in the provenance metadata), the repository revision
and the output of sessionInfo().
Continue
The diagnostic accuracy tutorial
covers classifier metrics and cutoffs; the exact tests tutorial covers Boschloo’s test
for two-by-two tables. ?result-contract describes the
result columns in full.
Acknowledgments and sources
This package is part of an inference suite (infercalcR, inferguidR and infereportR) that is deeply inspired by DescrTab2 [1], developed at the Institute of Medical Biometry (IMBI), Heidelberg University, by Jan Meis, Lukas Baumann, Maximilian Pilz and Lukas Sauer. Its test-choice tree, its validation statement and its practice of computing each test with an established implementation are the model for how the suite chooses, computes and reports statistical tests. DescrTab2 was archived on CRAN on 18 October 2025; its source code remains available on GitHub. The methods-consulting pages of the University of Zurich [2], with their decision tree, decision assistant, worked examples and explanations, inspired the guidance in the suite. Power and sample-size planning and the distribution calculator follow G*Power 3.1 [3, 4]. We thank the authors of these works for making them openly available. The package README cites all software the package uses.
- Meis J, Baumann L, Pilz M, Sauer L (2022). DescrTab2: Publication Quality Descriptive Statistics Tables. R package version 2.1.16. https://doi.org/10.32614/CRAN.package.DescrTab2. Source code, including the development version 2.1.26: https://github.com/imbi-heidelberg/DescrTab2.
- Universität Zürich, Methodenberatung. Datenanalyse mit SPSS: decision tree (Entscheidbaum), decision assistant (Entscheidassistent) and method pages with examples and explanations. https://www.methodenberatung.uzh.ch/de/datenanalyse_spss.html. Accessed 3 October 2026.
- Faul F, Erdfelder E, Lang A-G, Buchner A (2007). G*Power 3: A flexible statistical power analysis program for the social, behavioral, and biomedical sciences. Behavior Research Methods, 39(2), 175-191. https://doi.org/10.3758/BF03193146.
- Faul F, Erdfelder E, Buchner A, Lang A-G (2009). Statistical power analyses using G*Power 3.1: Tests for correlation and regression analyses. Behavior Research Methods, 41(4), 1149-1160. https://doi.org/10.3758/BRM.41.4.1149.
sessionInfo()
#> R version 4.6.2 (2026-10-29)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 26.04.1 LTS
#>
#> Matrix products: default
#> BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.32.so; LAPACK version 3.12.0
#>
#> locale:
#> [1] LC_CTYPE=C.UTF-8 LC_NUMERIC=C LC_TIME=C.UTF-8
#> [4] LC_COLLATE=C.UTF-8 LC_MONETARY=C.UTF-8 LC_MESSAGES=C.UTF-8
#> [7] LC_PAPER=C.UTF-8 LC_NAME=C LC_ADDRESS=C
#> [10] LC_TELEPHONE=C LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C
#>
#> time zone: UTC
#> tzcode source: system (glibc)
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] infercalcR_0.1.0
#>
#> loaded via a namespace (and not attached):
#> [1] digest_0.6.39 desc_1.4.3 R6_2.6.1 fastmap_1.2.0
#> [5] xfun_0.62 cachem_1.1.0 knitr_1.52 htmltools_0.5.9
#> [9] rmarkdown_2.32 lifecycle_1.0.5 cli_3.6.6 sass_0.4.10
#> [13] pkgdown_2.2.1 textshaping_1.0.5 jquerylib_0.1.4 systemfonts_1.3.2
#> [17] compiler_4.6.2 tools_4.6.2 ragg_1.5.2 bslib_0.12.0
#> [21] evaluate_1.0.5 yaml_2.3.12 brunnermunzel_2.0 otel_0.3.0
#> [25] jsonlite_2.0.0 rlang_1.3.0 fs_2.1.0