Hypothesis Testing with {statim}

Rationale

While t.test(x, alternative = "less") works and is simple, there are times you won’t remember the intent of the function call, e.g. which side of the inequality “less” refers to:

Is x hypothesized to be less than mu, or mu less than x?

statim takes a different approach on how statistical inference in R is done. For instance, statim has another way to write the null hypothesis itself as an algebraic expression instead (e.g. MU(x) < 120, RHO(x, y) == 0), so the code declares what the null hypothesis is, direction included, rather than which argument slot lives in.

statim is fully declarative, mainly through piped/grammar sementics like ggplot2, including on how you declare the estimation method in a statistical inference you want to perform. The question is, does that actually pay off?

This vignette is a showcase of statim which runs four questions from a real dataset through statim, base R, and rstatix, and is honest about where the extra syntax doesn’t earn its keep. These four are the straightforward cases; the harder one, where the syntax stops being optional, is covered in the Conclusion below.

Sample Problem

The dataset: 303 patients from the Cleveland Clinic. Originally collected by Robert Detrano, M.D., Ph.D., and hosted at the UCI Machine Learning Repository (Heart Disease dataset, Cleveland subset). This copy was redistributed by Daniel Bourke’s “zero-to-mastery-ml” repo, and is bundled with this package at inst/extdata/heart-disease.csv so the vignette builds without a network call.

The questions, framed as a cardiologist would ask them:

Is average resting blood pressure different from the clinical normal of 120 mmHg?
Do men and women differ in maximum heart rate achieved during stress testing?
Is there a linear relationship between age and maximum heart rate?
Is the proportion of male patients with fasting blood sugar above 120 mg/dL different from an assumed population baseline of 15%?

They are answered by: one-sample t-test, two-sample t-test, correlation test, one-sample proportion test, respectively.

Setup

statim’s author reaches for box day to day. It’s an R package that bring an another but better import system that forces every dependency to be declared explicitly, so a script’s imports double as its own dependency graph.

box::use(
    stats[t.test, cor.test, binom.test],
    statim[
        T_TEST, COR_TEST, P_TEST,
        # Current grammars
        define_model, prepare, via, state_null, conclude,
        # Multiple executions
        write_models, display,
        # Mappers
        x_by, rel, prop, on,
        # "Parameters" object callers
        MU, RHO, PI
    ],
    rstatix[t_test, cor_test, binom_test]
)

box::use(
    dplyr[keep_when = filter, mutate, glimpse],
    readr[read_csv]
)

Then import the CSV file from this package:

heart = read_csv(system.file("extdata", "heart-disease.csv", package = "statim")) |>
    mutate(
        sex = factor(sex, levels = c(0, 1), labels = c("Female", "Male")),
        fbs = factor(fbs, levels = c(0, 1), labels = c("Normal", "High"))
    )

 [1mRows:  [22m [34m303 [39m  [1mColumns:  [22m [34m14 [39m
 [36m── [39m  [1mColumn specification [22m  [36m──────────────────────────────────────────────────────── [39m
 [1mDelimiter: [22m ","
 [32mdbl [39m (14): age, sex, cp, trestbps, chol, fbs, restecg, thalach, exang, oldpea...

 [36mℹ [39m Use `spec()` to retrieve the full column specification for this data.
 [36mℹ [39m Specify the column types or set `show_col_types = FALSE` to quiet this message.

glimpse(heart)

Rows: 303
Columns: 14
$ age      <dbl> 63, 37, 41, 56, 57, 57, 56, 44, 52, 57, 54, 48, 49, 64, 58, 5…
$ sex      <fct> Male, Male, Female, Male, Female, Male, Female, Male, Male, M…
$ cp       <dbl> 3, 2, 1, 1, 0, 0, 1, 1, 2, 2, 0, 2, 1, 3, 3, 2, 2, 3, 0, 3, 0…
$ trestbps <dbl> 145, 130, 130, 120, 120, 140, 140, 120, 172, 150, 140, 130, 1…
$ chol     <dbl> 233, 250, 204, 236, 354, 192, 294, 263, 199, 168, 239, 275, 2…
$ fbs      <fct> High, Normal, Normal, Normal, Normal, Normal, Normal, Normal,…
$ restecg  <dbl> 0, 1, 0, 1, 1, 1, 0, 1, 1, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1…
$ thalach  <dbl> 150, 187, 172, 178, 163, 148, 153, 173, 162, 174, 160, 139, 1…
$ exang    <dbl> 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0…
$ oldpeak  <dbl> 2.3, 3.5, 1.4, 0.8, 0.6, 0.4, 1.3, 0.0, 0.5, 1.6, 1.2, 0.2, 0…
$ slope    <dbl> 0, 0, 2, 2, 2, 1, 1, 2, 2, 2, 2, 2, 2, 1, 2, 1, 2, 0, 2, 2, 1…
$ ca       <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 0…
$ thal     <dbl> 1, 2, 2, 2, 2, 1, 2, 3, 3, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3…
$ target   <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…

Here’s an overview for the columns to be used:

Column	Type	Description
`trestbps`	Continuous	Resting blood pressure (mm Hg)
`thalach`	Continuous	Maximum heart rate achieved
`age`	Continuous	Age in years
`sex`	Binary	Sex (0 = Female, 1 = Male)
`fbs`	Binary	Fasting blood sugar > 120 mg/dL (0 = No, 1 = Yes)

Question 1: Is resting blood pressure elevated?

H_0: \mu=120 \qquad H_1: \mu\neq120

This is the simplest case, so it’s worth showing every entry point once.

Codes

There are two layouts to perform one-sample t-test:

Using on()

heart |>
    define_model(on(trestbps)) |>
    prepare(T_TEST, .mu = 120) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : on 
Args : trestbps 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────
    term    estimate  true_mu  t_stat  p_val   
───────────────────────────────────────────────
  trestbps  131.624     120    11.537  <0.001  
───────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

────────────────────────────────
    term    lower_95  upper_95  
────────────────────────────────
  trestbps  129.641   133.607   
────────────────────────────────

T_TEST(on(trestbps), heart, .mu = 120)

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────
    term    estimate  true_mu  t_stat  p_val   
───────────────────────────────────────────────
  trestbps  131.624     120    11.537  <0.001  
───────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

────────────────────────────────
    term    lower_95  upper_95  
────────────────────────────────
  trestbps  129.641   133.607   
────────────────────────────────

Using formula syntax

heart |>
    define_model(trestbps ~ 1) |>
    prepare(T_TEST, .mu = 120) |>
    # update(.mu = 120) |> 
    conclude()


== Model ======================================================================= 

Variable Mapper : formula 
Args : trestbps ~ 1 
    left_var : 1 
    right_var : 0 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

─────────────────────────────────────────────────────────
  groups     type     est_type    est    t-stat   pval   
─────────────────────────────────────────────────────────
    1     one sample     mu     131.624  11.537  <0.001  
─────────────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────────────────────────
  groups     type     lower_95  upper_95  
──────────────────────────────────────────
    1     one sample  129.641   133.606   
──────────────────────────────────────────

T_TEST(trestbps ~ 1, heart, .mu = 120)

-- Summary ---------------------------------------------------------------------

─────────────────────────────────────────────────────────
  groups     type     est_type    est    t-stat   pval   
─────────────────────────────────────────────────────────
    1     one sample     mu     131.624  11.537  <0.001  
─────────────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────────────────────────
  groups     type     lower_95  upper_95  
──────────────────────────────────────────
    1     one sample  129.641   133.606   
──────────────────────────────────────────

t_test(heart, trestbps ~ 1, mu = 120)

# A tibble: 1 × 7
  .y.      group1 group2         n statistic    df        p
* <chr>    <chr>  <chr>      <int>     <dbl> <dbl>    <dbl>
1 trestbps 1      null model   303      11.5   302 9.34e-26

t.test(heart$trestbps, mu = 120)


    One Sample t-test

data:  heart$trestbps
t = 11.537, df = 302, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 120
95 percent confidence interval:
 129.6411 133.6065
sample estimates:
mean of x 
 131.6238

t.test(trestbps ~ 1, heart, mu = 120)


    One Sample t-test

data:  trestbps
t = 11.537, df = 302, p-value < 2.2e-16
alternative hypothesis: true mean is not equal to 120
95 percent confidence interval:
 129.6411 133.6065
sample estimates:
mean of x 
 131.6238

Verdict

All of the packages addresses the problem by performing one-sample t-test in a single line of code. However, what statim buys instead is legibility of intent: define_model(), prepare(), conclude() read like a sentence, and you can do more with statim’s main API, whereas t.test(x, mu = 120) reads like an API you have to already know. And also, ttest-on can utilize state_null(), it is not covered for the reason being the usage is not too different on declaring .mu.

Interpretation: 131.6 mmHg average, significantly above 120 (p < 0.001). This cohort runs elevated.

Question 2: Does max heart rate differ by sex?

H_0: \mu_{\text{thalach}\mid\text{sex=Female}}=\mu_{\text{thalach}\mid\text{sex=Male}} \qquad H_1: \mu_{\text{thalach}\mid\text{sex=Female}} \neq \mu_{\text{thalach}\mid\text{sex=Male}}

Codes

There are three layouts to perform two-sample t-test:

Using on()

This version requires via("two_sample") after prepare(T_TEST) to perform two-sample t-test. Since it requires via("two_sample"), you can’t use its eager form/one liner code.

female = heart$thalach[heart$sex == "Female"]
male = heart$thalach[heart$sex == "Male"]

# Requires via("two_sample")
# To perform two-sample t-test with `on()` layout
define_model(on(female, male)) |>
    prepare(T_TEST) |>
    via("two_sample") |>
    state_null(
        MU(female) == MU(male)
    ) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : on 
Args : female, male 

== T-Test · two_sample ========================================================= 

-- Summary ---------------------------------------------------------------------

────────────────────────────────────────────────────────
        group         estimate  t_stat    df     p_val  
────────────────────────────────────────────────────────
  1*female + -1*male   2.164    0.818   219.790  0.414  
────────────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────────────────────────
        group         lower_95  upper_95  
──────────────────────────────────────────
  1*female + -1*male   -3.050    7.378    
──────────────────────────────────────────

Note: Using state_null() is optional, unless the full null hypothesis expression is required.

Using x_by()

heart |>
    define_model(x_by(thalach, sex)) |>
    prepare(T_TEST) |>
    state_null(
        MU(thalach, sex == "Female") == MU(thalach, sex == "Male")
    ) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : x_by 
Args : thalach | sex 
    x_vars : 1 
    by_vars : 1 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────
  group  estimate  t_stat    df     p_val  
───────────────────────────────────────────
   sex    2.164    0.818   219.790  0.414  
───────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────
  group  lower_95  upper_95  
─────────────────────────────
   sex    -3.050    7.378    
─────────────────────────────

T_TEST(x_by(thalach, sex), heart)

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────
  group  estimate  t_stat    df     p_val  
───────────────────────────────────────────
   sex    -2.164   -0.818  219.790  0.414  
───────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────
  group  lower_95  upper_95  
─────────────────────────────
   sex    -7.378    3.050    
─────────────────────────────

Note: Using state_null() is optional, unless the full null hypothesis expression is required.

Using formula syntax

Currently, the <formula> layout doesn’t have translation for state_null().

heart |>
    define_model(thalach ~ sex) |>
    prepare(T_TEST) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : formula 
Args : thalach ~ sex 
    left_var : 1 
    right_var : 1 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

──────────────────────────────────────────────────────
  groups     type     est_type   est   t-stat  pval   
──────────────────────────────────────────────────────
   sex    two sample  mu_diff   2.164  0.818   0.414  
──────────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────────────────────────
  groups     type     lower_95  upper_95  
──────────────────────────────────────────
   sex    two sample   -3.051    7.378    
──────────────────────────────────────────

T_TEST(thalach ~ sex, heart)

-- Summary ---------------------------------------------------------------------

──────────────────────────────────────────────────────
  groups     type     est_type   est   t-stat  pval   
──────────────────────────────────────────────────────
   sex    two sample  mu_diff   2.164  0.818   0.414  
──────────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────────────────────────
  groups     type     lower_95  upper_95  
──────────────────────────────────────────
   sex    two sample   -3.051    7.378    
──────────────────────────────────────────

t_test(heart, thalach ~ sex)

# A tibble: 1 × 8
  .y.     group1 group2    n1    n2 statistic    df     p
* <chr>   <chr>  <chr>  <int> <int>     <dbl> <dbl> <dbl>
1 thalach Female Male      96   207     0.818  220. 0.414

Two forms, but the <formula> interface is usually more preferred than the regular vector one. For the regular vector, we can use the same vector from the statim section.

Regular Vector

t.test(female, male)


    Welch Two Sample t-test

data:  female and male
t = 0.8178, df = 219.79, p-value = 0.4144
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.050537  7.377831
sample estimates:
mean of x mean of y 
 151.1250  148.9614

Formula Syntax

t.test(thalach ~ sex, data = heart)


    Welch Two Sample t-test

data:  thalach by sex
t = 0.8178, df = 219.79, p-value = 0.4144
alternative hypothesis: true difference in means between group Female and group Male is not equal to 0
95 percent confidence interval:
 -3.050537  7.377831
sample estimates:
mean in group Female   mean in group Male 
            151.1250             148.9614

Verdict

This is where the two designs actually diverge, not just in syntax. Empirically, statim, rstatix, and base R treat “compare two groups” as something the <formula> already encodes, thalach ~ sex says everything. However, statim goes beyond that with x_by() and on(). `

x_by() says the same thing, but then lets you go further: this layout has state_null() translation, and you can write out MU(thalach, sex == "Female") == MU(thalach, sex == "Male") as an actual algebraic expression, group filters and all.
on(), on the other hand, has the same logic as x_by()’s except it treats the variables to be independent to each other, and its null hypothesis expression doesn’t use the <sex == "Male"> given argument.

That’s more typing for a question this simple. It stops being “more typing” and starts being “the only way to say what you mean” the moment the hypothesis isn’t a straight group comparison anymore.

Interpretation: no difference (p = 0.414). Sex isn’t doing any explanatory work here.

Question 3: Does max heart rate fall with age?

H_0: \rho_{\text{thalach, age}}=0 \qquad H_1: \rho_{\text{thalach, age}} \neq 0

Code

There are two layouts to perform correlation test:

Using rel()

heart |>
    define_model(rel(thalach, age)) |>
    prepare(COR_TEST) |>
    state_null(
        RHO(thalach, age) == 0
    ) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : rel 
Args : thalach ; age 
    x_vars : 1 
    resp_vars : 1 

== Correlation Test ============================================================ 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────────
      pair       estimate  statistic  df   p_val   
───────────────────────────────────────────────────
  age ~ thalach   -0.398    -7.539    301  <0.001  
───────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────────────
      pair       lower_95  upper_95  
─────────────────────────────────────
  age ~ thalach   -0.489    -0.299   
─────────────────────────────────────

COR_TEST(rel(thalach, age), heart)

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────────
      pair       estimate  statistic  df   p_val   
───────────────────────────────────────────────────
  age ~ thalach   -0.398    -7.539    301  <0.001  
───────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────────────
      pair       lower_95  upper_95  
─────────────────────────────────────
  age ~ thalach   -0.489    -0.299   
─────────────────────────────────────

Note: Using state_null() is optional, unless the full null hypothesis expression is required.

Using formula syntax

heart |>
    define_model(age ~ thalach) |>
    prepare(COR_TEST) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : formula 
Args : age ~ thalach 
    left_var : 1 
    right_var : 1 

== Correlation Test ============================================================ 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────────
      pair       estimate  statistic  df   p_val   
───────────────────────────────────────────────────
  age ~ thalach   -0.398    -7.539    301  <0.001  
───────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────────────
      pair       lower_95  upper_95  
─────────────────────────────────────
  age ~ thalach   -0.489    -0.299   
─────────────────────────────────────

COR_TEST(age ~ thalach, heart)

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────────
      pair       estimate  statistic  df   p_val   
───────────────────────────────────────────────────
  age ~ thalach   -0.398    -7.539    301  <0.001  
───────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────────────
      pair       lower_95  upper_95  
─────────────────────────────────────
  age ~ thalach   -0.489    -0.299   
─────────────────────────────────────

rstatix::cor_test(heart, age, thalach)

# A tibble: 1 × 9
  var1  var2      cor statistic    df        p conf.low conf.high method 
  <chr> <chr>   <dbl>     <dbl> <int>    <dbl>    <dbl>     <dbl> <chr>  
1 age   thalach  -0.4     -7.54   301 5.63e-13   -0.489    -0.299 Pearson

cor.test(~ thalach + age, heart)


    Pearson's product-moment correlation

data:  thalach and age
t = -7.5386, df = 301, p-value = 5.628e-13
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.4892312 -0.2992831
sample estimates:
       cor 
-0.3985219

Verdict

Small warning: base R’s <formula> here is ~ thalach + age, not thalach ~ age. There’s no dependent variable in a correlation, so the usual left/right convention doesn’t mean anything, but it’s easy to assume it does and misread the formula. {rel(thalach, age)} sidesteps the ambiguity by just naming both variables. rstatix takes some roundabout by borrowing dplyr::select()-style column picking, and the selected variable always in pairwise combination.

Interpretation: a real, moderate negative correlation (r = -0.398, p < 0.001). Heart rate ceiling drops with age, as expected.

Question 4: Is high fasting blood sugar unusually common in men?

H_0: \pi=0.15 \qquad H_1: \pi\neq0.15

Prepare the data first

males = keep_when(heart, sex == "Male")
n_high_fbs = sum(males$fbs == "High")
n_males = nrow(males)

Code

By default, P_TEST() performs a binomial test.

define_model(prop(n_high_fbs, n_males)) |>
    prepare(P_TEST) |>
    state_null(
        PI() == 0.15
    ) |>
    conclude()


== Model ======================================================================= 

Variable Mapper : prop 
Args : 33 / 207 
    x : 33 
    n : 207 

== Proportion Test ============================================================= 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────
  x    n   true_p  estimate  statistic  p_val  
───────────────────────────────────────────────
  33  207  0.150    0.159       33      0.697  
───────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────
  lower_95  upper_95  
──────────────────────
   0.112     0.216    
──────────────────────

P_TEST(prop(n_high_fbs, n_males), .p = 0.15)

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────────
  x    n   true_p  estimate  statistic  p_val  
───────────────────────────────────────────────
  33  207  0.150    0.159       33      0.697  
───────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

──────────────────────
  lower_95  upper_95  
──────────────────────
   0.112     0.216    
──────────────────────

Note: Using state_null() is optional, unless the full null hypothesis expression is required.

binom_test(n_high_fbs, n_males, p = 0.15)

# A tibble: 1 × 6
      n estimate conf.low conf.high     p p.signif
* <int>    <dbl>    <dbl>     <dbl> <dbl> <chr>   
1   207    0.159    0.112     0.217 0.697 ns

binom.test(n_high_fbs, n_males, p = 0.15)


    Exact binomial test

data:  n_high_fbs and n_males
number of successes = 33, number of trials = 207, p-value = 0.6969
alternative hypothesis: true probability of success is not equal to 0.15
95 percent confidence interval:
 0.1123500 0.2165365
sample estimates:
probability of success 
             0.1594203

Verdict

No <formula> anywhere in this section. A proportion test is just two numbers, x and n, and all three packages treat it that way. The only real question is where those two numbers live: as positional arguments (x, n, p =) in base R and rstatix, or wrapped in prop() so the pipeline keeps the same shape it had in Questions 1-3. Neither is more correct; statim’s version is only worth it if you’re already committed to the pipeline for other reasons.

Interpretation: 33 of 207 men (15.9%), not distinguishable from the 15% baseline (p = 0.697).

Multiple Executions

Every question above ran one test through one layout. write_models() lets you batch several layouts of the same test into a single pipeline, instead of writing out define_model() |> prepare() |> conclude() once per layout:

out =
    heart |>
    write_models(
        mod1 = on(trestbps),
        mod2 = x_by(thalach, sex),
        mod3 = thalach ~ sex
    ) |>
    prepare(T_TEST) |>
    conclude()

out


── 3 models · T-Test ─────────────────────────────────────────────────────────── 

mod1 : <cld_exec>
mod2 : <cld_exec>
mod3 : <cld_exec>

Use display() to inspect individual results.

Each name becomes its own lazy model behind the scenes, prepare(T_TEST) attaches the same test to all three at once, and conclude() runs each independently and hands back a <multi_exec>: tidy() or display() on it to inspect the individual results, as covered in Execution and Retrieval of Outputs.

display(out, 2)


1. mod1

== Model ======================================================================= 

Variable Mapper : on 
Args : trestbps 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

────────────────────────────────────────────────
    term    estimate  true_mu  t_stat   p_val   
────────────────────────────────────────────────
  trestbps  131.624      0     130.639  <0.001  
────────────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

────────────────────────────────
    term    lower_95  upper_95  
────────────────────────────────
  trestbps  129.641   133.607   
────────────────────────────────



2. mod2

== Model ======================================================================= 

Variable Mapper : x_by 
Args : thalach | sex 
    x_vars : 1 
    by_vars : 1 

== T-Test ====================================================================== 

-- Summary ---------------------------------------------------------------------

───────────────────────────────────────────
  group  estimate  t_stat    df     p_val  
───────────────────────────────────────────
   sex    -2.164   -0.818  219.790  0.414  
───────────────────────────────────────────


-- Confidence Interval ---------------------------------------------------------

─────────────────────────────
  group  lower_95  upper_95  
─────────────────────────────
   sex    -7.378    3.050    
─────────────────────────────

On a side note, prepare()’s ... forwards into the spec on a write_models() batch too, the same way it already does on a single <def_var>. So .mu = 120 below reaches every model in the batch whose fn actually has a .mu formal:

heart |>
    write_models(
        mod1 = on(trestbps),
        mod3 = thalach ~ sex
    ) |>
    prepare(T_TEST, .mu = 120) |>
    conclude()


── 2 models · T-Test ─────────────────────────────────────────────────────────── 

mod1 : <cld_exec>
mod3 : <cld_exec>

Use display() to inspect individual results.

mod1 picks it up because its one-sample fn declares .mu; mod3’s two-sample fn doesn’t, so the same argument is simply ignored for that model rather than erroring — inject_and_run() only pulls from all_args what a given fn’s formals actually ask for, whether the model came from a batch or a single-model pipeline.

Two things need extra care once a batch mixes layouts like this, rather than repeating the same layout across models.

`state_null()` doesn’t span mixed layouts

Not every layout parses a stated hypothesis the same way. Recall from Question 2: x_by(thalach, sex)’s claim_parser reads group filters straight out of the claim (MU(thalach, sex == "Female")), on() has no such filter at all, and the <formula> layout has no claim_parser translation of its own yet. A single state_null() call attached to a write_models() batch would have to mean the same thing across all three shapes at once. Therefore, right now, state_null() isn’t wired up for <multi_lazy> objects.

In practice, a write_models() batch is for comparing layouts (or variables) against a test’s own default hypothesis, not for stating one numeric claim across a batch of differently-shaped models. If a specific hypothesis needs testing, stick to one define_model() |> prepare() |> state_null() |> conclude() pipeline per layout, the way Questions 1-3 do it above.

`via()` applies the same call to every model in the batch

via() on a <multi_lazy> doesn’t ask which model it’s talking to — the method name and any arguments you pass are applied to every model in the batch. That’s only safe when every layout present shares the same estimation-method vocabulary. on()’s one-sample t-test and x_by()’s two-sample t-test don’t: on() needs via("two_sample") to compare two vectors (see Question 2’s first tab), while x_by() is already two-sample by default and has no reason to register that variant at all. Mixing the two in one batch and calling via() once means the variant name has to be registered, with matching arguments, for every model type present:

heart |>
    write_models(
        mod1 = on(trestbps),
        mod2 = x_by(thalach, sex)
    ) |>
    prepare(T_TEST) |>
    via("two_sample") |>   # registered for on(); x_by() is two-sample already and may not have it
    conclude()

 [1m [33mError [39m in `method(via, list(statim::test_lazy, class_character))`: [22m
 [1m [22m [33m! [39m No variant  [34m"two_sample" [39m registered for model type  [34m"x_by" [39m.
 [36mℹ [39m Available variants:  [34m"contrast" [39m,  [34m"multi" [39m,  [34m"boot" [39m, and  [34m"permute" [39m.

The safer default is to keep a batch to layouts that share the same estimation method — either all left at base, or all recalibrated to a variant every model type in the batch actually registers — and fall back to separate single-model pipelines the moment the variants diverge.

Conclusion

statim makes sure the simplicity from base R and other packages like rstatix must exists, otherwise the steeper learning curve will get you. One-liner codes exist, the piped/grammar syntax is inherited to make sure the spirits of ggplot2 and dplyr exist on statim space.

What this vignette can’t show you is the fifth question: a custom contrast, a hypothesis that isn’t a straight equality, a test family that doesn’t have a tidy formula shorthand yet. That’s the actual bet statim is making: a consistent define_model() |> prepare() |> state_null() |> conclude() shape you can extend once, rather than learn a new argument convention per test. Four questions above are pretty straightforward and aren’t the right test of that bet. Judge it on the harder one. If you wanna know what’s beyond testing the equality in the null hypothesis, there’s a dedicated example for that.

A practical guide using real clinical data

Rationale

Sample Problem

Setup

Question 1: Is resting blood pressure elevated?

Codes

Verdict

Question 2: Does max heart rate differ by sex?

Codes

Verdict

Question 3: Does max heart rate fall with age?

Code

Verdict

Question 4: Is high fasting blood sugar unusually common in men?

Code

Verdict

Multiple Executions

state_null() doesn’t span mixed layouts

via() applies the same call to every model in the batch

Conclusion

`state_null()` doesn’t span mixed layouts

`via()` applies the same call to every model in the batch