Details of Computation

This page gathers the computational details behind the reported quantities: the confidence interval for the geometric mean ratio, the within-subject variability estimates that feed the reference-scaling criteria, and the ANOVA p-values. For the model definitions themselves and the arguments that select them, see Models.

Confidence intervals

The confidence interval for the GMR is the primary inferential output of the analysis: in the standard case, bioequivalence is concluded from whether it lies within the acceptance limits (see Regulatory Criteria).

The log-GMR is the formulation contrast δ̂ with standard error SE; the reported interval is the TOST (1 − 2α) interval

\[\hat{\delta} \pm t_{\nu,\,1-\alpha}\, \mathrm{SE},\]

with α = (1 − level_tost)/2, exponentiated to the ratio scale (so the default level_tost = 0.9 gives α = 0.05 and a 90% interval for the GMR). The denominator degrees of freedom ν depend on the model:

• Parallel / paired (t-tests): the t-test degrees of freedom, which are nₜ + nᵣ − 2 for the equal-variance two-sample test (Welch–Satterthwaite when homogeneity = false), and n − 1 for the paired test.
• Fixed-effects linear model: the model's residual degrees of freedom.
• Mixed models: the Satterthwaite approximation (detailed under p-values).

The nonparametric interval

The nonparametric paradigm does not use the interval above. It fits no model, so there is no contrast δ̂, no standard error, and no degrees-of-freedom count.

The idea. The parametric interval starts from an estimate and adds a margin to it, the margin built from a fitted standard error and an assumed sampling distribution. A rank interval is assembled the other way around, out of the data themselves. Write down every pairwise comparison the design supports, sort them, and the interval is a pair of order statistics from that sorted list: count in Cα − 1 places from each end. The point estimate is the middle of the same list.

This is a confidence interval because a null value far enough from the data in either direction makes the rank test reject, and the interval is exactly the set of null values it does not reject. Since the test reads only ranks, how far to count in depends on the sample size and not on the endpoint's distribution, which is what makes the interval distribution-free. That property holds on both routes. The normal approximation approximates a null distribution that is itself free of the endpoint's distribution, so what it costs is accuracy of the level, not distribution-freeness; ties are what erodes the property, and they do so on the exact and approximate routes alike. Both points are taken up under Exact or approximate below.

The pairwise set. There are three of these, one per route, rather than one per procedure: the parallel and cross-sequence routes share the rank sum machinery but hand it different inputs. Throughout, values are on the log scale. Write Dᵢ for subject i's within-subject difference, test minus reference, and Xᵢᵀ, Xⱼᴿ for a subject's value under one formulation on the parallel route. Write ñ for the number of non-zero differences, nₜ and nᵣ for the sizes of the two formulation groups, and n₁ and n₂ for the numbers of subjects in each of the two sequences supplying a difference.

The last row needs explanation. The construction is defined on Walsh averages, as the signed rank route is, but restricted to the pairs that draw one subject from each sequence, since those are the pairs in which the period contrast cancels. It is reached by handing the rank sum machinery one sequence's differences against the negated differences of the other, which produces the sums Dᵢ + Dⱼ rather than the averages. The estimate and both bounds are therefore halved on the way out, which converts the sums into the averages the construction is defined on.

Finding Cα. The point estimate is the median of the set, which is the Hodges–Lehmann estimate. The interval runs between order statistics of the same set, from the Cα-th smallest to the Cα-th largest, which covers 1 − 2·P(W ≤ Cα − 1) for W the procedure's rank statistic under the null. So Cα is the only quantity left to determine, and it is taken as the largest index whose coverage still reaches the requested level. It comes either from the exact null distribution or, under the normal approximation, from

\[C_\alpha = \left\lceil \mu - z_{1-\alpha}\,\sigma - \tfrac{1}{2} \right\rceil,\]

where α = (1 − level_tost)/2 as above, so z₁₋α is the same normal quantile the parametric interval would use, μ and σ are the null mean and standard deviation of W, σ carries the tie correction, and the ½ is a continuity correction. Textbook recipes for this interval often omit both corrections, and a Cα computed that way can differ from this package's by one, which moves each bound to a neighboring order statistic; the two agree whenever those neighboring values coincide, as they often do on a tied endpoint. Both endpoints are then exponentiated to the ratio scale, as in the parametric case, after the halving on the cross-sequence route.

Two consequences with no counterpart above. Because Cα is an integer, only finitely many coverages are attainable, so the achieved level is the nearest attainable one at or above the nominal level: rank intervals are conservative, and increasingly so as the sample shrinks or ties multiply. It also bounds what a sample can deliver: the widest interval the construction admits runs between the two extremes of the set, and its coverage, 1 − 2/2^ñ for the signed rank and 1 − 2/C(n₁+n₂, n₁) for a rank sum, is the most any interval of that construction can carry. A requested level above it is refused with an ArgumentError naming the coverage the sample could reach, rather than the widest interval being reported under a label it does not meet. At the default level that needs at least five non-zero differences for the signed rank, and at least twenty pairs for a rank sum. That ceiling is the exact null's, and the exact test reaches it. The normal approximation, which a tied sample takes past the size-only range, reaches a little less: its Cα is ⌈m/2 − z·σ − ½⌉ for m pairwise estimates and the tie-corrected σ, so it admits an interval only while the level is below 2Φ((m/2 − ½)/σ) − 1. A request between the two bounds is refused in the same way, the message naming the approximation. And the interval is not symmetric about its point estimate, so the Hodges–Lehmann estimate is not the midpoint of the interval printed beside it.

Zeros. On the signed rank route, zero differences are discarded before any of this, as the signed rank statistic itself discards them, so the estimate, the interval, and the test all describe the same sample. That is why the count above is ñ rather than the number of subjects. The two rank sum routes discard nothing, having no notion of a zero to drop.

Exact or approximate

Whether Cα comes from the exact null distribution or from the normal approximation is not a choice the caller makes. Each procedure decides it from the sample it is given, by size first and ties second, and the rules are part of the specification rather than an implementation detail, so they are stated here in full. They are this package's rule and not a universal one: other implementations switch between exact and approximate at different sizes, and some ignore ties in the decision, so a sample can be exact under one rule and approximate under another.

Write a = Σₜ (t³ − t) for the tie adjustment, summed over the counts t of tied values at each tied position; a = 0 exactly when no two ranked values are equal. Note that each procedure ranks a different thing, so each counts ties in a different place.

Signed rank. Let ñ be the number of non-zero within-subject differences, the zeros having already been discarded, and let a be the tie adjustment over the absolute values of those ñ differences. The exact procedure is used when

\[\tilde n \le 15 \quad\text{or}\quad (\tilde n \le 50 \text{ and } a = 0),\]

and the normal approximation otherwise. Two differences equal in magnitude but opposite in sign are a tie for this purpose, since it is the absolute values that are ranked.

Rank sum, on the parallel route and on the cross-sequence route alike. Let N be the combined size of the two groups handed to the procedure: nₜ + nᵣ on the parallel route, and n₁ + n₂ on the cross-sequence route. Let a be the tie adjustment over those N pooled values. The exact procedure is used when

\[N \le 10 \quad\text{or}\quad (N \le 50 \text{ and } a = 0),\]

and the normal approximation otherwise. Note that on the cross-sequence route N counts the subjects supplying a within-subject difference, not the n₁ · n₂ pairs formed from them.

Three consequences follow. A tied endpoint such as Tmax does not force the approximation: a small enough sample lands in the size-only range, where ties are not consulted at all. Between the two thresholds, ties are what decides, and a single tie is enough to send a sample of any size in that band to the approximation. And (exact) on the Model line describes the null distribution, not the coverage: on a tied sample the exact interval is still built by inverting the tie-free distribution, so its achieved coverage remains approximate.

Within-subject variability

The within-subject standard deviation σw of a formulation, on the log scale, is a key quantity estimated in almost every analysis: it is what the reported coefficient of variation (CV) summarizes, and, on the reference-scaling criteria (the FDA RSABE and EMA ABEL), it is additionally what scales the acceptance limits. The coefficient of variation and the log-scale standard deviation convert via geocv2sigma and sigma2geocv:

\[\sigma_w = \sqrt{\ln\!\left(1 + \mathrm{CV}^{2}\right)}, \qquad \mathrm{CV} = \sqrt{e^{\sigma_w^{2}} - 1},\]

where CV is a proportion (the summary reports it as a percentage, 100·CV).

A CV is reported in almost every analysis, but how it is obtained (and whether it is genuinely within-subject) depends on the design:

So most analyses read the CV straight off the model already fitted for the formulation comparison. Only the cases that need per-formulation within-subject variances (replicated crossover designs and multi-period parallel designs) call on a separate auxiliary model, which is also what reference scaling requires.

The auxiliary model. This is a secondary linear regression fitted only to recover σw, distinct from the formulation-comparison model that produces the GMR and from that model's own variance components. Fitted over the subjects with repeated observations of a formulation, it regresses the log endpoint on design terms (subject, period, plus formulation and sequence depending on the variant), and takes its residual standard deviation as σw; this yields the reference-scaling quantities σ̂wᵣ (and σ̂wₜ) used by the criteria. Both pumas_be and run_be fit it automatically when it is needed, so you do not normally call the functions below yourself. They are exposed for inspection.

• The within_subject_variability function returns σw and its degrees of freedom k from the auxiliary linear model. On a partially replicated design a formulation that is not replicated has no estimable σw; passing allownonreplicated = true returns missing for it instead of raising an error.
• The within_subject_variability_ratio function returns the ratio σ̂wₜ/σ̂wᵣ and its upper confidence bound (used by the FDA NTID criterion).

This auxiliary model is one of four WSVModel variants, selected automatically but overridable through the wsv_model argument of run_be, as well as of within_subject_variability directly; pumas_be does not accept the argument, as the criterion determines the choice there (EMA_HighlyVariable sets additive_sequence_wsv_model). All four share period + id and differ in how formulation/sequence enter:

p-values

The ANOVA p-values are a diagnostic accompaniment to the confidence interval, not part of the bioequivalence decision. The parallel and paired t-test paradigms report no ANOVA table. For the linear and mixed crossover models, each term (formulation, sequence, and period) is tested with a Type III Wald F-test on the subvector β̂_R of its coefficients, with covariance block V_R extracted from V = vcov(model):

\[F = \tfrac{1}{q}\, \hat{\beta}_R^{\top} V_R^{-1} \hat{\beta}_R, \qquad p = \Pr\!\left(F_{q,\nu} > F\right),\]

where q = |R| is the numerator degrees of freedom of that term. The numerator statistic is computed identically for the fixed-effects and mixed models (the same β̂_R and V_R = vcov(model) form); the two differ only in the denominator degrees of freedom ν.

Denominator degrees of freedom ν, by model.

• Fixed-effects linear model. Here ν is the residual degrees of freedom: a single value, shared by the GMR confidence interval and all three ANOVA F-tests.
• Mixed models. Here ν is the Satterthwaite approximation with a small-sample, observed-information adjustment, and is generally non-integer. It is computed separately for each contrast, so it is not a single shared value: the GMR confidence interval and the formulation F-test use the same single-coefficient ν (the mixed model is only used for two-formulation designs, so q = 1 for formulation), while the sequence and period F-tests each use their own Satterthwaite ν.

The Degrees of Freedom row of the per-formulation block reports this same ν: the residual degrees of freedom for the linear model, the formulation-contrast Satterthwaite value for the mixed models. It is therefore equal to the denominator of the formulation F-test, but not in general to those of the sequence and period tests. Note that the row is taken from the confidence-interval computation rather than from stat_tests_pvalues; the two agree because both rest on the same single contrast.

Under the FDA criteria a second Degrees of Freedom row appears, in the Reference Scaling Analysis block. That one is unrelated to the above: it is the degrees of freedom d of σ̂wᵣ entering Howe's bound (see Within-subject variability), not a denominator for any test reported here.