Package 'matrixCorr'

Summary Accessor for Correlation Summaries

Description

Summary Accessor for Correlation Summaries

Usage

## S3 method for class 'summary.corr_result'
x[[i, ...]]

Arguments

x	A summary.corr_result object.
i	A column name or summary metadata key.
...	Unused.

Value

A summary column (if present) or summary metadata entry.

---

Summary Accessor for Correlation Summaries

Description

Summary Accessor for Correlation Summaries

Usage

## S3 method for class 'summary.corr_result'
x$name

Arguments

x	A summary.corr_result object.
name	A column name or summary metadata key.

Value

A summary column (if present) or summary metadata entry.

---

Edge-list data-frame view

Description

Edge-list data-frame view

Usage

## S3 method for class 'corr_edge_list'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

x	A corr_edge_list object.
row.names	Ignored.
optional	Ignored.
...	Unused.

Value

A data frame with columns row, col, value.

---

Bland-Altman statistics with confidence intervals

Description

Computes Bland-Altman mean difference and limits of agreement (LoA) between two numeric measurement vectors, including t-based confidence intervals for the mean difference and for each LoA using 'C++' backend. If group2 is omitted and group1 is a numeric matrix or data frame with at least two numeric columns, ba() computes all pairwise contrasts across methods and returns a pairwise Bland-Altman matrix object.

Note: Lin's concordance correlation coefficient (CCC) is a complementary, single-number summary of agreement (precision + accuracy). It is useful for quick screening or reporting an overall CI, but may miss systematic or magnitude-dependent bias; consider reporting CCC alongside Bland-Altman.

Usage

ba(
  group1,
  group2,
  loa_multiplier = 1.96,
  mode = 1L,
  conf_level = 0.95,
  n_threads = getOption("matrixCorr.threads", 1L),
  verbose = FALSE
)

## S3 method for class 'ba'
print(
  x,
  digits = 3,
  ci_digits = 3,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ba'
summary(object, digits = 3, ci_digits = 3, ...)

## S3 method for class 'summary.ba'
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ba'
plot(
  x,
  title = "Bland-Altman Plot",
  subtitle = NULL,
  point_alpha = 0.7,
  point_size = 2.2,
  line_size = 0.8,
  shade_ci = TRUE,
  shade_alpha = 0.08,
  smoother = c("none", "loess", "lm"),
  symmetrize_y = TRUE,
  show_value = TRUE,
  ...
)

## S3 method for class 'ba_matrix'
print(
  x,
  digits = 3,
  ci_digits = 3,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  style = c("pairs", "matrices"),
  ...
)

## S3 method for class 'ba_matrix'
summary(object, digits = 3, ci_digits = 3, ...)

## S3 method for class 'summary.ba_matrix'
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ba_matrix'
plot(
  x,
  pairs = NULL,
  against = NULL,
  facet_scales = c("free_y", "fixed"),
  title = "Bland-Altman (pairwise)",
  point_alpha = 0.6,
  point_size = 1.8,
  line_size = 0.7,
  shade_ci = TRUE,
  shade_alpha = 0.08,
  smoother = c("none", "loess", "lm"),
  show_value = TRUE,
  ...
)

Arguments

group1	Numeric vector of paired measurements, or a numeric matrix/data frame with at least two numeric columns when group2 is omitted.
group2	Optional numeric vector of paired measurements. If omitted, group1 must be a numeric matrix or data frame and all unordered numeric column pairs are analysed.
loa_multiplier	Positive scalar; the multiple of the standard deviation used to define the LoA (default 1.96 for nominal 95\ intervals always use $t_{n-1,\,1-\alpha/2}$ regardless of this choice.
mode	Integer; 1 uses group1 - group2, 2 uses group2 - group1.
conf_level	Confidence level for CIs (default 0.95).
n_threads	Integer $\geq 1$. Number of OpenMP threads. Defaults to getOption("matrixCorr.threads", 1L).
verbose	Logical; if TRUE, prints how many OpenMP threads are used.
x	A "ba_matrix" object.
digits	Number of digits for estimates (default 3).
ci_digits	Number of digits for CI bounds (default 3).
n	Optional row threshold for compact preview output.
topn	Optional number of leading/trailing rows to show when truncated.
max_vars	Optional maximum number of visible columns; NULL derives this from console width.
width	Optional display width; defaults to getOption("width").
show_ci	One of "yes" or "no".
...	Passed to ggplot2::theme() (ggplot path) or plot().
object	A "ba" object.
title	Plot title.
subtitle	Optional subtitle. If NULL, shows n and LoA summary.
point_alpha	Point transparency.
point_size	Point size.
line_size	Line width for mean/LoA.
shade_ci	Logical; if TRUE, draw shaded CI bands instead of 6 dashed lines.
shade_alpha	Transparency of CI bands.
smoother	One of "none", "loess", "lm" to visualize proportional bias.
symmetrize_y	Logical; if TRUE, y-axis centered at mean difference with symmetric limits.
show_value	Logical; included for a consistent plotting interface. Bland-Altman plots do not overlay numeric cell values, so this argument currently has no effect.
style	Show the pairwise result as "pairs" or "matrices".
pairs	Optional character vector of pair labels to display.
against	Optional single method name; if supplied, only contrasts involving that method are plotted.
facet_scales	Either "free_y" (default) or "fixed".

Details

Given paired measurements $(x_i, y_i)$, Bland-Altman analysis uses $d_i = x_i - y_i$ (or $y_i - x_i$ if mode = 2) and $m_i = (x_i + y_i)/2$. The mean difference $\bar d$ estimates bias. The limits of agreement (LoA) are $\bar d \pm z \cdot s_d$, where $s_d$ is the sample standard deviation of $d_i$ and $z$ (argument loa_multiplier) is typically 1.96 for nominal 95% LoA.

When group2 is omitted and group1 is a wide numeric matrix or data frame, the same two-method calculation is applied to every unordered column pair. The returned pairwise matrices follow the requested mode: with mode = 1, upper-triangle entries represent row minus column; with mode = 2, upper-triangle entries represent column minus row.

Confidence intervals use Student's $t$ distribution with $n-1$ degrees of freedom, with

• Mean-difference CI given by $\bar d \pm t_{n-1,\,1-\alpha/2}\, s_d/\sqrt{n}$; and

• LoA CI given by $(\bar d \pm z\, s_d) \;\pm\; t_{n-1,\,1-\alpha/2}\, s_d\,\sqrt{3/n}$.

Assumptions include approximately normal differences and roughly constant variability across the measurement range; if differences increase with magnitude, consider a transformation before analysis. Missing values are removed pairwise (rows with an NA in either input are dropped before calling the C++ backend).

Probability of agreement, available through prob_agree, is a tolerance-based companion to Bland-Altman analysis. Bland-Altman reports bias and limits of agreement for paired differences. prob_agree() instead uses the sampling distribution of estimated differences to quantify the probability that two estimated quantities or curves agree within a user-specified practical tolerance.

Value

If group1 and group2 are both supplied, an object of class "ba" (list) with elements:

• means, diffs: numeric vectors

• groups: data.frame used after NA removal

• n_obs: integer, number of complete pairs used.

• based.on: compatibility alias for n_obs.

• lower.limit, mean.diffs, upper.limit

• lines: named numeric vector (lower, mean, upper)

• CI.lines: named numeric vector for CIs of those lines

• loa_multiplier, critical.diff

If group2 is omitted and group1 is a wide numeric matrix or data frame, an object of class "ba_matrix" with pairwise components:

• bias: pairwise matrix of Bland-Altman mean differences.

• sd_loa: pairwise matrix of SDs of differences.

• loa_lower, loa_upper: pairwise matrices of LoA endpoints.

• width: pairwise matrix of LoA widths.

• n: integer pairwise matrix of complete-case counts.

• mean_ci_low, mean_ci_high: pairwise matrices of CI bounds for the mean difference.

• loa_lower_ci_low, loa_lower_ci_high: pairwise matrices of CI bounds for the lower LoA.

• loa_upper_ci_low, loa_upper_ci_high: pairwise matrices of CI bounds for the upper LoA.

• methods: character vector of analysed method names.

• loa_multiplier, mode: calculation settings reused for every pair.

Author(s)

Thiago de Paula Oliveira

References

Bland JM, Altman DG (1986). Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet, 307-310.

Bland JM, Altman DG (1999). Measuring agreement in method comparison studies. Statistical Methods in Medical Research, 8(2), 135-160.

See Also

print.ba, plot.ba, ccc, prob_agree, ccc_rm_ustat, ccc_rm_reml

Examples

set.seed(1)
x <- rnorm(100, 100, 10)
y <- x + rnorm(100, 0, 8)
fit_ba <- ba(x, y)
print(fit_ba)
estimate(fit_ba)
tidy(fit_ba)
confint(fit_ba)
plot(fit_ba)

# Pairwise Bland-Altman across 3 methods
set.seed(7)
wide3 <- data.frame(
  ref = rnorm(80, 100, 8),
  m2  = rnorm(80, 101, 8),
  m3  = rnorm(80,  99, 9)
)
fit_ba3 <- ba(wide3)
print(fit_ba3)
summary(fit_ba3)
tidy(fit_ba3)
plot(fit_ba3)

---

Bland-Altman for repeated measurements

Description

Repeated-measures Bland-Altman (BA) analysis for method comparison based on a mixed-effects model fitted to subject-time matched paired differences. The fitted model includes a subject-specific random intercept and, optionally, an AR(1) residual correlation structure within subject.

The function accepts either exactly two methods or $\ge 3$ methods. With exactly two methods it returns a single fitted BA object. With $\ge 3$ methods it fits the same model to every unordered method pair and returns pairwise matrices of results.

Required variables

• response: numeric measurements.

• subject: subject identifier.

• method: method label with at least two distinct levels.

• time: replicate/time key used to form within-subject pairs.

For any analysed pair of methods, only records where both methods are present for the same subject and integer-coerced time contribute to the fit. Rows with missing values in any required field are excluded for that analysed pair.

Usage

ba_rm(
  data = NULL,
  response,
  subject,
  method,
  time,
  conf_level = 0.95,
  n_threads = getOption("matrixCorr.threads", 1L),
  verbose = FALSE,
  loa_multiplier = 1.96,
  include_slope = FALSE,
  use_ar1 = FALSE,
  ar1_rho = NA_real_,
  max_iter = 200L,
  tol = 1e-06
)

## S3 method for class 'ba_repeated'
print(
  x,
  digits = 3,
  ci_digits = 3,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ba_repeated_matrix'
print(
  x,
  digits = 3,
  ci_digits = 3,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  style = c("pairs", "matrices"),
  ...
)

## S3 method for class 'ba_repeated'
plot(
  x,
  title = "Bland-Altman (repeated measurements)",
  subtitle = NULL,
  point_alpha = 0.7,
  point_size = 2.2,
  line_size = 0.8,
  shade_ci = TRUE,
  shade_alpha = 0.08,
  smoother = c("none", "loess", "lm"),
  symmetrize_y = TRUE,
  show_points = TRUE,
  show_value = TRUE,
  ...
)

## S3 method for class 'ba_repeated_matrix'
plot(
  x,
  pairs = NULL,
  against = NULL,
  facet_scales = c("free_y", "fixed"),
  title = "Bland-Altman (repeated, pairwise)",
  point_alpha = 0.6,
  point_size = 1.8,
  line_size = 0.7,
  shade_ci = TRUE,
  shade_alpha = 0.08,
  smoother = c("none", "loess", "lm"),
  show_points = TRUE,
  show_value = TRUE,
  ...
)

Arguments

data	Optional data frame-like object. If supplied, response, subject, method, and time may be column names. Objects not already inheriting from data.frame are first coerced with as.data.frame().
response	Numeric response vector, or a single character string naming the response column in data.
subject	Subject identifier vector (integer, numeric, or factor), or a single character string naming the subject column in data.
method	Method label vector (character, factor, integer, or numeric), or a single character string naming the method column in data. At least two distinct method levels are required.
time	Replicate/time index vector (integer or numeric), or a single character string naming the time column in data. Values are coerced to integer before pairing and before AR(1) contiguity checks.
conf_level	Confidence level for Wald confidence intervals for the reported bias and both LoA endpoints. Must lie in ⁠(0, 1)⁠. Default 0.95.
n_threads	Integer $\geq 1$. Number of OpenMP threads. Defaults to getOption("matrixCorr.threads", 1L).
verbose	Logical. If TRUE, print progress for the pairwise $\ge 3$-method path. It has no material effect in the exactly-two-method path.
loa_multiplier	Positive scalar giving the SD multiplier used to form the limits of agreement. Default is 1.96. In the exported ba_rm() interface this default is fixed and does not depend on conf_level.
include_slope	Logical. If TRUE, the model includes the paired mean as a fixed effect and estimates a proportional-bias slope. The reported BA centre remains the fitted mean difference at the centred reference paired mean used internally by the backend; the returned LoA remain horizontal bands and are not regression-adjusted curves.
use_ar1	Logical. If TRUE, request an AR(1) residual structure within subject over contiguous integer time blocks.
ar1_rho	Optional AR(1) parameter. Must satisfy abs(ar1_rho) < 0.999 when supplied. If NA and use_ar1 = TRUE, the backend estimates rho separately for each analysed pair.
max_iter	Maximum number of EM/GLS iterations used by the backend.
tol	Convergence tolerance for the backend EM/GLS iterations.
x	A "ba_repeated_matrix" object.
digits	Number of digits for estimates (default 3).
ci_digits	Number of digits for CI bounds (default 3).
n	Optional row threshold for compact preview output.
topn	Optional number of leading/trailing rows to show when truncated.
max_vars	Optional maximum number of visible columns; NULL derives this from console width.
width	Optional display width; defaults to getOption("width").
show_ci	One of "yes" or "no".
...	Additional theme adjustments passed to ggplot2::theme(...) (e.g., plot.title.position = "plot", axis.title.x = element_text(size=11)).
style	Show as pairs or matrix format?
title	Plot title (character scalar). Defaults to "Bland-Altman (repeated measurements)" for two methods and "Bland-Altman (repeated, pairwise)" for the faceted matrix plot.
subtitle	Optional subtitle (character scalar). If NULL, a compact summary is shown using the fitted object.
point_alpha	Numeric in ⁠[0, 1]⁠. Transparency for scatter points drawn at (pair mean, pair difference) when point data are available. Passed to ggplot2::geom_point(alpha = ...). Default 0.7.
point_size	Positive numeric. Size of scatter points; passed to ggplot2::geom_point(size = ...). Default 2.2.
line_size	Positive numeric. Line width for horizontal bands (bias and both LoA) and, when requested, the proportional-bias line. Passed to ggplot2::geom_hline(linewidth = ...) (and geom_abline). Default 0.8.
shade_ci	Logical. If TRUE and confidence intervals are available in the object (CI.lines for two methods; ⁠*_ci_*⁠ matrices for the pairwise case), semi-transparent rectangles are drawn to indicate CI bands for the bias and each LoA. If FALSE, dashed horizontal CI lines are drawn instead. Has no effect if CIs are not present. Default TRUE.
shade_alpha	Numeric in ⁠[0, 1]⁠. Opacity of the CI shading rectangles when shade_ci = TRUE. Passed to ggplot2::annotate("rect", alpha = ...). Default 0.08.
smoother	One of "none", "loess", or "lm". Adds an overlaid trend for differences vs means when points are drawn, to visualise proportional bias. "lm" fits a straight line with no SE ribbon; "loess" draws a locally-smoothed curve (span 0.9) with no SE ribbon; "none" draws no smoother. Ignored if show_points = FALSE or if no point data are available.
symmetrize_y	Logical (two-method plot only). If TRUE, the y-axis is centred at the estimated bias and expanded symmetrically to cover all elements used to compute the range (bands, CIs, and points if shown). Default TRUE.
show_points	Logical. If TRUE, per-pair points are drawn for the exactly-two-method path from the stored means and diffs vectors. Pairwise matrix plots draw points reconstructed from the canonical long-form columns stored by ba_rm(). If FALSE or if point data are unavailable, only the bands (and optional CI indicators) are drawn. Default TRUE.
show_value	Logical; included for a consistent plotting interface. Repeated-measures Bland-Altman plots do not overlay numeric cell values, so this argument currently has no effect.
pairs	(Faceted pairwise plot only.) Optional character vector of labels specifying which method contrasts to display. Labels must match the "row - column" convention used by print()/summary() (e.g., "B - A"). Defaults to all upper-triangle pairs.
against	(Faceted pairwise plot only.) Optional single method name. If supplied, facets are restricted to contrasts of the chosen method against all others. Ignored when pairs is provided.
facet_scales	(Faceted pairwise plot only.) Either "free_y" (default) to allow each facet its own y-axis limits, or "fixed" for a common scale across facets. Passed to ggplot2::facet_wrap(scales = ...).

Details

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:

1. fitting the corresponding iid model;

2. computing a moments-based lag-1 estimate from detrended residuals within contiguous blocks, used only as a seed; and

3. 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:

• the sample SD;

• the IQR-based scale $\mathrm{IQR}(m)/1.349$;

• the MAD-based scale $1.4826\,\mathrm{MAD}(m)$.

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.

Value

Either a "ba_repeated" object (exactly two methods) or a "ba_repeated_matrix" object ($\ge 3$ methods).

If exactly two methods are supplied, the returned "ba_repeated" object is a list with components:

• means: numeric vector of paired means $(y_1 + y_2)/2$ used for plotting helpers.

• diffs: numeric vector of paired differences $y_2 - y_1$ used for plotting helpers.

• n_obs: integer number of complete subject-time pairs used.

• based.on: compatibility alias for n_obs.

• mean.diffs: scalar model-based BA centre. When include_slope = FALSE, this is the fitted intercept of the paired- difference model. When include_slope = TRUE, this is the fitted mean difference at the centred reference paired mean used internally by the backend.

• lower.limit, upper.limit: scalar limits of agreement, computed as $\mu_0 \pm \texttt{loa\_multiplier}\sqrt{\sigma_u^2 + \sigma_e^2}$.

• lines: named numeric vector with entries lower, mean, and upper.

• CI.lines: named numeric vector containing Wald confidence interval bounds for the bias and both LoA endpoints: mean.diff.ci.lower, mean.diff.ci.upper, lower.limit.ci.lower, lower.limit.ci.upper, upper.limit.ci.lower, upper.limit.ci.upper.

• loa_multiplier: scalar LoA multiplier actually used.

• critical.diff: scalar LoA half-width $\texttt{loa\_multiplier} \times \texttt{sd\_loa}$.

• include_slope: logical, copied from the call.

• beta_slope: proportional-bias slope on the original paired-mean scale when include_slope = TRUE; otherwise NA.

• sigma2_subject: estimated variance of the subject-level random intercept on paired differences.

• sigma2_resid: estimated residual variance on paired differences.

• use_ar1: logical, copied from the call.

• residual_model: either "ar1" or "iid", indicating the final residual structure actually used.

• ar1_rho: AR(1) correlation actually used in the final fit when residual_model == "ar1"; otherwise NA.

• ar1_estimated: logical indicating whether ar1_rho was estimated internally (TRUE) or supplied by the user (FALSE) when the final residual model is AR(1); otherwise NA.

The confidence level is stored as attr(x, "conf.level").

If $\ge 3$ methods are supplied, the returned "ba_repeated_matrix" object is a list with components:

• bias: numeric $m \times m$ matrix of model-based BA centres. For indices $(j, k)$ with $j < k$, bias[j, k] estimates $\mathrm{method}_k - \mathrm{method}_j$. Thus the matrix orientation is column minus row, not row minus column. The diagonal is NA.

• sd_loa: numeric $m \times m$ matrix of LoA SDs, $\sqrt{\sigma_u^2 + \sigma_e^2}$. This matrix is symmetric.

• loa_lower, loa_upper: numeric $m \times m$ matrices of LoA endpoints corresponding to bias. These satisfy $\texttt{loa\_lower}[j,k] = -\texttt{loa\_upper}[k,j]$ and $\texttt{loa\_upper}[j,k] = -\texttt{loa\_lower}[k,j]$.

• width: numeric $m \times m$ matrix of LoA widths, loa_upper - loa_lower. This matrix is symmetric.

• n: integer $m \times m$ matrix giving the number of complete subject-time pairs used for each analysed contrast. Pairs with fewer than two complete matches are left as NA in the estimate matrices.

• mean_ci_low, mean_ci_high: numeric $m \times m$ matrices of Wald confidence interval bounds for bias.

• loa_lower_ci_low, loa_lower_ci_high: numeric $m \times m$ matrices of Wald confidence interval bounds for the lower LoA.

• loa_upper_ci_low, loa_upper_ci_high: numeric $m \times m$ matrices of Wald confidence interval bounds for the upper LoA.

• slope: optional numeric $m \times m$ matrix of proportional-bias slopes on the original paired-mean scale when include_slope = TRUE; otherwise NULL. This matrix is antisymmetric in sign because each fitted contrast is reversed across the transpose.

• methods: character vector of method levels defining matrix row and column order.

• loa_multiplier: scalar LoA multiplier actually used.

• conf_level: scalar confidence level used for the reported Wald intervals.

• use_ar1: logical, copied from the call.

• ar1_rho: scalar equal to the user-supplied common ar1_rho when use_ar1 = TRUE and a value was supplied; otherwise NA. This field does not store the per-pair estimated AR(1) parameters.

• residual_model: character $m \times m$ matrix whose entries are "ar1", "iid", or NA, indicating the final residual structure used for each pair.

• sigma2_subject: numeric $m \times m$ matrix of estimated subject-level random-intercept variances.

• sigma2_resid: numeric $m \times m$ matrix of estimated residual variances.

• ar1_rho_pair: optional numeric $m \times m$ matrix giving the AR(1) correlation actually used for each pair when the final residual model is "ar1"; otherwise NA for that entry. Present only when use_ar1 = TRUE.

• ar1_estimated: optional logical $m \times m$ matrix indicating whether the pair-specific ar1_rho_pair was estimated internally (TRUE) or supplied by the user (FALSE) for entries whose final residual model is "ar1"; otherwise NA. Present only when use_ar1 = TRUE.

• data_long: canonical long-form data frame with columns .response, .subject, .method, and .time retained for pairwise plotting helpers.

• mapping: named list mapping response, subject, method, and time to the canonical stored columns in data_long.

Author(s)

Thiago de Paula Oliveira

Examples

# -------- Simulate repeated-measures data --------
set.seed(1)

# design (no AR)
# subjects
S   <- 30L
# replicates per subject
Tm  <- 15L
subj <- rep(seq_len(S), each = Tm)
time <- rep(seq_len(Tm), times = S)

# subject signal centered at 0 so BA "bias" won't be driven by the mean level
mu_s  <- rnorm(S, mean = 0, sd = 8)
# constant within subject across replicates
true  <- mu_s[subj]

# common noise (no AR, i.i.d.)
sd_e <- 2
e0   <- rnorm(length(true), 0, sd_e)

# --- Methods ---
# M1: signal + noise
y1 <- true + e0

# M2: same precision as M1; here identical so M3 can be
#     almost perfectly the inverse of both M1 and M2
y2 <- y1 + rnorm(length(true), 0, 0.01)

# M3: perfect inverse of M1 and M2
y3 <- -y1   # = -(true + e0)

# M4: unrelated to all others (pure noise, different scale)
y4 <- rnorm(length(true), 3, 6)

data <- rbind(
  data.frame(y = y1, subject = subj, method = "M1", time = time),
  data.frame(y = y2, subject = subj, method = "M2", time = time),
  data.frame(y = y3, subject = subj, method = "M3", time = time),
  data.frame(y = y4, subject = subj, method = "M4", time = time)
)
data$method <- factor(data$method, levels = c("M1","M2","M3","M4"))

# quick sanity checks
with(data, {
  Y <- split(y, method)
  round(cor(cbind(M1 = Y$M1, M2 = Y$M2, M3 = Y$M3, M4 = Y$M4)), 3)
})

# Run BA (no AR)
ba4 <- ba_rm(
  data = data,
  response = "y", subject = "subject", method = "method", time = "time",
  loa_multiplier = 1.96, conf_level = 0.95,
  include_slope = FALSE, use_ar1 = FALSE
)
summary(ba4)
estimate(ba4)
tidy(ba4)
confint(ba4)
plot(ba4)

# -------- Simulate repeated-measures with AR(1) data --------
set.seed(123)
S <- 40L                      # subjects
Tm <- 50L                     # replicates per subject
methods <- c("A","B","C")     # N = 3 methods
rho <- 0.4                    # AR(1) within-subject across time

ar1_sim <- function(n, rho, sd = 1) {
  z <- rnorm(n)
  e <- numeric(n)
  e[1] <- z[1] * sd
  if (n > 1) for (t in 2:n) e[t] <- rho * e[t-1] + sqrt(1 - rho^2) * z[t] * sd
  e
}

# Subject baseline + time trend (latent "true" signal)
subj <- rep(seq_len(S), each = Tm)
time <- rep(seq_len(Tm), times = S)
# subject effects
mu_s  <- rnorm(S, 50, 7)
trend <- rep(seq_len(Tm) - mean(seq_len(Tm)), times = S) * 0.8
true  <- mu_s[subj] + trend

# Method-specific biases (B has +1.5 constant; C has slight proportional bias)
bias  <- c(A = 0, B = 1.5, C = -0.5)
# proportional component on "true"
prop  <- c(A = 0.00, B = 0.00, C = 0.10)

# Build long data: for each method, add AR(1) noise within subject over time
make_method <- function(meth, sd = 3) {
  e <- unlist(lapply(split(seq_along(time), subj),
                     function(ix) ar1_sim(length(ix), rho, sd)))
  y <- true * (1 + prop[meth]) + bias[meth] + e
  data.frame(y = y, subject = subj, method = meth, time = time,
             check.names = FALSE)
}

data <- do.call(rbind, lapply(methods, make_method))
data$method <- factor(data$method, levels = methods)

# -------- Repeated BA (pairwise matrix) ---------------------
baN <- ba_rm(
  response = data$y, subject = data$subject, method = data$method, time = data$time,
  loa_multiplier = 1.96, conf_level = 0.95,
  include_slope = FALSE,         # estimate proportional bias per pair
  use_ar1 = TRUE, ar1_rho = rho
)

# Matrices (row - column orientation)
print(baN)
summary(baN)
tidy(baN)

# Faceted BA scatter by pair
plot(baN, smoother = "lm", facet_scales = "free_y")

# -------- Two-method AR(1) path (A vs B only) ------------------------------
data_AB <- subset(data, method %in% c("A","B"))
baAB <- ba_rm(
  response = data_AB$y, subject = data_AB$subject,
  method = droplevels(data_AB$method), time = data_AB$time,
  include_slope = FALSE, use_ar1 = TRUE, ar1_rho = 0.4
)
print(baAB)
plot(baAB)

---

Biweight mid-correlation (bicor)

Description

Computes pairwise biweight mid-correlations for numeric data. Bicor is a robust, Pearson-like correlation that down-weights outliers and heavy-tailed observations. Optional large-sample confidence intervals are available as a derived feature.

Usage

bicor(
  data,
  na_method = c("error", "pairwise", "complete"),
  ci = FALSE,
  conf_level = 0.95,
  n_threads = getOption("matrixCorr.threads", 1L),
  output = c("matrix", "sparse", "edge_list"),
  threshold = 0,
  diag = TRUE,
  c_const = 9,
  max_p_outliers = 1,
  pearson_fallback = c("hybrid", "none", "all"),
  mad_consistent = FALSE,
  w = NULL,
  sparse_threshold = NULL
)

diag.bicor(x, ...)

## S3 method for class 'bicor'
print(
  x,
  digits = 4,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  ci_digits = 3,
  show_ci = NULL,
  na_print = "NA",
  ...
)

## S3 method for class 'bicor'
plot(
  x,
  title = "Biweight mid-correlation heatmap",
  reorder = c("none", "hclust"),
  triangle = c("full", "lower", "upper"),
  low_color = "indianred1",
  mid_color = "white",
  high_color = "steelblue1",
  value_text_size = 3,
  ci_text_size = 2.5,
  show_value = TRUE,
  na_fill = "grey90",
  ...
)

## S3 method for class 'bicor'
summary(
  object,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  ci_digits = 3,
  p_digits = 4,
  show_ci = NULL,
  ...
)

## S3 method for class 'summary.bicor'
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

Arguments

data	A numeric matrix or a data frame containing numeric columns. Factors, logicals and common time classes are dropped in the data-frame path. Missing values are not allowed unless na_method = "pairwise".
na_method	Character scalar controlling missing-data handling. "error" rejects missing, NaN, and infinite values. "pairwise" recomputes each estimate on its own pairwise complete-case overlap. This is permissive, but different matrix entries may be based on different rows. "complete" performs listwise deletion once across the retained numeric columns and then computes the estimator on the common complete sample.
ci	Logical (default FALSE). If TRUE, attach approximate Fisher-z confidence intervals for the off-diagonal biweight mid-correlations. This confidence interval is provided as an additional large-sample approximation.
conf_level	Confidence level used when ci = TRUE. Default is 0.95.
n_threads	Integer $\geq 1$. Number of OpenMP threads. Defaults to getOption("matrixCorr.threads", 1L).
output	Output representation for the computed estimates. "matrix" (default): full dense matrix; best when you need matrix algebra, dense heatmaps, or full compatibility with existing code. "sparse": sparse matrix from Matrix containing only retained entries; best when many values are dropped by thresholding. "edge_list": long-form data frame with columns row, col, value; convenient for filtering, joins, and network-style workflows.
threshold	Non-negative absolute-value filter for non-matrix outputs: keep entries with abs(value) >= threshold. Use threshold > 0 when you want only stronger associations (typically with output = "sparse" or "edge_list"). Keep threshold = 0 to retain all values. Must be 0 when output = "matrix".
diag	Logical; whether to include diagonal entries in "sparse" and "edge_list" outputs.
c_const	Positive numeric. Tukey biweight tuning constant applied to the raw MAD; default 9 (Langfelder & Horvath's convention).
max_p_outliers	Numeric in (0, 1]. Optional cap on the maximum proportion of outliers on each side; if < 1, side-specific rescaling maps those quantiles to |u|=1. Use 1 to disable.
pearson_fallback	Character scalar indicating the fallback policy. One of: "hybrid" (default): if a column has MAD = 0, that column uses Pearson standardisation, yielding a hybrid correlation. "none": return NA if a column has MAD = 0 or becomes degenerate after weighting. "all": force ordinary Pearson for all columns.
mad_consistent	Logical; if TRUE, use the normal-consistent MAD (MAD_raw * 1.4826) in the bicor weights. Default FALSE to match Langfelder & Horvath (2012).
w	Optional non-negative numeric vector of length nrow(data) giving row weights. When supplied, weighted medians/MADs are used and Tukey weights are multiplied by w before normalisation.
sparse_threshold	Optional numeric $\geq 0$. If supplied, sets entries with |r| < sparse_threshold to 0 and returns a sparse "ddiMatrix" (requires Matrix).
x	An object of class summary.bicor.
...	Additional arguments passed to ggplot2::theme() or other layers.
digits	Integer; number of decimal places used for the matrix.
n	Optional row threshold for compact preview output.
topn	Optional number of leading/trailing rows to show when truncated.
max_vars	Optional maximum number of visible columns; NULL derives this from console width.
width	Optional display width; defaults to getOption("width").
ci_digits	Integer; digits for bicor confidence limits in the pairwise summary.
show_ci	One of "yes" or "no".
na_print	Character; how to display missing values.
title	Plot title. Default is "Biweight mid-correlation heatmap".
reorder	Character; one of "none" (default) or "hclust". If "hclust", variables are reordered by complete-linkage clustering on the distance $d = 1 - r$, after replacing NA by 0 for clustering purposes only.
triangle	One of "full" (default), "lower", or "upper" to display the full matrix or a single triangle.
low_color, mid_color, high_color	Colours for the gradient in scale_fill_gradient2. Defaults are "indianred1", "white", "steelblue1".
value_text_size	Numeric; font size for cell labels. Set to NULL to suppress labels (recommended for large matrices).
ci_text_size	Text size for confidence-interval labels in the heatmap.
show_value	Logical; if TRUE (default), overlay numeric values on the heatmap tiles.
na_fill	Fill colour for NA cells. Default "grey90".
object	An object of class bicor.
p_digits	Integer; digits for bicor p-values in the pairwise summary.

Details

For a column $x = (x_a)_{a=1}^m$, let $\mathrm{med}(x)$ be the median and $\mathrm{MAD}(x) = \mathrm{med}(|x - \mathrm{med}(x)|)$ the (raw) median absolute deviation. If mad_consistent = TRUE, the consistent scale $\mathrm{MAD}^\star(x) = 1.4826\,\mathrm{MAD}(x)$ is used. With tuning constant $c>0$, define

$u_a = \frac{x_a - \mathrm{med}(x)}{c\,\mathrm{MAD}^{(\star)}(x)}.$

The Tukey biweight gives per-observation weights

$w_a = (1 - u_a^2)^2\,\mathbf{1}\{|u_a| < 1\}.$

Robust standardisation of a column is

$\tilde x_a = \frac{(x_a - \mathrm{med}(x))\,w_a}{ \sqrt{\sum_{b=1}^m \big[(x_b - \mathrm{med}(x))\,w_b\big]^2}}.$

For two columns $x,y$, the biweight mid-correlation is

$\mathrm{bicor}(x,y) = \sum_{a=1}^m \tilde x_a\,\tilde y_a \in [-1,1].$

Capping the maximum proportion of outliers (max_p_outliers). If max_p_outliers < 1, let $q_L = Q_x(\text{max\_p\_outliers})$ and $q_U = Q_x(1-\text{max\_p\_outliers})$ be the lower/upper quantiles of $x$. If the corresponding $|u|$ at either quantile exceeds 1, $u$ is rescaled separately on the negative and positive sides so that those quantiles land at $|u|=1$. This guarantees that all observations between the two quantiles receive positive weight. Note the bound applies per side, so up to $2\,\text{max\_p\_outliers}$ of observations can be treated as outliers overall.

Fallback when for zero MAD / degeneracy (pearson_fallback). If a column has $\mathrm{MAD}=0$ or the robust denominator becomes zero, the following rules apply:

• "none" when correlations involving that column are NA (diagonal remains 1).

• "hybrid" when only the affected column switches to Pearson standardisation $\bar x_a = (x_a - \overline{x}) / \sqrt{\sum_b (x_b - \overline{x})^2}$, yielding the hybrid correlation

  $\mathrm{bicor}_{\mathrm{hyb}}(x,y) = \sum_a \bar x_a\,\tilde y_a,$

  with the other column still robust-standardised.

• "all" when all columns use ordinary Pearson standardisation; the result equals stats::cor(..., method="pearson") when the NA policy matches.

Handling missing values (na_method).

• "error" (default): inputs must be finite; this yields a symmetric, positive semidefinite (PSD) matrix since $R = \tilde X^\top \tilde X$.

• "pairwise": each $R_{jk}$ is computed on the intersection of rows where both columns are finite. Pairs with fewer than 5 overlapping rows return NA (guarding against instability). Pairwise deletion can break PSD, as in the Pearson case.

Row weights (w). When w is supplied (non-negative, length $m$), the weighted median $\mathrm{med}_w(x)$ and weighted MAD $\mathrm{MAD}_w(x) = \mathrm{med}_w(|x - \mathrm{med}_w(x)|)$ are used to form $u$. The Tukey weights are then multiplied by the observation weights prior to normalisation:

$\tilde x_a = \frac{(x_a - \mathrm{med}_w(x))\,w_a\,w^{(\mathrm{obs})}_a}{ \sqrt{\sum_b \big[(x_b - \mathrm{med}_w(x))\,w_b\,w^{(\mathrm{obs})}_b\big]^2}},$

where $w^{(\mathrm{obs})}_a \ge 0$ are the user-supplied row weights and $w_a$ are the Tukey biweights built from the weighted median/MAD. Weighted pairwise behaves analogously on each column pair's overlap.

MAD choice (mad_consistent). Setting mad_consistent = TRUE multiplies the raw MAD by 1.4826 inside $u$. Equivalently, it uses an effective tuning constant $c^\star = c \times 1.4826$. The default FALSE reproduces the convention in Langfelder & Horvath (2012).

Optional sparsification (sparse_threshold). If provided, entries with $|r| < \text{sparse\_threshold}$ are set to 0 and the result is returned as a "ddiMatrix" (diagonal is forced to 1). This is a post-processing step that does not alter the per-pair estimates.

Computation and threads. Columns are robust-standardised in parallel and the matrix is formed as $R = \tilde X^\top \tilde X$. n_threads selects the number of OpenMP threads; by default it uses getOption("matrixCorr.threads", 1L).

Large-sample inference. For a pairwise estimate $r$ computed from $n_{jk}$ observed rows, the standard large-sample summaries use

$Z_{jk} = \frac{1}{2}\log\!\left(\frac{1 + r}{1 - r}\right)\sqrt{n_{jk} - 2}$

and

$T_{jk} = \sqrt{n_{jk} - 2}\;\frac{|r|}{\sqrt{1-r^2}}.$

The reported p-value is the two-sided Student-$t$ tail probability with $n_{jk}-2$ degrees of freedom. When ci = TRUE, the package also reports an approximate Fisher-z confidence interval obtained from

$z_{jk} = \operatorname{atanh}(r), \qquad \operatorname{SE}(z_{jk}) = \frac{1}{\sqrt{n_{jk}-3}},$

followed by back-transformation with tanh(). Confidence intervals are currently available only for dense, unweighted outputs.

Basic properties. $\mathrm{bicor}(a x + b,\; c y + d) = \mathrm{sign}(ac)\,\mathrm{bicor}(x,y)$. With no missing data (and with per-column hybrid/robust standardisation), the output is symmetric and PSD. As with Pearson, affine equivariance does not hold for the associated biweight midcovariance.

Value

A symmetric correlation matrix with class bicor (or a dgCMatrix if sparse_threshold is used), with attributes: method = "biweight_mid_correlation", description, and package = "matrixCorr". Downstream code should be prepared to handle either a dense numeric matrix or a sparse dgCMatrix. When ci = TRUE, the object also carries a ci attribute with elements est, lwr.ci, upr.ci, conf.level, and ci.method, together with an inference attribute containing the standard large-sample summary matrices estimate, statistic, p_value, Z, and n_obs. Pairwise complete-case counts are stored in attr(x, "diagnostics")$n_complete. Internally, all medians/MADs, Tukey weights, optional pairwise-NA handling, and OpenMP loops are implemented in the C++ helpers (bicor_*_cpp()), so the R wrapper mostly validates arguments and dispatches to the appropriate backend.

Invisibly returns x.

A ggplot object.

Author(s)

Thiago de Paula Oliveira

References

Langfelder, P. & Horvath, S. (2012). Fast R Functions for Robust Correlations and Hierarchical Clustering. Journal of Statistical Software, 46(11), 1-17. doi:10.18637/jss.v046.i11

Examples

set.seed(1)
X <- matrix(rnorm(2000 * 40), 2000, 40)
R <- bicor(X, c_const = 9, max_p_outliers = 1,
                       pearson_fallback = "hybrid")
print(attr(R, "method"))
summary(R)
R_ci <- bicor(X[, 1:5], ci = TRUE)
summary(R_ci)
estimate(R_ci)
tidy(R_ci)
ci(R_ci)
confint(R_ci)

# Interactive viewing (requires shiny)
if (interactive() && requireNamespace("shiny", quietly = TRUE)) {
  view_corr_shiny(R)
}

---

Biserial Correlation Between Continuous and Binary Variables

Description

Computes biserial correlations between continuous variables in data and binary variables in y. Both pairwise vector mode and rectangular matrix/data-frame mode are supported.

Usage

biserial(data, y, na_method = c("error", "pairwise"), ci = FALSE, p_value = FALSE,
  conf_level = 0.95, ...)

## S3 method for class 'biserial_corr'
print(
  x,
  digits = 4,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'biserial_corr'
plot(
  x,
  title = "Biserial correlation heatmap",
  low_color = "indianred1",
  high_color = "steelblue1",
  mid_color = "white",
  value_text_size = 4,
  ci_text_size = 3,
  show_value = TRUE,
  ...
)

## S3 method for class 'biserial_corr'
summary(
  object,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  ci_digits = 3,
  p_digits = 4,
  show_ci = NULL,
  ...
)

## S3 method for class 'summary.biserial_corr'
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

Arguments

data	A numeric vector, matrix, or data frame containing continuous variables.
y	A binary vector, matrix, or data frame. In data-frame mode, only two-level columns are retained.
na_method	Character scalar controlling missing-data handling. "error" rejects missing values. "pairwise" uses pairwise complete cases.
ci	Logical (default FALSE). If TRUE, attach approximate large-sample confidence intervals derived from a Fisher $z$-transformation of the biserial estimate.
p_value	Logical (default FALSE). If TRUE, attach model-based large-sample p-values, test statistics, and degrees of freedom for each biserial estimate.
conf_level	Confidence level used when ci = TRUE. Default is 0.95.
...	Additional arguments passed to print().
x	An object of class summary.biserial_corr.
digits	Integer; number of decimal places to print.
n	Optional row threshold for compact preview output.
topn	Optional number of leading/trailing rows to show when truncated.
max_vars	Optional maximum number of visible columns; NULL derives this from console width.
width	Optional display width; defaults to getOption("width").
show_ci	One of "yes" or "no".
title	Plot title. Default is "Biserial correlation heatmap".
low_color	Color for the minimum correlation.
high_color	Color for the maximum correlation.
mid_color	Color for zero correlation.
value_text_size	Font size used in tile labels.
ci_text_size	Text size for confidence intervals in the heatmap.
show_value	Logical; if TRUE (default), overlay numeric values on the heatmap tiles.
object	An object of class biserial_corr.
ci_digits	Integer; digits for biserial confidence limits in the pairwise summary.
p_digits	Integer; digits for biserial p-values in the pairwise summary.

Details

The biserial correlation is the special two-category case of the polyserial model. It assumes that a binary variable $Y$ arises by thresholding an unobserved standard-normal variable $Z$ that is jointly normal with a continuous variable $X$. Writing $p = P(Y = 1)$ and $q = 1-p$, let $z_p = \Phi^{-1}(p)$ and $\phi(z_p)$ be the standard-normal density evaluated at $z_p$. If $\bar x_1$ and $\bar x_0$ denote the sample means of $X$ in the two observed groups and $s_x$ is the sample standard deviation of $X$, the usual biserial estimator is

$r_b = \frac{\bar x_1 - \bar x_0}{s_x} \frac{pq}{\phi(z_p)}.$

This is exactly the estimator implemented in the underlying C++ kernel.

Assumptions. The biserial coefficient is appropriate when the observed binary variable is viewed as a thresholded version of an unobserved continuous latent variable that is jointly normal with the observed continuous variable. The optional p-values and confidence intervals adopt this latent-normal interpretation together with the usual large-sample approximations used for correlation coefficients. These inferential quantities are therefore model-based and should not be interpreted as distribution-free summaries.

Inference. When p_value = TRUE, the package reports the large-sample $t$-statistic

$t = r_b \sqrt{\frac{n - 2}{1 - r_b^2}},$

referenced to a Student $t$-distribution with $n - 2$ degrees of freedom. When ci = TRUE, the package forms an approximate Fisher $z$-interval by transforming $r_b$ with $z = \operatorname{atanh}(r_b)$, using standard error $1 / \sqrt{n - 3}$, and mapping the limits back with $\tanh(\cdot)$. The CI is therefore an internal large-sample extension and is only computed when explicitly requested.

In vector mode a single biserial correlation is returned. In matrix/data-frame mode, every numeric column of data is paired with every binary column of y, producing a rectangular matrix of continuous-by-binary biserial correlations.

Unlike the point-biserial correlation, which is just Pearson correlation on a 0/1 coding of the binary variable, the biserial coefficient explicitly assumes an underlying latent normal threshold model and rescales the mean difference accordingly.

Computational complexity. If data has $p_x$ continuous columns and y has $p_y$ binary columns, the matrix path computes $p_x p_y$ closed-form estimates with negligible extra memory beyond the output matrix.

Value

If both data and y are vectors, a numeric scalar. Otherwise a numeric matrix of class biserial_corr with rows corresponding to the continuous variables in data and columns to the binary variables in y. Matrix outputs carry attributes method, description, and package = "matrixCorr". When p_value = TRUE, the object also carries an inference attribute with matrices estimate, statistic, parameter, p_value, and n_obs. When ci = TRUE, it additionally carries a ci attribute with matrices lwr.ci and upr.ci, plus attr(x, "conf.level"). Scalar outputs keep the same point estimate and gain the same metadata only when inference is requested.

Author(s)

Thiago de Paula Oliveira

References

Olsson, U., Drasgow, F., & Dorans, N. J. (1982). The polyserial correlation coefficient. Psychometrika, 47(3), 337-347.

Fisher, R. A. (1921). On the probable error of a coefficient of correlation deduced from a small sample. Metron, 1, 3-32.

Examples

set.seed(126)
n <- 1000
Sigma <- matrix(c(
  1.00, 0.35, 0.50, 0.25,
  0.35, 1.00, 0.30, 0.55,
  0.50, 0.30, 1.00, 0.40,
  0.25, 0.55, 0.40, 1.00
), 4, 4, byrow = TRUE)

Z <- mnormt::rmnorm(n = n, mean = rep(0, 4), varcov = Sigma)
X <- data.frame(x1 = Z[, 1], x2 = Z[, 2])
Y <- data.frame(
  g1 = Z[, 3] > stats::qnorm(0.65),
  g2 = Z[, 4] > stats::qnorm(0.55)
)

bs <- biserial(X, Y, ci = TRUE, p_value = TRUE)
print(bs, digits = 3)
summary(bs)
estimate(bs)
tidy(bs)
ci(bs)
confint(bs)
plot(bs)

---

Pairwise Lin's concordance correlation coefficient

Description

Computes all pairwise Lin's Concordance Correlation Coefficients (CCC) from the numeric columns of a matrix or data frame. CCC measures both precision (Pearson correlation) and accuracy (closeness to the 45-degree line). This function is backed by a high-performance 'C++' implementation.

Lin's CCC quantifies the concordance between a new test/measurement and a gold-standard for the same variable. Like a correlation, CCC ranges from -1 to 1 with perfect agreement at 1, and it cannot exceed the absolute value of the Pearson correlation between variables. It can be legitimately computed even with small samples (e.g., 10 observations), and results are often similar to intraclass correlation coefficients. CCC provides a single summary of agreement, but it may not capture systematic bias; a Bland-Altman plot (differences vs. means) is recommended to visualize bias, proportional trends, and heteroscedasticity (see ba).

Usage

ccc(
  data,
  ci = FALSE,
  conf_level = 0.95,
  na_method = c("error", "complete", "pairwise"),
  n_threads = getOption("matrixCorr.threads", 1L),
  output = c("matrix", "sparse", "edge_list"),
  threshold = 0,
  diag = TRUE,
  verbose = FALSE
)

## S3 method for class 'ccc'
print(
  x,
  digits = 4,
  ci_digits = 4,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ccc'
summary(
  object,
  digits = 4,
  ci_digits = 2,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'summary.ccc'
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)

## S3 method for class 'ccc'
plot(
  x,
  title = "Lin's Concordance Correlation Heatmap",
  low_color = "indianred1",
  high_color = "steelblue1",
  mid_color = "white",
  value_text_size = 4,
  ci_text_size = 3,
  show_value = TRUE,
  ...
)

Arguments

data	A numeric matrix or data frame with at least two numeric columns. Non-numeric columns will be ignored.
ci	Logical; if TRUE, return lower and upper confidence bounds
conf_level	Confidence level for CI, default = 0.95
na_method	Character scalar controlling missing-data handling. "error" rejects missing, NaN, and infinite values. "complete" removes rows incomplete in any retained numeric column before computing all CCC estimates. "pairwise" recomputes each method pair on its own complete finite overlap.
n_threads	Integer $\geq 1$. Number of OpenMP threads. Defaults to getOption("matrixCorr.threads", 1L).
output	Output representation for the computed estimates. "matrix" (default): full dense matrix; best when you need matrix algebra, dense heatmaps, or full compatibility with existing code. "sparse": sparse matrix from Matrix containing only retained entries; best when many values are dropped by thresholding. "edge_list": long-form data frame with columns row, col, value; convenient for filtering, joins, and network-style workflows.
threshold	Non-negative absolute-value filter for non-matrix outputs: keep entries with abs(value) >= threshold. Use threshold > 0 when you want only stronger associations (typically with output = "sparse" or "edge_list"). Keep threshold = 0 to retain all values. Must be 0 when output = "matrix".
diag	Logical; whether to include diagonal entries in "sparse" and "edge_list" outputs.
verbose	Logical; if TRUE, prints how many threads are used
x	An object of class "ccc" (either a matrix or a list with CIs).
digits	Integer; decimals for CCC estimates (default 4).
ci_digits	Integer; decimals for CI bounds (default 2).
n	Optional row threshold for compact preview output.
topn	Optional number of leading/trailing rows to show when truncated.
max_vars	Optional maximum number of visible columns; NULL derives this from console width.
width	Optional display width; defaults to getOption("width").
show_ci	One of "yes" or "no".
...	Passed to ggplot2::theme().
object	A "ccc" or "ccc_ci" object to summarize.
title	Title for the plot.
low_color	Color for low CCC values.
high_color	Color for high CCC values.
mid_color	Color for mid CCC values.
value_text_size	Text size for CCC values in the heatmap.
ci_text_size	Text size for confidence intervals.
show_value	Logical; if TRUE (default), overlay numeric values on the heatmap tiles.

Details

Lin's CCC is defined as

$\rho_c \;=\; \frac{2\,\mathrm{cov}(X, Y)} {\sigma_X^2 + \sigma_Y^2 + (\mu_X - \mu_Y)^2},$

where $\mu_X,\mu_Y$ are the means, $\sigma_X^2,\sigma_Y^2$ the variances, and $\mathrm{cov}(X,Y)$ the covariance. Equivalently,

$\rho_c \;=\; r \times C_b, \qquad r \;=\; \frac{\mathrm{cov}(X,Y)}{\sigma_X \sigma_Y}, \quad C_b \;=\; \frac{2 \sigma_X \sigma_Y} {\sigma_X^2 + \sigma_Y^2 + (\mu_X - \mu_Y)^2}.$

Hence $|\rho_c| \le |r| \le 1$, $\rho_c = r$ iff $\mu_X=\mu_Y$ and $\sigma_X=\sigma_Y$, and $\rho_c=1$ iff, in addition, $r=1$. CCC is symmetric in $(X,Y)$ and penalises both location and scale differences; unlike Pearson's $r$, it is not invariant to affine transformations that change means or variances.

When ci = TRUE, large-sample confidence intervals for $\rho_c$ are returned for each pair. The implementation uses Lin's delta-method standard error and then forms limits on a Fisher-z transformed CCC scale before mapping back to $[-1, 1]$. For speed, CIs are omitted when ci = FALSE.

If either variable has zero variance, $\rho_c$ is undefined and NA is returned for that pair (including the diagonal).

Missing-data handling follows the same contract as pearson_corr and spearman_rho. With na_method = "error" (default), missing and non-finite values are rejected. With "complete", rows incomplete in any retained numeric column are removed once before all pairwise CCC estimates are computed. With "pairwise", each method pair is computed on its own complete finite overlap, and n_complete diagnostics may vary by pair. Pairwise CIs require at least three complete observations for that pair.

Negative CCC estimates are fully supported. Matrix output preserves the signed estimate, while thresholded "sparse" and "edge_list" outputs retain entries according to abs(CCC) >= threshold.

Probability of agreement, available through prob_agree, answers a different question from CCC. CCC is a coefficient-based summary of concordance between paired measurements. prob_agree() follows Stevens and Anderson-Cook (2017) and estimates the probability that two estimated quantities or curves differ by no more than a user-specified practical tolerance.

Value

A symmetric numeric matrix with class "ccc" and attributes:

• method: The method used ("Lin's concordance")

• description: Description string

If ci = FALSE, returns matrix of class "ccc". If ci = TRUE, returns a list with elements: est, lwr.ci, upr.ci.

For summary.ccc, a data frame with columns item1, item2, estimate, and (optionally) lwr, upr, plus n_complete when available.

Author(s)

Thiago de Paula Oliveira

References

Lin L (1989). A concordance correlation coefficient to evaluate reproducibility. Biometrics 45: 255-268.

Lin L (2000). A note on the concordance correlation coefficient. Biometrics 56: 324-325.

Bland J, Altman D (1986). Statistical methods for assessing agreement between two methods of clinical measurement. The Lancet 327: 307-310.

See Also

print.ccc, plot.ccc, ba, prob_agree, and cia. ccc() answers the question "How well do two paired measurements agree overall, accounting for both correlation and mean/scale bias?". In contrast, cia() answers "Are two or more methods interchangeable at the individual level relative to within-method replicate disagreement?".

For repeated measurements look at ccc_rm_reml, ccc_rm_ustat or ba_rm

Examples

# Example with multivariate normal data
Sigma <- matrix(c(1, 0.5, 0.3,
                  0.5, 1, 0.4,
                  0.3, 0.4, 1), nrow = 3)
mu <- c(0, 0, 0)
set.seed(123)
mat_mvn <- MASS::mvrnorm(n = 100, mu = mu, Sigma = Sigma)
result_mvn <- ccc(mat_mvn, ci = TRUE)
print(result_mvn)
summary(result_mvn)
estimate(result_mvn)
tidy(result_mvn)
ci(result_mvn)
confint(result_mvn)
plot(result_mvn)

# Interactive viewing (requires shiny)
if (interactive() && requireNamespace("shiny", quietly = TRUE)) {
  view_corr_shiny(result_mvn)
}

---

Poisson GLMM concordance correlation for count agreement

Description

Computes pairwise generalized concordance correlation coefficients (CCC) for count outcomes measured repeatedly by two or more methods, observers, or devices. The current implementation fits Poisson-log generalized linear mixed models and reports total, inter-method, and intra-method agreement summaries from the fitted mean and variance components.

Usage

ccc_glmm(
  data,
  response,
  subject,
  method,
  replicate = NULL,
  family = "poisson",
  link = "log",
  overdispersion = c("none", "pearson"),
  include_subject_method = FALSE,
  ci = FALSE,
  conf_level = 0.95,
  max_iter = 1000,
  tol = 1e-08,
  n_threads = getOption("matrixCorr.threads", 1L),
  verbose = FALSE
)

Arguments

data	A data frame containing the measurements.
response	Character. Name of the non-negative integer-like count column.
subject	Character. Name of the subject identifier column.
method	Character. Name of the method/observer column.
replicate	Optional character. Name of the replicate column. If NULL, one replicate per subject-method cell is assumed.
family	Character. Currently only "poisson" is implemented.
link	Character. Currently only "log" is implemented.
overdispersion	Character. "none" fixes $\phi = 1$; "pearson" estimates $\phi$ from post-fit Pearson residuals.
include_subject_method	Logical. If TRUE, include an additional subject-by-method random intercept variance component.
ci	Logical. If TRUE, compute delta-method confidence intervals for the total, inter-method, and intra-method CCC estimates.
conf_level	Confidence level for delta-method confidence intervals.
max_iter	Positive integer. Maximum optimiser iterations.
tol	Positive numeric convergence tolerance.
n_threads	Integer $\geq 1$. Number of OpenMP threads to make available to compiled code.
verbose	Logical. If TRUE, emit progress messages.

Details

The fitted model for a pair of methods is

$Y_{ijl} \mid \alpha_i, \gamma_{ij} \sim \mathrm{Poisson}(\mu_{ijl}), \qquad \log(\mu_{ijl}) = \eta_j + \alpha_i + \gamma_{ij},$

where $i$ indexes subjects, $j$ indexes methods, and $l$ indexes replicate readings. The subject effect $\alpha_i$ captures between-subject heterogeneity. When include_subject_method = TRUE, $\gamma_{ij}$ captures subject-specific method departures; otherwise it is fixed at zero. The fixed method difference contributes to disagreement through sigma2_method.

The model is fitted by marginal maximum likelihood using Gauss-Hermite quadrature. The random-intercept model uses 40 quadrature points. The subject-by-method model uses tensor-product quadrature with fewer points per dimension because it integrates over a three-dimensional subject block. The reported CCC quantities are then computed from the fitted Poisson-log mean and variance components.

The main matrix value, rho_ccc, is the total agreement coefficient. Use it as the primary overall count-agreement summary when each individual reading is the unit of inference. It penalizes lack of between-subject signal, systematic method bias, subject-by-method disagreement, and Poisson residual variation.

rho_ccc_inter is the inter-method agreement coefficient for averages over replicated readings. It is useful when decisions are based on the mean of $m$ replicate readings per subject-method cell. Replication reduces only the residual count variation term; it does not dilute systematic method disagreement or between-subject variation.

rho_ccc_intra_method1 and rho_ccc_intra_method2 are method-specific repeatability coefficients. Use them to diagnose whether one method is internally more repeatable than the other on the count scale. These are not direct method-comparison coefficients; they describe within-method reproducibility after accounting for the fitted mean and random effects.

precision isolates the share of non-systematic variation attributable to subject ranking/heterogeneity, while accuracy is the ratio rho_ccc / precision. Low accuracy with reasonable precision indicates that disagreement is driven mainly by method bias or extra method-specific variation rather than poor subject discrimination.

overdispersion = "none" fixes $\phi = 1$, the Poisson model. overdispersion = "pearson" replaces the residual term by a Pearson dispersion estimate. Use the Pearson adjustment as a sensitivity analysis when counts appear more variable than the Poisson model allows; it should generally reduce CCC when extra-Poisson variation is present.

Confidence intervals

When ci = TRUE, large-sample confidence intervals are computed for rho_ccc, rho_ccc_inter, rho_ccc_intra_method1, and rho_ccc_intra_method2. The implementation uses a delta-method standard error based on the fitted GLMM parameter vector and the inverse Hessian of the marginal negative log-likelihood:

$\mathrm{Var}\{g(\hat\theta)\} \approx \nabla g(\hat\theta)^\top \mathrm{Var}(\hat\theta) \nabla g(\hat\theta),$

where $g(\theta)$ is the relevant CCC function. Gradients are evaluated numerically by central finite differences. The reported point estimates are not replaced by $g(\hat\theta)$; they remain the values from the standard point-estimate path.

Intervals are formed on a Fisher-Z transformed CCC scale, $z = 0.5\log\{(1+\rho)/(1-\rho)\}$, and then back-transformed to the CCC scale. This is the same broad strategy used elsewhere for CCC-style intervals because it is usually more stable than a raw Wald interval near the boundaries. CI limits are not forcibly truncated to $[0, 1]$.

For overdispersion = "pearson", $\phi$ is a post-fit Pearson dispersion estimate rather than a likelihood parameter. Delta-method confidence intervals therefore treat $\phi$ as fixed and a warning is issued.

This function currently implements the Poisson-log count-data case. The family and link arguments are present for API stability and future extensions, but only family = "poisson" and link = "log" are currently supported. The returned estimates are variance-component CCCs constrained to $[0, 1]$; they are not Lin's raw moment CCC and should not be expected to produce negative values.

Value

A symmetric numeric matrix of class c("ccc_glmm", "ccc") containing rho_ccc. Additional pairwise matrices are stored as attributes:

• rho_ccc_inter: agreement for replicated method averages.

• rho_ccc_intra_method1, rho_ccc_intra_method2: method-specific repeatability coefficients.

• sigma2_subject, sigma2_method, sigma2_subject_method: fitted variance/disagreement components.

• phi, mu, precision, accuracy: fitted count-scale diagnostics used in the CCC decomposition.

• beta0, beta_method, nll, logLik, convergence_code: fitting diagnostics.

• n_obs, n_subjects, m_reps: design diagnostics.

• when ci = TRUE, standard error and confidence-limit matrices for rho_ccc, rho_ccc_inter, rho_ccc_intra_method1, and rho_ccc_intra_method2.

References

Carrasco JL (2010). A generalized concordance correlation coefficient based on the variance components generalized linear mixed models for overdispersed count data. Biometrics.

Examples

# Example 1: overall agreement for two Poisson count methods.
# Here the methods have similar means, so rho_ccc is mainly driven by
# between-subject count heterogeneity versus Poisson residual noise.
set.seed(1)
df1 <- expand.grid(
  subject = factor(seq_len(12)),
  method = factor(c("A", "B")),
  replicate = factor(seq_len(2))
)
subject_eff <- rnorm(12, 0, 0.5)
df1$eta <- 1.1 + subject_eff[as.integer(df1$subject)]
df1$y <- rpois(nrow(df1), exp(df1$eta))

fit1 <- ccc_glmm(df1, "y", "subject", "method", replicate = "replicate",
ci = TRUE)
fit1
summary(fit1)
estimate(fit1)
tidy(fit1)
ci(fit1)
confint(fit1)


# Example 2: method bias lowers total agreement.
# This tests whether a systematic method shift is reflected in rho_ccc
# and the accuracy component.
set.seed(2)
df2 <- expand.grid(
  subject = factor(seq_len(12)),
  method = factor(c("A", "B")),
  replicate = factor(seq_len(2))
)
subject_eff <- rnorm(12, 0, 0.6)
df2$eta <- 1.0 + subject_eff[as.integer(df2$subject)] +
  ifelse(df2$method == "B", 0.6, 0)
df2$y <- rpois(nrow(df2), exp(df2$eta))
fit2 <- ccc_glmm(df2, "y", "subject", "method", replicate = "replicate")
summary(fit2)

# Example 3: subject-by-method variation.
# This tests whether individual subjects respond differently by method.
# The subject-method component is useful when disagreement is not explained
# by a single fixed method bias.
set.seed(3)
df3 <- expand.grid(
  subject = factor(seq_len(10)),
  method = factor(c("A", "B")),
  replicate = factor(seq_len(2))
)
subject_eff <- rnorm(10, 0, 0.4)
subject_method_eff <- matrix(rnorm(20, 0, 0.25), nrow = 10)
method_id <- as.integer(df3$method)
df3$eta <- 1.1 + subject_eff[as.integer(df3$subject)] +
  subject_method_eff[cbind(as.integer(df3$subject), method_id)]
df3$y <- rpois(nrow(df3), exp(df3$eta))

fit3 <- ccc_glmm(
  df3, "y", "subject", "method",
  replicate = "replicate",
  include_subject_method = TRUE
)
attr(fit3, "sigma2_subject_method")


# Example 4: four methods.
# This tests pairwise agreement across several count methods and helps
# identify which method pairs have the strongest total CCC.
set.seed(4)
df4 <- expand.grid(
  subject = factor(seq_len(10)),
  method = factor(c("A", "B", "C", "D")),
  replicate = factor(seq_len(2))
)
subject_eff <- rnorm(10, 0, 0.5)
method_bias <- c(A = 0, B = 0.1, C = 0.4, D = -0.2)
df4$eta <- 1.0 + subject_eff[as.integer(df4$subject)] +
  method_bias[as.character(df4$method)]
df4$y <- rpois(nrow(df4), exp(df4$eta))
fit4 <- ccc_glmm(df4, "y", "subject", "method", replicate = "replicate")
fit4
summary(fit4, n = 3)

---

Concordance Correlation via REML (Linear Mixed-Effects Model)

Description

Compute Lin's Concordance Correlation Coefficient (CCC) from a linear mixed-effects model fitted by REML. The fixed-effects part can include method and/or time (optionally their interaction), with a subject-specific random intercept to capture between-subject variation. Large $n \times n$ inversions are avoided by solving small per-subject systems.

Assumption: time levels are treated as regular, equally spaced visits indexed by their order within subject. The AR(1) residual model is in discrete time on the visit index (not calendar time). NA time codes break the serial run. Gaps in the factor levels are ignored (adjacent observed visits are treated as lag-1).

Usage

ccc_rm_reml(
  data,
  response,
  subject,
  method = NULL,
  time = NULL,
  ci = FALSE,
  conf_level = 0.95,
  n_threads = getOption("matrixCorr.threads", 1L),
  ci_mode = c("auto", "raw", "logit"),
  verbose = FALSE,
  digits = 4,
  use_message = TRUE,
  interaction = FALSE,
  max_iter = 100,
  tol = 1e-06,
  Dmat = NULL,
  Dmat_type = c("time-avg", "typical-visit", "weighted-avg", "weighted-sq"),
  Dmat_weights = NULL,
  Dmat_rescale = TRUE,
  ar = c("none", "ar1"),
  ar_rho = NA_real_,
  slope = c("none", "subject", "method", "custom"),
  slope_var = NULL,
  slope_Z = NULL,
  drop_zero_cols = TRUE,
  vc_select = c("auto", "none"),
  vc_alpha = 0.05,
  vc_test_order = c("subj_time", "subj_method"),
  include_subj_method = NULL,
  include_subj_time = NULL,
  sb_zero_tol = 1e-10
)

Arguments

data	A data frame.
response	Character. Response variable name.
subject	Character. Subject ID variable name.
method	Character or NULL. Optional column name of method factor (added to fixed effects).
time	Character or NULL. Optional column name of time factor (added to fixed effects).
ci	Logical. If TRUE, return a CI container; limits are computed by a large-sample delta method for CCC (see CIs note below).
conf_level	Numeric in $(0,1)$. Confidence level when ci = TRUE (default 0.95).
n_threads	Integer $\geq 1$. Number of OpenMP threads to use for computation. Defaults to getOption("matrixCorr.threads", 1L).
ci_mode	Character scalar; one of c("auto","raw","logit"). Controls how confidence intervals are computed when ci = TRUE. If "raw", a Wald CI is formed on the CCC scale and truncated to [0,1]. If "logit", a Wald CI is computed on the $\mathrm{logit}(\mathrm{CCC})$ scale and back-transformed to the original scale (often more stable near 0 or 1). If "auto" (default), the method is chosen per estimate based on simple diagnostics (e.g., proximity to the [0,1] boundary / numerical stability), typically preferring "logit" near the boundaries and "raw" otherwise.
verbose	Logical. If TRUE, prints a structured summary of the fitted variance components and $S_B$ for each fit. Default FALSE.
digits	Integer $(\ge 0)$. Number of decimal places to use in the printed summary when verbose = TRUE. Default 4.
use_message	Logical. When verbose = TRUE, verbose summaries are emitted with cli messages.
interaction	Logical. Include method:time interaction? (default FALSE).
max_iter	Integer. Maximum iterations for variance-component updates (default 100).
tol	Numeric. Convergence tolerance on parameter change (default 1e-6).
Dmat	Optional $n_t \times n_t$ numeric matrix to weight/aggregate time-specific fixed biases in the $S_B$ quadratic form. If supplied, it is used (after optional mass rescaling; see Dmat_rescale) whenever at least two present time levels exist; otherwise it is ignored. If Dmat is NULL, a canonical kernel $D_m$ is constructed from Dmat_type and Dmat_weights (see below). Dmat should be symmetric positive semidefinite; small asymmetries are symmetrized internally.
Dmat_type	Character, one of c("time-avg","typical-visit", "weighted-avg","weighted-sq"). Only used when Dmat = NULL. It selects the aggregation target for time-specific fixed biases in $S_B$. Options are: "time-avg": square of the time-averaged bias, $D_m=(1/n_t)\,11^\top$. "typical-visit": average of squared per-time biases, $D_m=I_{n_t}$. "weighted-avg": square of a weighted average, $D_m=n_t\,w\,w^\top$ with $\sum w=1$. "weighted-sq": weighted average of squared biases, $D_m=n_t\,\mathrm{diag}(w)$ with $\sum w=1$. Pick "time-avg" for CCC targeting the time-averaged measurement; pick "typical-visit" for CCC targeting a randomly sampled visit (typical occasion). Default "time-avg".
Dmat_weights	Optional numeric weights $w$ used when Dmat_type %in% c("weighted-avg","weighted-sq"). Must be nonnegative and finite. If names(w) are provided, they should match the full time levels in data; they are aligned to the present time subset per fit. If unnamed, the length must equal the number of present time levels. In all cases $w$ is internally normalized to sum to 1.
Dmat_rescale	Logical. When TRUE (default), the supplied/built $D_m$ is rescaled to satisfy the simple mass rule $1^\top D_m 1 = n_t$. This keeps the $S_B$ denominator invariant and harmonizes with the $\kappa$-shrinkage used for variance terms.
ar	Character. Residual correlation structure: "none" (iid) or "ar1" for subject-level AR(1) correlation within contiguous time runs. Default c("none").
ar_rho	Numeric in $(-0.999,\,0.999)$ or NA. If ar = "ar1" and ar_rho is finite, it is treated as fixed. If ar = "ar1" and ar_rho = NA, $\rho$ is estimated by profiling a 1-D objective (REML when available; an approximation otherwise). Default NA_real_.
slope	Character. Optional extra random-effect design $Z$. With "subject" a single random slope is added (one column in $Z$); with "method" one column per method level is added; with "custom" you provide slope_Z directly. Default c("none","subject","method","custom").
slope_var	For slope %in% c("subject","method"), a character string giving the name of a column in data used as the slope regressor (e.g., centered time). It is looked up inside data; do not pass the vector itself. NAs are treated as zeros in $Z$.
slope_Z	For slope = "custom", a numeric matrix with $n$ rows (same order as data) providing the full extra random-effect design $Z$. Each column of slope_Z has its own variance component $\sigma_{Z,j}^2$; columns are treated as uncorrelated (diagonal block in $G$). Ignored otherwise.
drop_zero_cols	Logical. When slope = "method", drop all-zero columns of $Z$ after subsetting (useful in pairwise fits). Default TRUE.
vc_select	Character scalar; one of c("auto","none"). Controls how the subject by method $\sigma^2_{A\times M}$ ("subj_method") and subject by time $\sigma^2_{A\times T}$ ("subj_time") variance components are included. If "auto" (default), the function performs boundary-aware REML likelihood-ratio tests (LRTs; null on the boundary at zero with a half-$\chi^2_1$ reference) to decide whether to retain each component, in the order given by vc_test_order. If "none", no testing is done and inclusion is taken from include_subj_method/include_subj_time (or, if NULL, from the mere presence of the corresponding factor in the design). In pairwise fits, the decision is made independently for each method pair.
vc_alpha	Numeric scalar in $(0,1)$; default 0.05. Per-component significance level for the boundary-aware REML LRTs used when vc_select = "auto". The tests are one-sided for variance components on the boundary and are not multiplicity-adjusted.
vc_test_order	Character vector (length 2) with a permutation of c("subj_time","subj_method"); default c("subj_time","subj_method"). Specifies the order in which the two variance components are tested when vc_select = "auto". The component tested first may be dropped before testing the second. If a factor is absent in the design (e.g., no time factor so "subj_time" is undefined), the corresponding test is skipped.
include_subj_method, include_subj_time	Logical scalars or NULL. When vc_select = "none", these control whether the $\sigma^2_{A\times M}$ ("subj_method") and $\sigma^2_{A\times T}$ ("subj_time") random effects are included (TRUE) or excluded (FALSE) in the model. If NULL (default), inclusion defaults to the presence of the corresponding factor in the data (i.e., at least two method/time levels). When vc_select = "auto", these arguments are ignored (automatic selection is used instead).
sb_zero_tol	Non-negative numeric scalar; default 1e-10. Numerical threshold for the fixed-effect dispersion term $S_B$. After computing $\widehat{S_B}$ and its delta-method variance, if $\widehat{S_B} \le$ sb_zero_tol or non-finite, the procedure treats $S_B$ as fixed at zero in the delta step. It sets $d_{S_B}=0$ and $\mathrm{Var}(\widehat{S_B})=0$, preventing numerical blow-ups of SE(CCC) when $\widehat{S_B}\to 0$ and the fixed-effects variance is ill-conditioned for the contrast. This stabilises inference in rare boundary cases; it has no effect when $\widehat{S_B}$ is comfortably above the threshold.

Details

For measurement $y_{ij}$ on subject $i$ under fixed levels (method, time), we fit

$y = X\beta + Zu + \varepsilon,\qquad u \sim N(0,\,G),\ \varepsilon \sim N(0,\,R).$

Notation: $m$ subjects, $n=\sum_i n_i$ total rows; $nm$ method levels; $nt$ time levels; $q_Z$ extra random-slope columns (if any); $r=1+nm+nt$ (or $1+nm+nt+q_Z$ with slopes). Here $Z$ is the subject-structured random-effects design and $G$ is block-diagonal at the subject level with the following per-subject parameterisation. Specifically,

• one random intercept with variance $\sigma_A^2$;

• optionally, method deviations (one column per method level) with a common variance $\sigma_{A\times M}^2$ and zero covariances across levels (i.e., multiple of an identity);

• optionally, time deviations (one column per time level) with a common variance $\sigma_{A\times T}^2$ and zero covariances across levels;

• optionally, an extra random effect aligned with $Z$ (random slope), where each column has its own variance $\sigma_{Z,j}^2$ and columns are uncorrelated.

The fixed-effects design is ~ 1 + method + time and, if interaction=TRUE, + method:time.

Residual correlation $R$ (regular, equally spaced time). Write $R_i=\sigma_E^2\,C_i(\rho)$. With ar="none", $C_i=I$. With ar="ar1", within-subject residuals follow a discrete AR(1) process along the visit index after sorting by increasing time level. Ties retain input order, and any NA time code breaks the series so each contiguous block of non-NA times forms a run. The correlation between adjacent observed visits in a run is $\rho$; we do not use calendar-time gaps. Internally we work with the precision of the AR(1) correlation: for a run of length $L\ge 2$, the tridiagonal inverse has

$(C^{-1})_{11}=(C^{-1})_{LL}=\frac{1}{1-\rho^2},\quad (C^{-1})_{tt}=\frac{1+\rho^2}{1-\rho^2}\ (2\le t\le L-1),\quad (C^{-1})_{t,t+1}=(C^{-1})_{t+1,t}=\frac{-\rho}{1-\rho^2}.$

The working inverse is $\,R_i^{-1}=\sigma_E^{-2}\,C_i(\rho)^{-1}$.

Per-subject Woodbury system. For subject $i$ with $n_i$ rows, define the per-subject random-effects design $U_i$ (columns: intercept, method indicators, time indicators; dimension $\,r=1+nm+nt\,$). The core never forms $V_i = R_i + U_i G U_i^\top$ explicitly. Instead,

$M_i \;=\; G^{-1} \;+\; U_i^\top R_i^{-1} U_i,$

and accumulates GLS blocks via rank-$r$ corrections using $\,V_i^{-1} = R_i^{-1}-R_i^{-1}U_i M_i^{-1}U_i^\top R_i^{-1}\,$:

$X^\top V^{-1} X \;=\; \sum_i \Big[ X_i^\top R_i^{-1}X_i \;-\; (X_i^\top R_i^{-1}U_i)\, M_i^{-1}\,(U_i^\top R_i^{-1}X_i) \Big],$

$X^\top V^{-1} y \;=\; \sum_i \Big[ X_i^\top R_i^{-1}y_i \;-\; (X_i^\top R_i^{-1}U_i)\,M_i^{-1}\, (U_i^\top R_i^{-1}y_i) \Big].$

Because $G^{-1}$ is diagonal with positive entries, each $M_i$ is symmetric positive definite; solves/inversions use symmetric-PD routines with a small diagonal ridge and a pseudo-inverse if needed.

Random-slope $Z$. Besides $U_i$, the function can include an extra design $Z_i$.

• slope="subject": $Z$ has one column (the regressor in slope_var); $Z_{i}$ is the subject-$i$ block, with its own variance $\sigma_{Z,1}^2$.

• slope="method": $Z$ has one column per method level; row $t$ uses the slope regressor if its method equals level $\ell$, otherwise 0; all-zero columns can be dropped via drop_zero_cols=TRUE after subsetting. Each column has its own variance $\sigma_{Z,\ell}^2$.

• slope="custom": $Z$ is provided fully via slope_Z. Each column is an independent random effect with its own variance $\sigma_{Z,j}^2$; cross-covariances among columns are set to 0.

Computations simply augment $\tilde U_i=[U_i\ Z_i]$ and the corresponding inverse-variance block. The EM updates then include, for each column $j=1,\ldots,q_Z$,

$\sigma_{Z,j}^{2\,(new)} \;=\; \frac{1}{m} \sum_{i=1}^m \Big( b_{i,\text{extra},j}^2 + (M_i^{-1})_{\text{extra},jj} \Big) \quad (\text{if } q_Z>0).$

Interpretation: the $\sigma_{Z,j}^2$ represent additional within-subject variability explained by the slope regressor(s) in column $j$ and are not part of the CCC denominator (agreement across methods/time).

EM-style variance-component updates. With current $\hat\beta$, form residuals $r_i = y_i - X_i\hat\beta$. The BLUPs and conditional covariances are

$b_i \;=\; M_i^{-1}\,(U_i^\top R_i^{-1} r_i), \qquad \mathrm{Var}(b_i\mid y) \;=\; M_i^{-1}.$

Let $e_i=r_i-U_i b_i$. Expected squares then yield closed-form updates:

$\sigma_A^{2\,(new)} \;=\; \frac{1}{m}\sum_i \Big( b_{i,0}^2 + (M_i^{-1})_{00} \Big),$

$\sigma_{A\times M}^{2\,(new)} \;=\; \frac{1}{m\,nm} \sum_i \sum_{\ell=1}^{nm}\!\Big( b_{i,\ell}^2 + (M_i^{-1})_{\ell\ell} \Big) \quad (\text{if } nm>0),$

$\sigma_{A\times T}^{2\,(new)} \;=\; \frac{1}{m\,nt} \sum_i \sum_{t=1}^{nt} \Big( b_{i,t}^2 + (M_i^{-1})_{tt} \Big) \quad (\text{if } nt>0),$

$\sigma_E^{2\,(new)} \;=\; \frac{1}{n} \sum_i \Big( e_i^\top C_i(\rho)^{-1} e_i \;+\; \mathrm{tr}\!\big(M_i^{-1}U_i^\top C_i(\rho)^{-1} U_i\big) \Big),$

together with the per-column update for $\sigma_{Z,j}^2$ given above. Iterate until the $\ell_1$ change across components is $<$tol or max_iter is reached.

Fixed-effect dispersion $S_B$: choosing the time-kernel $D_m$.

Let $d = L^\top \hat\beta$ stack the within-time, pairwise method differences, grouped by time as $d=(d_1^\top,\ldots,d_{n_t}^\top)^\top$ with $d_t \in \mathbb{R}^{P}$ and $P = n_m(n_m-1)$. The symmetric positive semidefinite kernel $D_m \succeq 0$ selects which functional of the bias profile $t \mapsto d_t$ is targeted by $S_B$. Internally, the code rescales any supplied/built $D_m$ to satisfy $1^\top D_m 1 = n_t$ for stability and comparability.

• Dmat_type = "time-avg" (square of the time-averaged bias). Let

  $w \;=\; \frac{1}{n_t}\,\mathbf{1}_{n_t}, \qquad D_m \;\propto\; I_P \otimes (w\,w^\top),$

  so that

  $d^\top D_m d \;\propto\; \sum_{p=1}^{P} \left( \frac{1}{n_t}\sum_{t=1}^{n_t} d_{t,p} \right)^{\!2}.$

  Methods have equal $\textit{time-averaged}$ means within subject, i.e. $\sum_{t=1}^{n_t} d_{t,p}/n_t = 0$ for all $p$. Appropriate when decisions depend on an average over time and opposite-signed biases are allowed to cancel.

• Dmat_type = "typical-visit" (average of squared per-time biases). With equal visit probability, take

  $D_m \;\propto\; I_P \otimes \mathrm{diag}\!\Big(\tfrac{1}{n_t},\ldots,\tfrac{1}{n_t}\Big),$

  yielding

  $d^\top D_m d \;\propto\; \frac{1}{n_t} \sum_{t=1}^{n_t}\sum_{p=1}^{P} d_{t,p}^{\,2}.$

  Methods agree on a $\textit{typical}$ occasion drawn uniformly from the visit set. Use when each visit matters on its own; alternating signs $d_{t,p}$ do not cancel.

• Dmat_type = "weighted-avg" (square of a weighted time average). For user weights $a=(a_1,\ldots,a_{n_t})^\top$ with $a_t \ge 0$, set

  $w \;=\; \frac{a}{\sum_{t=1}^{n_t} a_t}, \qquad D_m \;\propto\; I_P \otimes (w\,w^\top),$

  so that

  $d^\top D_m d \;\propto\; \sum_{p=1}^{P} \left( \sum_{t=1}^{n_t} w_t\, d_{t,p} \right)^{\!2}.$

  Methods have equal $\textit{weighted time-averaged}$ means, i.e. $\sum_{t=1}^{n_t} w_t\, d_{t,p} = 0$ for all $p$. Use when some visits (e.g., baseline/harvest) are a priori more influential; opposite-signed biases may cancel according to $w$.

• Dmat_type = "weighted-sq" (weighted average of squared per-time biases). With the same weights $w$, take

  $D_m \;\propto\; I_P \otimes \mathrm{diag}(w_1,\ldots,w_{n_t}),$

  giving

  $d^\top D_m d \;\propto\; \sum_{t=1}^{n_t} w_t \sum_{p=1}^{P} d_{t,p}^{\,2}.$

  Methods agree at visits sampled with probabilities $\{w_t\}$, counting each visit's discrepancy on its own. Use when per-visit agreement is required but some visits should be emphasised more than others.

Time-averaging for CCC (regular visits). The reported CCC targets agreement of the time-averaged measurements per method within subject by default (Dmat_type="time-avg"). Averaging over $T$ non-NA visits shrinks time-varying components by

$\kappa_T^{(g)} \;=\; 1/T, \qquad \kappa_T^{(e)} \;=\; \{T + 2\sum_{k=1}^{T-1}(T-k)\rho^k\}/T^2,$

with $\kappa_T^{(e)}=1/T$ when residuals are i.i.d. With unbalanced $T$, the implementation averages the per-(subject,method) $\kappa$ values across the pairs contributing to CCC and then clamps $\kappa_T^{(e)}$ to $[10^{-12},\,1]$ for numerical stability. Choosing Dmat_type="typical-visit" makes $S_B$ match the interpretation of a randomly sampled occasion instead.

Concordance correlation coefficient. The CCC used is

$\mathrm{CCC} \;=\; \frac{\sigma_A^2 + \kappa_T^{(g)}\,\sigma_{A\times T}^2} {\sigma_A^2 + \sigma_{A\times M}^2 + \kappa_T^{(g)}\,\sigma_{A\times T}^2 + S_B + \kappa_T^{(e)}\,\sigma_E^2}.$

Special cases: with no method factor, $S_B=\sigma_{A\times M}^2=0$; with a single time level, $\sigma_{A\times T}^2=0$ (no $\kappa$-shrinkage). When $T=1$ or $\rho=0$, both $\kappa$-factors equal 1. The extra random-effect variances $\{\sigma_{Z,j}^2\}$ (if used) are not included.

CIs / SEs (delta method for CCC). Let

$\theta \;=\; \big(\sigma_A^2,\ \sigma_{A\times M}^2,\ \sigma_{A\times T}^2,\ \sigma_E^2,\ S_B\big)^\top,$

and write $\mathrm{CCC}(\theta)=N/D$ with $N=\sigma_A^2+\kappa_T^{(g)}\sigma_{A\times T}^2$ and $D=\sigma_A^2+\sigma_{A\times M}^2+\kappa_T^{(g)}\sigma_{A\times T}^2+S_B+\kappa_T^{(e)}\sigma_E^2$. The gradient components are

$\frac{\partial\,\mathrm{CCC}}{\partial \sigma_A^2} \;=\; \frac{\sigma_{A\times M}^2 + S_B + \kappa_T^{(e)}\sigma_E^2}{D^2},$

$\frac{\partial\,\mathrm{CCC}}{\partial \sigma_{A\times M}^2} \;=\; -\,\frac{N}{D^2}, \qquad \frac{\partial\,\mathrm{CCC}}{\partial \sigma_{A\times T}^2} \;=\; \frac{\kappa_T^{(g)}\big(\sigma_{A\times M}^2 + S_B + \kappa_T^{(e)}\sigma_E^2\big)}{D^2},$

$\frac{\partial\,\mathrm{CCC}}{\partial \sigma_E^2} \;=\; -\,\frac{\kappa_T^{(e)}\,N}{D^2}, \qquad \frac{\partial\,\mathrm{CCC}}{\partial S_B} \;=\; -\,\frac{N}{D^2}.$

Estimating $\mathrm{Var}(\hat\theta)$. The EM updates write each variance component as an average of per-subject quantities. For subject $i$,

$t_{A,i} \;=\; b_{i,0}^2 + (M_i^{-1})_{00},\qquad t_{M,i} \;=\; \frac{1}{nm}\sum_{\ell=1}^{nm} \Big(b_{i,\ell}^2 + (M_i^{-1})_{\ell\ell}\Big),$

$t_{T,i} \;=\; \frac{1}{nt}\sum_{j=1}^{nt} \Big(b_{i,j}^2 + (M_i^{-1})_{jj}\Big),\qquad s_i \;=\; \frac{e_i^\top C_i(\rho)^{-1} e_i + \mathrm{tr}\!\big(M_i^{-1}U_i^\top C_i(\rho)^{-1} U_i\big)}{n_i},$

where $b_i = M_i^{-1}(U_i^\top R_i^{-1} r_i)$ and $M_i = G^{-1} + U_i^\top R_i^{-1} U_i$. With $m$ subjects, form the empirical covariance of the stacked subject vectors and scale by $m$ to approximate the covariance of the means:

$\widehat{\mathrm{Cov}}\!\left( \begin{bmatrix} t_{A,\cdot} \\ t_{M,\cdot} \\ t_{T,\cdot} \end{bmatrix} \right) \;\approx\; \frac{1}{m}\, \mathrm{Cov}_i\!\left( \begin{bmatrix} t_{A,i} \\ t_{M,i} \\ t_{T,i} \end{bmatrix}\right).$

(Drop rows/columns as needed when nm==0 or nt==0.)

The residual variance estimator is a weighted mean $\hat\sigma_E^2=\sum_i w_i s_i$ with $w_i=n_i/n$. Its variance is approximated by the variance of a weighted mean of independent terms,

$\widehat{\mathrm{Var}}(\hat\sigma_E^2) \;\approx\; \Big(\sum_i w_i^2\Big)\,\widehat{\mathrm{Var}}(s_i),$

where $\widehat{\mathrm{Var}}(s_i)$ is the sample variance across subjects. The method-dispersion term uses the quadratic-form delta already computed for $S_B$:

$\widehat{\mathrm{Var}}(S_B) \;=\; \frac{2\,\mathrm{tr}\!\big((A_{\!fix}\,\mathrm{Var}(\hat\beta))^2\big) \;+\; 4\,\hat\beta^\top A_{\!fix}\,\mathrm{Var}(\hat\beta)\, A_{\!fix}\,\hat\beta} {\big[nm\,(nm-1)\,\max(nt,1)\big]^2},$

with $A_{\!fix}=L\,D_m\,L^\top$.

Putting it together. Assemble $\widehat{\mathrm{Var}}(\hat\theta)$ by combining the $(\sigma_A^2,\sigma_{A\times M}^2,\sigma_{A\times T}^2)$ covariance block from the subject-level empirical covariance, add the $\widehat{\mathrm{Var}}(\hat\sigma_E^2)$ and $\widehat{\mathrm{Var}}(S_B)$ terms on the diagonal, and ignore cross-covariances across these blocks (a standard large-sample simplification). Then

$\widehat{\mathrm{se}}\{\widehat{\mathrm{CCC}}\} \;=\; \sqrt{\,\nabla \mathrm{CCC}(\hat\theta)^\top\, \widehat{\mathrm{Var}}(\hat\theta)\, \nabla \mathrm{CCC}(\hat\theta)\,}.$

A two-sided $(1-\alpha)$ normal CI is

$\widehat{\mathrm{CCC}} \;\pm\; z_{1-\alpha/2}\, \widehat{\mathrm{se}}\{\widehat{\mathrm{CCC}}\},$

truncated to $[0,1]$ in the output for convenience. When $S_B$ is truncated at 0 or samples are very small/imbalanced, the normal CI can be mildly anti-conservative near the boundary; a logit transform for CCC or a subject-level (cluster) bootstrap can be used for sensitivity analysis.

Choosing $\rho$ for AR(1). When ar="ar1" and ar_rho = NA, $\rho$ is estimated by profiling the REML log-likelihood at $(\hat\beta,\hat G,\hat\sigma_E^2)$. With very few visits per subject, $\rho$ can be weakly identified; consider sensitivity checks over a plausible range.

Notes on stability and performance

All per-subject solves are $\,r\times r$ with $r=1+nm+nt+q_Z$, so cost scales with the number of subjects and the fixed-effects dimension rather than the total number of observations. Solvers use symmetric-PD paths with a small diagonal ridge and pseudo-inverse, which helps for very small/unbalanced subsets and near-boundary estimates. For AR(1), observations are ordered by time within subject; NA time codes break the run, and gaps between factor levels are treated as regular steps (elapsed time is not used).

Heteroscedastic slopes across $Z$ columns are supported. Each $Z$ column has its own variance component $\sigma_{Z,j}^2$, but cross-covariances among $Z$ columns are set to zero (diagonal block). Column rescaling changes the implied prior on $b_{i,\text{extra}}$ but does not introduce correlations.

Threading and BLAS guards

The C++ backend uses OpenMP loops while also forcing vendor BLAS libraries to run single-threaded so that overall CPU usage stays predictable. This guard is applied to OpenBLAS, Apple's Accelerate, and Intel MKL when their runtime controls are available. You can opt out manually by setting MATRIXCORR_DISABLE_BLAS_GUARD=1 in the environment before loading matrixCorr.

Model overview

Internally, the call is routed to ccc_lmm_reml_pairwise(), which fits one repeated-measures mixed model per pair of methods. Each model includes:

• subject random intercepts (always)

• optional subject-by-method (⁠sigma^2_{A \times M}⁠) and subject-by-time (⁠sigma^2_{A \times T}⁠) variance components

• optional random slopes specified via slope/slope_var/slope_Z

• residual structure ar = "none" (iid) or ar = "ar1"

D-matrix options (Dmat_type, Dmat, Dmat_weights) control how time averaging operates when translating variance components into CCC summaries.

Author(s)

Thiago de Paula Oliveira

References

Lin L (1989). A concordance correlation coefficient to evaluate reproducibility. Biometrics, 45: 255-268.

Lin L (2000). A note on the concordance correlation coefficient. Biometrics, 56: 324-325.

Carrasco, J. L. et al. (2013). Estimation of the concordance correlation coefficient for repeated measures using SAS and R. Computer Methods and Programs in Biomedicine, 109(3), 293-304. doi:10.1016/j.cmpb.2012.09.002

King et al. (2007). A Class of Repeated Measures Concordance Correlation Coefficients. Journal of Biopharmaceutical Statistics, 17(4). doi:10.1080/10543400701329455

See Also

build_L_Dm_Z_cpp for constructing $L$/$D_m$/$Z$; ccc_rm_ustat for a U-statistic alternative; and cccrm for a reference approach via nlme.

Examples

# ====================================================================
# 1) Subject x METHOD variance present, no time
#     y_{i,m} = mu + b_m + u_i + w_{i,m} + e_{i,m}
#     with u_i ~ N(0, s_A^2), w_{i,m} ~ N(0, s_{AxM}^2)
# ====================================================================
set.seed(102)
n_subj <- 60
n_time <- 8

id     <- factor(rep(seq_len(n_subj), each = 2 * n_time))
time   <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))
method <- factor(rep(rep(c("A","B"),    each  = n_time), times = n_subj))

sigA  <- 0.6   # subject
sigAM <- 0.3   # subject x method
sigAT <- 0.5   # subject x time
sigE  <- 0.4   # residual
# Expected estimate S_B = 0.2^2 = 0.04
biasB <- 0.2   # fixed method bias

# random effects
u_i <- rnorm(n_subj, 0, sqrt(sigA))
u   <- u_i[as.integer(id)]

sm      <- interaction(id, method, drop = TRUE)
w_im_lv <- rnorm(nlevels(sm), 0, sqrt(sigAM))
w_im    <- w_im_lv[as.integer(sm)]

st      <- interaction(id, time, drop = TRUE)
g_it_lv <- rnorm(nlevels(st), 0, sqrt(sigAT))
g_it    <- g_it_lv[as.integer(st)]

# residuals & response
e <- rnorm(length(id), 0, sqrt(sigE))
y <- (method == "B") * biasB + u + w_im + g_it + e

dat_both <- data.frame(y, id, method, time)

# Both sigma2_subject_method and sigma2_subject_time are identifiable here
fit_both <- ccc_rm_reml(dat_both, "y", "id", method = "method", time = "time",
                         vc_select = "auto", ci = TRUE, verbose = TRUE)
summary(fit_both)
estimate(fit_both)
tidy(fit_both)
ci(fit_both)
confint(fit_both)
plot(fit_both)

# Interactive viewing (requires shiny)
if (interactive() && requireNamespace("shiny", quietly = TRUE)) {
  view_corr_shiny(fit_both)
}

# ====================================================================
# 2) Subject x TIME variance present (sag > 0) with two methods
#     y_{i,m,t} = mu + b_m + u_i + g_{i,t} + e_{i,m,t}
#     where g_{i,t} ~ N(0, s_{AxT}^2) shared across methods at time t
# ====================================================================
set.seed(202)
n_subj <- 60; n_time <- 14
id     <- factor(rep(seq_len(n_subj), each = 2 * n_time))
method <- factor(rep(rep(c("A","B"), each = n_time), times = n_subj))
time   <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))

sigA  <- 0.7
sigAT <- 0.5
sigE  <- 0.5
biasB <- 0.25

u   <- rnorm(n_subj, 0, sqrt(sigA))[as.integer(id)]
gIT <- rnorm(n_subj * n_time, 0, sqrt(sigAT))
g   <- gIT[ (as.integer(id) - 1L) * n_time + as.integer(time) ]
y   <- (method == "B") * biasB + u + g + rnorm(length(id), 0, sqrt(sigE))
dat_sag <- data.frame(y, id, method, time)

# sigma_AT should be retained; sigma_AM may be dropped (since w_{i,m}=0)
fit_sag <- ccc_rm_reml(dat_sag, "y", "id", method = "method", time = "time",
                        vc_select = "auto", verbose = TRUE)
summary(fit_sag)
plot(fit_sag)

# ====================================================================
# 3) BOTH components present: sab > 0 and sag > 0 (2 methods x T times)
#     y_{i,m,t} = mu + b_m + u_i + w_{i,m} + g_{i,t} + e_{i,m,t}
# ====================================================================
set.seed(303)
n_subj <- 60; n_time <- 4
id     <- factor(rep(seq_len(n_subj), each = 2 * n_time))
method <- factor(rep(rep(c("A","B"), each = n_time), times = n_subj))
time   <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))

sigA  <- 0.8
sigAM <- 0.3
sigAT <- 0.4
sigE  <- 0.5
biasB <- 0.2

u   <- rnorm(n_subj, 0, sqrt(sigA))[as.integer(id)]
# (subject, method) random deviations: repeat per (i,m) across its times
wIM <- rnorm(n_subj * 2, 0, sqrt(sigAM))
w   <- wIM[ (as.integer(id) - 1L) * 2 + as.integer(method) ]
# (subject, time) random deviations: shared across methods at time t
gIT <- rnorm(n_subj * n_time, 0, sqrt(sigAT))
g   <- gIT[ (as.integer(id) - 1L) * n_time + as.integer(time) ]
y   <- (method == "B") * biasB + u + w + g + rnorm(length(id), 0, sqrt(sigE))
dat_both <- data.frame(y, id, method, time)

fit_both <- ccc_rm_reml(dat_both, "y", "id", method = "method", time = "time",
                         vc_select = "auto", verbose = TRUE, ci = TRUE)
summary(fit_both)
plot(fit_both)

# If you want to force-include both VCs (skip testing):
fit_both_forced <-
 ccc_rm_reml(dat_both, "y", "id", method = "method", time = "time",
              vc_select = "none", include_subj_method  = TRUE,
              include_subj_time  = TRUE, verbose = TRUE)
summary(fit_both_forced)
plot(fit_both_forced)

# ====================================================================
# 4) D_m choices: time-averaged (default) vs typical visit
# ====================================================================
# Time-average
ccc_rm_reml(dat_both, "y", "id", method = "method", time = "time",
             vc_select = "none", include_subj_method  = TRUE,
             include_subj_time  = TRUE, Dmat_type = "time-avg")

# Typical visit
ccc_rm_reml(dat_both, "y", "id", method = "method", time = "time",
             vc_select = "none", include_subj_method  = TRUE,
             include_subj_time  = TRUE, Dmat_type = "typical-visit")

# ====================================================================
# 5) AR(1) residual correlation with fixed rho (larger example)
# ====================================================================

set.seed(10)
n_subj   <- 40
n_time   <- 10
methods  <- c("A", "B", "C", "D")
nm       <- length(methods)
id     <- factor(rep(seq_len(n_subj), each = n_time * nm))
method <- factor(rep(rep(methods, each = n_time), times = n_subj),
                 levels = methods)
time   <- factor(rep(rep(seq_len(n_time), times = nm), times = n_subj))

beta0    <- 0
beta_t   <- 0.2
bias_met <- c(A = 0.00, B = 0.30, C = -0.15, D = 0.05)
sigA     <- 1.0
rho_true <- 0.6
sigE     <- 0.7

t_num <- as.integer(time)
t_c   <- t_num - mean(seq_len(n_time))
mu    <- beta0 + beta_t * t_c + bias_met[as.character(method)]

u_subj <- rnorm(n_subj, 0, sqrt(sigA))
u      <- u_subj[as.integer(id)]

e <- numeric(length(id))
for (s in seq_len(n_subj)) {
  for (m in methods) {
    idx <- which(id == levels(id)[s] & method == m)
    e[idx] <- stats::arima.sim(list(ar = rho_true), n = n_time, sd = sigE)
  }
}
y <- mu + u + e
dat_ar4 <- data.frame(y = y, id = id, method = method, time = time)

fit4 <- ccc_rm_reml(dat_ar4,
             response = "y", subject = "id", method = "method", time = "time",
             ar = "ar1", ar_rho = 0.6, verbose = TRUE)
fit4
summary(fit4)
plot(fit4)


# ====================================================================
# 6) Random slope variants (subject, method, custom Z)
# ====================================================================

## By SUBJECT
set.seed(2)
n_subj <- 60; n_time <- 4
id  <- factor(rep(seq_len(n_subj), each = 2 * n_time))
tim <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))
method <- factor(rep(rep(c("A","B"), each = n_time), times = n_subj))
subj <- as.integer(id)
slope_i <- rnorm(n_subj, 0, 0.15)
slope_vec <- slope_i[subj]
base <- rnorm(n_subj, 0, 1.0)[subj]
tnum <- as.integer(tim)
y <- base + 0.3*(method=="B") + slope_vec*(tnum - mean(seq_len(n_time))) +
     rnorm(length(id), 0, 0.5)
dat_s <- data.frame(y, id, method, time = tim)
dat_s$t_num <- as.integer(dat_s$time)
dat_s$t_c   <- ave(dat_s$t_num, dat_s$id, FUN = function(v) v - mean(v))
ccc_rm_reml(dat_s, "y", "id", method = "method", time = "time",
             slope = "subject", slope_var = "t_c", verbose = TRUE)

## By METHOD
set.seed(3)
n_subj <- 60; n_time <- 4
id  <- factor(rep(seq_len(n_subj), each = 2 * n_time))
tim <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))
method <- factor(rep(rep(c("A","B"), each = n_time), times = n_subj))
slope_m <- ifelse(method=="B", 0.25, 0.00)
base <- rnorm(n_subj, 0, 1.0)[as.integer(id)]
tnum <- as.integer(tim)
y <- base + 0.3*(method=="B") + slope_m*(tnum - mean(seq_len(n_time))) +
     rnorm(length(id), 0, 0.5)
dat_m <- data.frame(y, id, method, time = tim)
dat_m$t_num <- as.integer(dat_m$time)
dat_m$t_c   <- ave(dat_m$t_num, dat_m$id, FUN = function(v) v - mean(v))
ccc_rm_reml(dat_m, "y", "id", method = "method", time = "time",
             slope = "method", slope_var = "t_c", verbose = TRUE)

## SUBJECT + METHOD random slopes (custom Z)
set.seed(4)
n_subj <- 50; n_time <- 4
id  <- factor(rep(seq_len(n_subj), each = 2 * n_time))
tim <- factor(rep(rep(seq_len(n_time), times = 2), times = n_subj))
method <- factor(rep(rep(c("A","B"), each = n_time), times = n_subj))
subj <- as.integer(id)
slope_subj <- rnorm(n_subj, 0, 0.12)[subj]
slope_B    <- ifelse(method=="B", 0.18, 0.00)
tnum <- as.integer(tim)
base <- rnorm(n_subj, 0, 1.0)[subj]
y <- base + 0.3*(method=="B") +
     (slope_subj + slope_B) * (tnum - mean(seq_len(n_time))) +
     rnorm(length(id), 0, 0.5)
dat_bothRS <- data.frame(y, id, method, time = tim)
dat_bothRS$t_num <- as.integer(dat_bothRS$time)
dat_bothRS$t_c   <- ave(dat_bothRS$t_num, dat_bothRS$id, FUN = function(v) v - mean(v))
MM <- model.matrix(~ 0 + method, data = dat_bothRS)
Z_custom <- cbind(
  subj_slope = dat_bothRS$t_c,
  MM * dat_bothRS$t_c
)
ccc_rm_reml(dat_bothRS, "y", "id", method = "method", time = "time",
             slope = "custom", slope_Z = Z_custom, verbose = TRUE)

---