API reference

The public API, in the groups the R package's reference index uses and in the same order, so a name is found in the same place on either site. Import every name here from the top-level package, for example from pilotr import simulate, power.

The R package covers more ground. Calibration, sweeps, the generated-analysis emitters and the no-code app have no Python counterpart yet, so the groups they occupy on that side do not appear below.

Design specifications

Load and validate the portable JSON specification that everything else reads. The specification is the same file in either language, and SPEC_VERSION is the version of that format this release understands.

pilotr.simulate.load_spec

load_spec(path, validate=True)

Load a JSON design specification from a file.

The specification is validated by default, because several ways of getting one wrong produce plausible data and no error at all: a mistyped coefficient key resolves to no column and so silently sets that effect to zero, and a response parameter left over from another family is ignored. Validation also refuses a specification declaring a spec_version newer than this implementation understands, and never reads such a file in part.

The file is read as R's jsonlite reads it, so that one file means one design in both twins: as UTF-8 with or without a byte-order mark, whatever the locale; with a one-element array read as its value where the schema has a single value, and [] read as an empty coefficients, slopes or correlations object, both as jsonlite::write_json() writes them by default. A file that repeats a key within one object is refused whatever validate says. See "Reading a specification" in spec/SPEC.md.

Parameters:

Name	Type	Description	Default
path	str or PathLike	Path to a JSON design-specification file, as a string or a pathlib.Path.	required
validate	bool	Whether to validate the specification after reading it. True (the default) applies validate_spec strictly, and False skips validation. Any other value validates leniently, with strict=False, so that an unknown field draws a warning, not an error. None, 0 and "lenient" all do so, as in the R twin.	True

Returns:

Type	Description
dict	The parsed specification, ready for simulate or power.

Raises:

Type	Description
ValueError	If the file repeats a key within one object, or, when validating, if the specification is invalid.

pilotr.validate.validate_spec

validate_spec(spec, strict: bool = True)

Validate a design specification, returning it so the call can be chained.

Checks the specification against the portable schema and against the cross-field rules the schema cannot express, and checks that its declared spec_version is one this implementation understands. Called by load_spec by default.

Names are checked together as well as one by one. Two columns with one name, a contrast column defined by two factors or named like another column, an interaction whose analysis column already exists, a level listed twice, a correlation that pairs a term with itself or gives one pair twice, a factor both between and within a unit, and a random.item entry in a design without items each validated and then changed the data in silence. The rules are set out under Names in the specification.

Several factors may be between one unit only in a specification declaring version 0.4 or later, which allocates them jointly, one block of units per combination of their levels. Up to version 0.3, pilotr assigned the levels of each between factor to blocks of its unit on its own, so two between factors over one unit fell into the same or overlapping blocks. A 2 x 2 between-subjects design over 40 subjects gave cells of 20, 0, 0 and 20, and its effects could not be estimated apart, so a specification declaring an earlier version is refused for them. Every version reads the same design written as one between factor whose levels are the cells, as the specification shows under Worked encodings.

From version 0.4, a within factor that sets counterbalance shows each subject each item at one level, rotating the levels over subjects and items as a Latin square does. It needs an item unit, and it is refused with per_subject, since every subject has to see every item. A warning says when the subjects or the items are not a multiple of the number of lists, which is the number of combinations of the counterbalanced factors' levels, so that the lists are unequal. The simulated data name each subject's list in a column list, so no other column of such a design may take that name.

Version negotiation covers the other direction. A specification that uses a feature introduced in 0.3 must declare 0.3 or later, and one that uses a feature introduced in 0.4 must declare 0.4. A specification declaring a version newer than SPEC_VERSION is refused outright, and one with no spec_version is treated as 0.2.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A parsed design specification, or the path to a JSON file holding one. A path is read as load_spec reads it and then validated, as the R twin's validate_spec() reads a path.	required
strict	bool	Whether an unrecognised field is an error (the default) or a warning. Pass False to load a specification carrying private annotations, accepting that a misspelled field will then be ignored in silence.	True

Returns:

Type	Description
dict	The specification as given, or as read from the file.

Raises:

Type	Description
ValueError	If the specification is invalid. All problems found are reported together. A file that repeats a key within one object is refused as load_spec refuses it.

Warns:

Type	Description
UserWarning	For an unrecognised field when strict is false, and in either mode for a within factor whose vary_within leaves out a unit of the design, which pilotr crosses with every unit anyway. From spec_version 0.4 the latter is an error. Also in either mode for counterbalanced factors whose lists are unequal, because the subjects or the items are not a multiple of their number.

pilotr.validate.SPEC_VERSION module-attribute

SPEC_VERSION = '0.4'

pilotr.examples.pilotr_example

pilotr_example(name: str | None = None) -> list[str] | str

List the bundled example specifications, or return the path to one.

pilotr ships one ready-to-run specification per design family, as JSON. These are the same files that drive the R twin and the no-code app.

Parameters:

Name	Type	Description	Default
name	str	The base name of an example, with or without the .json extension, for example "between_2group_gaussian". When None (the default), the available example names are returned in place of a path.	None

Returns:

Type	Description
list of str, or str	When name is None, the available example names. Otherwise, the path to that example's JSON file, ready to pass to :func:pilotr.simulate.load_spec.

Simulation

Draw a data set from a specification, returned together with the design information an analysis needs. model_data and model_formula give the analysis model the specification implies, the one power_mixed fits.

pilotr.simulate.simulate

simulate(spec, validate=True) -> Dataset

Simulate a data set from a design specification.

Build a linear predictor from the fixed effect sizes (categorical contrasts, continuous predictors and their interactions) plus the crossed by-subject and by-item random intercepts and slopes, then map it through the chosen response family.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A design specification, either an already-parsed dict or the path to a JSON spec file, as a string or a pathlib.Path. The format is documented at https://github.com/pablobernabeu/pilotr/blob/main/spec/SPEC.md.	required
validate	bool	Whether to validate the specification first. The default True catches the errors that would otherwise pass silently, such as a mistyped coefficient key, which resolves to no column and so sets that effect to zero. Validation costs a few milliseconds, so replicate loops validate once and then pass False. Any other value validates leniently, as in load_spec.	True

Returns:

Type	Description
Dataset	A table with one row per observation: a subject column, an optional item column, any grouping, factor, and continuous-predictor columns, and the response column named by spec["response"]["name"]. A design with a counterbalanced factor also has a list column after item, holding the list each subject follows, counted from 1.

pilotr.simulate.Dataset

A simulated data set: named columns and a list of row dicts.

Returned by simulate. Lightweight and dependency-free; build a pandas DataFrame from dataset.rows when you need one.

Attributes:

Name	Type	Description
columns	list of str	Column names, in order.
rows	list of dict	One dict per observation, keyed by column name.

column

column(name)

Return the values in column name as a list.

head

head(n=6)

Return the first n rows (default 6) as a list of dicts.

to_csv

to_csv(path)

Write the data set to path as CSV: a header row then one row per observation.

The file is UTF-8 whatever the locale, with LF line endings on every platform. A field holding a comma, a double quote or a newline is quoted, so a label such as "low, fast" stays one field. A float is written with repr, which reads back as the same double.

Parameters:

Name	Type	Description	Default
path	str or PathLike	The file to write, which is replaced if it exists.	required

pilotr.autoformula.model_data

model_data(spec, dataset)

The modelling data for a simulated data set, as the R twin's model_data() builds it.

Three kinds of column are added to the data, in this order. .y is the analysis response: log(y) for the lognormal family, log(y - shift) for shifted_lognormal and the response itself for every other family. Each factor's contrast columns follow, read from its labels. Every factor's level index is taken before any contrast column is written, since a contrast column may carry its own factor's name and then replaces the labels. Last comes one product column per interaction key among the coefficients and the random slopes, a_b for "a:b", the product of its parts taken from left to right.

A column that already exists is replaced where it stands, and a new one goes at the end, as in the R twin. The contrast and product columns equal the R twin's bit for bit, given the same data. .y passes through the C library's log(), which can differ in the last place between R and Python on Windows.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A design specification (dict or path to a JSON file). It is validated first.	required
dataset	Dataset or DataFrame	A data set simulated from spec, as simulate returns it, or a DataFrame with its columns. A DataFrame is copied, so the one given is left as it was.	required

Returns:

Type	Description
DataFrame	The data with the .y, contrast and product columns that model_formula reads.

Raises:

Type	Description
TypeError	If dataset is neither a Dataset nor a DataFrame.

Notes

Requires pandas, from the mixed extra.

pilotr.autoformula.model_formula

model_formula(spec, structure='maximal')

The lme4 formula a specification implies, as the R twin's model_formula() writes it.

The formula's response is .y, the column model_data adds. Each coefficient is a fixed term, an interaction "a:b" as its product column a_b. Each grouping factor has a random intercept and the random slopes the specification declares.

The three structures are the usual steps of a rule that simplifies a model when its fit is singular or fails to converge (Barr et al., 2013, doi:10.1016/j.jml.2012.11.001; Matuschek et al., 2017, doi:10.1016/j.jml.2017.01.001). "maximal" has every random slope the specification declares, correlated with the intercept where the specification correlates them, and uncorrelated (||) where it does not. "zcp", the zero-correlation-parameter model, keeps those slopes and writes || for every grouping factor with a slope. A grouping factor without slopes keeps (1 | g), since lme4 refuses (1 || g). "intercepts" keeps a random intercept for each grouping factor and drops every slope. The R package's power_mixed() takes these names as the rungs of such a rule, which this package does not simulate.

The string is the R twin's, character for character, which tools/parity/formula_cross.py checks on every shipped example and parity case. power_mixed fits the "maximal" model, with each random effect as an independent variance component.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A design specification (dict or path to a JSON file). It is validated first.	required
structure	('maximal', 'zcp', 'intercepts')	The random-effects structure (default "maximal").	"maximal"

Returns:

Type	Description
str	The formula, such as ".y ~ cond + (1 + cond | subject) + (1 + cond | item)".

Raises:

Type	Description
ValueError	If structure is not one of the three names, or if the specification has neither a coefficient nor a random effect, so that the formula would have no terms. The R twin refuses such a specification in the same words.

Power and design analysis

Estimate power by simulation, at one sample size or across a range of them, solve the resulting curve for the value that meets a target, and write the sample-size justification.

pilotr.power.power

power(spec, n_sims=1000, alpha=0.05, workers=1, alternative='two.sided')

Simulation-based power and design analysis for a two-group Gaussian design.

Repeatedly simulate from the specification, apply a two-sample t-test, and report power together with the Type S (sign) and Type M (magnitude) errors of Gelman and Carlin (2014), computed over the significant replicates.

The t-test treats every row as an independent observation, which is valid only when each subject contributes one row and no rows share a cluster. power therefore takes a Gaussian design with exactly one two-level factor between subjects and no item unit, within factor or grouping factor besides subject. Continuous predictors and by-subject random effects are allowed, since with one row per subject they vary independently from row to row. Any other design is refused with a message naming what makes its rows correlated. The R package's power_mixed() fits the model that such a design implies. Take 30 subjects crossed with 20 items, with by-subject and residual standard deviations of 1 and a by-item one of 0.3. With no true effect, a t-test of every row is significant in about half of all replicates.

Each level of the factor needs at least two subjects, so a design needs at least 4. The subjects are split into two blocks in level order, and the second level holds N // 2. With one subject at a level, the pooled variance rests on the other level alone, and with one at each it has no degree of freedom at all.

A directional prediction calls for a one-sided test, which alternative selects. The estimate is the mean of the second level minus that of the first, whose true value is true_effect. "greater" tests whether it is positive and "less" whether it is negative, by scipy.stats.ttest_ind with the same alternative. The one-sided p-value is half the two-sided one when the estimate lies in the tested direction. A one-sided test at alpha therefore rejects in exactly the replicates in which the two-sided test at 2 * alpha rejects with the estimate in that direction. Type S and Type M are taken over the replicates significant in the tested direction. When that is the direction of the true effect, Type S is 0 by construction, and Type M is the exaggeration among the replicates that would support the prediction.

Results from earlier versions were two-sided. A one-sided power at alpha, in the direction of the true effect, can be read off a two-sided result at 2 * alpha: it is power * (1 - type_s). Type M is recovered only when that result's Type S is 0, since the wrong-signed replicates cannot otherwise be told apart from the others.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A two-group Gaussian design specification with one row per subject (dict or path to a JSON file).	required
n_sims	int	Number of Monte Carlo replicates (default 1000).	1000
alpha	float	Significance level of the test (default 0.05), two-sided unless alternative says otherwise.	0.05
workers	int	Number of local worker processes over which to spread the replicates (default 1, serial). The replicate seeds are derived once from the specification's seed, so any worker count returns results identical to a serial run.	1
alternative	('two.sided', 'greater', 'less')	The test's alternative hypothesis (default "two.sided"), spelt as in the R twin. "greater" and "less" give a one-sided test of a directional prediction.	"two.sided"

Returns:

Type

Description

dict

Keys: n_sims, alpha, alternative, power, power_mcse (its Monte Carlo standard error), power_lo and power_hi (its 95% Wilson interval), n_significant, true_effect, mean_estimate, type_s, type_s_lo and type_s_hi (its 95% Wilson interval), type_m and type_m_mcse (its Monte Carlo standard error). Both design-analysis quantities are nan when no replicate reached significance and when the true effect is zero, as in a null condition: neither is defined without a true value to compare against, and Type M divides by it. Their interval and standard error are nan then too, and type_m_mcse is also nan when a single replicate reached significance, since a standard deviation needs two. A specification with no coefficient for the factor's contrast has a true effect of 0. Only the significant replicates inform Type S and Type M, and they are often few. Type S comes with a Wilson interval, since its standard error is zero whenever no sign is wrong. Type M comes with its standard error, the standard deviation of the exaggeration ratios over the square root of their number. Each is computed as in the R twin's power_design(), operation for operation, so the same counts give the same doubles.

Raises:

Type	Description
NotImplementedError	If the design is not Gaussian, does not have exactly one two-level between factor, or has rows that are correlated because of an item unit, a within factor or a grouping factor besides subject.
ValueError	If alternative is not "two.sided", "greater" or "less", or if a level of the between factor holds fewer than two subjects.

Notes

Requires scipy (imported lazily, in each worker process when parallel); install the power or dev extra.

pilotr.power.power_curve

power_curve(spec, subject_ns, n_sims=1000, alpha=0.05, workers=1, seeds='independent', alternative='two.sided')

Power curve over sample size for a two-group Gaussian design.

Sweep the number of subjects and compute power at each grid point.

Each grid point draws its own replicates. Point j of a sweep over J sizes takes replicate_seeds(seed, J)[j] as its seed, the rule of the R twin's sweep_spec(). Its replicates take their seeds from that one, so no two points share their random numbers. solve_curve fits the points as independent binomial counts, and its interval keeps its nominal coverage only when they are. The rule is prefix-stable: appending sizes leaves the earlier points' results as they were. A point's seed follows its position in subject_ns, though, so reordering the sizes, or inserting one, changes the results at every point that moves. seeds="shared" gives every point the specification's own seed, as up to 0.3.0. The points then share their noise. Over 1,000 emulated curves for each of twelve settings on a tight grid of sample sizes, the nominal 95% interval of solve_curve covered the true size in 87.8% of solves with shared seeds and in 96.5% with independent ones.

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A two-group Gaussian design specification with one row per subject, as power requires.	required
subject_ns	iterable of int	Subject counts to evaluate: any whole numbers, such as a range, a list or a NumPy array from np.arange or np.linspace.	required
n_sims	int	Replicates per grid point (default 1000).	1000
alpha	float	Significance level (default 0.05), two-sided unless alternative says otherwise.	0.05
workers	int	Number of local worker processes over which to spread the replicates at each grid point (default 1, serial). One process pool is started and reused across the whole sweep, and any worker count returns results identical to a serial run.	1
seeds	('independent', 'shared')	How the grid points are seeded (default "independent"). See above.	"independent"
alternative	('two.sided', 'greater', 'less')	The alternative hypothesis of the test at every grid point (default "two.sided"), as in power.	"two.sided"

Returns:

Type

Description

list of dict

One dict per grid point, with keys n_subject, power, power_mcse, power_lo, power_hi, type_s, type_s_lo, type_s_hi, type_m, type_m_mcse, n_sims and n_significant, each as power returns it. n_subject is a Python int, whatever type the grid held, so the records serialise as JSON. power_lo and power_hi give the band to draw around a plotted curve. The two counts are what solve_curve weights the curve by, since a power estimate over 1000 replicates should count for more than one over 50.

Raises:

Type	Description
ValueError	If seeds is neither "independent" nor "shared", or alternative is not "two.sided", "greater" or "less", or a size in subject_ns is below 4, which leaves a level of the between factor with fewer than two subjects (see power).

See Also

solve_curve, target_n : solve this curve for the sample size that meets a target power.

pilotr.power.power_mixed

power_mixed(spec, n_sims=50, alpha=0.05, workers=1, alternative='two.sided', focal=None)

Mixed-effects simulation-based power in Python, for the model the specification implies, via statsmodels MixedLM.

The model comes from the specification, as in the R twin: model_formula writes it and model_data builds its columns. Every coefficient is a fixed effect, an interaction "a:b" entering as its product column a_b. Every grouping factor that the specification gives random effects contributes the random intercept and the random slopes it declares, an extra grouping factor such as a site included, and nothing else. Each random effect is an independent variance component, 0 + C(g) for an intercept and 0 + C(g):term for a slope in statsmodels' notation, fitted over a single group so that crossed factors are fitted together. The response is analysed as log(y) for the lognormal family and as log(y - shift) for shifted_lognormal, whose coefficients are on the log scale. Every other family is analysed on its own scale. A replicate whose response rounds to 0, or to the shift, has no logarithm and returns no estimate, as in the R twin.

Up to 0.3, the model was the same whatever the specification declared, yv ~ cc with by-subject and by-item intercept and slope components, and it tested the first contrast of a single within factor crossed with items, refusing every other design. Any design that declares a random effect is now fitted, and one without random effects is refused, since MixedLM then has no variance component to fit. The R package's power_mixed() fits such a design with lm() or glm(), and power analyses the two-group Gaussian design.

Where the specification correlates a grouping factor's intercept and slopes, the independent components leave the correlations out, and one warning per call names the grouping factors concerned. The R package's power_mixed() estimates them, with | in its formula.

focal chooses the effects to test, as in the R twin. By default it is every coefficient, with the true value the specification gives it. A list of names leaves the true values unknown, which leaves Type S and Type M undefined, and a dict from names to values gives the true values to compare with, which may differ from those simulated. A name is a column of the model, so "a_b" for the interaction "a:b", or "(Intercept)". A name that the model lacks draws a warning in the R twin's words, and its results are empty. A coefficient the specification does not list is not in the model, so a null effect is planned with a coefficient of 0.

For bernoulli, poisson, ordinal and beta, the response's own scale is not the one the coefficients are written on, which is the logit scale, or the log scale for poisson. The estimate of a linear model is then a difference on the response scale, such as a difference in the probability of a correct response, and comparing it with the true coefficient says nothing about exaggeration. For these families, power_mixed warns and returns mean_estimate and type_m as nan for every effect. type_s compares signs alone, which the link preserves for a main effect. It does not preserve the sign of an interaction, since a bounded response compresses differences near the ends of its range, so the type_s of an interaction is nan too, as in the R twin. Power is the rate at which the linear model's test rejects. The R package's power_mixed() fits bernoulli and poisson responses with lme4::glmer() on their link scale, and its generate_design_analysis() writes a Bayesian design analysis for any family.

Every replicate is analysed with this one model, singular or not. A preregistered simplification rule, which moves to a simpler model whenever a fit is singular, has power and a Type I error rate of its own. It is simulated only in the R package, whose power_mixed() and precision_design() take a ladder of models as their formula.

This function has no joint test of several coefficients, such as the two contrasts of a three-level factor or the two coefficients of its interaction with a two-level factor. The power of such a multi-degree-of-freedom test differs from that of either coefficient. The R package's power_mixed() plans it when its focal names sets of coefficients, as in list(interaction = c("A1_B", "A2_B")).

This is pilotr's own simulation loop over the portable design specification, not a wrapper around an existing power package. It covers territory pioneered by simr (Green and MacLeod, 2016, doi:10.1111/2041-210x.12504) and mixedpower (Kumle, Vo and Draschkow, 2021, doi:10.3758/s13428-021-01546-0); pilotr differs in being driven by the portable cross-language spec, in reporting Type S and Type M errors and in built-in parallelisation via workers (default 1, serial), which returns results identical to a serial run for any worker count.

Each replicate is fitted by REML with Powell's method, and L-BFGS takes over if Powell's does not converge. statsmodels' default chain of optimisers stopped short of the REML optimum in most crossed fits, which earlier versions described as statsmodels overstating random-slope variance. Over the first 12 replicates of the crossed_mixed_rt example, the standard error of the effect now matches the one from lme4's uncorrelated model, (1 + cond || subject) + (1 + cond || item), to four decimal places in 11. The same 6 replicates are significant in both. A replicate whose Powell fit raises is fitted again with the default chain and counts as a fit with a warning. The random-effects design is held in dense matrices when it is small and in sparse ones from 100,000 entries, rows times random coefficients, which is quicker for a large design and reaches the same optimum.

The fits take time, and more of it as the model grows. On a Windows laptop, with NumPy held to one thread, a replicate of the 12 x 8 design of the Python guide took about 0.3 s and one of crossed_mixed_rt (30 subjects, 24 items) about 2 s. One of reading_time_continuous, with 50 subjects, 40 items, seven coefficients and four variance components, took about 7 s. workers spreads the replicates over processes.

The p-values are Wald z tests. A Wald z treats the estimate over its standard error as normal and so ignores the uncertainty in the variance components. With few subjects or items, it rejects more readily than a test with Satterthwaite's degrees of freedom (Luke, 2017, doi:10.3758/s13428-016-0809-y), which is the test of the R package's power_mixed(). In small designs, the power here therefore tends to exceed R's on the same data. On the 12 x 8 design of the documentation, 6 of 12 replicates are significant here and 4 of 12 in R. With 8 subjects, 6 items, random intercepts alone and no true effect, the test here rejected the null hypothesis in 0.070 of 300 replicates (Monte Carlo standard error 0.015). On the same data, the R package's Satterthwaite test of the same model rejected in 0.057 (0.013). tools/calibration/wald_null_rate.py measures both.

A directional prediction calls for a one-sided test, which alternative selects for every focal effect. "greater" tests whether the coefficient is positive and "less" whether it is negative. The one-sided p-value is the normal tail on the tested side of the same Wald z, so it is half the two-sided p-value when the estimate lies in the tested direction. A one-sided test at alpha therefore rejects in exactly the replicates in which the two-sided test at 2 * alpha rejects with the estimate in that direction. Type S and Type M are taken over the replicates significant in the tested direction, and backend names the direction. A two-sided result at 2 * alpha from an earlier version gives the one-sided power at alpha, in the direction of the true effect, as power * (1 - type_s).

Parameters:

Name	Type	Description	Default
spec	(dict, str or PathLike)	A design specification (dict or path to a JSON file) that declares at least one random effect.	required
n_sims	int	Number of Monte Carlo replicates (default 50, smaller than power's 1000 because each replicate fits a mixed model).	50
alpha	float	Significance level of the test (default 0.05), two-sided unless alternative says otherwise.	0.05
workers	int	Number of local worker processes over which to spread the replicates (default 1, serial). The replicate seeds are derived once from the specification's seed, so any worker count returns results identical to a serial run.	1
alternative	('two.sided', 'greater', 'less')	The alternative hypothesis of every focal effect's test (default "two.sided"), as in power.	"two.sided"
focal	None, str, list of str or dict	The effects to test (default None, every coefficient, with the specification's true values). Names alone leave the true values unknown, and a dict from names to numbers gives them.	None

Returns:

Type

Description

dict

Keys: backend (the estimator and its test), n_sims, alpha, the fit counts n_attempted, n_returned, n_converged, n_singular and n_warning, and effects, a dict with one entry per focal effect, in the order of focal. Each entry is a dict with keys power, power_mcse, power_lo, power_hi, n_significant, true_effect, mean_estimate, type_s, type_s_lo, type_s_hi, type_m and type_m_mcse, as in power. pandas.DataFrame.from_dict(result["effects"], orient="index") makes a table of them. A design with a single within factor whose first contrast is a focal effect also has that effect's twelve quantities at the top level, where earlier versions returned them. Those keys are deprecated and will be removed in the next release, and reading one raises a FutureWarning that names its place in effects. The fit counts mean what they mean in the R twin. n_attempted is n_sims, and n_returned counts the replicates that returned a fit. Each effect's power is the proportion of significant results among the fits that give it a finite estimate and p-value, not among n_sims. n_singular counts fits with a variance component whose standard deviation is below 1e-4 of the residual one, the test of lme4's isSingular(). n_warning counts fits whose optimisers did not converge, whose Powell fit raised or whose Hessian was not positive definite. A component at zero can fail statsmodels' check of the Hessian by itself, so that notice counts only against a fit that is not singular. n_converged counts the fits that are neither singular nor warned. statsmodels' notice that a fit may be on the boundary counts towards none of them, since it comes with every fit that has a variance component below 0.01. Singular fits and fits with a warning stay in power, as in the R twin. power_mcse, power_lo and power_hi are taken over the same fits. type_s and type_m are nan when no replicate reached significance and when the true effect is zero or unknown. mean_estimate and type_m are also nan for a bernoulli, poisson, ordinal or beta response, whose true effect is on a link scale, and so is the type_s of an interaction. The interval of type_s and the standard error of type_m are nan wherever their rate is, and type_m_mcse also when a single replicate reached significance.

Raises:

Type	Description
ValueError	If alternative is not "two.sided", "greater" or "less", if focal is none of its forms or names an effect twice, or if there is no effect to test, as for a specification without coefficients and no focal.
NotImplementedError	If the specification declares no random effect.

Warns:

Type	Description
UserWarning	Once per call, when the specification correlates the random effects of a grouping factor, which the independent components leave out.
UserWarning	Once per call, for a bernoulli, poisson, ordinal or beta response, saying that the linear model is fitted on the response's own scale and that the mean estimate and Type M are withheld. The R twin raises the same words for the families it fits in that way.
UserWarning	When a focal effect never appeared in a fitted model, in the R twin's words.
FutureWarning	When a deprecated top-level key of the result is read.

Notes

Requires the mixed extra (statsmodels, pandas and scipy, imported lazily in each worker process when parallel).

pilotr.solve.solve_curve

solve_curve(curve, target, x=None, y=None, n=None, effect=None, transform='sqrt', level=0.95)

Solve a simulated design curve for the value that meets a target.

Take the curve a sweep has already produced, fit the decision rate against the swept value, and solve for the value at which the rate meets a target. The solved value comes with a confidence interval, because a point read off a simulated curve without one repeats the overconfidence that design analysis exists to expose.

The input is the list of records power_curve returns, used as it stands. The swept value is taken from the leading column, which is where the curve functions put it, and the rate from power or p_meaningful, whichever the curve carries. Each rate is a proportion over a known number of replicates, and that count, read from n_returned, n_converged or n_sims, weights the fit: a rate over 200 replicates should count for more than a rate over 20.

The fit is a binomial regression with a probit link, and the solved value is the swept value at which the fitted rate equals target. The probit is chosen because power is a normal tail probability: under the normal approximation to a two-group comparison, the probit of power is linear in the square root of the sample size, so the model has the shape a design analysis already implies. Measured against R's stats::power.t.test across twelve combinations of effect size and target power on each of three grid shapes, at 400 replicates a point, the solved sample size fell within 2.6% of the analytic answer on average. A logit fitted the same way did as well on that draw, at 2.5%, and worse on another, at 3.5% against 2.2%. Curves this noisy do not separate the two links, and the probit is chosen for its shape.

The interval is the delta-method interval of R's MASS::dose.p, computed on the scale named by transform and mapped back, so it is symmetric on that scale and asymmetric on the natural one. That asymmetry is the honest shape: at the top of a power curve a given change in rate costs far more sample size than the same change lower down. Where the two-parameter model does not describe the curve, the interval is widened by the heterogeneity factor of probit analysis, Pearson's chi-square over its degrees of freedom (Finney, 1971), reported as dispersion. It is floored at 1, so a well-fitting curve is left alone and a badly-fitting one cannot report a narrower interval than its own residuals justify.

The default transform of "sqrt" suits a sample-size axis, where a rate rises with the square root of the sample size. Sweep something else, an effect size or a random-effect standard deviation, and "identity" is usually right.

Nothing here extrapolates. A curve whose rates do not straddle the target is refused, with the range it did cover reported, and so is a fit that solves outside the swept range. A curve whose fitted slope cannot be told from zero is refused too: the crossing is then compatible with any value at all, and an interval that said otherwise would be false.

A curve over subjects in a crossed design need not rise towards 1. When an effect varies from item to item, as a by-item random slope says it does, power at a fixed number of items levels off below 1 (Westfall, Kenny and Judd, 2014; Judd, Westfall and Kenny, 2017). Adding subjects cannot lift it past that ceiling. The probit model rises towards 1 and cannot follow such a plateau. A curve that levels off near the target can then be solved far from the true crossing. One that levels off just below the target never crosses it, but Monte Carlo noise can lift a simulated rate over the target, and the curve can then be solved as though it crossed. To compare designs that trade subjects for items, solve curves over the number of items, one at each of a few subject counts. Power over items levels off in the same way, at a ceiling set by the by-subject slope and the number of subjects. The subject counts must be large enough to put that ceiling well above the target.

What the interval covers is the Monte Carlo uncertainty of the fit, not the gap between the fitted shape and the true curve, and it treats the points of the curve as independent. power_curve now draws each point's replicates from a seed of its own. Over 1,000 simulated curves for each combination of four effect sizes and three target powers, at 400 replicates a point, the nominal 95% interval covered the analytic sample size in 96.5%, 95.9% and 96.2% of solves on the three grid shapes. For an effect size swept at a fixed sample size, it covered the analytic effect in 95.7% of solves. When every point reused the specification's seed, as up to 0.3.0, the figures were 87.8%, 94.0%, 93.5% and 70.3%. A curve assembled by hand, from separate runs that each started from the same seed, shares its noise in the same way. In the calibration's 36 checks, the solved size sat within 2.6% of the analytic answer on average and within 7.7% at worst. Nearly all of that error is the Monte Carlo noise the interval is describing, and it falls with the square root of the replicate count. Raise the count far enough and the interval narrows onto a fitted shape that is still slightly the wrong shape, so replicates alone do not make a solved size arbitrarily accurate. The remedies are more replicates, a finer grid or a design with a closed form to check against. A finer grid helps only while its points are independent. With shared seeds, its points share more of their noise, and the tight grid, whose points lie closest together, covered 87.8% where the coarse one covered 94.0%.

Parameters:

Name	Type	Description	Default
curve	sequence of dict, or data frame	A curve, as returned by power_curve: one record per swept value, holding the swept value, a decision rate, and the number of replicates behind it. A pandas data frame, or any table with a to_dict("records") method, is read as its records.	required
target	float	The decision rate to solve for, strictly between 0 and 1.	required
x	str	Name of the column holding the swept value. None, the default, takes the leading column.	None
y	str	Name of the column holding the decision rate. None, the default, takes power or p_meaningful, whichever is present.	None
n	(str, float or sequence)	The number of replicates behind each rate, either the name of a column or a numeric value, a NumPy number included. None, the default, takes n_returned, n_converged or n_sims, whichever is present.	None
effect	str	Which focal effect to solve for, when the curve holds more than one. Matched against the effect or param column. None, the default, uses every record, which is correct only when the curve holds one effect.	None
transform	str	The scale the swept value is fitted on: "sqrt" (the default, for a sample size), "identity" or "log".	'sqrt'
level	float	Confidence level for the reported interval (default 0.95).	0.95

Returns:

Type

Description

dict

Keys: value (the solved swept value), lo and hi (its confidence bounds), level, target, se (the delta-method standard error on the fitted scale, the scale on which the interval is symmetric), dispersion (the heterogeneity factor applied, 1 where the model fits), x and y (the columns used), transform, intercept and slope (the fitted coefficients), n_points (the number of curve points the fit used), and x_min and x_max (the swept range). A bound is allowed to fall outside that range. When one does, the sweep was too narrow to pin the value down and should be widened. A dispersion well above 1 says the curve is not the shape the model assumes, so the solve deserves a wider grid or more replicates before it deserves any trust.

Raises:

Type	Description
ValueError	If target or level is not strictly between 0 and 1, if the curve carries fewer than three usable points, if its rate does not vary, if its rates do not straddle target within the swept range, if the fitted slope cannot be told from zero, or if the solve lands outside the swept range.

References

Fieller, E. C. (1954). Some problems in interval estimation. Journal of the Royal Statistical Society: Series B, 16(2), 175-185. doi:10.1111/j.2517-6161.1954.tb00159.x

Finney, D. J. (1971). Probit analysis (3rd ed.). Cambridge University Press.

Judd, C. M., Westfall, J. and Kenny, D. A. (2017). Experiments with more than one random factor: Designs, analytic models, and statistical power. Annual Review of Psychology, 68, 601-625. doi:10.1146/annurev-psych-122414-033702

Westfall, J., Kenny, D. A. and Judd, C. M. (2014). Statistical power and optimal design in experiments in which samples of participants respond to samples of stimuli. Journal of Experimental Psychology: General, 143(5), 2020-2045. doi:10.1037/xge0000014

Examples:

>>> curve = [{"n_subject": 40, "power": 0.38, "n_sims": 100},
...          {"n_subject": 100, "power": 0.74, "n_sims": 100},
...          {"n_subject": 160, "power": 0.89, "n_sims": 100}]
>>> round(solve_curve(curve, target=0.8)["value"])
120

pilotr.solve.target_n

target_n(curve, target=0.8, **kwargs)

Solve a power curve for the sample size that reaches a target power.

The sample-size case of solve_curve, and the number a power analysis is usually run to obtain. Takes the curve a sweep over sample size has produced and returns the size at which power reaches target, rounded up to a whole number of units alongside the exact solution.

Everything solve_curve does applies here, including its refusals: a curve that never reaches the target within the sizes it swept is refused outright, and the reported interval can extend past the largest size simulated, which means the sweep was too narrow to settle the question.

The whole-number fields round up rather than to nearest, because a design cannot recruit a fraction of a subject and rounding down would leave the study short of the target it was sized for.

Parameters:

Name	Type	Description	Default
curve	sequence of dict, or data frame	A power curve, as returned by power_curve, or a pandas data frame of one.	required
target	float	The power to reach (default 0.8, the convention this package's plots draw a line at).	0.8
**kwargs		Further arguments passed to solve_curve, such as effect to pick one focal effect out of a curve holding several, or level for the interval.	{}

Returns:

Type	Description
dict	The dict solve_curve returns, with n, n_lo and n_hi added: value, lo and hi rounded up to whole numbers.

Examples:

>>> curve = [{"n_subject": 40, "power": 0.38, "n_sims": 100},
...          {"n_subject": 100, "power": 0.74, "n_sims": 100},
...          {"n_subject": 160, "power": 0.89, "n_sims": 100}]
>>> target_n(curve)["n"]
120

pilotr.report.power_report

power_report(x, solved=None, spec_file=None, justification=None, null_result=None)

Write a sample-size justification in Markdown from a power result.

Turns the result of power or power_mixed, with the target_n result that solved a curve for the sample size, into a sample-size justification, ready for a preregistration, a registered report or a grant application.

The sections follow the recommendations of Lakens (2022, Table 4, doi:10.1525/collabra.33267) for reporting an a-priori power analysis. The primary analysis comes first, with its model, test, alpha and sidedness, then the desired power and the sample size solved for it, with its interval. Each focal effect follows, with its true value, the justification given for that value and its power. Last comes the power when the null hypothesis is true. The simulation is then accounted for as Morris, White and Crowther (2019, doi:10.1002/sim.8086) recommend for simulation studies: the replicates simulated and those that returned an estimate, the fits that were singular or carried a warning, and the Monte Carlo standard error of every rate. The report ends with what a reader needs to reproduce the run: the software and its versions, the specification's seed and an MD5 checksum of the specification file's bytes. The versions are those loaded when the report is written, so write it in the session that ran the analysis.

Each rate's Monte Carlo standard error and 95% Wilson interval are computed in the report from the number of significant replicates and the number of replicates behind the rate, with the functions that compute them in the results. Type S has its Wilson interval from the number of wrong signs among the significant replicates. The R twin's power_report() writes the same text, byte for byte, from the same counts, which tools/parity/report_cross.py in the source repository checks.

The rates are those of the design x simulated. A sample size solved from a curve need not be that design's, so a justification usually pairs the solved size with a run of power_mixed at that size. spec_file names the specification of that run, and the report lists its units. A result carries neither its specification nor its seed, so without spec_file the report says that both are missing.

null_result is the result of the same analysis with every focal effect set to zero. Its rejection rate is the Type I error rate of the test, which for a Wald z with few subjects or items can exceed alpha.

The R twin also reports a precision_design() result, on the interval and the region of practical equivalence, which has no Python counterpart.

Parameters:

Name	Type	Description	Default
x	dict	The result of power or power_mixed at one design. A curve goes to target_n, and its result to solved. A power_mixed result is read from its effects where it has them.	required
solved	dict	The result of target_n, or of solve_curve for an axis other than sample size, which gives the desired power and the value solved for it.	None
spec_file	str or PathLike	The path of the specification file the analysis read. The report gives the MD5 checksum of its bytes, its seed and its units.	None
justification	str or dict	The justification of each focal effect's size: one string for every focal effect, or a dict from focal effect names to strings. An effect without one is reported as having none.	None
null_result	dict	The result of the same analysis with every focal effect set to zero. A result whose focal effects are not zero is refused.	None

Returns:

Type	Description
str	The Markdown, ending in a newline. The same arguments always give the same text.

Raises:

Type	Description
ValueError	If x or null_result is not a result at one design, if solved is not the result of target_n or solve_curve, if justification is neither a string nor a dict of strings naming focal effects, if null_result has a focal effect that is not zero, or if spec_file names no file.

Examples:

>>> from pilotr import power, power_report
>>> spec = {"name": "d", "seed": 1, "units": {"subject": {"n": 64}},
...         "factors": [{"name": "group", "levels": ["a", "b"],
...                      "contrasts": {"effect": [-0.5, 0.5]}, "between": "subject"}],
...         "fixed": {"intercept": 100, "coefficients": {"effect": 5}},
...         "response": {"family": "gaussian", "name": "score", "sigma": 10}}
>>> report = power_report(power(spec, n_sims=50),
...                       justification="Half a standard deviation.")
>>> print(report)

Reproducibility

The shared random-number stream and the numerical primitives that make a Python run and an R run agree, bit for bit wherever the two share a maths library.

pilotr.core.RNG

Shared cross-language random-number generator.

The combined linear congruential generator of L'Ecuyer (1988), implemented identically in R and Python so that a given seed yields the same stream in both. All intermediate products stay below 2**53, so the arithmetic is exact in IEEE-754 doubles. The draw-order contract is documented in the specification at https://github.com/pablobernabeu/pilotr/blob/main/spec/SPEC.md.

Parameters:

Name	Type	Description	Default
seed	int	Seed for the generator (coerced to a non-negative integer).	required

uniform

uniform() -> float

Return one draw from the standard uniform distribution on (0, 1).

normal

normal() -> float

Return one draw from the standard normal distribution.

normals

normals(k: int) -> list[float]

Return a list of k standard-normal draws.

pilotr.core.replicate_seeds

replicate_seeds(base, n: int) -> list[int]

The seeds pilotr's replicate loops give to their replicates, from a specification's seed.

Until 0.3 the rule was base + i, and consecutive seeds are not independent streams in this generator. Seeding is affine in the seed, setting s1 to 1 + (seed mod 2147483562) and s2 from s1, and each component recurrence is linear, so the k-th draws of consecutive seeds differ by very nearly the same amount mod 1. For the first draw that step is -0.0091: seeds 1000 to 1005 begin with 0.7027, 0.6937, 0.6846, 0.6755, 0.6665 and 0.6574. Draws in arithmetic progression mod 1 with step d correlate 1 - 6|d|(1 - |d|), so over 2,000 replicates the first draw of replicate i correlated 0.95 with the first draw of i + 1. A longer warm-up changes d and leaves the progression in place: with 50 warm-up draws in place of 10, the step is -0.32 and the correlation -0.30.

An arithmetic scramble does not fix that, because the seeding rule is affine in the seed: adding a Weyl increment and applying a Lehmer step left the first draws correlated at -0.27. Any affine rule for the seeds keeps the draws in progression, as Matsumoto et al. (2007) describe for the initialisation of other linear generators. Drawing the seeds from the shared generator does work, since successive outputs of the combined generator are what it exists to make look independent. Duplicates are skipped, so no two replicates are handed the same seed and silently produce identical data, and the skipping is deterministic, so this matches the R implementation exactly.

Parameters:

Name	Type	Description	Default
base	int	The specification's seed.	required
n	int	How many replicate seeds to return.	required

Returns:

Type	Description
list of int	n distinct seeds.

References

Matsumoto, M., Wada, I., Kuramoto, A., & Ashihara, H. (2007). Common defects in initialization of pseudorandom number generators. ACM Transactions on Modeling and Computer Simulation, 17(4), Article 15. https://doi.org/10.1145/1276927.1276928

pilotr.core.as241

as241(p: float) -> float

Inverse standard-normal CDF (quantile function), Wichura's (1988) Algorithm AS 241.

The PPND16 routine underlying R's qnorm, accurate to full double precision. The tail branch, which about 15% of uniform draws reach, calls log(). Where Python and R do not share a maths library, the two can therefore return different doubles for the same p. On Windows, 14 of 1,000,000 draws differed, by 1 to 3 units in the last place.

Parameters:

Name	Type	Description	Default
p	float	A probability in the open interval (0, 1).	required

Returns:

Type	Description
float	The standard-normal quantile for p.

pilotr.core.inv_logit

inv_logit(x: float) -> float

Logistic (inverse-logit) link, 1 / (1 + exp(-x)).

The two branches keep the exponent negative on either side of zero, so a large |x| underflows towards 0 or 1 where a single expression would overflow math.exp. The R twin, .inv_logit, branches at the same value, so both languages evaluate the same expression for a given x.

Parameters:

Name	Type	Description	Default
x	float	A value on the logit scale.	required

Returns:

Type	Description
float	The corresponding probability, in (0, 1).