# -------- 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 <- bland_altman_repeated(
data = data,
response = "y", subject = "subject", method = "method", time = "time",
two = 1.96, conf_level = 0.95,
include_slope = FALSE, use_ar1 = FALSE
)
summary(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 <- bland_altman_repeated(
response = data$y, subject = data$subject, method = data$method, time = data$time,
two = 1.96, conf_level = 0.95,
include_slope = FALSE, # estimate proportional bias per pair
use_ar1 = TRUE # model AR(1) within-subject
)
# Matrices (row - column orientation)
print(baN)
summary(baN)
# Faceted BA scatter by pair
plot(baN, smoother = "lm", facet_scales = "free_y")
# -------- Two-method path (A vs B only) -----------------------------------
data_AB <- subset(data, method %in% c("A","B"))
baAB <- bland_altman_repeated(
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)
Run the code above in your browser using DataLab