Anytime-Valid OLS

Sequentially valid coefficient inference from the OLS summary interface

1 Why Anytime-Valid?

Ordinary linear regression summaries are fixed-sample summaries. They are calibrated for a sample size chosen before looking at the data. Anytime-valid inference is useful when the same analysis might be checked repeatedly as observations arrive, because its p-values and confidence intervals are calibrated for optional stopping.

crabbymetrics exposes this through the usual OLS.summary(...) method rather than through a separate model class:

out = model.summary(vcov="vanilla", anytime_valid=True, g=g_star)

The coefficient estimates are still the ordinary least-squares estimates. The extra fields in the summary are the anytime-valid p-values, confidence intervals, and omnibus F-test p-value.

2 Data-Generating Example

This mirrors the structure of the avlm README example: three centered covariates, a binary treatment, and treatment-by-covariate interactions. The coefficient on trt is a baseline treatment effect, while the interaction coefficients allow that effect to vary with the covariates.

Show code

import numpy as np
from html import escape
from IPython.display import HTML, display

import crabbymetrics as cm

np.set_printoptions(precision=4, suppress=True)


def html_table(headers, rows):
    parts = ["<table>", "<thead>", "<tr>"]
    parts.extend(f"<th>{escape(str(header))}</th>" for header in headers)
    parts.extend(["</tr>", "</thead>", "<tbody>"])
    for row in rows:
        parts.append("<tr>")
        parts.extend(f"<td>{cell}</td>" for cell in row)
        parts.append("</tr>")
    parts.extend(["</tbody>", "</table>"])
    return "".join(parts)


rng = np.random.default_rng(1)
n = 100
x = rng.normal(size=(n, 3))
x = x - x.mean(axis=0)
trt = rng.choice([0.0, 1.0], size=n)
noise = rng.normal(size=n)

y = (
    1.0
    + 1.4 * x[:, 2]
    + 2.3 * trt
    + 2.0 * x[:, 0] * trt
    + 3.0 * x[:, 1] * trt
    + noise
)

design = np.column_stack([x, trt, x * trt[:, None]])
coef_names = ["Intercept", "X1", "X2", "X3", "trt", "X1:trt", "X2:trt", "X3:trt"]

3 Fit OLS Once

The fitted model is just cm.OLS(). The optimal_g(...) helper chooses the mixture precision parameter that minimizes the anytime-valid confidence-interval width at this target sample size and significance level.

The second argument to optimal_g is the total coefficient count, including the intercept fitted internally by OLS.

Show code

model = cm.OLS()
model.fit(design, y)

alpha = 0.05
g_star = cm.optimal_g(n, design.shape[1] + 1, alpha)

plain = model.summary(vcov="vanilla")
av = model.summary(vcov="vanilla", anytime_valid=True, g=g_star, level=1 - alpha)

print("g_star:", round(g_star, 4))

g_star: 11.5815

The usual summary fields remain available. For example, the first three slope coefficients are unchanged when anytime-valid inference is requested:

Show code

print("plain slope estimates:", np.round(plain["coef"][:3], 4))
print("AV slope estimates:   ", np.round(av["coef"][:3], 4))

plain slope estimates: [0.0257 0.1183 1.4786]
AV slope estimates:    [0.0257 0.1183 1.4786]

4 Anytime-Valid Coefficient Table

The anytime-valid summary adds full-parameter fields:

estimate: intercept followed by slopes
std_error: matching standard errors
t_value: estimate divided by standard error
p_value: anytime-valid coefficient p-values
confint: anytime-valid confidence intervals
f_statistic and f_p_value: omnibus test for all non-intercept coefficients

Show code

rows = []
for i, name in enumerate(coef_names):
    rows.append(
        [
            name,
            f"{av['estimate'][i]:.4f}",
            f"{av['std_error'][i]:.4f}",
            f"{av['p_value'][i]:.4g}",
            f"[{av['confint'][i, 0]:.4f}, {av['confint'][i, 1]:.4f}]",
        ]
    )

display(
    HTML(
        html_table(
            ["Term", "Estimate", "Std. Error", "AV p-value", "95% AV interval"],
            rows,
        )
    )
)

Term	Estimate	Std. Error	AV p-value	95% AV interval
Intercept	1.1679	0.1193	1.432e-12	[0.7975, 1.5383]
X1	0.0257	0.1391	1	[-0.4062, 0.4575]
X2	0.1183	0.1520	1	[-0.3535, 0.5900]
X3	1.4786	0.1390	5.278e-14	[1.0472, 1.9100]
trt	2.0934	0.1843	3.403e-15	[1.5215, 2.6653]
X1:trt	1.8064	0.1937	8.725e-12	[1.2053, 2.4075]
X2:trt	2.9676	0.2098	1.892e-19	[2.3163, 3.6189]
X3:trt	0.1927	0.2082	1	[-0.4536, 0.8390]

The omnibus anytime-valid F-test is reported in the same dictionary:

Show code

print("AV F-statistic:", round(av["f_statistic"], 4))
print("AV F-test p-value:", f"{av['f_p_value']:.4g}")

AV F-statistic: 131.9361
AV F-test p-value: 1.467e-33

5 Robust Variant

The anytime-valid option composes with the ordinary summary(vcov=...) covariance interface. Here is the same fit with HC1 standard errors:

Show code

av_hc1 = model.summary(vcov="hc1", anytime_valid=True, g=g_star, level=0.95)

rows = []
for i, name in enumerate(coef_names):
    rows.append(
        [
            name,
            f"{av_hc1['estimate'][i]:.4f}",
            f"{av_hc1['std_error'][i]:.4f}",
            f"{av_hc1['p_value'][i]:.4g}",
            f"[{av_hc1['confint'][i, 0]:.4f}, {av_hc1['confint'][i, 1]:.4f}]",
        ]
    )

display(
    HTML(
        html_table(
            ["Term", "Estimate", "HC1 Std. Error", "AV p-value", "95% AV interval"],
            rows,
        )
    )
)

Term	Estimate	HC1 Std. Error	AV p-value	95% AV interval
Intercept	1.1679	0.1148	3.133e-13	[0.8117, 1.5241]
X1	0.0257	0.1077	1	[-0.3085, 0.3599]
X2	0.1183	0.1214	1	[-0.2584, 0.4950]
X3	1.4786	0.1483	7.033e-13	[1.0182, 1.9390]
trt	2.0934	0.1879	7.739e-15	[1.5103, 2.6766]
X1:trt	1.8064	0.2120	2.146e-10	[1.1484, 2.4643]
X2:trt	2.9676	0.2080	1.245e-19	[2.3221, 3.6131]
X3:trt	0.1927	0.2094	1	[-0.4573, 0.8428]

6 Compact Helper

cm.av(...) is a small convenience wrapper around the same OLS.summary(...) call. It returns the summary dictionary directly, not a new model object.

Show code

helper = cm.av(model, g=g_star, vcov="vanilla")
print(np.allclose(helper["p_value"], av["p_value"]))
print(np.allclose(helper["confint"], av["confint"]))

True
True