For a selected pair of methods \((a, b)\), the backend first forms complete
within-subject pairs at matched subject and integer-coerced
time. Let
$$
d_{it} = y_{itb} - y_{ita},
\qquad
m_{it} = \frac{y_{ita} + y_{itb}}{2},
$$
where \(d_{it}\) is the paired difference and \(m_{it}\) is the paired
mean for subject \(i\) at time/replicate \(t\). Only complete
subject-time matches contribute to that pairwise fit.
If multiple rows are present for the same subject-time-method
combination within an analysed pair, the backend keeps the last encountered
value for that combination when forming the pair. The function therefore
implicitly assumes at most one observation per subject-time-method
cell for each analysed contrast.
The fitted model for each analysed pair is
$$
d_{it} = \beta_0 + \beta_1 x_{it} + u_i + \varepsilon_{it},
$$
where \(x_{it} = m_{it}\) if include_slope = TRUE and the term is omitted
otherwise; \(u_i \sim \mathcal{N}(0, \sigma_u^2)\) is a subject-specific
random intercept; and the within-subject residual vector satisfies
\(\mathrm{Cov}(\varepsilon_i) = \sigma_e^2 R_i\).
When use_ar1 = FALSE, \(R_i = I\). When use_ar1 = TRUE, the backend
works with the residual precision matrix \(C_i = R_i^{-1}\) over
contiguous time blocks within subject and uses \(\sigma_e^2 C_i^{-1}\) as
the residual covariance.
AR(1) residual structure
Within each subject, paired observations are ordered by integer-coerced
time. AR(1) correlation is applied only over strictly contiguous runs
satisfying \(t_{k+1} = t_k + 1\). Gaps break the run. Negative times, and
any isolated positions not belonging to a contiguous run, are treated as
independent singletons.
For a contiguous run of length \(L\) and correlation parameter \(\rho\),
the block precision matrix is
$$
C
=
\frac{1}{1-\rho^2}
\begin{bmatrix}
1 & -\rho & & & \\
-\rho & 1+\rho^2 & -\rho & & \\
& \ddots & \ddots & \ddots & \\
& & -\rho & 1+\rho^2 & -\rho \\
& & & -\rho & 1
\end{bmatrix},
$$
with a very small ridge added to the diagonal for numerical stability.
If use_ar1 = TRUE and ar1_rho is supplied, that value is used after
validation and clipping to the admissible numerical range handled by the
backend.
If use_ar1 = TRUE and ar1_rho = NA, the backend estimates rho
separately for each analysed pair by:
fitting the corresponding iid model;
computing a moments-based lag-1 estimate from detrended residuals
within contiguous blocks, used only as a seed; and
refining that seed by a short profile search over rho using the
profiled REML log-likelihood.
In the exported ba_rm() wrapper, if an AR(1) fit for a given analysed pair
fails specifically because the backend EM/GLS routine did not converge to
admissible finite variance-component estimates, the wrapper retries that pair
with iid residuals. If the iid refit succeeds, the final reported residual
model for that pair is "iid" and a warning is issued. Other AR(1) failures
are not simplified and are propagated as errors.
Internal centring and scaling for the proportional-bias slope
When include_slope = TRUE, the paired mean regressor is centred and scaled
internally before fitting. Let \(\bar m\) be the mean of the observed paired
means. The backend chooses a scaling denominator from:
It uses the first of these that is not judged near-zero relative to the
largest finite positive candidate scale, under a threshold proportional to
\(\sqrt{\epsilon_{\mathrm{mach}}}\). If all candidate scales are treated as
near-zero, the fit stops with an error because the proportional-bias slope is
not estimable on the observed paired-mean scale.
The returned beta_slope is back-transformed to the original paired-mean
scale. The returned BA centre is the fitted mean difference at the centred
reference paired mean \(\bar m\), not the original-scale intercept
coefficient.
Estimation
The backend uses a stabilised EM/GLS scheme.
Conditional on current variance components, the fixed effects are updated by
GLS using the marginal precision of the paired differences after integrating
out the random subject intercept. The resulting fixed-effect covariance used
in the confidence-interval calculations is the GLS covariance
$$
\mathrm{Var}(\hat\beta \mid \hat\theta)
=
\left( \sum_i X_i^\top V_i^{-1} X_i \right)^{-1}.
$$
Given updated fixed effects, the variance components are refreshed by EM using
the conditional moments of the subject random intercept and the residual
quadratic forms. Variance updates are ratio-damped and clipped to admissible
ranges for numerical stability.
Reported BA centre and limits of agreement
The reported BA centre is always model-based.
When include_slope = FALSE, it is the fitted intercept of the paired-
difference mixed model.
When include_slope = TRUE, it is the fitted mean difference at the centred
reference paired mean used internally by the backend.
The reported limits of agreement are
$$
\mu_0 \pm \texttt{loa\_multiplier}\sqrt{\sigma_u^2 + \sigma_e^2},
$$
where \(\mu_0\) is the reported model-based BA centre. These LoA are for a
single new paired difference from a random subject under the fitted model.
Under the implemented parameterisation, AR(1) correlation affects the
off-diagonal within-subject covariance structure and therefore the estimation
of the model parameters and their uncertainty, but not the marginal variance
of a single paired difference. Consequently rho does not appear explicitly
in the LoA point-estimate formula.
Confidence intervals
The backend returns Wald confidence intervals for the reported BA centre and
for both LoA endpoints.
These intervals combine:
the conditional GLS uncertainty in the fixed effects at the fitted
covariance parameters; and
a delta-method propagation of covariance-parameter uncertainty from
the observed information matrix of the profiled REML log-likelihood.
The covariance-parameter vector is profiled on transformed scales:
log-variances for \(\sigma_u^2\) and \(\sigma_e^2\), and, when rho is
estimated internally under AR(1), a transformed correlation parameter
mapped back by \(\rho = 0.95\tanh(z)\).
Numerical central finite differences are used to approximate both the
observed Hessian of the profiled REML log-likelihood and the gradients of the
reported derived quantities. The resulting variances are combined and the
final intervals are formed with the normal quantile corresponding to
conf_level.
Exactly two methods versus \(\ge 3\) methods
With exactly two methods, at least two complete subject-time pairs are
required; otherwise the function errors.
With \(\ge 3\) methods, the function analyses every unordered pair of method
levels. For a given pair with fewer than two complete subject-time matches,
that contrast is skipped and the corresponding matrix entries remain NA.
For a fitted contrast between methods in matrix positions \((j, k)\) with
\(j < k\), the stored orientation is:
$$
\texttt{bias}[j,k] \approx \mathrm{method}_k - \mathrm{method}_j.
$$
Hence the transposed entry changes sign, while sd_loa and width are
symmetric.
Identifiability and safeguards
Separate estimation of the residual and subject-level variance components
requires sufficient complete within-subject replication after pairing. If the
paired data are not adequate to separate these components, the fit stops with
an identifiability error.
If the model is conceptually estimable but no finite positive pooled
within-subject variance can be formed during initialisation, the backend uses
\(0.5 \times v_{\mathrm{ref}}\) only as a temporary positive starting value
for the EM routine and records a warning string in the backend output. The
exported wrapper does not otherwise modify the final estimates.
If the EM/GLS routine fails to reach admissible finite variance-component
estimates, the backend throws an explicit convergence error rather than
returning fallback estimates.