cia_rm: Repeated-Measures Coefficient of Individual Agreement

Description

Computes pairwise repeated-measures coefficients of individual agreement (CIA) from long-format matched repeated-measures data using the categorical condition ANOVA formulation of Haber, Gao, and Barnhart (2010).

Usage

cia_rm(
  data,
  response,
  subject,
  method = NULL,
  time = NULL,
  ci = FALSE,
  conf_level = 0.95,
  n_threads = getOption("matrixCorr.threads", 1L),
  estimator = c("vc_constrained", "mom_unconstrained"),
  inference = c("delta", "bootstrap", "none"),
  verbose = FALSE,
  digits = 4,
  use_message = TRUE,
  homogeneous = FALSE,
  B = 1000L,
  seed = NULL,
  ...
)
# S3 method for cia_rm
summary(
  object,
  digits = 4,
  ci_digits = 3,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for cia_rm
print(
  x,
  digits = 4,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for summary.cia_rm
print(
  x,
  digits = NULL,
  n = NULL,
  topn = NULL,
  max_vars = NULL,
  width = NULL,
  show_ci = NULL,
  ...
)
# S3 method for cia_rm
plot(
  x,
  title = NULL,
  facet_by_pair = FALSE,
  facet_scales = c("fixed", "free_y"),
  show_ci = NULL,
  show_common = TRUE,
  point_size = 2.2,
  line_size = 0.7,
  ci_alpha = 0.16,
  ci_linewidth = 0.45,
  ...
)

Value

If ci = FALSE, the result is a list of class c("cia_rm", "cia"). If ci = TRUE, the result is a list of class c("cia_rm", "cia_ci", "cia").

Main components:

est: condition-specific pairwise CIA estimates.
common: pairwise overall summary CIA estimates from the reduced homogeneous model.
condition: the ordered condition/time labels used in the output.
se: standard errors for est when ci = TRUE.
lwr.ci and upr.ci: lower and upper confidence limits for est when ci = TRUE.
common.se: standard errors for common when ci = TRUE.
common.lwr.ci and common.upr.ci: lower and upper confidence limits for common when ci = TRUE.

When there are three or more methods, additional overall multi-method components are returned:

overall: condition-specific overall CIA across all methods.
overall.common: the overall summary CIA across all methods when homogeneous = TRUE.
overall.se: standard errors for overall when ci = TRUE.
overall.lwr.ci and overall.upr.ci: lower and upper confidence limits for overall when ci = TRUE.
overall.common.se: standard errors for overall.common when ci = TRUE.
overall.common.lwr.ci and overall.common.upr.ci: lower and upper confidence limits for overall.common when ci = TRUE.

The object carries pairwise ANOVA diagnostics as attributes, including sigma2_error, sigma2_subject_method, repeatability, homogeneity_F, homogeneity_p, df_method_time, df_error, n_obs, n_subjects, n_methods, n_times, time_levels, and method_levels. Here homogeneity_F is the method-time interaction test statistic for each method pair, and homogeneity_p is the corresponding p-value for the null hypothesis of homogeneous agreement across conditions.

Arguments

data: A data frame in long format.
response: Character scalar naming the numeric response column.
subject: Character scalar naming the subject identifier column.
method: Character scalar naming the method column. Required.
time: Character scalar naming the repeated condition column. Required.
ci: Logical; if TRUE, attach confidence intervals.
conf_level: Confidence level for intervals when ci = TRUE.
n_threads: Positive integer thread hint passed to the C++ backend.
estimator: One of "vc_constrained" or "mom_unconstrained". "vc_constrained" is the default bounded variance-component estimator. It constrains the estimated subject-method variance component to be non-negative before converting it back to CIA. "mom_unconstrained" is the raw method-of-moments CIA ratio and may exceed 1.
inference: One of "delta", "bootstrap", or "none". The delta method is the default when ci = TRUE.
verbose: Logical; if TRUE, emit progress messages for CI work.
digits: Integer print precision carried on the returned object.
use_message: Logical; if TRUE, use package messaging helpers when verbose = TRUE.
homogeneous: Logical; if TRUE, also compute the pooled/common CIA from the reduced model without method-time interaction.
B: Number of subject-bootstrap resamples when ci = TRUE.
seed: Optional positive integer seed for reproducible bootstrap resampling.
...: Additional arguments passed to ggplot2::theme().
object: A cia_rm object.
ci_digits: Integer; number of digits for confidence interval bounds.
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.
width: Optional display width.
show_ci: Logical; if TRUE and confidence intervals are present, show them on the plot. By default this is inferred from the object class.
x: A summary.cia_rm object.
title: Optional plot title.
facet_by_pair: Logical; if TRUE, draw one panel per method pair.
facet_scales: Passed to ggplot2::facet_wrap(scales = ...) when facet_by_pair = TRUE. One of "fixed" or "free_y".
show_common: Logical; if TRUE and common estimates are available, draw dashed horizontal reference lines for the pooled/common CIA.
point_size: Numeric point size passed to ggplot2::geom_point().
line_size: Numeric line width passed to ggplot2::geom_line().
ci_alpha: Numeric alpha used for confidence ribbons in continuous plots.
ci_linewidth: Numeric line width used for confidence-interval error bars in categorical plots.

Author

Thiago de Paula Oliveira

Details

cia_rm() is for matched repeated measurements under conditions such as raters, visits, laboratories, treatments, or time points. It is not a technical-replicate estimator. If the same subject has true technical replicates within each subject-method cell, use cia() instead.

Let Y_ijk denote the measurement on subject i, with method j, under condition k, and consider the model $$ Y_{ijk} = \mu + \alpha_i + \beta_j + \gamma_k + (\alpha\beta)_{ij} + (\alpha\gamma)_{ik} + (\beta\gamma)_{jk} + e_{ijk}. $$ Here subject is random, method and time/condition are fixed, subject-method and subject-time are random interactions, method-time is a fixed interaction, and e_ijk is the residual.

For each method pair and condition k, the repeated-measures CIA is $$ \psi_k = \frac{2\sigma_e^2}{d_k^2 + 2\sigma_{\alpha\beta}^2 + 2\sigma_e^2}, $$ where d_k is the condition-specific mean difference between the two methods, $\sigma_{\alpha\beta}^2$ is the subject-method variance component, and $\sigma_e^2$ is the residual variance component.

The function fits this estimator separately to every method pair. The homogeneity of agreement across conditions is tested by the method-time interaction: $$ F = \frac{MS_{\mathrm{method}\times\mathrm{time}}}{MS_{\mathrm{error}}}. $$ The returned homogeneity_F attribute stores this test statistic for each method pair, and homogeneity_p stores the corresponding upper-tail F-test p-value using degrees of freedom df_method_time and df_error. Larger homogeneity_F values and smaller homogeneity_p values indicate stronger evidence that agreement changes across conditions, meaning that the condition-specific CIA curves should be interpreted directly rather than relying on a single pooled/common coefficient. Conversely, a small homogeneity_F and a non-small homogeneity_p indicate that the data do not show strong evidence of method-by-condition heterogeneity, so a common estimate may be a reasonable summary if it is scientifically meaningful.

When homogeneous = TRUE, the function also reports a common or pooled CIA for each method pair from the reduced model that pools the method-time sum of squares into the residual term. That common estimate is meaningful only when agreement is reasonably homogeneous across conditions.

Confidence intervals can be computed by a delta-method normal approximation or by a subject-level percentile bootstrap. The delta-method interval is the default. For each method pair and condition k, let $\hat\psi_k$ be the ANOVA estimator, let $g_k$ be its numerical gradient with respect to the subject-level moment vector, and let $\Sigma$ be the empirical covariance matrix of that moment vector divided by the number of subjects. The delta-method standard error is $$ \widehat{\mathrm{se}}(\hat\psi_k) = \sqrt{\max\left(g_k^\top \Sigma g_k, 0\right)}. $$ The raw normal interval is $$ \hat\psi_k \pm z_{1-\alpha/2}\,\widehat{\mathrm{se}}(\hat\psi_k). $$ Under the bootstrap option, the interval is given by the empirical percentile limits from subject-resampled estimates. The argument estimator controls whether the reported result is the literal method-of-moments ratio or a bounded variance-component variant. Under estimator = "vc_constrained", the estimated subject-method variance component is constrained to be non-negative before converting it back to CIA, and the reported interval limits are also truncated to the CIA parameter space, $$ \mathrm{lwr} = \max(0, \mathrm{lwr}_{\mathrm{raw}}), \qquad \mathrm{upr} = \min(1, \mathrm{upr}_{\mathrm{raw}}). $$ This keeps the reported interval inside the CIA parameter space [0, 1]. Use estimator = "mom_unconstrained" to inspect the literal raw method-of-moments estimator and the corresponding unbounded interval on the estimator scale.

The current implementation requires exactly one observation in every subject-method-time cell. It therefore targets the balanced repeated-measures ANOVA setting from the cited paper and returns an error otherwise.

References

Haber M, Gao J, Barnhart HX. (2010). Evaluation of agreement between measurement methods from data with matched repeated measurements via the coefficient of individual agreement. Journal of Data Science, 8, 457-469.

Barnhart HX, Kosinski AS, Haber M. (2007). Assessing individual agreement. Journal of Biopharmaceutical Statistics, 17(4), 697-719.

Examples

Run this code

# Example 1
set.seed(1)
dat_rater <- expand.grid(
  id = factor(sprintf("s%02d", 1:8)),
  method = factor(c("Device_A", "Device_B")),
  time = factor(c("Rater_1", "Rater_2", "Rater_3")),
  KEEP.OUT.ATTRS = FALSE
)
subj_eff <- rnorm(nlevels(dat_rater$id), sd = 1)[dat_rater$id]
method_shift <- c(Device_A = 0, Device_B = 0.25)[dat_rater$method]
rater_shift <- c(Rater_1 = 0, Rater_2 = 0.15, Rater_3 = -0.10)[dat_rater$time]
interaction_shift <- ifelse(
  dat_rater$method == "Device_B" & dat_rater$time == "Rater_3",
  0.12, 0
)
dat_rater$y <- subj_eff + method_shift + rater_shift +
  interaction_shift + rnorm(nrow(dat_rater), sd = 0.25)

fit_rater <- cia_rm(
  dat_rater,
  response = "y",
  subject = "id",
  method = "method",
  time = "time",
  homogeneous = TRUE
)
print(fit_rater)
summary(fit_rater)
estimate(fit_rater)
tidy(fit_rater)
plot(fit_rater)

# Example 2
set.seed(2)
dat_time <- expand.grid(
  id = factor(sprintf("s%02d", 1:10)),
  method = factor(c("Assay_A", "Assay_B", "Assay_C", "Assay_D")),
  time = factor(
    c("baseline", "week2", "month1", "month2", "month3"),
    levels = c("baseline", "week2", "month1", "month2", "month3")
  ),
  KEEP.OUT.ATTRS = FALSE
)
subj_eff <- rnorm(nlevels(dat_time$id), sd = 0.9)[dat_time$id]
method_shift <- c(Assay_A = 0, Assay_B = 0.20, Assay_C = -0.10, Assay_D = 0.35)[dat_time$method]
time_shift <- c(
  baseline = 0, week2 = 0.10, month1 = 0.22, month2 = 0.32, month3 = 0.42
)[dat_time$time]
interaction_shift <- ifelse(dat_time$method == "Assay_B" & dat_time$time == "month3", 0.10, 0) +
  ifelse(dat_time$method == "Assay_D" & dat_time$time == "month2", -0.08, 0)
dat_time$y <- subj_eff + method_shift + time_shift +
  interaction_shift + rnorm(nrow(dat_time), sd = 0.30)

fit_time <- cia_rm(
  dat_time,
  response = "y",
  subject = "id",
  method = "method",
  time = "time",
  ci = TRUE
)
print(fit_time)
plot(fit_time)

# Example 3
set.seed(3)
dat_days <- expand.grid(
  id = factor(sprintf("s%02d", 1:12)),
  method = factor(c("Sensor_A", "Sensor_B", "Sensor_C")),
  time = 1:15,
  KEEP.OUT.ATTRS = FALSE
)
subj_eff <- rnorm(nlevels(dat_days$id), sd = 0.8)[dat_days$id]
method_shift <- c(Sensor_A = 0, Sensor_B = 0.15, Sensor_C = -0.08)[dat_days$method]
day_trend <- 0.05 * dat_days$time
interaction_shift <- ifelse(dat_days$method == "Sensor_B", 0.01 * dat_days$time, 0) +
  ifelse(dat_days$method == "Sensor_C", -0.005 * dat_days$time, 0)
dat_days$y <- subj_eff + method_shift + day_trend +
  interaction_shift + rnorm(nrow(dat_days), sd = 0.22)

fit_days <- cia_rm(
  dat_days,
  response = "y",
  subject = "id",
  method = "method",
  time = "time",
  ci = TRUE
)
plot(fit_days, facet_by_pair = TRUE)

Run the code above in your browser using DataLab