Nothing
R/tablong.R
In R4VN: Health Data Analysis and Publication-Ready Reporting
Defines functions
.r4vn_long_contrasts_table
.r4vn_long_tests_table
.r4vn_long_descriptive_rows
.r4vn_long_show_html
.r4vn_long_html_document
.r4vn_long_html_table
.r4vn_long_build_rows
.r4vn_long_change_name
.r4vn_long_effect_name
.r4vn_long_observed_cell
.r4vn_long_summary_count
.r4vn_long_summary_binary
.r4vn_long_summary_continuous
.r4vn_long_contrast_rows
.r4vn_long_contrast
.r4vn_long_grid_row
.r4vn_long_fit_one
.r4vn_long_extract_beta_vcov
.r4vn_long_term_columns
.r4vn_long_wald
.r4vn_long_fixed_rank
.r4vn_long_safe_scalar
.r4vn_long_glm_robust_vcov
.r4vn_long_prepare_wide
.r4vn_long_prepare_long_common
.r4vn_long_time_source
.r4vn_long_common_wide_label
.r4vn_long_event
.r4vn_long_apply_covariates
.r4vn_long_numeric
.r4vn_long_factor
.r4vn_long_adjusted_meta
.r4vn_long_meta
.r4vn_long_parse_symbol
.r4vn_long_parse_spec_text
.r4vn_long_fmt_ci
.r4vn_long_fmt_p
.r4vn_long_fmt
.r4vn_long_label
.r4vn_long_levels
.r4vn_long_escape
# ============================================================================
# R4VN::tablong()
# Longitudinal / repeated-measures publication table
#
# Design principles
# - Same culture as tab(): data = NULL uses active data; vars() declares outcomes.
# - Accepts long and wide data; wide data are reshaped internally only.
# - Continuous repeated outcomes: random-intercept model via recommended nlme;
# automatic base-R cluster-robust fallback if nlme is unavailable.
# - Binary/count repeated outcomes: marginal regression with subject-clustered
# robust sandwich variance calculated internally by R4VN.
# - rr = TRUE / pr = TRUE: modified Poisson with robust variance; geepack is
# optional only for explicitly requested AR(1) GEE.
# - Repeated cross-sectional data are detected automatically when IDs do not
# repeat across time (or no ID is supplied).
# - Returns c("r4vn_tablong", "r4vn_tab") so tabexport() works unchanged.
# ============================================================================
.r4vn_long_escape <- function(x) {
x <- as.character(x)
x[is.na(x)] <- ""
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
gsub("'", "'", x, fixed = TRUE)
}
.r4vn_long_levels <- function(x) {
z <- x[!is.na(x)]
if (!length(z)) return(character())
if (is.factor(x)) {
lev <- levels(x)
return(lev[lev %in% as.character(z)])
}
if (is.logical(x)) return(as.character(c(FALSE, TRUE)[c(FALSE, TRUE) %in% z]))
if (is.numeric(x)) return(as.character(sort(unique(z))))
unique(as.character(z))
}
.r4vn_long_label <- function(x, fallback) {
lab <- attr(x, "label", exact = TRUE)
if (is.null(lab) || !length(lab) || is.na(lab[1L]) || !nzchar(as.character(lab[1L]))) {
return(fallback)
}
as.character(lab[1L])
}
.r4vn_long_fmt <- function(x, digits = 2) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
formatC(x, format = "f", digits = digits)
}
.r4vn_long_fmt_p <- function(p, digits = 3) {
if (!length(p) || is.na(p) || !is.finite(p)) return("")
limit <- 10^(-digits)
if (p < limit) return(paste0("<", formatC(limit, format = "f", digits = digits)))
formatC(p, format = "f", digits = digits)
}
.r4vn_long_fmt_ci <- function(est, low, high, digits = 2) {
if (any(!is.finite(c(est, low, high)))) return("")
paste0(
formatC(est, format = "f", digits = digits),
" (",
formatC(low, format = "f", digits = digits),
", ",
formatC(high, format = "f", digits = digits),
")"
)
}
.r4vn_long_parse_spec_text <- function(text, allow_continuous = TRUE) {
out <- list(
variable = text,
type = "categorical",
reference_index = 1L,
continuous = FALSE,
specification = text
)
if (grepl("^b[1-9][0-9]*\\.", text)) {
out$reference_index <- as.integer(sub("^b([1-9][0-9]*)\\..*$", "\\1", text))
out$variable <- sub("^b[1-9][0-9]*\\.", "", text)
} else if (startsWith(text, "c.")) {
if (!allow_continuous) {
stop("The `c.` prefix is not allowed here.", call. = FALSE)
}
out$type <- "mean"
out$continuous <- TRUE
out$reference_index <- NA_integer_
out$variable <- sub("^c\\.", "", text)
} else if (startsWith(text, "q.")) {
out$type <- "median"
out$reference_index <- NA_integer_
out$variable <- sub("^q\\.", "", text)
} else if (startsWith(text, "f.")) {
out$type <- "full"
out$reference_index <- NA_integer_
out$variable <- sub("^f\\.", "", text)
}
if (!nzchar(out$variable)) stop("A prefix must be followed by a variable name.", call. = FALSE)
out
}
.r4vn_long_parse_symbol <- function(expr, role, allow_continuous = TRUE) {
if (!is.symbol(expr)) {
stop("`", role, "` must be an unquoted variable name in long data.", call. = FALSE)
}
.r4vn_long_parse_spec_text(as.character(expr), allow_continuous = allow_continuous)
}
.r4vn_long_meta <- function(x, argument = "vars", data = NULL) {
if (!inherits(x, "r4vn_vars")) {
stop("`", argument, "` must be created using `vars()`.", call. = FALSE)
}
required <- c("variable", "type", "reference_index")
if (!all(required %in% names(x))) {
stop("`", argument, "` is not a valid current R4VN `vars()` object.", call. = FALSE)
}
if (is.data.frame(data)) {
resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
if (is.function(resolver)) {
x <- resolver(x, data = data, default_type = "auto", strict = TRUE)
}
}
x
}
.r4vn_long_adjusted_meta <- function(expr, missing_arg, env, data = NULL) {
if (isTRUE(missing_arg) || identical(expr, quote(NULL))) return(NULL)
value <- eval(expr, envir = env)
.r4vn_long_meta(value, "adjusted", data = data)
}
.r4vn_long_factor <- function(x, reference_index = 1L, reference_value = NULL) {
display_levels <- .r4vn_long_levels(x)
if (length(display_levels) < 1L) {
return(list(
x = factor(x),
display_levels = character(),
model_levels = character(),
reference = NA_character_
))
}
reference <- NULL
if (!is.null(reference_value) && length(reference_value) && !is.na(reference_value[1L])) {
candidate <- as.character(reference_value[1L])
hit <- match(candidate, display_levels)
if (is.na(hit)) {
stop("Reference level `", candidate, "` was not found.", call. = FALSE)
}
reference <- display_levels[hit]
} else {
if (is.na(reference_index) || reference_index < 1L || reference_index > length(display_levels)) {
reference_index <- 1L
}
reference <- display_levels[reference_index]
}
model_levels <- c(reference, setdiff(display_levels, reference))
list(
x = factor(as.character(x), levels = model_levels),
display_levels = display_levels,
model_levels = model_levels,
reference = reference
)
}
.r4vn_long_numeric <- function(x, variable) {
y <- suppressWarnings(as.numeric(x))
if (all(is.na(y)) && any(!is.na(x))) {
stop("Variable `", variable, "` cannot be converted to numeric.", call. = FALSE)
}
y
}
.r4vn_long_apply_covariates <- function(data, meta) {
if (is.null(meta) || !nrow(meta)) {
return(list(data = data, names = character(), map = data.frame()))
}
names_internal <- character(nrow(meta))
map <- meta
for (i in seq_len(nrow(meta))) {
variable <- meta$variable[i]
if (!variable %in% names(data)) {
stop("Adjusted variable `", variable, "` was not found in `data`.", call. = FALSE)
}
internal <- paste0(".z", i)
names_internal[i] <- internal
if (identical(meta$type[i], "categorical")) {
f <- .r4vn_long_factor(
data[[variable]],
reference_index = meta$reference_index[i]
)
data[[internal]] <- f$x
} else {
data[[internal]] <- .r4vn_long_numeric(data[[variable]], variable)
}
}
map$internal <- names_internal
list(data = data, names = names_internal, map = map)
}
.r4vn_long_event <- function(x, event, outcome_name) {
observed <- .r4vn_long_levels(x)
if (length(observed) != 2L) {
stop(
"Binary outcome `", outcome_name, "` must have exactly two observed levels; found ",
length(observed), ".",
call. = FALSE
)
}
chosen <- NULL
if (!is.null(event)) {
if (length(event) > 1L && !is.null(names(event)) && outcome_name %in% names(event)) {
chosen <- as.character(event[[outcome_name]])
} else if (length(event) == 1L) {
chosen <- as.character(event[1L])
}
}
if (is.null(chosen) || is.na(chosen) || !nzchar(chosen)) chosen <- observed[length(observed)]
if (!chosen %in% observed) {
stop("Event `", chosen, "` was not found in outcome `", outcome_name, "`.", call. = FALSE)
}
chosen
}
.r4vn_long_common_wide_label <- function(data, variables) {
labs <- vapply(variables, function(v) .r4vn_long_label(data[[v]], ""), character(1))
labs <- unique(labs[nzchar(labs)])
if (length(labs) == 1L) return(labs)
stripped <- sub("([_.]?(baseline|base|pre|before|post|after|month|m|visit|t)?[_.]?[0-9]+)$",
"", variables, ignore.case = TRUE)
stripped <- sub("([_.](baseline|base|pre|before|post|after))$", "",
stripped, ignore.case = TRUE)
stripped <- unique(stripped[nzchar(stripped)])
if (length(stripped) == 1L) return(stripped)
variables[1L]
}
.r4vn_long_time_source <- function(expr, data, meta_n, env) {
if (identical(expr, quote(NULL))) {
if (meta_n > 1L) return(list(mode = "wide", labels = NULL, spec = NULL))
stop("`time` is required for long data.", call. = FALSE)
}
if (is.symbol(expr)) {
spec <- .r4vn_long_parse_spec_text(as.character(expr), allow_continuous = TRUE)
if (spec$variable %in% names(data)) {
return(list(mode = "long", labels = NULL, spec = spec))
}
value <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (!is.null(value) && length(value) == meta_n && meta_n > 1L) {
return(list(mode = "wide", labels = as.character(value), spec = NULL))
}
stop(
"`time` variable `", spec$variable,
"` was not found in `data`. For wide data, use a vector such as ",
'`time = c("Baseline", "Month 3", "Month 6")`.',
call. = FALSE
)
}
value <- eval(expr, envir = env)
if (meta_n <= 1L) {
stop("A vector of time labels is only used with wide data containing multiple repeated outcome variables.", call. = FALSE)
}
if (length(value) != meta_n) {
stop("The number of `time` labels must equal the number of repeated variables in `vars()`.", call. = FALSE)
}
list(mode = "wide", labels = as.character(value), spec = NULL)
}
.r4vn_long_prepare_long_common <- function(data, time_spec, id_expr, id_missing,
by_expr, by_missing, adjusted_meta,
ref, env) {
time_name <- time_spec$variable
if (!time_name %in% names(data)) stop("`time` variable was not found.", call. = FALSE)
out <- data
if (isTRUE(time_spec$continuous)) {
out$.time <- .r4vn_long_numeric(out[[time_name]], time_name)
time_display <- sort(unique(out$.time[is.finite(out$.time)]))
time_reference <- if (length(time_display)) min(time_display) else NA_real_
time_continuous <- TRUE
} else {
tf <- .r4vn_long_factor(
out[[time_name]],
reference_index = time_spec$reference_index,
reference_value = ref
)
out$.time <- tf$x
time_display <- tf$display_levels
time_reference <- tf$reference
time_continuous <- FALSE
}
id_name <- NULL
if (isTRUE(id_missing) || identical(id_expr, quote(NULL))) {
out$.id <- seq_len(nrow(out))
} else {
id_spec <- .r4vn_long_parse_symbol(id_expr, "id", allow_continuous = FALSE)
id_name <- id_spec$variable
if (!id_name %in% names(out)) stop("`id` variable `", id_name, "` was not found.", call. = FALSE)
out$.id <- out[[id_name]]
}
by_name <- NULL
by_display <- NULL
by_reference <- NULL
if (!isTRUE(by_missing) && !identical(by_expr, quote(NULL))) {
by_spec <- .r4vn_long_parse_symbol(by_expr, "by", allow_continuous = FALSE)
by_name <- by_spec$variable
if (!by_name %in% names(out)) stop("`by` variable `", by_name, "` was not found.", call. = FALSE)
bf <- .r4vn_long_factor(out[[by_name]], reference_index = by_spec$reference_index)
if (length(bf$display_levels) < 2L) stop("`by` must contain at least two observed groups.", call. = FALSE)
out$.by <- bf$x
by_display <- bf$display_levels
by_reference <- bf$reference
}
cov <- .r4vn_long_apply_covariates(out, adjusted_meta)
out <- cov$data
repeated <- FALSE
if (!is.null(id_name)) {
valid <- !is.na(out$.id) & !is.na(out$.time)
if (any(valid)) {
key <- split(as.character(out$.time[valid]), as.character(out$.id[valid]))
repeated <- any(vapply(key, function(z) length(unique(z)) > 1L, logical(1)))
}
}
list(
data = out,
time_name = time_name,
time_display = time_display,
time_reference = time_reference,
time_continuous = time_continuous,
id_name = id_name,
repeated = repeated,
by_name = by_name,
by_display = by_display,
by_reference = by_reference,
covariates = cov$names,
covariate_map = cov$map
)
}
.r4vn_long_prepare_wide <- function(data, meta, time_labels, id_expr, id_missing,
by_expr, by_missing, adjusted_meta,
ref, exposure_expr, exposure_missing, env) {
variables <- meta$variable
missing_vars <- setdiff(variables, names(data))
if (length(missing_vars)) {
stop("Variables not found in `data`: ", paste(missing_vars, collapse = ", "), ".", call. = FALSE)
}
types <- unique(meta$type)
if (length(types) > 1L) {
continuous_types <- c("mean", "median", "full")
if (!all(types %in% continuous_types)) {
stop("All repeated variables in wide data must describe the same outcome type.", call. = FALSE)
}
}
if (is.null(time_labels)) {
time_labels <- variables
}
if (anyDuplicated(time_labels)) stop("Time labels must be unique.", call. = FALSE)
id_name <- NULL
base_id <- seq_len(nrow(data))
if (!isTRUE(id_missing) && !identical(id_expr, quote(NULL))) {
id_spec <- .r4vn_long_parse_symbol(id_expr, "id", allow_continuous = FALSE)
id_name <- id_spec$variable
if (!id_name %in% names(data)) stop("`id` variable `", id_name, "` was not found.", call. = FALSE)
base_id <- data[[id_name]]
}
by_name <- NULL
by_spec <- NULL
if (!isTRUE(by_missing) && !identical(by_expr, quote(NULL))) {
by_spec <- .r4vn_long_parse_symbol(by_expr, "by", allow_continuous = FALSE)
by_name <- by_spec$variable
if (!by_name %in% names(data)) stop("`by` variable `", by_name, "` was not found.", call. = FALSE)
}
exposure_mode <- "none"
exposure_vars <- NULL
exposure_name <- NULL
if (!isTRUE(exposure_missing) && !identical(exposure_expr, quote(NULL))) {
exposure_value <- tryCatch(eval(exposure_expr, envir = env), error = function(e) NULL)
if (inherits(exposure_value, "r4vn_vars")) {
exposure_vars <- exposure_value$variable
if (length(exposure_vars) != length(variables)) {
stop("Wide `exposure = vars(...)` must contain one exposure variable per repeated outcome variable.", call. = FALSE)
}
if (any(!exposure_vars %in% names(data))) stop("Some exposure variables were not found in `data`.", call. = FALSE)
exposure_mode <- "wide"
} else if (is.symbol(exposure_expr)) {
exposure_name <- as.character(exposure_expr)
if (!exposure_name %in% names(data)) stop("`exposure` variable was not found.", call. = FALSE)
exposure_mode <- "single"
} else {
stop("In wide data, `exposure` must be one variable or `vars(...)` with one variable per time point.", call. = FALSE)
}
}
pieces <- vector("list", length(variables))
for (j in seq_along(variables)) {
piece <- data
piece$.id <- base_id
piece$.time_source <- time_labels[j]
if (all(meta$type %in% c("mean", "median", "full"))) {
piece$.outcome <- .r4vn_long_numeric(data[[variables[j]]], variables[j])
} else {
piece$.outcome <- as.character(data[[variables[j]]])
}
if (identical(exposure_mode, "wide")) piece$.exposure <- data[[exposure_vars[j]]]
if (identical(exposure_mode, "single")) piece$.exposure <- data[[exposure_name]]
pieces[[j]] <- piece
}
out <- do.call(rbind, pieces)
rownames(out) <- NULL
tf <- .r4vn_long_factor(
out$.time_source,
reference_index = 1L,
reference_value = ref
)
out$.time <- tf$x
by_display <- NULL
by_reference <- NULL
if (!is.null(by_name)) {
bf <- .r4vn_long_factor(out[[by_name]], reference_index = by_spec$reference_index)
if (length(bf$display_levels) < 2L) stop("`by` must contain at least two observed groups.", call. = FALSE)
out$.by <- bf$x
by_display <- bf$display_levels
by_reference <- bf$reference
}
cov <- .r4vn_long_apply_covariates(out, adjusted_meta)
out <- cov$data
list(
data = out,
outcome_name = .r4vn_long_common_wide_label(data, variables),
time_name = NULL,
time_display = tf$display_levels,
time_reference = tf$reference,
time_continuous = FALSE,
id_name = id_name,
repeated = TRUE,
by_name = by_name,
by_display = by_display,
by_reference = by_reference,
covariates = cov$names,
covariate_map = cov$map,
exposure_mode = exposure_mode
)
}
.r4vn_long_glm_robust_vcov <- function(fit, cluster = NULL) {
X <- tryCatch(stats::model.matrix(fit), error = function(e) NULL)
if (is.null(X) || !nrow(X) || !ncol(X)) return(NULL)
if (inherits(fit, "glm")) {
wr <- tryCatch(stats::residuals(fit, type = "working"), error = function(e) NULL)
ww <- fit$weights
if (is.null(wr) || is.null(ww) || length(wr) != nrow(X) || length(ww) != nrow(X)) return(NULL)
} else {
wr <- tryCatch(stats::residuals(fit), error = function(e) NULL)
ww <- fit$weights
if (is.null(ww)) ww <- rep(1, nrow(X))
if (is.null(wr) || length(wr) != nrow(X) || length(ww) != nrow(X)) return(NULL)
}
score <- X * as.vector(wr * ww)
information <- crossprod(X, X * as.vector(ww))
bread <- tryCatch(
solve(information),
error = function(e) tryCatch(qr.solve(information), error = function(e2) NULL)
)
if (is.null(bread)) return(NULL)
correction <- 1
if (!is.null(cluster)) {
cluster <- as.character(cluster)
if (length(cluster) != nrow(X)) return(NULL)
cluster[is.na(cluster)] <- "<NA>"
U <- rowsum(score, group = cluster, reorder = FALSE)
meat <- crossprod(U)
G <- nrow(U)
N <- nrow(X)
P <- qr(X)$rank
if (G > 1L && N > P) {
correction <- (G / (G - 1)) * ((N - 1) / (N - P))
}
} else {
meat <- crossprod(score)
N <- nrow(X)
P <- qr(X)$rank
if (N > P) correction <- N / (N - P)
}
V <- correction * bread %*% meat %*% bread
dimnames(V) <- list(colnames(X), colnames(X))
V
}
.r4vn_long_safe_scalar <- function(x) {
if (is.null(x) || !length(x)) return(NA_real_)
z <- suppressWarnings(as.numeric(x[1L]))
if (!length(z) || is.na(z) || !is.finite(z)) NA_real_ else z
}
.r4vn_long_fixed_rank <- function(formula, data) {
mf <- tryCatch(
stats::model.frame(formula, data = data, na.action = stats::na.omit),
error = function(e) NULL
)
if (is.null(mf) || !nrow(mf)) {
return(list(full_rank = TRUE, rank = NA_integer_, columns = NA_integer_))
}
mm <- tryCatch(stats::model.matrix(formula, data = mf), error = function(e) NULL)
if (is.null(mm) || !ncol(mm)) {
return(list(full_rank = TRUE, rank = NA_integer_, columns = NA_integer_))
}
q <- qr(mm)
list(full_rank = q$rank == ncol(mm), rank = q$rank, columns = ncol(mm))
}
.r4vn_long_wald <- function(beta, V, terms) {
terms <- intersect(terms, names(beta))
if (!length(terms)) return(NA_real_)
b <- beta[terms]
VV <- V[terms, terms, drop = FALSE]
good <- is.finite(b) & is.finite(diag(VV)) & diag(VV) > 0
b <- b[good]
VV <- VV[good, good, drop = FALSE]
if (!length(b)) return(NA_real_)
inv <- tryCatch(solve(VV), error = function(e) tryCatch(qr.solve(VV), error = function(e2) NULL))
if (is.null(inv)) return(NA_real_)
df <- qr(VV)$rank
if (df < 1L) return(NA_real_)
stat <- as.numeric(t(b) %*% inv %*% b)
if (!is.finite(stat)) return(NA_real_)
stats::pchisq(stat, df = df, lower.tail = FALSE)
}
.r4vn_long_term_columns <- function(fit, fixed_formula, term_label) {
mf <- tryCatch(stats::model.frame(fixed_formula, data = fit$model), error = function(e) NULL)
if (is.null(mf)) {
mf <- tryCatch(fit$model, error = function(e) NULL)
}
if (is.null(mf)) return(character())
mm <- tryCatch(stats::model.matrix(fixed_formula, data = mf), error = function(e) NULL)
if (is.null(mm)) return(character())
assign <- attr(mm, "assign")
labels <- attr(stats::terms(fixed_formula), "term.labels")
position <- which(labels == term_label)
if (!length(position) && identical(term_label, ".time:.by")) {
position <- which(labels %in% c(".time:.by", ".by:.time"))
}
if (!length(position)) return(character())
colnames(mm)[assign %in% position]
}
.r4vn_long_extract_beta_vcov <- function(fit, engine, robust = FALSE, cluster = NULL) {
if (identical(engine, "mixed")) {
if (inherits(fit, "lme")) {
beta <- nlme::fixef(fit)
V <- as.matrix(stats::vcov(fit))
} else {
beta <- stats::coef(fit)
V <- as.matrix(stats::vcov(fit))
}
} else {
beta <- stats::coef(fit)
V <- if (isTRUE(robust)) .r4vn_long_glm_robust_vcov(fit, cluster = cluster) else as.matrix(stats::vcov(fit))
}
list(beta = beta, V = V)
}
.r4vn_long_fit_one <- function(data, outcome_type, effect_type, repeated,
gee, ar1, slope, time_continuous,
covariates, exposure = FALSE) {
has_by <- ".by" %in% names(data)
z_rhs <- if (length(covariates)) paste(covariates, collapse = " + ") else ""
offset_rhs <- if (isTRUE(exposure)) "offset(log(.exposure))" else ""
join_rhs <- function(parts) {
parts <- parts[nzchar(parts)]
if (!length(parts)) "1" else paste(parts, collapse = " + ")
}
main_time <- ".time"
main_by <- if (has_by) ".by" else ""
full_core <- if (has_by) ".time * .by" else ".time"
full_rhs <- join_rhs(c(full_core, z_rhs, offset_rhs))
add_rhs <- join_rhs(c(main_time, main_by, z_rhs, offset_rhs))
time_reduced_rhs <- join_rhs(c(main_by, z_rhs, offset_rhs))
by_reduced_rhs <- join_rhs(c(main_time, z_rhs, offset_rhs))
null_rhs <- join_rhs(c(z_rhs, offset_rhs))
fixed_full <- stats::as.formula(paste(".outcome ~", full_rhs))
fixed_add <- stats::as.formula(paste(".outcome ~", add_rhs))
fixed_time_reduced <- stats::as.formula(paste(".outcome ~", time_reduced_rhs))
fixed_by_reduced <- stats::as.formula(paste(".outcome ~", by_reduced_rhs))
fixed_null <- stats::as.formula(paste(".outcome ~", null_rhs))
engine <- NULL
link <- "identity"
family <- NULL
robust <- FALSE
fits <- list()
correlation <- "independent"
package_note <- NULL
if (identical(outcome_type, "continuous")) {
family <- stats::gaussian()
link <- "identity"
} else if (identical(outcome_type, "binary")) {
if (effect_type %in% c("RR", "PR")) {
family <- stats::poisson(link = "log")
link <- "log"
} else {
family <- stats::binomial(link = "logit")
link <- "logit"
}
} else if (identical(outcome_type, "count")) {
family <- stats::poisson(link = "log")
link <- "log"
} else {
stop("Unsupported outcome type.", call. = FALSE)
}
# Advanced AR(1) GEE remains available when geepack is installed. Routine
# longitudinal RR/PR, binary and count analyses do not require geepack.
request_gee <- isTRUE(gee) || effect_type %in% c("RR", "PR")
use_geepack <- isTRUE(repeated) && isTRUE(ar1) && isTRUE(request_gee) &&
requireNamespace("geepack", quietly = TRUE)
if (isTRUE(repeated) && isTRUE(ar1) && isTRUE(request_gee) && !use_geepack) {
warning(
"`ar1 = TRUE` requires optional package `geepack`; using working-independence cluster-robust inference instead.",
call. = FALSE
)
package_note <- "AR(1) was requested but geepack was unavailable; working-independence cluster-robust inference was used."
}
if (isTRUE(repeated) && use_geepack) {
time_order <- if (is.factor(data$.time)) as.numeric(data$.time) else data$.time
data <- data[order(as.character(data$.id), time_order, na.last = TRUE), , drop = FALSE]
rownames(data) <- NULL
engine <- "gee"
correlation <- "AR(1)"
fit_fun <- function(fixed) {
geepack::geeglm(
formula = fixed,
data = data,
id = data$.id,
family = family,
corstr = "ar1",
std.err = "san.se",
na.action = stats::na.omit
)
}
fits$full <- fit_fun(fixed_full)
if (has_by) fits$add <- fit_fun(fixed_add)
beta_vcov <- list(
beta = stats::coef(fits$full),
V = as.matrix(stats::vcov(fits$full))
)
if (has_by) {
add_beta <- stats::coef(fits$add)
add_V <- as.matrix(stats::vcov(fits$add))
time_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".time")
group_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".by")
interaction_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time:.by")
p_time <- .r4vn_long_wald(add_beta, add_V, time_cols)
p_group <- .r4vn_long_wald(add_beta, add_V, group_cols)
p_interaction <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, interaction_cols)
} else {
time_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time")
p_time <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, time_cols)
p_group <- NA_real_
p_interaction <- NA_real_
}
singular <- FALSE
df_contrast <- Inf
} else if (isTRUE(repeated) && identical(outcome_type, "continuous") && !isTRUE(gee) &&
requireNamespace("nlme", quietly = TRUE)) {
# nlme is an R recommended package and is included in ordinary R installs.
# Before fitting, detect perfect collinearity in the fixed-effects design.
# This avoids exposing cryptic nlme errors such as
# "Singularity in backsolve at level 0, block 1" to R4VN users.
rank_check <- .r4vn_long_fixed_rank(fixed_full, data)
if (!isTRUE(rank_check$full_rank)) {
stop(
"The adjusted longitudinal model cannot be estimated because the fixed-effects predictors are perfectly collinear. ",
"This commonly occurs when an adjustment variable is identical to, or completely determined by, the group/time variables. ",
"Remove or recode the redundant adjustment variable and run tablong() again.",
call. = FALSE
)
}
engine <- "mixed"
correlation <- "random intercept"
fit_fun <- function(fixed) {
random_formula <- if (isTRUE(slope) && isTRUE(time_continuous)) {
stats::as.formula("~ 1 + .time | .id")
} else {
stats::as.formula("~ 1 | .id")
}
nlme::lme(
fixed = fixed,
random = random_formula,
data = data,
method = "ML",
na.action = stats::na.omit,
control = nlme::lmeControl(returnObject = TRUE)
)
}
fits$full <- fit_fun(fixed_full)
if (has_by) {
fits$add <- fit_fun(fixed_add)
fits$time_reduced <- fit_fun(fixed_time_reduced)
fits$by_reduced <- fit_fun(fixed_by_reduced)
} else {
fits$null <- fit_fun(fixed_null)
}
extract_p_compare <- function(a, b) {
z <- tryCatch(stats::anova(a, b), error = function(e) NULL)
if (is.null(z) || nrow(z) < 2L) return(NA_real_)
pcol <- grep("p-value|Pr\\(", names(z), value = TRUE, ignore.case = TRUE)
if (!length(pcol)) return(NA_real_)
.r4vn_long_safe_scalar(z[[pcol[1L]]][2L])
}
if (has_by) {
p_time <- extract_p_compare(fits$time_reduced, fits$add)
p_group <- extract_p_compare(fits$by_reduced, fits$add)
p_interaction <- extract_p_compare(fits$add, fits$full)
} else {
p_time <- extract_p_compare(fits$null, fits$full)
p_group <- NA_real_
p_interaction <- NA_real_
}
beta_vcov <- .r4vn_long_extract_beta_vcov(fits$full, engine)
singular <- FALSE
df_contrast <- Inf
} else {
# Base-R marginal models. For repeated subjects, cluster-robust sandwich
# variance is calculated internally by R4VN, so sandwich/geepack/lme4 are
# not required for routine binary, RR/PR, count, or Gaussian GEE analyses.
repeated_robust <- isTRUE(repeated)
robust <- repeated_robust || (identical(outcome_type, "binary") && effect_type %in% c("RR", "PR"))
engine <- if (repeated_robust) "cluster_robust" else "independent"
correlation <- if (repeated_robust) "working independence" else "independent"
if (isTRUE(repeated) && identical(outcome_type, "continuous") && !isTRUE(gee) &&
!requireNamespace("nlme", quietly = TRUE)) {
warning(
"Recommended package `nlme` is unavailable; using a marginal linear model with subject-clustered robust standard errors.",
call. = FALSE
)
package_note <- "nlme was unavailable; a marginal linear model with subject-clustered robust standard errors was used."
}
if (isTRUE(slope) && identical(outcome_type, "continuous") && !identical(engine, "mixed")) {
warning("`slope = TRUE` requires the recommended package `nlme`; the random slope was not fitted.", call. = FALSE)
}
fit_fun <- function(fixed) {
if (identical(outcome_type, "continuous")) {
stats::lm(fixed, data = data, na.action = stats::na.omit)
} else {
stats::glm(fixed, data = data, family = family, na.action = stats::na.omit)
}
}
fits$full <- fit_fun(fixed_full)
if (has_by) fits$add <- fit_fun(fixed_add)
cluster <- if (repeated_robust) data$.id else NULL
full_V <- if (robust) .r4vn_long_glm_robust_vcov(fits$full, cluster = cluster) else as.matrix(stats::vcov(fits$full))
if (is.null(full_V)) full_V <- as.matrix(stats::vcov(fits$full))
beta_vcov <- list(beta = stats::coef(fits$full), V = full_V)
if (has_by) {
add_V <- if (robust) .r4vn_long_glm_robust_vcov(fits$add, cluster = cluster) else as.matrix(stats::vcov(fits$add))
if (is.null(add_V)) add_V <- as.matrix(stats::vcov(fits$add))
add_beta <- stats::coef(fits$add)
time_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".time")
group_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".by")
interaction_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time:.by")
p_time <- .r4vn_long_wald(add_beta, add_V, time_cols)
p_group <- .r4vn_long_wald(add_beta, add_V, group_cols)
p_interaction <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, interaction_cols)
} else {
time_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time")
p_time <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, time_cols)
p_group <- NA_real_
p_interaction <- NA_real_
}
singular <- FALSE
if (repeated_robust) {
G <- length(unique(data$.id[!is.na(data$.id)]))
df_contrast <- if (G > 1L) G - 1L else Inf
} else {
df_contrast <- if (identical(outcome_type, "continuous")) stats::df.residual(fits$full) else Inf
}
}
list(
fit = fits$full,
fits = fits,
engine = engine,
family = family,
link = link,
fixed_formula = fixed_full,
beta = beta_vcov$beta,
V = beta_vcov$V,
p_time = p_time,
p_group = p_group,
p_interaction = p_interaction,
singular = singular,
df_contrast = df_contrast,
correlation = correlation,
package_note = package_note
)
}
.r4vn_long_grid_row <- function(data, fixed_formula, time_value, by_value = NULL) {
nd <- data.frame(.outcome = 0)
if (is.factor(data$.time)) {
nd$.time <- factor(as.character(time_value), levels = levels(data$.time))
} else {
nd$.time <- as.numeric(time_value)
}
if (".by" %in% names(data)) {
nd$.by <- factor(as.character(by_value), levels = levels(data$.by))
}
z_names <- grep("^\\.z[0-9]+$", names(data), value = TRUE)
for (z in z_names) {
if (is.factor(data[[z]])) {
nd[[z]] <- factor(levels(data[[z]])[1L], levels = levels(data[[z]]))
} else {
value <- mean(data[[z]], na.rm = TRUE)
if (!is.finite(value)) value <- 0
nd[[z]] <- value
}
}
if (".exposure" %in% names(data)) nd$.exposure <- 1
mm <- stats::model.matrix(fixed_formula, data = nd)
mm
}
.r4vn_long_contrast <- function(engine_fit, L, level = 0.95) {
beta <- engine_fit$beta
V <- engine_fit$V
common <- intersect(names(beta), colnames(L))
if (!length(common)) {
return(c(estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
}
L2 <- as.numeric(L[, common, drop = FALSE])
names(L2) <- common
b <- beta[common]
VV <- V[common, common, drop = FALSE]
estimate_link <- sum(L2 * b)
variance <- as.numeric(t(L2) %*% VV %*% L2)
if (!is.finite(variance) || variance < 0) {
return(c(estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
}
se <- sqrt(variance)
alpha <- 1 - level
df <- engine_fit$df_contrast
critical <- if (is.finite(df)) stats::qt(1 - alpha / 2, df = df) else stats::qnorm(1 - alpha / 2)
lower_link <- estimate_link - critical * se
upper_link <- estimate_link + critical * se
statistic <- if (se > 0) estimate_link / se else NA_real_
p <- if (!is.finite(statistic)) NA_real_ else {
if (is.finite(df)) 2 * stats::pt(-abs(statistic), df = df)
else 2 * stats::pnorm(-abs(statistic))
}
if (engine_fit$link %in% c("log", "logit")) {
c(
estimate = exp(estimate_link),
lower = exp(lower_link),
upper = exp(upper_link),
p = p
)
} else {
c(
estimate = estimate_link,
lower = lower_link,
upper = upper_link,
p = p
)
}
}
.r4vn_long_contrast_rows <- function(data, fit, time_display, time_reference,
by_display, time_continuous, change,
pairwise, level, adjust) {
has_by <- ".by" %in% names(data)
out <- list()
index <- 0L
add <- function(type, time, group1 = "", group2 = "", comparison = "", result) {
index <<- index + 1L
out[[index]] <<- data.frame(
type = type,
time = as.character(time),
group1 = as.character(group1),
group2 = as.character(group2),
comparison = as.character(comparison),
estimate = unname(result["estimate"]),
lower = unname(result["lower"]),
upper = unname(result["upper"]),
p = unname(result["p"]),
stringsAsFactors = FALSE
)
}
xrow <- function(time, group = NULL) {
.r4vn_long_grid_row(data, fit$fixed_formula, time, group)
}
if (isTRUE(time_continuous)) {
groups <- if (has_by) by_display else ""
base_time <- if (length(time_display)) min(as.numeric(time_display)) else 0
if (has_by && length(by_display) == 2L) {
for (tt in time_display) {
X1 <- xrow(tt, by_display[1L])
X2 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, X2 - X1, level)
add("between", tt, by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
}
for (g in groups) {
X0 <- xrow(base_time, if (has_by) g else NULL)
X1 <- xrow(base_time + 1, if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, X1 - X0, level)
add("slope", "Per 1 time unit", if (has_by) g else "", "", "Per 1 time unit", result)
}
if (has_by && length(by_display) == 2L) {
X00 <- xrow(base_time, by_display[1L])
X01 <- xrow(base_time + 1, by_display[1L])
X10 <- xrow(base_time, by_display[2L])
X11 <- xrow(base_time + 1, by_display[2L])
result <- .r4vn_long_contrast(fit, (X11 - X10) - (X01 - X00), level)
add("slope_difference", "Per 1 time unit", by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
} else {
if (has_by && length(by_display) == 2L) {
for (tt in time_display) {
X1 <- xrow(tt, by_display[1L])
X2 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, X2 - X1, level)
add("between", tt, by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
}
if (isTRUE(change)) {
follow <- setdiff(time_display, as.character(time_reference))
groups <- if (has_by) by_display else ""
for (tt in follow) {
for (g in groups) {
X0 <- xrow(time_reference, if (has_by) g else NULL)
X1 <- xrow(tt, if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, X1 - X0, level)
add("change", tt, if (has_by) g else "", "",
paste(tt, "vs", time_reference), result)
}
if (has_by && length(by_display) == 2L) {
X00 <- xrow(time_reference, by_display[1L])
X01 <- xrow(tt, by_display[1L])
X10 <- xrow(time_reference, by_display[2L])
X11 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, (X11 - X10) - (X01 - X00), level)
add("change_difference", tt, by_display[1L], by_display[2L],
paste0("(", by_display[2L], " change) vs (", by_display[1L], " change)"),
result)
}
}
}
if (isTRUE(pairwise)) {
if (length(time_display) > 1L) {
tp <- utils::combn(time_display, 2L, simplify = FALSE)
groups <- if (has_by) by_display else ""
for (pair in tp) {
for (g in groups) {
Xa <- xrow(pair[1L], if (has_by) g else NULL)
Xb <- xrow(pair[2L], if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, Xb - Xa, level)
add("pairwise_time", pair[2L], if (has_by) g else "", "",
paste(pair[2L], "vs", pair[1L]), result)
}
}
}
if (has_by && length(by_display) > 1L) {
gp <- utils::combn(by_display, 2L, simplify = FALSE)
for (tt in time_display) {
for (pair in gp) {
Xa <- xrow(tt, pair[1L])
Xb <- xrow(tt, pair[2L])
result <- .r4vn_long_contrast(fit, Xb - Xa, level)
add("pairwise_group", tt, pair[1L], pair[2L],
paste(pair[2L], "vs", pair[1L]), result)
}
}
}
}
}
if (!length(out)) {
return(data.frame(
type = character(), time = character(), group1 = character(),
group2 = character(), comparison = character(),
estimate = numeric(), lower = numeric(), upper = numeric(), p = numeric(),
p_adjusted = numeric(), stringsAsFactors = FALSE
))
}
result <- do.call(rbind, out)
if (!identical(adjust, "none")) {
result$p_adjusted <- stats::p.adjust(result$p, method = adjust)
} else {
result$p_adjusted <- result$p
}
result
}
.r4vn_long_summary_continuous <- function(x, type, digits, show_n) {
y <- suppressWarnings(as.numeric(x))
y <- y[is.finite(y)]
n <- length(y)
if (!n) return("")
if (identical(type, "median")) {
q <- stats::quantile(y, c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 7)
text <- paste0(
.r4vn_long_fmt(q[2L], digits), " (",
.r4vn_long_fmt(q[1L], digits), ", ",
.r4vn_long_fmt(q[3L], digits), ")"
)
} else if (identical(type, "full")) {
q <- stats::quantile(y, c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 7)
text <- paste0(
.r4vn_long_fmt(mean(y), digits), " (", .r4vn_long_fmt(stats::sd(y), digits), "); ",
.r4vn_long_fmt(q[2L], digits), " [", .r4vn_long_fmt(q[1L], digits), ", ",
.r4vn_long_fmt(q[3L], digits), "]; ",
.r4vn_long_fmt(min(y), digits), "-", .r4vn_long_fmt(max(y), digits)
)
} else {
text <- paste0(
.r4vn_long_fmt(mean(y), digits),
" (",
.r4vn_long_fmt(stats::sd(y), digits),
")"
)
}
if (isTRUE(show_n)) text <- paste0(text, " [n=", n, "]")
text
}
.r4vn_long_summary_binary <- function(x, event, digits) {
keep <- !is.na(x)
n <- sum(keep)
if (!n) return("")
cases <- sum(as.character(x[keep]) == as.character(event))
pct <- 100 * cases / n
paste0(cases, "/", n, " (", .r4vn_long_fmt(pct, digits), "%)")
}
.r4vn_long_summary_count <- function(x, digits, show_n) {
y <- suppressWarnings(as.numeric(x))
y <- y[is.finite(y)]
n <- length(y)
if (!n) return("")
text <- paste0(.r4vn_long_fmt(mean(y), digits), " (", .r4vn_long_fmt(stats::sd(y), digits), ")")
if (isTRUE(show_n)) text <- paste0(text, " [n=", n, "]")
text
}
.r4vn_long_observed_cell <- function(data, time_value, group_value, outcome_type,
summary_type, event, digits, show_n,
time_continuous) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(time_value)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(time_value)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(group_value)
}
x <- if (".outcome_display" %in% names(data)) {
data$.outcome_display[mask]
} else {
data$.outcome[mask]
}
if (identical(outcome_type, "continuous")) {
return(.r4vn_long_summary_continuous(x, summary_type, digits, show_n))
}
if (identical(outcome_type, "binary")) {
return(.r4vn_long_summary_binary(x, event, digits))
}
.r4vn_long_summary_count(x, digits, show_n)
}
.r4vn_long_effect_name <- function(outcome_type, effect_type) {
if (identical(outcome_type, "continuous")) return("Difference (95% CI)")
if (identical(outcome_type, "count")) return("IRR (95% CI)")
paste0(effect_type, " (95% CI)")
}
.r4vn_long_change_name <- function(outcome_type, effect_type) {
if (identical(outcome_type, "continuous")) return("Change (95% CI)")
if (identical(outcome_type, "count")) return("IRR (95% CI)")
paste0(effect_type, " (95% CI)")
}
.r4vn_long_build_rows <- function(data, outcome_label, outcome_type, summary_type,
event, effect_type, time_display, time_reference,
by_display, time_continuous, contrasts, tests,
change, digits, effect_digits, p_digits,
show_n) {
has_by <- ".by" %in% names(data)
two_groups <- has_by && length(by_display) == 2L
effect_header <- .r4vn_long_effect_name(outcome_type, effect_type)
rows <- list()
rid <- 0L
add_row <- function(values) {
rid <<- rid + 1L
rows[[rid]] <<- values
}
if (!has_by) {
for (tt in time_display) {
cell <- .r4vn_long_observed_cell(
data, tt, NULL, outcome_type, summary_type, event,
digits, show_n, time_continuous
)
if (isTRUE(time_continuous)) {
effect_text <- ""
p_text <- ""
} else {
effect_text <- if (as.character(tt) == as.character(time_reference)) {
"Ref"
} else {
hit <- contrasts[
contrasts$type == "change" &
contrasts$time == as.character(tt),
, drop = FALSE
]
if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
}
hitp <- contrasts[
contrasts$type == "change" &
contrasts$time == as.character(tt),
, drop = FALSE
]
p_text <- if (nrow(hitp)) .r4vn_long_fmt_p(hitp$p_adjusted[1L], p_digits) else ""
}
add_row(c(
Outcome = outcome_label,
Time = as.character(tt),
Summary = cell,
Effect = effect_text,
`p-value` = p_text
))
}
if (isTRUE(time_continuous)) {
hit <- contrasts[contrasts$type == "slope", , drop = FALSE]
add_row(c(
Outcome = "",
Time = "Change per 1 time unit",
Summary = "",
Effect = if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else "",
`p-value` = if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
))
}
add_row(c(
Outcome = "",
Time = "Overall time",
Summary = "",
Effect = "",
`p-value` = .r4vn_long_fmt_p(tests$p_time, p_digits)
))
} else {
for (tt in time_display) {
values <- c(Outcome = outcome_label, Time = as.character(tt))
for (g in by_display) {
values <- c(
values,
stats::setNames(
.r4vn_long_observed_cell(
data, tt, g, outcome_type, summary_type, event,
digits, show_n, time_continuous
),
g
)
)
}
if (two_groups) {
hit <- contrasts[contrasts$type == "between" & contrasts$time == as.character(tt), , drop = FALSE]
if (nrow(hit)) {
values <- c(
values,
stats::setNames(
.r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits),
effect_header
),
`p-value` = .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits)
)
} else {
values <- c(values, stats::setNames("", effect_header), `p-value` = "")
}
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
if (isTRUE(time_continuous)) {
slope_rows <- contrasts[contrasts$type == "slope", , drop = FALSE]
if (nrow(slope_rows)) {
values <- c(Outcome = "", Time = "Change per 1 time unit")
for (g in by_display) {
hit <- slope_rows[slope_rows$group1 == g, , drop = FALSE]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
values <- c(values, stats::setNames(txt, g))
}
if (two_groups) {
hit <- contrasts[contrasts$type == "slope_difference", , drop = FALSE]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
pv <- if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
values <- c(values, stats::setNames(txt, effect_header), `p-value` = pv)
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
} else if (isTRUE(change)) {
follow <- setdiff(as.character(time_display), as.character(time_reference))
for (tt in follow) {
values <- c(Outcome = "", Time = paste0("Change: ", tt, " vs ", time_reference))
for (g in by_display) {
hit <- contrasts[
contrasts$type == "change" &
contrasts$time == tt &
contrasts$group1 == g,
, drop = FALSE
]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
values <- c(values, stats::setNames(txt, g))
}
if (two_groups) {
hit <- contrasts[
contrasts$type == "change_difference" &
contrasts$time == tt,
, drop = FALSE
]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
pv <- if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
values <- c(values, stats::setNames(txt, effect_header), `p-value` = pv)
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
}
test_rows <- list(
c("Overall time", tests$p_time),
c("Overall group", tests$p_group),
c("Time x group", tests$p_interaction)
)
for (z in test_rows) {
values <- c(Outcome = "", Time = z[1L])
for (g in by_display) values <- c(values, stats::setNames("", g))
if (two_groups) values <- c(values, stats::setNames("", effect_header))
values <- c(values, `p-value` = .r4vn_long_fmt_p(as.numeric(z[2L]), p_digits))
add_row(values)
}
}
out <- as.data.frame(do.call(rbind, rows), stringsAsFactors = FALSE, check.names = FALSE)
rownames(out) <- NULL
if ("Effect" %in% names(out)) names(out)[names(out) == "Effect"] <- .r4vn_long_change_name(outcome_type, effect_type)
out
}
.r4vn_long_html_table <- function(data, title, notes, bold_p, p_bold) {
cols <- names(data)
header <- paste0("<th>", .r4vn_long_escape(cols), "</th>", collapse = "")
body <- character(nrow(data))
last_outcome <- NULL
for (i in seq_len(nrow(data))) {
cells <- as.character(data[i, , drop = TRUE])
outcome <- cells[1L]
if (nzchar(outcome) && identical(outcome, last_outcome)) cells[1L] <- ""
if (nzchar(outcome)) last_outcome <- outcome
td <- character(length(cells))
for (j in seq_along(cells)) {
value <- .r4vn_long_escape(cells[j])
class_attr <- ""
if (j == 1L && nzchar(cells[j])) {
value <- paste0('<span class="variable-name">', value, "</span>")
}
if (identical(cols[j], "p-value") && isTRUE(bold_p)) {
raw <- suppressWarnings(as.numeric(sub("^<", "", cells[j])))
if (is.finite(raw) && raw < p_bold) value <- paste0("<strong>", value, "</strong>")
}
td[j] <- paste0("<td", class_attr, ">", value, "</td>")
}
row_class <- if (grepl("^(Overall|Time x group)", cells[2L])) ' class="diagnostic-row"' else ""
body[i] <- paste0("<tr", row_class, ">", paste0(td, collapse = ""), "</tr>")
}
note_html <- if (length(notes)) {
paste0('<div class="model-note">', .r4vn_long_escape(notes), "</div>", collapse = "")
} else ""
paste0(
'<section class="r4vn-table">',
if (!is.null(title) && nzchar(title)) paste0('<div class="table-title">', .r4vn_long_escape(title), "</div>") else "",
"<table><thead><tr>", header, "</tr></thead><tbody>",
paste(body, collapse = ""),
"</tbody></table>",
note_html,
"</section>"
)
}
.r4vn_long_html_document <- function(table_html) {
css <- paste0(
"body{font-family:Arial,'Helvetica Neue',sans-serif;background:#fff;color:#111;margin:18px}",
".table-title{font-size:18px;font-weight:700;margin:0 0 8px}.section-title{font-size:15px;font-weight:700;margin:18px 0 6px}",
".r4vn-table{margin-bottom:28px;overflow-x:auto}",
".r4vn-table table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #111;border-bottom:2px solid #111}",
".r4vn-table th{padding:5px 9px;border-bottom:1.5px solid #111;font-weight:700;white-space:nowrap;text-align:center;background:#fff}",
".r4vn-table th:first-child{text-align:left;min-width:220px}",
".r4vn-table td{padding:4px 9px;vertical-align:top;border:0;white-space:nowrap;text-align:center}",
".r4vn-table td:first-child{text-align:left}",
".variable-name{font-weight:700}",
".diagnostic-row td{border-top:1px solid #aaa;font-style:italic}",
".model-note{font-size:12px;color:#333;margin-top:6px;line-height:1.35}",
".r4vn-plot{margin:18px 0 28px}.plot-wrap{max-width:1100px;overflow-x:auto}.plot-wrap svg{display:block;max-width:100%;height:auto;background:#fff}"
)
paste0(
'<!DOCTYPE html><html><head><meta charset="UTF-8">',
'<meta name="viewport" content="width=device-width,initial-scale=1">',
"<style>", css, "</style></head><body>",
table_html,
"</body></html>"
)
}
.r4vn_long_show_html <- function(file) {
if (!interactive()) return(invisible(file))
viewer <- getOption("viewer")
normalized <- normalizePath(file, winslash = "/", mustWork = FALSE)
if (is.function(viewer)) viewer(normalized) else utils::browseURL(normalized)
invisible(file)
}
.r4vn_long_descriptive_rows <- function(data, outcome_label, outcome_type,
summary_type, event, time_display,
by_display, time_continuous,
digits, show_n) {
groups <- if (".by" %in% names(data)) by_display else "Overall"
rows <- list()
k <- 0L
for (tt in time_display) {
for (gg in groups) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(tt)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(tt)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(gg)
}
# Work with row indices inside the current time x group cell. `observed`
# is cell-sized, whereas `mask` is data-sized; combining them directly
# would recycle vectors and can produce incorrect subject counts/warnings.
idx <- which(mask)
x <- if (".outcome_display" %in% names(data)) data$.outcome_display[idx] else data$.outcome[idx]
observed <- !is.na(x)
if (is.numeric(x)) observed <- observed & is.finite(x)
observed_idx <- idx[observed]
k <- k + 1L
rows[[k]] <- data.frame(
Outcome = outcome_label,
Time = as.character(tt),
Group = as.character(gg),
N = sum(observed),
Missing = sum(!observed),
Subjects = if (length(observed_idx)) {
length(unique(data$.id[observed_idx][!is.na(data$.id[observed_idx])]))
} else {
0L
},
Summary = .r4vn_long_observed_cell(
data, tt, if (".by" %in% names(data)) gg else NULL,
outcome_type, summary_type, event, digits, show_n, time_continuous
),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_tests_table <- function(results, digits = 3L) {
if (!length(results)) return(data.frame())
rows <- lapply(results, function(z) {
data.frame(
Outcome = z$label,
Test = c("Overall time", "Overall group", "Time x group"),
p = c(z$tests$p_time, z$tests$p_group, z$tests$p_interaction),
stringsAsFactors = FALSE,
check.names = FALSE
)
})
out <- do.call(rbind, rows)
out <- out[is.finite(out$p), , drop = FALSE]
rownames(out) <- NULL
if (!nrow(out)) return(out)
out$`p-value` <- vapply(out$p, .r4vn_long_fmt_p, character(1), digits = digits)
out$p <- NULL
out
}
.r4vn_long_contrasts_table <- function(results, effect_digits = 2L, p_digits = 3L) {
if (!length(results)) return(data.frame())
rows <- list()
k <- 0L
for (z in results) {
cc <- z$contrasts
if (is.null(cc) || !nrow(cc)) next
for (i in seq_len(nrow(cc))) {
k <- k + 1L
type_label <- switch(
as.character(cc$type[i]),
between = "Between groups",
change = "Change from reference time",
change_difference = "Difference in change",
pairwise_time = "Pairwise time comparison",
pairwise_group = "Pairwise group comparison",
slope = "Slope",
slope_difference = "Difference in slopes",
as.character(cc$type[i])
)
desc <- as.character(cc$comparison[i])
if (!nzchar(desc)) desc <- as.character(cc$type[i])
if (as.character(cc$type[i]) %in% c("between", "pairwise_group") && nzchar(as.character(cc$time[i]))) {
desc <- paste0(desc, " at ", cc$time[i])
}
if (as.character(cc$type[i]) %in% c("change", "pairwise_time") && nzchar(as.character(cc$group1[i]))) {
desc <- paste0(desc, " in ", cc$group1[i])
}
rows[[k]] <- data.frame(
Outcome = z$label,
Type = type_label,
Contrast = desc,
Effect = z$effect,
`Estimate (95% CI)` = .r4vn_long_fmt_ci(cc$estimate[i], cc$lower[i], cc$upper[i], effect_digits),
`p-value` = .r4vn_long_fmt_p(cc$p_adjusted[i], p_digits),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_diagnostics_one <- function(fit, model_data, outcome_label,
outcome_type, time_display,
time_continuous, repeated) {
model <- fit$fit
n_obs <- .r4vn_long_safe_scalar(tryCatch(stats::nobs(model), error = function(e) nrow(model_data)))
if (!is.finite(n_obs)) n_obs <- nrow(model_data)
subject_ok <- !is.na(model_data$.id)
n_subjects <- length(unique(model_data$.id[subject_ok]))
complete_subjects <- NA_integer_
if (isTRUE(repeated) && !isTRUE(time_continuous) && length(time_display)) {
z <- split(as.character(model_data$.time[subject_ok]), as.character(model_data$.id[subject_ok]))
complete_subjects <- sum(vapply(z, function(x) length(unique(x)) == length(time_display), logical(1)))
}
convergence <- "OK"
if (identical(fit$engine, "mixed")) {
if (inherits(model, "lme")) convergence <- "OK"
} else if (identical(fit$engine, "gee")) {
err <- .r4vn_long_safe_scalar(tryCatch(model$geese$error, error = function(e) NA_real_))
if (is.finite(err) && err != 0) convergence <- paste0("GEE error code ", err)
} else if (inherits(model, "glm") && isFALSE(model$converged)) {
convergence <- "Model did not converge"
}
# AIC/BIC are not reported for GEE or cluster-robust marginal inference.
# In particular, geeglm may return numeric(0), which previously triggered
# `if (is.finite(aic))` with "argument is of length zero".
if (fit$engine %in% c("gee", "cluster_robust")) {
aic <- NA_real_
bic <- NA_real_
} else {
aic <- .r4vn_long_safe_scalar(tryCatch(stats::AIC(model), error = function(e) NA_real_))
bic <- .r4vn_long_safe_scalar(tryCatch(stats::BIC(model), error = function(e) NA_real_))
}
icc <- NA_real_
if (identical(fit$engine, "mixed") && identical(outcome_type, "continuous")) {
if (inherits(model, "lme")) {
vc <- tryCatch(nlme::VarCorr(model), error = function(e) NULL)
if (!is.null(vc)) {
vals <- suppressWarnings(as.numeric(vc[, "Variance"]))
rn <- rownames(vc)
vi <- vals[grepl("Intercept", rn, fixed = TRUE)]
vr <- vals[grepl("Residual", rn, fixed = TRUE)]
if (length(vi) && length(vr) && is.finite(vi[1L] + vr[length(vr)]) && vi[1L] + vr[length(vr)] > 0) {
icc <- vi[1L] / (vi[1L] + vr[length(vr)])
}
}
}
}
engine_label <- switch(
fit$engine,
mixed = "Linear mixed model",
gee = "GEE",
cluster_robust = "Marginal robust regression",
independent = "Independent regression",
fit$engine
)
data.frame(
Outcome = outcome_label,
Type = outcome_type,
Engine = engine_label,
Correlation = if (!is.null(fit$correlation)) fit$correlation else "",
Observations = as.integer(n_obs),
Subjects = as.integer(n_subjects),
`Complete subjects` = complete_subjects,
Singular = if (identical(fit$engine, "mixed")) isTRUE(fit$singular) else NA,
Convergence = convergence,
AIC = if (is.finite(aic)) round(aic, 2) else NA_real_,
BIC = if (is.finite(bic)) round(bic, 2) else NA_real_,
ICC = if (is.finite(icc)) round(icc, 3) else NA_real_,
stringsAsFactors = FALSE,
check.names = FALSE
)
}
.r4vn_long_interpretation <- function(results, alpha = 0.05, p_digits = 3L) {
if (!length(results)) return(data.frame())
rows <- list()
k <- 0L
add <- function(outcome, section, text) {
k <<- k + 1L
rows[[k]] <<- data.frame(
Outcome = outcome, Section = section, Interpretation = text,
stringsAsFactors = FALSE, check.names = FALSE
)
}
for (z in results) {
pt <- z$tests$p_time
pg <- z$tests$p_group
pi <- z$tests$p_interaction
if (is.finite(pi)) {
add(
z$label, "Time x group",
if (pi < alpha) {
paste0("The change over time differs between groups (p=", .r4vn_long_fmt_p(pi, p_digits),
"). Interpret the time-specific and change contrasts rather than the main effects alone.")
} else {
paste0("There is no statistical evidence that the temporal pattern differs between groups (p=",
.r4vn_long_fmt_p(pi, p_digits), ").")
}
)
}
if (is.finite(pt)) {
add(
z$label, "Time",
paste0(if (pt < alpha) "There is evidence of an overall time effect" else "There is no statistical evidence of an overall time effect",
" (p=", .r4vn_long_fmt_p(pt, p_digits), ").")
)
}
if (is.finite(pg)) {
add(
z$label, "Group",
paste0(if (pg < alpha) "There is evidence of an overall group effect" else "There is no statistical evidence of an overall group effect",
" (p=", .r4vn_long_fmt_p(pg, p_digits), ").")
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_plot_rows <- function(data, outcome_label, outcome_type, event,
time_display, by_display, time_continuous,
level = 0.95) {
groups <- if (".by" %in% names(data)) by_display else "Overall"
alpha <- 1 - level
rows <- list()
k <- 0L
for (tt in time_display) {
for (gg in groups) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(tt)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(tt)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(gg)
}
x_display <- if (".outcome_display" %in% names(data)) data$.outcome_display[mask] else data$.outcome[mask]
estimate <- lower <- upper <- NA_real_
n <- sum(!is.na(x_display))
measure <- "Observed value"
if (identical(outcome_type, "continuous")) {
y <- suppressWarnings(as.numeric(x_display))
y <- y[is.finite(y)]
n <- length(y)
if (n) {
estimate <- mean(y)
se <- if (n > 1L) stats::sd(y) / sqrt(n) else NA_real_
crit <- if (n > 1L) stats::qt(1 - alpha / 2, df = n - 1L) else NA_real_
if (is.finite(se) && is.finite(crit)) {
lower <- estimate - crit * se
upper <- estimate + crit * se
}
}
measure <- "Mean"
} else if (identical(outcome_type, "binary")) {
keep <- !is.na(x_display)
n <- sum(keep)
cases <- sum(as.character(x_display[keep]) == as.character(event))
if (n) {
estimate <- 100 * cases / n
ci <- tryCatch(suppressWarnings(stats::prop.test(cases, n, conf.level = level, correct = FALSE)$conf.int),
error = function(e) c(NA_real_, NA_real_))
lower <- 100 * ci[1L]
upper <- 100 * ci[2L]
}
measure <- paste0("Event percentage (", event, ")")
} else {
y <- suppressWarnings(as.numeric(x_display))
keep <- is.finite(y)
n <- sum(keep)
if (".exposure" %in% names(data)) {
e <- suppressWarnings(as.numeric(data$.exposure[mask]))
good <- keep & is.finite(e) & e > 0
events <- sum(y[good])
pt <- sum(e[good])
if (pt > 0) {
estimate <- events / pt
lower <- if (events > 0) stats::qchisq(alpha / 2, 2 * events) / (2 * pt) else 0
upper <- stats::qchisq(1 - alpha / 2, 2 * (events + 1)) / (2 * pt)
}
measure <- "Incidence rate per 1 person-time"
} else {
yy <- y[keep]
if (length(yy)) {
estimate <- mean(yy)
se <- if (length(yy) > 1L) stats::sd(yy) / sqrt(length(yy)) else NA_real_
crit <- if (length(yy) > 1L) stats::qt(1 - alpha / 2, df = length(yy) - 1L) else NA_real_
if (is.finite(se) && is.finite(crit)) {
lower <- max(0, estimate - crit * se)
upper <- estimate + crit * se
}
}
measure <- "Mean count"
}
}
k <- k + 1L
rows[[k]] <- data.frame(
outcome = outcome_label,
time = if (isTRUE(time_continuous)) as.numeric(tt) else as.character(tt),
time_order = match(as.character(tt), as.character(time_display)),
group = as.character(gg),
estimate = estimate,
lower = lower,
upper = upper,
n = n,
measure = measure,
time_continuous = isTRUE(time_continuous),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_build_plot <- function(plot_data, ci = TRUE, plot_args = list()) {
if (is.null(plot_data) || !nrow(plot_data)) return(NULL)
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
if (length(plot_args) && is.null(names(plot_args))) stop("`plot_args` must be a named list.", call. = FALSE)
use_ci <- if (!is.null(plot_args$ci)) isTRUE(plot_args$ci) else isTRUE(ci)
structure(
list(data = plot_data, ci = use_ci, args = plot_args),
class = "r4vn_long_plot"
)
}
.r4vn_long_draw_plot <- function(plot_data, ci = TRUE, plot_args = list()) {
if (is.null(plot_data) || !nrow(plot_data)) return(invisible(NULL))
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
title <- if (!is.null(plot_args$title)) as.character(plot_args$title)[1L] else "Longitudinal profile"
xlab <- if (!is.null(plot_args$xlab)) as.character(plot_args$xlab)[1L] else "Time"
ylab <- if (!is.null(plot_args$ylab)) as.character(plot_args$ylab)[1L] else "Observed estimate (95% CI)"
line_width <- if (!is.null(plot_args$line_width)) as.numeric(plot_args$line_width)[1L] else 1.8
point_size <- if (!is.null(plot_args$point_size)) as.numeric(plot_args$point_size)[1L] else 1.05
base_size <- if (!is.null(plot_args$base_size)) as.numeric(plot_args$base_size)[1L] else 11
legend_position <- if (!is.null(plot_args$legend_position)) as.character(plot_args$legend_position)[1L] else "bottom"
font_family <- if (!is.null(plot_args$font_family)) as.character(plot_args$font_family)[1L] else "sans"
if (!nzchar(font_family)) font_family <- "sans"
use_ci <- if (!is.null(plot_args$ci)) isTRUE(plot_args$ci) else isTRUE(ci)
grid <- if (!is.null(plot_args$grid)) isTRUE(plot_args$grid) else TRUE
outcomes <- unique(as.character(plot_data$outcome))
n_out <- length(outcomes)
if (!n_out) return(invisible(NULL))
ncol_panels <- if (n_out <= 1L) 1L else if (n_out <= 4L) 2L else 3L
nrow_panels <- ceiling(n_out / ncol_panels)
old <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(old), add = TRUE)
graphics::par(
mfrow = c(nrow_panels, ncol_panels),
mar = c(4.3, 4.4, 3.2, 1.2),
oma = c(0, 0, if (n_out > 1L) 1.7 else 0, 0),
family = font_family,
cex = max(0.65, base_size / 11)
)
for (outcome in outcomes) {
pd <- plot_data[as.character(plot_data$outcome) == outcome, , drop = FALSE]
continuous_x <- all(pd$time_continuous)
if (continuous_x) {
x_all <- suppressWarnings(as.numeric(pd$time))
x_levels <- sort(unique(x_all[is.finite(x_all)]))
x_at <- x_levels
x_labels <- format(x_levels, trim = TRUE)
} else {
ord <- order(pd$time_order)
x_levels <- unique(as.character(pd$time[ord]))
x_all <- match(as.character(pd$time), x_levels)
x_at <- seq_along(x_levels)
x_labels <- x_levels
}
finite_y <- is.finite(pd$estimate)
y_values <- pd$estimate[finite_y]
if (use_ci) {
y_values <- c(y_values, pd$lower[is.finite(pd$lower)], pd$upper[is.finite(pd$upper)])
}
if (!length(y_values)) y_values <- c(0, 1)
yr <- range(y_values, finite = TRUE)
if (!all(is.finite(yr))) yr <- c(0, 1)
if (diff(yr) == 0) yr <- yr + c(-0.5, 0.5)
pad <- 0.07 * diff(yr)
ylim <- yr + c(-pad, pad)
if (continuous_x) {
xr <- range(x_all, finite = TRUE)
if (!all(is.finite(xr))) xr <- c(0, 1)
if (diff(xr) == 0) xr <- xr + c(-0.5, 0.5)
xlim <- xr
} else {
xlim <- c(0.6, max(1.4, length(x_levels) + 0.4))
}
panel_title <- if (n_out > 1L) outcome else title
graphics::plot(
NA_real_, NA_real_, type = "n", xlim = xlim, ylim = ylim,
xaxt = "n", xlab = xlab, ylab = ylab, main = panel_title,
bty = "l", las = 1
)
graphics::axis(1, at = x_at, labels = x_labels)
if (isTRUE(grid)) {
graphics::abline(h = graphics::axTicks(2), col = "grey90", lty = 1, lwd = 0.7)
}
groups <- unique(as.character(pd$group))
if (!length(groups)) groups <- "Overall"
colors <- plot_args$colors
if (is.null(colors) || !length(colors)) {
colors <- if (length(groups) == 1L) "black" else grDevices::hcl.colors(length(groups), palette = "Dark 3")
}
colors <- rep(colors, length.out = length(groups))
pch <- plot_args$point_shapes
if (is.null(pch) || !length(pch)) pch <- c(16, 17, 15, 18, 3, 4, 8)
pch <- rep(pch, length.out = length(groups))
lty <- plot_args$line_types
if (is.null(lty) || !length(lty)) lty <- seq_along(groups)
lty <- rep(lty, length.out = length(groups))
for (j in seq_along(groups)) {
gg <- groups[j]
z <- pd[as.character(pd$group) == gg, , drop = FALSE]
z <- z[order(z$time_order), , drop = FALSE]
xx <- if (continuous_x) suppressWarnings(as.numeric(z$time)) else match(as.character(z$time), x_levels)
good <- is.finite(xx) & is.finite(z$estimate)
if (sum(good) >= 1L) {
graphics::lines(xx[good], z$estimate[good], col = colors[j], lty = lty[j], lwd = line_width)
graphics::points(xx[good], z$estimate[good], col = colors[j], pch = pch[j], cex = point_size)
}
if (use_ci) {
ci_good <- is.finite(xx) & is.finite(z$lower) & is.finite(z$upper)
if (any(ci_good)) {
graphics::arrows(
xx[ci_good], z$lower[ci_good], xx[ci_good], z$upper[ci_good],
angle = 90, code = 3, length = 0.035,
col = colors[j], lwd = max(0.8, line_width * 0.65)
)
}
}
}
if (!(length(groups) == 1L && identical(groups, "Overall")) && !identical(legend_position, "none")) {
pos <- switch(
tolower(legend_position),
right = "topright", left = "topleft", top = "top", bottom = "bottom",
topright = "topright", topleft = "topleft", bottomright = "bottomright",
bottomleft = "bottomleft", "bottom"
)
graphics::legend(
pos, legend = groups, col = colors, lty = lty, pch = pch,
bty = "n", horiz = identical(pos, "bottom"), cex = 0.85
)
}
}
if (n_out > 1L && nzchar(title)) graphics::mtext(title, outer = TRUE, side = 3, line = 0.3, font = 2)
invisible(NULL)
}
.r4vn_long_base64 <- function(x) {
if (!length(x)) return("")
alphabet <- strsplit("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/", "", fixed = TRUE)[[1L]]
n <- length(x)
pad <- (3L - n %% 3L) %% 3L
vals <- as.integer(x)
if (pad) vals <- c(vals, rep(0L, pad))
m <- matrix(vals, ncol = 3L, byrow = TRUE)
i1 <- bitwShiftR(m[, 1L], 2L)
i2 <- bitwOr(bitwShiftL(bitwAnd(m[, 1L], 3L), 4L), bitwShiftR(m[, 2L], 4L))
i3 <- bitwOr(bitwShiftL(bitwAnd(m[, 2L], 15L), 2L), bitwShiftR(m[, 3L], 6L))
i4 <- bitwAnd(m[, 3L], 63L)
out <- as.vector(rbind(alphabet[i1 + 1L], alphabet[i2 + 1L], alphabet[i3 + 1L], alphabet[i4 + 1L]))
if (pad >= 1L) out[length(out)] <- "="
if (pad == 2L) out[length(out) - 1L] <- "="
paste0(out, collapse = "")
}
.r4vn_long_plot_html <- function(graph) {
if (is.null(graph) || !inherits(graph, "r4vn_long_plot")) return("")
pd <- graph$data
if (is.null(pd) || !nrow(pd)) return("")
n_out <- length(unique(as.character(pd$outcome)))
width <- if (!is.null(graph$args$viewer_width)) as.numeric(graph$args$viewer_width)[1L] else 9.2
height <- if (!is.null(graph$args$viewer_height)) as.numeric(graph$args$viewer_height)[1L] else {
if (n_out <= 1L) 5.2 else 4.4 * ceiling(n_out / if (n_out <= 4L) 2 else 3)
}
if (!is.finite(width) || width <= 0) width <- 9.2
if (!is.finite(height) || height <= 0) height <- 5.2
# Prefer vector SVG in Viewer. If the local R build cannot create SVG,
# fall back to an inline PNG data URI; both paths use only base/recommended R.
tf <- tempfile(pattern = "r4vn-tablong-plot-", fileext = ".svg")
opened <- FALSE
ok <- tryCatch({
svg_family <- if (!is.null(graph$args$font_family)) as.character(graph$args$font_family)[1L] else "sans"
if (!nzchar(svg_family)) svg_family <- "sans"
grDevices::svg(tf, width = width, height = height, pointsize = 11, family = svg_family)
opened <- TRUE
.r4vn_long_draw_plot(pd, ci = graph$ci, plot_args = graph$args)
grDevices::dev.off()
opened <- FALSE
TRUE
}, error = function(e) FALSE)
if (opened) try(grDevices::dev.off(), silent = TRUE)
plot_markup <- ""
if (ok && file.exists(tf)) {
lines <- readLines(tf, warn = FALSE, encoding = "UTF-8")
start <- grep("<svg", lines, fixed = TRUE)[1L]
if (length(start) && !is.na(start)) {
plot_markup <- paste(lines[start:length(lines)], collapse = "\n")
}
}
unlink(tf)
if (!nzchar(plot_markup)) {
png_file <- tempfile(pattern = "r4vn-tablong-plot-", fileext = ".png")
png_open <- FALSE
png_ok <- tryCatch({
grDevices::png(
png_file,
width = max(900L, as.integer(width * 120)),
height = max(560L, as.integer(height * 120)),
res = 120,
type = if (.Platform$OS.type == "windows") "windows" else "cairo"
)
png_open <- TRUE
.r4vn_long_draw_plot(pd, ci = graph$ci, plot_args = graph$args)
grDevices::dev.off()
png_open <- FALSE
TRUE
}, error = function(e) FALSE)
if (png_open) try(grDevices::dev.off(), silent = TRUE)
if (png_ok && file.exists(png_file)) {
raw <- readBin(png_file, what = "raw", n = file.info(png_file)$size)
plot_markup <- paste0(
'<img alt="Longitudinal plot" src="data:image/png;base64,',
.r4vn_long_base64(raw), '" />'
)
}
unlink(png_file)
}
if (!nzchar(plot_markup)) return("")
paste0(
'<section class="r4vn-plot">',
'<div class="section-title">Longitudinal plot</div>',
'<div class="plot-wrap">', plot_markup, '</div>',
'</section>'
)
}
.r4vn_long_html_simple_table <- function(data, title) {
if (is.null(data) || !is.data.frame(data) || !nrow(data)) return("")
x <- data
# Viewer readability: a one/few-row table with many columns is easier to
# inspect after transposition. Keep the returned R object unchanged; this is
# presentation-only and follows the same compact Viewer rule used elsewhere
# in R4VN.
if (ncol(x) >= 8L && nrow(x) <= 4L) {
labels <- if ("Outcome" %in% names(x)) as.character(x$Outcome) else paste0("Result ", seq_len(nrow(x)))
labels <- make.unique(ifelse(is.na(labels) | !nzchar(labels), paste0("Result ", seq_along(labels)), labels))
payload <- if ("Outcome" %in% names(x)) x[, setdiff(names(x), "Outcome"), drop = FALSE] else x
tx <- data.frame(Metric = names(payload), stringsAsFactors = FALSE, check.names = FALSE)
for (i in seq_len(nrow(payload))) {
tx[[labels[i]]] <- vapply(payload, function(v) as.character(v[i]), character(1))
}
x <- tx
}
for (nm in names(x)) {
if (is.numeric(x[[nm]])) {
x[[nm]] <- ifelse(is.na(x[[nm]]), "", format(x[[nm]], trim = TRUE, scientific = FALSE))
} else {
x[[nm]][is.na(x[[nm]])] <- ""
}
}
header <- paste0("<th>", .r4vn_long_escape(names(x)), "</th>", collapse = "")
body <- vapply(seq_len(nrow(x)), function(i) {
vals <- vapply(as.list(x[i, , drop = FALSE]), function(z) as.character(z[1L]), character(1))
paste0("<tr>", paste0("<td>", .r4vn_long_escape(vals), "</td>", collapse = ""), "</tr>")
}, character(1))
paste0(
'<section class="r4vn-table secondary-table">',
'<div class="section-title">', .r4vn_long_escape(title), '</div>',
'<table><thead><tr>', header, '</tr></thead><tbody>',
paste0(body, collapse = ""), '</tbody></table></section>'
)
}
#' Longitudinal and Repeated-Measures Analysis
#'
#' Performs publication-ready longitudinal or repeated-measures analysis from
#' either long or wide data. `tablong()` is designed for the usual biomedical
#' workflow: describe each time point, test overall time and group effects,
#' test the time-by-group interaction, estimate clinically interpretable
#' contrasts, retain fitted models for advanced use, and optionally create a
#' longitudinal profile plot and a cautious interpretation table.
#'
#' The interface follows the R4VN principle of keeping routine analysis simple.
#' In most studies the essential call is only `vars()`, `time`, `id`, and
#' optionally `by`. Repeated continuous outcomes use a random-intercept model
#' when R's recommended `nlme` package is available. Repeated binary/count
#' outcomes use marginal regression with subject-clustered robust standard errors
#' calculated internally by R4VN, so routine analyses need no extra package.
#'
#' @param data Optional data frame. If `NULL`, the active data selected by
#' `usedf()` are used.
#' @param vars Outcome specification created by `vars()`. In long data, several
#' outcomes can be analyzed in one call. Continuous outcomes use the `c.`
#' prefix, `q.` requests median (IQR) descriptive display, and `f.` requests a
#' fuller continuous summary. An unprefixed two-level variable is treated as
#' binary. In wide data, the variables in `vars()` are repeated measurements
#' of the same outcome.
#' @param time In long data, an unquoted time variable. Use `c.timevar` for a
#' continuous linear time effect or `b2.timevar`, `b3.timevar`, etc. to select
#' the categorical reference level. In wide data, provide display labels such
#' as `c("Baseline", "Month 3", "Month 6")`; when omitted, the repeated
#' variable names are used as time labels.
#' @param id Optional subject identifier. In wide data it is optional because
#' each source row represents one subject. In long data, repeated IDs trigger
#' longitudinal analysis with within-subject correlation. If no ID is supplied,
#' or IDs do not repeat across time, observations are analyzed as repeated
#' cross-sectional samples.
#' @param by Optional grouping variable, for example treatment group. Prefix
#' `b2.`, `b3.`, etc. selects the reference group.
#' @param ref Optional categorical reference time label. This overrides a `bN.`
#' prefix supplied in `time`.
#' @param event Event level for binary outcomes. A single value applies to every
#' binary outcome; a named vector can specify a different event for each
#' outcome, for example `event = c(controlled = "Yes", admitted = "Yes")`.
#' @param adjusted Optional adjustment variables created by `vars()`. Use the
#' same R4VN prefixes as elsewhere, for example `vars(c.age, sex, b2.site)`.
#' @param gee Logical. Request a population-average marginal model for repeated
#' data. By default R4VN uses working-independence regression with subject-
#' clustered robust sandwich standard errors calculated internally. No extra
#' package is required. Set `ar1 = TRUE` to request AR(1) GEE through optional
#' package `geepack` when it is installed.
#' @param ar1 Logical. Request an AR(1) working correlation for a marginal GEE.
#' This advanced option uses optional package `geepack`. If it is unavailable,
#' `tablong()` warns and falls back to working-independence cluster-robust
#' inference instead of stopping the analysis.
#' @param slope Logical. With a mixed model and continuous time
#' (`time = c.month`), include a subject-specific random linear time slope in
#' addition to the random intercept.
#' @param or,rr,pr Logical effect switches for binary outcomes. Odds ratio is
#' the default when none is selected. `rr = TRUE` reports risk ratios and
#' `pr = TRUE` reports prevalence ratios using modified Poisson regression with
#' robust variance; repeated subjects use subject-clustered robust variance.
#' No external sandwich/GEE package is needed unless `ar1 = TRUE` is requested.
#' Only one switch may be `TRUE`.
#' @param count Logical. Treat numeric outcomes as non-negative counts and fit
#' Poisson models. Count outcomes report incidence-rate ratios (IRR).
#' @param exposure Optional positive exposure/person-time variable for count
#' models. In wide data it may also be `vars(exp0, exp1, ...)`, with one
#' exposure variable per repeated count variable.
#' @param change Logical. For categorical time, report change from the reference
#' time. With exactly two groups, also report the difference in change, i.e.
#' the usual difference-in-differences contrast. Default `TRUE`.
#' @param pairwise Logical. Calculate all available time and group pairwise
#' contrasts and retain them in `$contrasts` and `$contrasts_table`. The main
#' publication table remains compact. Default `FALSE`.
#' @param adjust Multiplicity adjustment applied to contrast p-values. Any method
#' accepted by `p.adjust()` may be used, including `"none"`, `"holm"`,
#' `"bonferroni"`, and `"BH"`.
#' @param missing Logical. Append cell-specific `n` to continuous/count summary
#' cells. Binary cells always show event/total. Detailed observed and missing
#' counts are always available in `$descriptive`.
#' @param level Confidence level, default 0.95.
#' @param digit Decimal places for descriptive summaries.
#' @param p_digit Decimal places for p-values.
#' @param effect_digit Decimal places for model effects and confidence intervals.
#' @param bold_p Logical. Bold p-values below `p_bold` in the HTML Viewer.
#' @param p_bold Threshold used by `bold_p`.
#' @param diagnostics Logical. Include the compact model-diagnostics table in the
#' HTML Viewer. Diagnostics are always retained in `$diagnostics`; default
#' `FALSE` keeps the primary Viewer concise.
#' @param diagnosis Logical singular alias for `diagnostics`, provided for consistency with other R4VN regression commands. Default `FALSE`. When explicitly supplied, it overrides `diagnostics`; omitting it preserves backward-compatible use of `diagnostics`.
#' @param interpretation Logical. Add a cautious deterministic interpretation
#' table. The default is `FALSE`. The interpretation emphasizes the
#' time-by-group interaction when present and does not replace substantive or
#' clinical interpretation by the researcher.
#' @param plot Logical. Create an observed longitudinal profile plot with 95%
#' confidence intervals, include the same plot directly in the HTML Viewer,
#' display it in the Plot pane, and store its specification in `$graph` and
#' `$plots$trajectory`. Plotting uses base R graphics; `ggplot2` is not required.
#' Default `FALSE`.
#' @param plot_args Named list controlling the profile plot. Supported entries
#' include `title`, `xlab`, `ylab`, `ci`, `line_width`, `point_size`,
#' `base_size`, `legend_position`, `font_family`, `colors`, `point_shapes`,
#' `line_types`, and `grid`. Generic `sans` is the default font for reliable
#' display in RStudio Viewer, browsers, Windows, macOS, and Linux.
#' @param name Logical. Display the original variable name after its variable
#' label in the main table.
#' @param title Optional table title.
#' @param file Optional HTML file path. When omitted, a temporary HTML file is
#' created. This file is the formatted Viewer report, not a replacement for
#' `tabexport()`.
#' @param raw Logical retained for backward compatibility. Raw models, tests,
#' contrasts, standardized long data, and reporting tables are always retained
#' in the returned object.
#' @param show Logical. Open the formatted HTML report in the Viewer/browser.
#' Default `TRUE`.
#'
#' @details
#' **Data format.** If `time` names a column in `data`, input is treated as long.
#' If `vars()` contains multiple repeated variables and `time` is a vector of
#' labels (or omitted), input is treated as wide and is reshaped internally.
#' The original data frame is never modified.
#'
#' **Continuous outcomes.** Repeated subjects use a random-intercept linear
#' mixed model through R's recommended `nlme` package when available.
#' `slope = TRUE` adds a random linear time slope when time is continuous. If
#' `nlme` is unavailable, R4VN falls back to a marginal linear model with
#' subject-clustered robust standard errors. Repeated cross-sectional data use
#' ordinary linear models. `gee = TRUE` explicitly requests the marginal model.
#'
#' **Binary outcomes.** The default effect is an odds ratio from logistic
#' regression. For repeated subjects, R4VN calculates subject-clustered robust
#' standard errors internally. `rr = TRUE` and `pr = TRUE` use modified Poisson
#' regression with robust variance and report RR or PR. This avoids requiring
#' `lme4`, `sandwich`, or `geepack` for routine binary longitudinal analysis.
#'
#' **Count outcomes.** `count = TRUE` fits a Poisson model and reports IRR.
#' Repeated subjects use subject-clustered robust standard errors calculated
#' internally. Supplying `exposure` adds `offset(log(exposure))` and the observed
#' plot is an incidence rate per one person-time unit.
#'
#' **Categorical time.** The model includes time, group when supplied, and the
#' time-by-group interaction. The main table reports observed summaries at each
#' time. With two groups it also reports the between-group effect at each time,
#' within-group change from the reference time, and the difference in change.
#' The omnibus `Time x group` p-value is the formal test that temporal changes
#' differ between groups.
#'
#' **Continuous time.** `time = c.month` estimates change per one time unit. With
#' a group variable, group-specific slopes and their difference are returned.
#'
#' **More than two groups.** The main table remains intentionally compact and
#' shows omnibus tests. Set `pairwise = TRUE` to obtain all model-based pairwise
#' comparisons in `$contrasts_table`; use `adjust` to control multiplicity.
#'
#' **Descriptive prefixes.** `q.` and `f.` affect the observed descriptive
#' summary only. Inferential effects remain based on the selected mean model;
#' they do not fit median regression.
#'
#' **Missing values.** Each model uses observations complete for that outcome,
#' time, ID/group, exposure if required, and adjustment variables. Descriptive
#' counts are retained separately so that missingness and attrition can be
#' reviewed before publication.
#'
#' **Package dependencies.** Routine `tablong()` analyses and plots are designed
#' to run with base/recommended R only. `nlme` is used for continuous random-
#' effects models and is bundled with standard R installations. `geepack` is
#' optional and used only when an AR(1) GEE is explicitly requested. `lme4`,
#' `sandwich`, `broom`, and `ggplot2` are not required by `tablong()`.
#'
#' **Returned reporting contract.** For programmatic reuse, the object includes
#' a flat main table plus descriptive, omnibus-test, contrast, diagnostics,
#' interpretation, plot, model, and metadata components. The object also
#' inherits from `r4vn_tab`, so existing `tabexport()` workflows continue to
#' work.
#'
#' @return Invisibly returns an object of class
#' `c("r4vn_tablong", "r4vn_tab")`. Important components are:
#' `$data` (main publication table), `$descriptive`, `$tests` and
#' `$tests_table`, `$contrasts` and `$contrasts_table`, `$diagnostics`,
#' `$interpretation`, `$tables`, `$graph`, `$plots`, `$models`, `$long_data`,
#' `$metadata`, `$html`, `$file`, `$results`, and `$call`.
#'
#' @seealso `vars()`, `tab()`, `tabexport()`, `usedf()`
#' @family R4VN tables
#' @export
#'
#' @examples
#' # -------------------------------------------------------------------------
#' # 1. Repeated cross-sectional continuous outcome: no optional package needed
#' # -------------------------------------------------------------------------
#' set.seed(11)
#' d <- data.frame(
#' period = factor(rep(c("Before", "After"), each = 80),
#' levels = c("Before", "After")),
#' group = factor(rep(rep(c("Control", "Intervention"), each = 40), 2)),
#' age = rnorm(160, 45, 10)
#' )
#' d$score <- 50 + 2 * (d$period == "After") +
#' 5 * (d$group == "Intervention") +
#' 4 * (d$period == "After" & d$group == "Intervention") +
#' 0.15 * d$age + rnorm(160, 0, 8)
#'
#' z <- tablong(
#' d, vars = vars(c.score), time = period, by = group,
#' adjusted = vars(c.age), show = FALSE
#' )
#' z$data
#' z$tests_table
#' z$contrasts_table
#'
#' \donttest{
#' # -----------------------------------------------------------------------
#' # 2. The same analysis with interpretation and a publication profile plot
#' # -----------------------------------------------------------------------
#' z2 <- tablong(
#' d, vars = vars(c.score), time = period, by = group,
#' adjusted = vars(c.age), interpretation = TRUE,
#' plot = TRUE, show = FALSE
#' )
#' z2$interpretation
#' z2$diagnostics
#' plot(z2)
#'
#' # -----------------------------------------------------------------------
#' # 3. No comparison group: change over time only
#' # -----------------------------------------------------------------------
#' z_time <- tablong(
#' d, vars = vars(c.score), time = period,
#' adjusted = vars(c.age), show = FALSE
#' )
#' z_time$tests_table
#'
#' # -----------------------------------------------------------------------
#' # 4. Change the reference time by value or by bN. prefix
#' # -----------------------------------------------------------------------
#' z_ref1 <- tablong(d, vars = vars(c.score), time = period,
#' by = group, ref = "After", show = FALSE)
#' z_ref2 <- tablong(d, vars = vars(c.score), time = b2.period,
#' by = group, show = FALSE)
#'
#' # -----------------------------------------------------------------------
#' # 5. Median/IQR or full descriptive display, while inference remains a
#' # mean model
#' # -----------------------------------------------------------------------
#' z_median <- tablong(d, vars = vars(q.score), time = period,
#' by = group, show = FALSE)
#' z_full <- tablong(d, vars = vars(f.score), time = period,
#' by = group, missing = TRUE, show = FALSE)
#'
#' # -----------------------------------------------------------------------
#' # 6. Long repeated data: linear mixed model
#' # -----------------------------------------------------------------------
#' set.seed(12)
#' n_subject <- 60
#' dl <- expand.grid(
#' id = seq_len(n_subject),
#' visit = factor(c("Baseline", "Month 3", "Month 6"),
#' levels = c("Baseline", "Month 3", "Month 6"))
#' )
#' dl <- dl[order(dl$id, dl$visit), ]
#' trt <- factor(sample(c("Control", "Intervention"), n_subject, TRUE),
#' levels = c("Control", "Intervention"))
#' dl$treatment <- rep(trt, each = 3)
#' u <- rnorm(n_subject, 0, 6)
#' dl$sbp <- 140 + u[dl$id] - 3 * (dl$visit == "Month 3") -
#' 5 * (dl$visit == "Month 6") -
#' 4 * (dl$treatment == "Intervention" & dl$visit == "Month 6") +
#' rnorm(nrow(dl), 0, 5)
#'
#' mixed <- tablong(
#' dl, vars = vars(c.sbp), time = visit, id = id,
#' by = treatment, diagnostics = TRUE, show = FALSE
#' )
#' mixed$models[[1]]
#' mixed$diagnostics
#'
#' # ---------------------------------------------------------------------
#' # 7. Wide repeated data: internally converted to long format
#' # ---------------------------------------------------------------------
#' dw <- data.frame(
#' id = seq_len(n_subject), treatment = trt,
#' sbp0 = rnorm(n_subject, 140, 10)
#' )
#' dw$sbp3 <- dw$sbp0 - 3 + rnorm(n_subject, 0, 4)
#' dw$sbp6 <- dw$sbp0 - 5 - 4 * (dw$treatment == "Intervention") +
#' rnorm(n_subject, 0, 4)
#' wide <- tablong(
#' dw, vars = vars(c.sbp0, c.sbp3, c.sbp6),
#' time = c("Baseline", "Month 3", "Month 6"),
#' id = id, by = treatment, show = FALSE
#' )
#' wide$input_format
#' head(wide$long_data)
#'
#' # ---------------------------------------------------------------------
#' # 8. Binary repeated outcome: marginal logistic model with clustered SE and OR
#' # ---------------------------------------------------------------------
#' p <- plogis(-1 + 0.4 * (dl$visit == "Month 6") +
#' 0.6 * (dl$treatment == "Intervention"))
#' dl$controlled <- factor(rbinom(nrow(dl), 1, p),
#' levels = c(0, 1), labels = c("No", "Yes"))
#' binary_or <- tablong(
#' dl, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", show = FALSE
#' )
#' binary_or$contrasts_table
#'
#' # ---------------------------------------------------------------------
#' # 9. Continuous time and random slope
#' # ---------------------------------------------------------------------
#' ds <- expand.grid(id = seq_len(50), month = c(0, 3, 6, 12))
#' ds <- ds[order(ds$id, ds$month), ]
#' ds$group <- factor(rep(sample(c("Control", "Intervention"), 50, TRUE), each = 4))
#' b0 <- rnorm(50, 0, 5)
#' b1 <- rnorm(50, 0, 0.15)
#' ds$score <- 50 + b0[ds$id] + (-0.3 + b1[ds$id]) * ds$month -
#' 0.25 * ds$month * (ds$group == "Intervention") + rnorm(nrow(ds), 0, 3)
#' slope_fit <- tablong(
#' ds, vars = vars(c.score), time = c.month, id = id,
#' by = group, slope = TRUE, show = FALSE
#' )
#' slope_fit$contrasts_table
#'
#' # ---------------------------------------------------------------------
#' # 10. Repeated count outcome with person-time offset -> IRR
#' # ---------------------------------------------------------------------
#' dc <- dl
#' dc$person_time <- runif(nrow(dc), 0.8, 1.2)
#' rate <- exp(0.2 + 0.2 * (dc$visit == "Month 6") -
#' 0.3 * (dc$treatment == "Intervention" & dc$visit == "Month 6"))
#' dc$events <- rpois(nrow(dc), rate * dc$person_time)
#' count_fit <- tablong(
#' dc, vars = vars(c.events), time = visit, id = id,
#' by = treatment, count = TRUE, exposure = person_time,
#' show = FALSE
#' )
#' count_fit$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 11. RR/PR without extra packages; optional AR(1) GEE
#' # -----------------------------------------------------------------------
#' set.seed(13)
#' n_subject <- 70
#' dg <- expand.grid(
#' id = seq_len(n_subject),
#' visit = factor(c("Baseline", "Month 6"),
#' levels = c("Baseline", "Month 6"))
#' )
#' dg <- dg[order(dg$id, dg$visit), ]
#' dg$treatment <- factor(
#' rep(sample(c("Control", "Intervention"), n_subject, TRUE), each = 2),
#' levels = c("Control", "Intervention")
#' )
#' p <- plogis(-1 + 0.3 * (dg$visit == "Month 6") +
#' 0.4 * (dg$treatment == "Intervention"))
#' dg$controlled <- factor(rbinom(nrow(dg), 1, p),
#' levels = c(0, 1), labels = c("No", "Yes"))
#'
#' fit_rr <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", rr = TRUE, show = FALSE
#' )
#' fit_pr <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", pr = TRUE,
#' show = FALSE
#' )
#' fit_rr$contrasts_table
#' fit_pr$contrasts_table
#'
#' # AR(1) is advanced and uses geepack only when explicitly requested.
#' if (requireNamespace("geepack", quietly = TRUE)) {
#' fit_pr_ar1 <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", pr = TRUE, gee = TRUE, ar1 = TRUE,
#' show = FALSE
#' )
#' }
#'
#' # -----------------------------------------------------------------------
#' # 12. Three or more groups: keep main table compact, request pairwise tests
#' # -----------------------------------------------------------------------
#' set.seed(14)
#' dm <- data.frame(
#' period = factor(rep(c("Baseline", "Follow-up"), each = 90),
#' levels = c("Baseline", "Follow-up")),
#' arm = factor(rep(rep(c("A", "B", "C"), each = 30), 2))
#' )
#' dm$score <- rnorm(nrow(dm), 50 + 2 * (dm$period == "Follow-up") +
#' 2 * (dm$arm == "B") + 4 * (dm$arm == "C"), 7)
#' multi_arm <- tablong(
#' dm, vars = vars(c.score), time = period, by = arm,
#' pairwise = TRUE, adjust = "holm", show = FALSE
#' )
#' multi_arm$tests_table
#' multi_arm$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 13. Several outcomes in one long-data analysis
#' # -----------------------------------------------------------------------
#' d$positive <- factor(
#' rbinom(nrow(d), 1, plogis(-1 + 0.5 * (d$period == "After"))),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#' multi_outcome <- tablong(
#' d, vars = vars(c.score, positive), time = period,
#' by = group, event = "Yes", show = FALSE
#' )
#' multi_outcome$data
#' multi_outcome$descriptive
#'
#' # -----------------------------------------------------------------------
#' # 14. Export remains compatible with ordinary R4VN table workflows
#' # -----------------------------------------------------------------------
#' export_data <- tabexport(z)
#' head(export_data)
#'
#' # -----------------------------------------------------------------------
#' # 15. Missing outcomes/attrition: inspect counts before publication
#' # -----------------------------------------------------------------------
#' d_missing <- d
#' d_missing$score[c(2, 7, 21, 100)] <- NA
#' miss_fit <- tablong(
#' d_missing, vars = vars(c.score), time = period, by = group,
#' missing = TRUE, diagnostics = TRUE, show = FALSE
#' )
#' miss_fit$descriptive
#' miss_fit$diagnostics
#'
#' # -----------------------------------------------------------------------
#' # 16. Several binary outcomes can use a named event vector
#' # -----------------------------------------------------------------------
#' d$admitted <- factor(
#' rbinom(nrow(d), 1, plogis(-1.4 + 0.4 * (d$period == "After"))),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#' binary_set <- tablong(
#' d, vars = vars(positive, admitted), time = period, by = group,
#' event = c(positive = "Yes", admitted = "Yes"), show = FALSE
#' )
#' binary_set$tests_table
#'
#' # -----------------------------------------------------------------------
#' # 17. Continuous Gaussian GEE with AR(1) working correlation
#' # -----------------------------------------------------------------------
#' if (requireNamespace("geepack", quietly = TRUE)) {
#' set.seed(17)
#' dg2 <- expand.grid(id = seq_len(60), month = c(0, 3, 6, 12))
#' dg2 <- dg2[order(dg2$id, dg2$month), ]
#' dg2$group <- factor(rep(sample(c("Control", "Intervention"), 60, TRUE), each = 4))
#' dg2$score <- 55 - 0.2 * dg2$month -
#' 0.15 * dg2$month * (dg2$group == "Intervention") + rnorm(nrow(dg2), 0, 5)
#' gee_cont <- tablong(
#' dg2, vars = vars(c.score), time = c.month, id = id, by = group,
#' gee = TRUE, ar1 = TRUE, show = FALSE
#' )
#' gee_cont$contrasts_table
#' }
#'
#' # -----------------------------------------------------------------------
#' # 18. Wide count data can supply one person-time variable per time point
#' # -----------------------------------------------------------------------
#' set.seed(18)
#' nw <- 50
#' wc <- data.frame(
#' id = seq_len(nw),
#' group = factor(sample(c("Control", "Intervention"), nw, TRUE)),
#' pt0 = runif(nw, 0.8, 1.2),
#' pt6 = runif(nw, 0.8, 1.2)
#' )
#' wc$event0 <- rpois(nw, 1.2 * wc$pt0)
#' wc$event6 <- rpois(nw,
#' exp(log(1.2) - 0.25 * (wc$group == "Intervention")) * wc$pt6)
#' wide_count <- tablong(
#' wc, vars = vars(c.event0, c.event6),
#' time = c("Baseline", "Month 6"), id = id, by = group,
#' count = TRUE, exposure = vars(pt0, pt6), show = FALSE
#' )
#' wide_count$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 19. Replot an existing result without refitting the statistical model
#' # -----------------------------------------------------------------------
#' plot(z, ci = FALSE, title = "Observed longitudinal profile",
#' base_size = 12, legend_position = "right")
#' }
tablong <- function(data = NULL, vars,
time = NULL, id = NULL, by = NULL,
ref = NULL, event = NULL, adjusted = NULL,
gee = FALSE, ar1 = FALSE, slope = FALSE,
or = FALSE, rr = FALSE, pr = FALSE,
count = FALSE, exposure = NULL,
change = TRUE, pairwise = FALSE,
adjust = "none", missing = FALSE,
level = 0.95, digit = 1, p_digit = 3, effect_digit = 2,
bold_p = TRUE, p_bold = 0.05,
diagnostics = FALSE, diagnosis = FALSE, interpretation = FALSE,
plot = FALSE, plot_args = list(),
name = FALSE, title = NULL, file = NULL,
raw = FALSE, show = TRUE) {
caller <- parent.frame()
data <- .r4vn_resolve_analysis_data(data)
meta <- .r4vn_long_meta(vars, "vars", data = data)
flags <- c(or = isTRUE(or), rr = isTRUE(rr), pr = isTRUE(pr))
if (sum(flags) > 1L) stop("Use only one of `or = TRUE`, `rr = TRUE`, or `pr = TRUE`.", call. = FALSE)
if (!is.logical(count) || length(count) != 1L || is.na(count)) stop("`count` must be TRUE or FALSE.", call. = FALSE)
if (isTRUE(count) && any(flags)) {
stop("`or`, `rr`, and `pr` are binary-outcome options and cannot be combined with `count = TRUE`.", call. = FALSE)
}
if (!is.logical(gee) || length(gee) != 1L || is.na(gee)) stop("`gee` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(ar1) || length(ar1) != 1L || is.na(ar1)) stop("`ar1` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(slope) || length(slope) != 1L || is.na(slope)) stop("`slope` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(change) || length(change) != 1L || is.na(change)) stop("`change` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(pairwise) || length(pairwise) != 1L || is.na(pairwise)) stop("`pairwise` must be TRUE or FALSE.", call. = FALSE)
if (!missing(diagnosis)) {
if (!is.logical(diagnosis) || length(diagnosis) != 1L || is.na(diagnosis)) stop("`diagnosis` must be TRUE or FALSE.", call. = FALSE)
diagnostics <- diagnosis
}
if (!is.logical(diagnostics) || length(diagnostics) != 1L || is.na(diagnostics)) stop("`diagnostics` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(interpretation) || length(interpretation) != 1L || is.na(interpretation)) stop("`interpretation` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(plot) || length(plot) != 1L || is.na(plot)) stop("`plot` must be TRUE or FALSE.", call. = FALSE)
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
if (length(plot_args) && is.null(names(plot_args))) stop("`plot_args` must be a named list.", call. = FALSE)
if (!is.numeric(level) || length(level) != 1L || !is.finite(level) || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
if (!adjust %in% stats::p.adjust.methods) stop("`adjust` must be one of: ", paste(stats::p.adjust.methods, collapse = ", "), ".", call. = FALSE)
adjusted_meta <- .r4vn_long_adjusted_meta(
substitute(adjusted),
missing(adjusted),
caller,
data = data
)
time_source <- .r4vn_long_time_source(
substitute(time),
data,
nrow(meta),
caller
)
exposure_expr <- substitute(exposure)
if (identical(time_source$mode, "wide")) {
prep <- .r4vn_long_prepare_wide(
data = data,
meta = meta,
time_labels = time_source$labels,
id_expr = substitute(id),
id_missing = missing(id),
by_expr = substitute(by),
by_missing = missing(by),
adjusted_meta = adjusted_meta,
ref = ref,
exposure_expr = exposure_expr,
exposure_missing = missing(exposure),
env = caller
)
input_format <- "wide"
} else {
prep <- .r4vn_long_prepare_long_common(
data = data,
time_spec = time_source$spec,
id_expr = substitute(id),
id_missing = missing(id),
by_expr = substitute(by),
by_missing = missing(by),
adjusted_meta = adjusted_meta,
ref = ref,
env = caller
)
input_format <- "long"
if (!missing(exposure) && !identical(exposure_expr, quote(NULL))) {
if (!is.symbol(exposure_expr)) stop("In long data, `exposure` must be one unquoted variable.", call. = FALSE)
exposure_name <- as.character(exposure_expr)
if (!exposure_name %in% names(prep$data)) stop("`exposure` variable was not found.", call. = FALSE)
prep$data$.exposure <- prep$data[[exposure_name]]
}
}
if (isTRUE(slope) && !isTRUE(prep$time_continuous)) {
stop("`slope = TRUE` is available only with continuous time, for example `time = c.month`.", call. = FALSE)
}
if (isTRUE(ar1) && !isTRUE(gee) && !isTRUE(rr) && !isTRUE(pr)) {
warning("`ar1 = TRUE` is used only with GEE. Set `gee = TRUE`, `rr = TRUE`, or `pr = TRUE`.", call. = FALSE)
}
design <- if (isTRUE(prep$repeated)) "longitudinal" else "repeated cross-sectional"
design_note <- if (identical(input_format, "wide")) {
"Wide data were converted to long format internally; the original data were not modified."
} else if (identical(design, "repeated cross-sectional")) {
"No subject ID was observed at more than one time point; observations across periods were treated as independent repeated cross-sectional samples."
} else {
"Repeated observations from the same subject were modeled using within-subject correlation."
}
if (identical(input_format, "long") && !isTRUE(prep$repeated) && is.null(prep$id_name)) {
design_note <- "No `id` was supplied; observations across time were treated as independent repeated cross-sectional samples."
}
results <- list()
all_rows <- list()
all_models <- list()
all_tests <- list()
all_contrasts <- list()
all_descriptive <- list()
all_diagnostics <- list()
all_plot_data <- list()
long_storage <- list()
notes <- character()
if (identical(input_format, "wide")) {
outcome_iterations <- 1L
} else {
outcome_iterations <- seq_len(nrow(meta))
}
for (ii in outcome_iterations) {
if (identical(input_format, "wide")) {
od <- prep$data
outcome_name <- prep$outcome_name
summary_type <- meta$type[1L]
original_variable <- paste(meta$variable, collapse = ", ")
} else {
outcome_name <- meta$variable[ii]
if (!outcome_name %in% names(prep$data)) stop("Outcome `", outcome_name, "` was not found in `data`.", call. = FALSE)
od <- prep$data
summary_type <- meta$type[ii]
original_variable <- outcome_name
if (summary_type %in% c("mean", "median", "full")) {
od$.outcome <- .r4vn_long_numeric(od[[outcome_name]], outcome_name)
} else {
od$.outcome <- od[[outcome_name]]
}
}
if (identical(input_format, "wide")) {
outcome_label <- prep$outcome_name
} else {
outcome_label <- .r4vn_long_label(data[[outcome_name]], outcome_name)
}
if (isTRUE(name) && !identical(outcome_label, original_variable)) {
outcome_label <- paste0(outcome_label, " (", original_variable, ")")
}
# Preserve the observed outcome exactly as supplied for descriptive
# summaries. The model copy below may be recoded (for example Yes/No -> 0/1).
od$.outcome_display <- od$.outcome
outcome_type <- NULL
event_level <- NULL
if (isTRUE(count)) {
numeric_outcome <- suppressWarnings(as.numeric(od$.outcome))
nonmissing <- numeric_outcome[is.finite(numeric_outcome)]
if (length(nonmissing) && (any(nonmissing < 0) || any(abs(nonmissing - round(nonmissing)) > 1e-8))) {
stop("`count = TRUE` requires non-negative integer outcomes.", call. = FALSE)
}
od$.outcome <- numeric_outcome
outcome_type <- "count"
effect_type <- "IRR"
if (".exposure" %in% names(od)) {
od$.exposure <- suppressWarnings(as.numeric(od$.exposure))
if (any(!is.na(od$.exposure) & (!is.finite(od$.exposure) | od$.exposure <= 0))) {
stop("`exposure` must be positive wherever it is observed.", call. = FALSE)
}
}
} else if (summary_type %in% c("mean", "median", "full")) {
od$.outcome <- suppressWarnings(as.numeric(od$.outcome))
outcome_type <- "continuous"
effect_type <- "Difference"
} else {
observed <- .r4vn_long_levels(od$.outcome)
if (length(observed) != 2L) {
stop(
"Outcome `", outcome_name, "` is categorical with ", length(observed),
" observed levels. `tablong()` currently supports continuous, binary, and count outcomes.",
call. = FALSE
)
}
event_level <- .r4vn_long_event(od$.outcome, event, outcome_name)
od$.outcome <- as.integer(as.character(od$.outcome) == as.character(event_level))
attr(od$.outcome, "event") <- event_level
outcome_type <- "binary"
effect_type <- if (isTRUE(rr)) "RR" else if (isTRUE(pr)) "PR" else "OR"
}
model_vars <- c(".outcome", ".time", ".id", prep$covariates)
if (!is.null(prep$by_name)) model_vars <- c(model_vars, ".by")
if (".exposure" %in% names(od)) model_vars <- c(model_vars, ".exposure")
model_vars <- unique(model_vars)
complete_base <- stats::complete.cases(od[, model_vars, drop = FALSE])
model_data <- od[complete_base, , drop = FALSE]
if (nrow(model_data) < 5L) stop("Too few complete observations for outcome `", outcome_name, "`.", call. = FALSE)
if (identical(outcome_type, "binary") && length(unique(model_data$.outcome)) < 2L) {
stop("Binary outcome `", outcome_name, "` has no variation after removing missing model data.", call. = FALSE)
}
fit <- .r4vn_long_fit_one(
data = model_data,
outcome_type = outcome_type,
effect_type = effect_type,
repeated = prep$repeated,
gee = gee,
ar1 = ar1,
slope = slope,
time_continuous = prep$time_continuous,
covariates = prep$covariates,
exposure = ".exposure" %in% names(model_data)
)
if (isTRUE(prep$time_continuous)) {
display_times <- prep$time_display
} else {
display_times <- prep$time_display
}
contrasts <- .r4vn_long_contrast_rows(
data = model_data,
fit = fit,
time_display = display_times,
time_reference = prep$time_reference,
by_display = prep$by_display,
time_continuous = prep$time_continuous,
change = change,
pairwise = pairwise,
level = level,
adjust = adjust
)
tests <- list(
p_time = fit$p_time,
p_group = fit$p_group,
p_interaction = fit$p_interaction
)
rows <- .r4vn_long_build_rows(
data = od,
outcome_label = outcome_label,
outcome_type = outcome_type,
summary_type = summary_type,
event = event_level,
effect_type = effect_type,
time_display = display_times,
time_reference = prep$time_reference,
by_display = prep$by_display,
time_continuous = prep$time_continuous,
contrasts = contrasts,
tests = tests,
change = change,
digits = digit,
effect_digits = effect_digit,
p_digits = p_digit,
show_n = missing
)
descriptive_rows <- .r4vn_long_descriptive_rows(
data = od, outcome_label = outcome_label, outcome_type = outcome_type,
summary_type = summary_type, event = event_level,
time_display = display_times, by_display = prep$by_display,
time_continuous = prep$time_continuous, digits = digit, show_n = missing
)
diagnostic_row <- .r4vn_long_diagnostics_one(
fit = fit, model_data = model_data, outcome_label = outcome_label,
outcome_type = outcome_type, time_display = display_times,
time_continuous = prep$time_continuous, repeated = prep$repeated
)
plot_rows <- .r4vn_long_plot_rows(
data = od, outcome_label = outcome_label, outcome_type = outcome_type,
event = event_level, time_display = display_times,
by_display = prep$by_display, time_continuous = prep$time_continuous,
level = level
)
if (isTRUE(fit$singular)) {
notes <- c(notes, paste0(outcome_label, ": the mixed model produced a singular random-effects fit; inspect `$models` before publication."))
}
if (!is.null(fit$package_note) && length(fit$package_note) && nzchar(fit$package_note[1L])) {
notes <- c(notes, paste0(outcome_label, ": ", fit$package_note[1L]))
}
if (summary_type %in% c("median", "full") && identical(outcome_type, "continuous")) {
notes <- c(
notes,
paste0(
outcome_label,
": `", if (identical(summary_type, "median")) "q." else "f.",
"` changes the observed descriptive summary only; inferential effects are mean-model effects."
)
)
}
key <- make.unique(c(names(all_models), outcome_name))[length(all_models) + 1L]
all_rows[[key]] <- rows
all_models[[key]] <- fit$fit
all_tests[[key]] <- tests
all_contrasts[[key]] <- contrasts
all_descriptive[[key]] <- descriptive_rows
all_diagnostics[[key]] <- diagnostic_row
all_plot_data[[key]] <- plot_rows
long_storage[[key]] <- od
results[[key]] <- list(
outcome = outcome_name,
label = outcome_label,
type = outcome_type,
event = event_level,
effect = effect_type,
model_engine = fit$engine,
model = fit$fit,
tests = tests,
contrasts = contrasts,
data = rows
)
}
# Multiple outcomes may use different effect scales (for example a mean
# difference for SBP and an OR for a binary outcome). Standardize only the
# effect-column heading when needed so all outcome blocks can share one table.
effect_columns <- unique(unlist(lapply(all_rows, function(z) {
grep("\\(95% CI\\)$", names(z), value = TRUE)
}), use.names = FALSE))
if (length(effect_columns) > 1L) {
all_rows <- lapply(all_rows, function(z) {
hit <- grep("\\(95% CI\\)$", names(z), value = TRUE)
if (length(hit) == 1L) names(z)[names(z) == hit] <- "Model effect (95% CI)"
z
})
scale_notes <- vapply(results, function(z) {
scale <- if (identical(z$type, "continuous")) "mean difference" else z$effect
paste0(z$label, ": model effect is ", scale, ".")
}, character(1))
notes <- c(notes, scale_notes)
}
flat <- do.call(rbind, all_rows)
rownames(flat) <- NULL
descriptive_table <- if (length(all_descriptive)) do.call(rbind, all_descriptive) else data.frame()
if (nrow(descriptive_table)) rownames(descriptive_table) <- NULL
diagnostics_table <- if (length(all_diagnostics)) do.call(rbind, all_diagnostics) else data.frame()
if (nrow(diagnostics_table)) rownames(diagnostics_table) <- NULL
plot_data <- if (length(all_plot_data)) do.call(rbind, all_plot_data) else data.frame()
if (nrow(plot_data)) rownames(plot_data) <- NULL
tests_table <- .r4vn_long_tests_table(results, digits = p_digit)
contrasts_table <- .r4vn_long_contrasts_table(
results, effect_digits = effect_digit, p_digits = p_digit
)
interpretation_table <- if (isTRUE(interpretation)) {
.r4vn_long_interpretation(results, alpha = 0.05, p_digits = p_digit)
} else NULL
graph <- if (isTRUE(plot)) {
.r4vn_long_build_plot(plot_data, ci = TRUE, plot_args = plot_args)
} else NULL
base_notes <- c(
design_note,
if (identical(design, "longitudinal") && any(vapply(results, function(z) identical(z$model_engine, "mixed"), logical(1)))) {
"Continuous repeated outcomes used a random-intercept linear mixed model from R's recommended `nlme` package; `slope = TRUE` adds a random linear time slope when time is continuous."
} else NULL,
if (any(vapply(results, function(z) identical(z$model_engine, "cluster_robust"), logical(1)))) {
"Repeated binary/count outcomes and marginal models used working-independence regression with subject-clustered robust sandwich standard errors calculated internally by R4VN."
} else NULL,
if (any(vapply(results, function(z) identical(z$model_engine, "gee"), logical(1)))) {
"AR(1) GEE used robust sandwich standard errors through the optional `geepack` package."
} else NULL,
if (any(vapply(results, function(z) z$effect %in% c("RR", "PR") && identical(z$model_engine, "independent"), logical(1)))) {
"RR/PR in independent samples used modified Poisson regression with robust sandwich standard errors calculated internally by R4VN."
} else NULL,
if (!is.null(prep$by_name)) {
"Time and group p-values are omnibus main-effect tests from the additive model; the time x group p-value tests the interaction in the full model."
} else {
"The overall time p-value tests whether the time effect is zero."
},
if (isTRUE(change) && !isTRUE(prep$time_continuous)) {
paste0("Change contrasts use ", prep$time_reference, " as the reference time.")
} else NULL,
if (!identical(adjust, "none")) {
paste0("Displayed contrast p-values use ", adjust, " multiplicity adjustment; confidence intervals are unadjusted.")
} else NULL
)
notes <- unique(c(base_notes, notes))
notes <- notes[nzchar(notes)]
if (is.null(title) || !length(title) || is.na(title[1L]) || !nzchar(as.character(title[1L]))) {
title <- "Longitudinal / repeated-measures analysis"
} else {
title <- as.character(title[1L])
}
table_html <- .r4vn_long_html_table(flat, title, notes, bold_p, p_bold)
plot_html <- if (isTRUE(plot) && !is.null(graph)) .r4vn_long_plot_html(graph) else ""
secondary_html <- paste0(
plot_html,
.r4vn_long_html_simple_table(tests_table, "Omnibus tests"),
if (isTRUE(diagnostics)) .r4vn_long_html_simple_table(diagnostics_table, "Model diagnostics") else "",
if (isTRUE(interpretation)) .r4vn_long_html_simple_table(interpretation_table, "Interpretation") else ""
)
html <- .r4vn_long_html_document(paste0(table_html, secondary_html))
if (is.null(file) || !length(file) || is.na(file[1L]) || !nzchar(as.character(file[1L]))) {
file <- tempfile(pattern = "r4vn-tablong-", fileext = ".html")
} else {
file <- path.expand(as.character(file[1L]))
if (!grepl("\\.html?$", file, ignore.case = TRUE)) file <- paste0(file, ".html")
dir.create(dirname(file), recursive = TRUE, showWarnings = FALSE)
}
writeLines(enc2utf8(html), file, useBytes = TRUE)
file <- normalizePath(file, winslash = "/", mustWork = TRUE)
if (isTRUE(show)) .r4vn_long_show_html(file)
if (isTRUE(plot) && !is.null(graph) && interactive()) {
.r4vn_long_draw_plot(graph$data, ci = graph$ci, plot_args = graph$args)
}
combined_long <- if (length(long_storage) == 1L) {
long_storage[[1L]]
} else {
do.call(rbind, lapply(names(long_storage), function(nm) {
z <- long_storage[[nm]]
z$.r4vn_outcome_name <- nm
z
}))
}
rownames(combined_long) <- NULL
tables <- list(
`Main table` = flat,
Descriptive = descriptive_table,
`Omnibus tests` = tests_table,
Contrasts = contrasts_table,
Diagnostics = diagnostics_table
)
if (isTRUE(interpretation)) tables$Interpretation <- interpretation_table
result <- list(
data = flat,
descriptive = descriptive_table,
estimates = list(contrasts = contrasts_table),
table_html = table_html,
html = html,
file = file,
models = all_models,
tests = all_tests,
tests_table = tests_table,
contrasts = all_contrasts,
contrasts_table = contrasts_table,
diagnostics = diagnostics_table,
interpretation = interpretation_table,
tables = tables,
plot_data = plot_data,
graph = graph,
plots = list(trajectory = graph),
long_data = combined_long,
results = results,
input_format = input_format,
design = design,
time = list(
variable = prep$time_name,
continuous = prep$time_continuous,
levels = prep$time_display,
reference = prep$time_reference
),
by = prep$by_name,
id = prep$id_name,
adjusted = adjusted_meta,
metadata = list(
input_format = input_format, design = design,
repeated = isTRUE(prep$repeated), time_variable = prep$time_name,
time_continuous = isTRUE(prep$time_continuous),
time_levels = prep$time_display, time_reference = prep$time_reference,
group_variable = prep$by_name, id_variable = prep$id_name,
change = isTRUE(change), pairwise = isTRUE(pairwise), adjust = adjust,
confidence_level = level, interpretation = isTRUE(interpretation),
plot = isTRUE(plot)
),
call = match.call()
)
class(result) <- c("r4vn_tablong", "r4vn_tab", "list")
invisible(result)
}
#' Print an R4VN Longitudinal Table
#'
#' @param x Object returned by `tablong()`.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @keywords internal
#' @method print r4vn_tablong
#' @export
print.r4vn_tablong <- function(x, ...) {
print(x$data, row.names = FALSE)
invisible(x)
}
#' Plot an R4VN Longitudinal Analysis
#'
#' Recreates the observed longitudinal profile stored by `tablong()`. This is
#' useful when the original analysis used `plot = FALSE` or when different plot
#' labels/sizing are wanted without refitting the statistical model.
#'
#' @param x Object returned by `tablong()`.
#' @param ci Show 95% confidence intervals. Default `TRUE`.
#' @param ... Named plot options accepted through `plot_args`, including
#' `title`, `xlab`, `ylab`, `line_width`, `point_size`, `base_size`, and
#' `legend_position`.
#' @return An R4VN plot specification, invisibly; the plot is drawn with base R graphics.
#' @method plot r4vn_tablong
#' @export
plot.r4vn_tablong <- function(x, ci = TRUE, ...) {
args <- list(...)
g <- .r4vn_long_build_plot(x$plot_data, ci = ci, plot_args = args)
if (!is.null(g)) .r4vn_long_draw_plot(g$data, ci = g$ci, plot_args = g$args)
invisible(g)
}
attr(tablong, "r4vn_version") <- "tablong-1.2.0-2026-09-02"
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R4VN documentation built on Sept. 30, 2026, 5:13 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions .r4vn_long_contrasts_table .r4vn_long_tests_table .r4vn_long_descriptive_rows .r4vn_long_show_html .r4vn_long_html_document .r4vn_long_html_table .r4vn_long_build_rows .r4vn_long_change_name .r4vn_long_effect_name .r4vn_long_observed_cell .r4vn_long_summary_count .r4vn_long_summary_binary .r4vn_long_summary_continuous .r4vn_long_contrast_rows .r4vn_long_contrast .r4vn_long_grid_row .r4vn_long_fit_one .r4vn_long_extract_beta_vcov .r4vn_long_term_columns .r4vn_long_wald .r4vn_long_fixed_rank .r4vn_long_safe_scalar .r4vn_long_glm_robust_vcov .r4vn_long_prepare_wide .r4vn_long_prepare_long_common .r4vn_long_time_source .r4vn_long_common_wide_label .r4vn_long_event .r4vn_long_apply_covariates .r4vn_long_numeric .r4vn_long_factor .r4vn_long_adjusted_meta .r4vn_long_meta .r4vn_long_parse_symbol .r4vn_long_parse_spec_text .r4vn_long_fmt_ci .r4vn_long_fmt_p .r4vn_long_fmt .r4vn_long_label .r4vn_long_levels .r4vn_long_escape
# ============================================================================
# R4VN::tablong()
# Longitudinal / repeated-measures publication table
#
# Design principles
# - Same culture as tab(): data = NULL uses active data; vars() declares outcomes.
# - Accepts long and wide data; wide data are reshaped internally only.
# - Continuous repeated outcomes: random-intercept model via recommended nlme;
# automatic base-R cluster-robust fallback if nlme is unavailable.
# - Binary/count repeated outcomes: marginal regression with subject-clustered
# robust sandwich variance calculated internally by R4VN.
# - rr = TRUE / pr = TRUE: modified Poisson with robust variance; geepack is
# optional only for explicitly requested AR(1) GEE.
# - Repeated cross-sectional data are detected automatically when IDs do not
# repeat across time (or no ID is supplied).
# - Returns c("r4vn_tablong", "r4vn_tab") so tabexport() works unchanged.
# ============================================================================
.r4vn_long_escape <- function(x) {
x <- as.character(x)
x[is.na(x)] <- ""
x <- gsub("&", "&", x, fixed = TRUE)
x <- gsub("<", "<", x, fixed = TRUE)
x <- gsub(">", ">", x, fixed = TRUE)
x <- gsub('"', """, x, fixed = TRUE)
gsub("'", "'", x, fixed = TRUE)
}
.r4vn_long_levels <- function(x) {
z <- x[!is.na(x)]
if (!length(z)) return(character())
if (is.factor(x)) {
lev <- levels(x)
return(lev[lev %in% as.character(z)])
}
if (is.logical(x)) return(as.character(c(FALSE, TRUE)[c(FALSE, TRUE) %in% z]))
if (is.numeric(x)) return(as.character(sort(unique(z))))
unique(as.character(z))
}
.r4vn_long_label <- function(x, fallback) {
lab <- attr(x, "label", exact = TRUE)
if (is.null(lab) || !length(lab) || is.na(lab[1L]) || !nzchar(as.character(lab[1L]))) {
return(fallback)
}
as.character(lab[1L])
}
.r4vn_long_fmt <- function(x, digits = 2) {
if (!length(x) || is.na(x) || !is.finite(x)) return("")
formatC(x, format = "f", digits = digits)
}
.r4vn_long_fmt_p <- function(p, digits = 3) {
if (!length(p) || is.na(p) || !is.finite(p)) return("")
limit <- 10^(-digits)
if (p < limit) return(paste0("<", formatC(limit, format = "f", digits = digits)))
formatC(p, format = "f", digits = digits)
}
.r4vn_long_fmt_ci <- function(est, low, high, digits = 2) {
if (any(!is.finite(c(est, low, high)))) return("")
paste0(
formatC(est, format = "f", digits = digits),
" (",
formatC(low, format = "f", digits = digits),
", ",
formatC(high, format = "f", digits = digits),
")"
)
}
.r4vn_long_parse_spec_text <- function(text, allow_continuous = TRUE) {
out <- list(
variable = text,
type = "categorical",
reference_index = 1L,
continuous = FALSE,
specification = text
)
if (grepl("^b[1-9][0-9]*\\.", text)) {
out$reference_index <- as.integer(sub("^b([1-9][0-9]*)\\..*$", "\\1", text))
out$variable <- sub("^b[1-9][0-9]*\\.", "", text)
} else if (startsWith(text, "c.")) {
if (!allow_continuous) {
stop("The `c.` prefix is not allowed here.", call. = FALSE)
}
out$type <- "mean"
out$continuous <- TRUE
out$reference_index <- NA_integer_
out$variable <- sub("^c\\.", "", text)
} else if (startsWith(text, "q.")) {
out$type <- "median"
out$reference_index <- NA_integer_
out$variable <- sub("^q\\.", "", text)
} else if (startsWith(text, "f.")) {
out$type <- "full"
out$reference_index <- NA_integer_
out$variable <- sub("^f\\.", "", text)
}
if (!nzchar(out$variable)) stop("A prefix must be followed by a variable name.", call. = FALSE)
out
}
.r4vn_long_parse_symbol <- function(expr, role, allow_continuous = TRUE) {
if (!is.symbol(expr)) {
stop("`", role, "` must be an unquoted variable name in long data.", call. = FALSE)
}
.r4vn_long_parse_spec_text(as.character(expr), allow_continuous = allow_continuous)
}
.r4vn_long_meta <- function(x, argument = "vars", data = NULL) {
if (!inherits(x, "r4vn_vars")) {
stop("`", argument, "` must be created using `vars()`.", call. = FALSE)
}
required <- c("variable", "type", "reference_index")
if (!all(required %in% names(x))) {
stop("`", argument, "` is not a valid current R4VN `vars()` object.", call. = FALSE)
}
if (is.data.frame(data)) {
resolver <- get0(".r4vn_resolve_vars", mode = "function", inherits = TRUE)
if (is.function(resolver)) {
x <- resolver(x, data = data, default_type = "auto", strict = TRUE)
}
}
x
}
.r4vn_long_adjusted_meta <- function(expr, missing_arg, env, data = NULL) {
if (isTRUE(missing_arg) || identical(expr, quote(NULL))) return(NULL)
value <- eval(expr, envir = env)
.r4vn_long_meta(value, "adjusted", data = data)
}
.r4vn_long_factor <- function(x, reference_index = 1L, reference_value = NULL) {
display_levels <- .r4vn_long_levels(x)
if (length(display_levels) < 1L) {
return(list(
x = factor(x),
display_levels = character(),
model_levels = character(),
reference = NA_character_
))
}
reference <- NULL
if (!is.null(reference_value) && length(reference_value) && !is.na(reference_value[1L])) {
candidate <- as.character(reference_value[1L])
hit <- match(candidate, display_levels)
if (is.na(hit)) {
stop("Reference level `", candidate, "` was not found.", call. = FALSE)
}
reference <- display_levels[hit]
} else {
if (is.na(reference_index) || reference_index < 1L || reference_index > length(display_levels)) {
reference_index <- 1L
}
reference <- display_levels[reference_index]
}
model_levels <- c(reference, setdiff(display_levels, reference))
list(
x = factor(as.character(x), levels = model_levels),
display_levels = display_levels,
model_levels = model_levels,
reference = reference
)
}
.r4vn_long_numeric <- function(x, variable) {
y <- suppressWarnings(as.numeric(x))
if (all(is.na(y)) && any(!is.na(x))) {
stop("Variable `", variable, "` cannot be converted to numeric.", call. = FALSE)
}
y
}
.r4vn_long_apply_covariates <- function(data, meta) {
if (is.null(meta) || !nrow(meta)) {
return(list(data = data, names = character(), map = data.frame()))
}
names_internal <- character(nrow(meta))
map <- meta
for (i in seq_len(nrow(meta))) {
variable <- meta$variable[i]
if (!variable %in% names(data)) {
stop("Adjusted variable `", variable, "` was not found in `data`.", call. = FALSE)
}
internal <- paste0(".z", i)
names_internal[i] <- internal
if (identical(meta$type[i], "categorical")) {
f <- .r4vn_long_factor(
data[[variable]],
reference_index = meta$reference_index[i]
)
data[[internal]] <- f$x
} else {
data[[internal]] <- .r4vn_long_numeric(data[[variable]], variable)
}
}
map$internal <- names_internal
list(data = data, names = names_internal, map = map)
}
.r4vn_long_event <- function(x, event, outcome_name) {
observed <- .r4vn_long_levels(x)
if (length(observed) != 2L) {
stop(
"Binary outcome `", outcome_name, "` must have exactly two observed levels; found ",
length(observed), ".",
call. = FALSE
)
}
chosen <- NULL
if (!is.null(event)) {
if (length(event) > 1L && !is.null(names(event)) && outcome_name %in% names(event)) {
chosen <- as.character(event[[outcome_name]])
} else if (length(event) == 1L) {
chosen <- as.character(event[1L])
}
}
if (is.null(chosen) || is.na(chosen) || !nzchar(chosen)) chosen <- observed[length(observed)]
if (!chosen %in% observed) {
stop("Event `", chosen, "` was not found in outcome `", outcome_name, "`.", call. = FALSE)
}
chosen
}
.r4vn_long_common_wide_label <- function(data, variables) {
labs <- vapply(variables, function(v) .r4vn_long_label(data[[v]], ""), character(1))
labs <- unique(labs[nzchar(labs)])
if (length(labs) == 1L) return(labs)
stripped <- sub("([_.]?(baseline|base|pre|before|post|after|month|m|visit|t)?[_.]?[0-9]+)$",
"", variables, ignore.case = TRUE)
stripped <- sub("([_.](baseline|base|pre|before|post|after))$", "",
stripped, ignore.case = TRUE)
stripped <- unique(stripped[nzchar(stripped)])
if (length(stripped) == 1L) return(stripped)
variables[1L]
}
.r4vn_long_time_source <- function(expr, data, meta_n, env) {
if (identical(expr, quote(NULL))) {
if (meta_n > 1L) return(list(mode = "wide", labels = NULL, spec = NULL))
stop("`time` is required for long data.", call. = FALSE)
}
if (is.symbol(expr)) {
spec <- .r4vn_long_parse_spec_text(as.character(expr), allow_continuous = TRUE)
if (spec$variable %in% names(data)) {
return(list(mode = "long", labels = NULL, spec = spec))
}
value <- tryCatch(eval(expr, envir = env), error = function(e) NULL)
if (!is.null(value) && length(value) == meta_n && meta_n > 1L) {
return(list(mode = "wide", labels = as.character(value), spec = NULL))
}
stop(
"`time` variable `", spec$variable,
"` was not found in `data`. For wide data, use a vector such as ",
'`time = c("Baseline", "Month 3", "Month 6")`.',
call. = FALSE
)
}
value <- eval(expr, envir = env)
if (meta_n <= 1L) {
stop("A vector of time labels is only used with wide data containing multiple repeated outcome variables.", call. = FALSE)
}
if (length(value) != meta_n) {
stop("The number of `time` labels must equal the number of repeated variables in `vars()`.", call. = FALSE)
}
list(mode = "wide", labels = as.character(value), spec = NULL)
}
.r4vn_long_prepare_long_common <- function(data, time_spec, id_expr, id_missing,
by_expr, by_missing, adjusted_meta,
ref, env) {
time_name <- time_spec$variable
if (!time_name %in% names(data)) stop("`time` variable was not found.", call. = FALSE)
out <- data
if (isTRUE(time_spec$continuous)) {
out$.time <- .r4vn_long_numeric(out[[time_name]], time_name)
time_display <- sort(unique(out$.time[is.finite(out$.time)]))
time_reference <- if (length(time_display)) min(time_display) else NA_real_
time_continuous <- TRUE
} else {
tf <- .r4vn_long_factor(
out[[time_name]],
reference_index = time_spec$reference_index,
reference_value = ref
)
out$.time <- tf$x
time_display <- tf$display_levels
time_reference <- tf$reference
time_continuous <- FALSE
}
id_name <- NULL
if (isTRUE(id_missing) || identical(id_expr, quote(NULL))) {
out$.id <- seq_len(nrow(out))
} else {
id_spec <- .r4vn_long_parse_symbol(id_expr, "id", allow_continuous = FALSE)
id_name <- id_spec$variable
if (!id_name %in% names(out)) stop("`id` variable `", id_name, "` was not found.", call. = FALSE)
out$.id <- out[[id_name]]
}
by_name <- NULL
by_display <- NULL
by_reference <- NULL
if (!isTRUE(by_missing) && !identical(by_expr, quote(NULL))) {
by_spec <- .r4vn_long_parse_symbol(by_expr, "by", allow_continuous = FALSE)
by_name <- by_spec$variable
if (!by_name %in% names(out)) stop("`by` variable `", by_name, "` was not found.", call. = FALSE)
bf <- .r4vn_long_factor(out[[by_name]], reference_index = by_spec$reference_index)
if (length(bf$display_levels) < 2L) stop("`by` must contain at least two observed groups.", call. = FALSE)
out$.by <- bf$x
by_display <- bf$display_levels
by_reference <- bf$reference
}
cov <- .r4vn_long_apply_covariates(out, adjusted_meta)
out <- cov$data
repeated <- FALSE
if (!is.null(id_name)) {
valid <- !is.na(out$.id) & !is.na(out$.time)
if (any(valid)) {
key <- split(as.character(out$.time[valid]), as.character(out$.id[valid]))
repeated <- any(vapply(key, function(z) length(unique(z)) > 1L, logical(1)))
}
}
list(
data = out,
time_name = time_name,
time_display = time_display,
time_reference = time_reference,
time_continuous = time_continuous,
id_name = id_name,
repeated = repeated,
by_name = by_name,
by_display = by_display,
by_reference = by_reference,
covariates = cov$names,
covariate_map = cov$map
)
}
.r4vn_long_prepare_wide <- function(data, meta, time_labels, id_expr, id_missing,
by_expr, by_missing, adjusted_meta,
ref, exposure_expr, exposure_missing, env) {
variables <- meta$variable
missing_vars <- setdiff(variables, names(data))
if (length(missing_vars)) {
stop("Variables not found in `data`: ", paste(missing_vars, collapse = ", "), ".", call. = FALSE)
}
types <- unique(meta$type)
if (length(types) > 1L) {
continuous_types <- c("mean", "median", "full")
if (!all(types %in% continuous_types)) {
stop("All repeated variables in wide data must describe the same outcome type.", call. = FALSE)
}
}
if (is.null(time_labels)) {
time_labels <- variables
}
if (anyDuplicated(time_labels)) stop("Time labels must be unique.", call. = FALSE)
id_name <- NULL
base_id <- seq_len(nrow(data))
if (!isTRUE(id_missing) && !identical(id_expr, quote(NULL))) {
id_spec <- .r4vn_long_parse_symbol(id_expr, "id", allow_continuous = FALSE)
id_name <- id_spec$variable
if (!id_name %in% names(data)) stop("`id` variable `", id_name, "` was not found.", call. = FALSE)
base_id <- data[[id_name]]
}
by_name <- NULL
by_spec <- NULL
if (!isTRUE(by_missing) && !identical(by_expr, quote(NULL))) {
by_spec <- .r4vn_long_parse_symbol(by_expr, "by", allow_continuous = FALSE)
by_name <- by_spec$variable
if (!by_name %in% names(data)) stop("`by` variable `", by_name, "` was not found.", call. = FALSE)
}
exposure_mode <- "none"
exposure_vars <- NULL
exposure_name <- NULL
if (!isTRUE(exposure_missing) && !identical(exposure_expr, quote(NULL))) {
exposure_value <- tryCatch(eval(exposure_expr, envir = env), error = function(e) NULL)
if (inherits(exposure_value, "r4vn_vars")) {
exposure_vars <- exposure_value$variable
if (length(exposure_vars) != length(variables)) {
stop("Wide `exposure = vars(...)` must contain one exposure variable per repeated outcome variable.", call. = FALSE)
}
if (any(!exposure_vars %in% names(data))) stop("Some exposure variables were not found in `data`.", call. = FALSE)
exposure_mode <- "wide"
} else if (is.symbol(exposure_expr)) {
exposure_name <- as.character(exposure_expr)
if (!exposure_name %in% names(data)) stop("`exposure` variable was not found.", call. = FALSE)
exposure_mode <- "single"
} else {
stop("In wide data, `exposure` must be one variable or `vars(...)` with one variable per time point.", call. = FALSE)
}
}
pieces <- vector("list", length(variables))
for (j in seq_along(variables)) {
piece <- data
piece$.id <- base_id
piece$.time_source <- time_labels[j]
if (all(meta$type %in% c("mean", "median", "full"))) {
piece$.outcome <- .r4vn_long_numeric(data[[variables[j]]], variables[j])
} else {
piece$.outcome <- as.character(data[[variables[j]]])
}
if (identical(exposure_mode, "wide")) piece$.exposure <- data[[exposure_vars[j]]]
if (identical(exposure_mode, "single")) piece$.exposure <- data[[exposure_name]]
pieces[[j]] <- piece
}
out <- do.call(rbind, pieces)
rownames(out) <- NULL
tf <- .r4vn_long_factor(
out$.time_source,
reference_index = 1L,
reference_value = ref
)
out$.time <- tf$x
by_display <- NULL
by_reference <- NULL
if (!is.null(by_name)) {
bf <- .r4vn_long_factor(out[[by_name]], reference_index = by_spec$reference_index)
if (length(bf$display_levels) < 2L) stop("`by` must contain at least two observed groups.", call. = FALSE)
out$.by <- bf$x
by_display <- bf$display_levels
by_reference <- bf$reference
}
cov <- .r4vn_long_apply_covariates(out, adjusted_meta)
out <- cov$data
list(
data = out,
outcome_name = .r4vn_long_common_wide_label(data, variables),
time_name = NULL,
time_display = tf$display_levels,
time_reference = tf$reference,
time_continuous = FALSE,
id_name = id_name,
repeated = TRUE,
by_name = by_name,
by_display = by_display,
by_reference = by_reference,
covariates = cov$names,
covariate_map = cov$map,
exposure_mode = exposure_mode
)
}
.r4vn_long_glm_robust_vcov <- function(fit, cluster = NULL) {
X <- tryCatch(stats::model.matrix(fit), error = function(e) NULL)
if (is.null(X) || !nrow(X) || !ncol(X)) return(NULL)
if (inherits(fit, "glm")) {
wr <- tryCatch(stats::residuals(fit, type = "working"), error = function(e) NULL)
ww <- fit$weights
if (is.null(wr) || is.null(ww) || length(wr) != nrow(X) || length(ww) != nrow(X)) return(NULL)
} else {
wr <- tryCatch(stats::residuals(fit), error = function(e) NULL)
ww <- fit$weights
if (is.null(ww)) ww <- rep(1, nrow(X))
if (is.null(wr) || length(wr) != nrow(X) || length(ww) != nrow(X)) return(NULL)
}
score <- X * as.vector(wr * ww)
information <- crossprod(X, X * as.vector(ww))
bread <- tryCatch(
solve(information),
error = function(e) tryCatch(qr.solve(information), error = function(e2) NULL)
)
if (is.null(bread)) return(NULL)
correction <- 1
if (!is.null(cluster)) {
cluster <- as.character(cluster)
if (length(cluster) != nrow(X)) return(NULL)
cluster[is.na(cluster)] <- "<NA>"
U <- rowsum(score, group = cluster, reorder = FALSE)
meat <- crossprod(U)
G <- nrow(U)
N <- nrow(X)
P <- qr(X)$rank
if (G > 1L && N > P) {
correction <- (G / (G - 1)) * ((N - 1) / (N - P))
}
} else {
meat <- crossprod(score)
N <- nrow(X)
P <- qr(X)$rank
if (N > P) correction <- N / (N - P)
}
V <- correction * bread %*% meat %*% bread
dimnames(V) <- list(colnames(X), colnames(X))
V
}
.r4vn_long_safe_scalar <- function(x) {
if (is.null(x) || !length(x)) return(NA_real_)
z <- suppressWarnings(as.numeric(x[1L]))
if (!length(z) || is.na(z) || !is.finite(z)) NA_real_ else z
}
.r4vn_long_fixed_rank <- function(formula, data) {
mf <- tryCatch(
stats::model.frame(formula, data = data, na.action = stats::na.omit),
error = function(e) NULL
)
if (is.null(mf) || !nrow(mf)) {
return(list(full_rank = TRUE, rank = NA_integer_, columns = NA_integer_))
}
mm <- tryCatch(stats::model.matrix(formula, data = mf), error = function(e) NULL)
if (is.null(mm) || !ncol(mm)) {
return(list(full_rank = TRUE, rank = NA_integer_, columns = NA_integer_))
}
q <- qr(mm)
list(full_rank = q$rank == ncol(mm), rank = q$rank, columns = ncol(mm))
}
.r4vn_long_wald <- function(beta, V, terms) {
terms <- intersect(terms, names(beta))
if (!length(terms)) return(NA_real_)
b <- beta[terms]
VV <- V[terms, terms, drop = FALSE]
good <- is.finite(b) & is.finite(diag(VV)) & diag(VV) > 0
b <- b[good]
VV <- VV[good, good, drop = FALSE]
if (!length(b)) return(NA_real_)
inv <- tryCatch(solve(VV), error = function(e) tryCatch(qr.solve(VV), error = function(e2) NULL))
if (is.null(inv)) return(NA_real_)
df <- qr(VV)$rank
if (df < 1L) return(NA_real_)
stat <- as.numeric(t(b) %*% inv %*% b)
if (!is.finite(stat)) return(NA_real_)
stats::pchisq(stat, df = df, lower.tail = FALSE)
}
.r4vn_long_term_columns <- function(fit, fixed_formula, term_label) {
mf <- tryCatch(stats::model.frame(fixed_formula, data = fit$model), error = function(e) NULL)
if (is.null(mf)) {
mf <- tryCatch(fit$model, error = function(e) NULL)
}
if (is.null(mf)) return(character())
mm <- tryCatch(stats::model.matrix(fixed_formula, data = mf), error = function(e) NULL)
if (is.null(mm)) return(character())
assign <- attr(mm, "assign")
labels <- attr(stats::terms(fixed_formula), "term.labels")
position <- which(labels == term_label)
if (!length(position) && identical(term_label, ".time:.by")) {
position <- which(labels %in% c(".time:.by", ".by:.time"))
}
if (!length(position)) return(character())
colnames(mm)[assign %in% position]
}
.r4vn_long_extract_beta_vcov <- function(fit, engine, robust = FALSE, cluster = NULL) {
if (identical(engine, "mixed")) {
if (inherits(fit, "lme")) {
beta <- nlme::fixef(fit)
V <- as.matrix(stats::vcov(fit))
} else {
beta <- stats::coef(fit)
V <- as.matrix(stats::vcov(fit))
}
} else {
beta <- stats::coef(fit)
V <- if (isTRUE(robust)) .r4vn_long_glm_robust_vcov(fit, cluster = cluster) else as.matrix(stats::vcov(fit))
}
list(beta = beta, V = V)
}
.r4vn_long_fit_one <- function(data, outcome_type, effect_type, repeated,
gee, ar1, slope, time_continuous,
covariates, exposure = FALSE) {
has_by <- ".by" %in% names(data)
z_rhs <- if (length(covariates)) paste(covariates, collapse = " + ") else ""
offset_rhs <- if (isTRUE(exposure)) "offset(log(.exposure))" else ""
join_rhs <- function(parts) {
parts <- parts[nzchar(parts)]
if (!length(parts)) "1" else paste(parts, collapse = " + ")
}
main_time <- ".time"
main_by <- if (has_by) ".by" else ""
full_core <- if (has_by) ".time * .by" else ".time"
full_rhs <- join_rhs(c(full_core, z_rhs, offset_rhs))
add_rhs <- join_rhs(c(main_time, main_by, z_rhs, offset_rhs))
time_reduced_rhs <- join_rhs(c(main_by, z_rhs, offset_rhs))
by_reduced_rhs <- join_rhs(c(main_time, z_rhs, offset_rhs))
null_rhs <- join_rhs(c(z_rhs, offset_rhs))
fixed_full <- stats::as.formula(paste(".outcome ~", full_rhs))
fixed_add <- stats::as.formula(paste(".outcome ~", add_rhs))
fixed_time_reduced <- stats::as.formula(paste(".outcome ~", time_reduced_rhs))
fixed_by_reduced <- stats::as.formula(paste(".outcome ~", by_reduced_rhs))
fixed_null <- stats::as.formula(paste(".outcome ~", null_rhs))
engine <- NULL
link <- "identity"
family <- NULL
robust <- FALSE
fits <- list()
correlation <- "independent"
package_note <- NULL
if (identical(outcome_type, "continuous")) {
family <- stats::gaussian()
link <- "identity"
} else if (identical(outcome_type, "binary")) {
if (effect_type %in% c("RR", "PR")) {
family <- stats::poisson(link = "log")
link <- "log"
} else {
family <- stats::binomial(link = "logit")
link <- "logit"
}
} else if (identical(outcome_type, "count")) {
family <- stats::poisson(link = "log")
link <- "log"
} else {
stop("Unsupported outcome type.", call. = FALSE)
}
# Advanced AR(1) GEE remains available when geepack is installed. Routine
# longitudinal RR/PR, binary and count analyses do not require geepack.
request_gee <- isTRUE(gee) || effect_type %in% c("RR", "PR")
use_geepack <- isTRUE(repeated) && isTRUE(ar1) && isTRUE(request_gee) &&
requireNamespace("geepack", quietly = TRUE)
if (isTRUE(repeated) && isTRUE(ar1) && isTRUE(request_gee) && !use_geepack) {
warning(
"`ar1 = TRUE` requires optional package `geepack`; using working-independence cluster-robust inference instead.",
call. = FALSE
)
package_note <- "AR(1) was requested but geepack was unavailable; working-independence cluster-robust inference was used."
}
if (isTRUE(repeated) && use_geepack) {
time_order <- if (is.factor(data$.time)) as.numeric(data$.time) else data$.time
data <- data[order(as.character(data$.id), time_order, na.last = TRUE), , drop = FALSE]
rownames(data) <- NULL
engine <- "gee"
correlation <- "AR(1)"
fit_fun <- function(fixed) {
geepack::geeglm(
formula = fixed,
data = data,
id = data$.id,
family = family,
corstr = "ar1",
std.err = "san.se",
na.action = stats::na.omit
)
}
fits$full <- fit_fun(fixed_full)
if (has_by) fits$add <- fit_fun(fixed_add)
beta_vcov <- list(
beta = stats::coef(fits$full),
V = as.matrix(stats::vcov(fits$full))
)
if (has_by) {
add_beta <- stats::coef(fits$add)
add_V <- as.matrix(stats::vcov(fits$add))
time_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".time")
group_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".by")
interaction_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time:.by")
p_time <- .r4vn_long_wald(add_beta, add_V, time_cols)
p_group <- .r4vn_long_wald(add_beta, add_V, group_cols)
p_interaction <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, interaction_cols)
} else {
time_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time")
p_time <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, time_cols)
p_group <- NA_real_
p_interaction <- NA_real_
}
singular <- FALSE
df_contrast <- Inf
} else if (isTRUE(repeated) && identical(outcome_type, "continuous") && !isTRUE(gee) &&
requireNamespace("nlme", quietly = TRUE)) {
# nlme is an R recommended package and is included in ordinary R installs.
# Before fitting, detect perfect collinearity in the fixed-effects design.
# This avoids exposing cryptic nlme errors such as
# "Singularity in backsolve at level 0, block 1" to R4VN users.
rank_check <- .r4vn_long_fixed_rank(fixed_full, data)
if (!isTRUE(rank_check$full_rank)) {
stop(
"The adjusted longitudinal model cannot be estimated because the fixed-effects predictors are perfectly collinear. ",
"This commonly occurs when an adjustment variable is identical to, or completely determined by, the group/time variables. ",
"Remove or recode the redundant adjustment variable and run tablong() again.",
call. = FALSE
)
}
engine <- "mixed"
correlation <- "random intercept"
fit_fun <- function(fixed) {
random_formula <- if (isTRUE(slope) && isTRUE(time_continuous)) {
stats::as.formula("~ 1 + .time | .id")
} else {
stats::as.formula("~ 1 | .id")
}
nlme::lme(
fixed = fixed,
random = random_formula,
data = data,
method = "ML",
na.action = stats::na.omit,
control = nlme::lmeControl(returnObject = TRUE)
)
}
fits$full <- fit_fun(fixed_full)
if (has_by) {
fits$add <- fit_fun(fixed_add)
fits$time_reduced <- fit_fun(fixed_time_reduced)
fits$by_reduced <- fit_fun(fixed_by_reduced)
} else {
fits$null <- fit_fun(fixed_null)
}
extract_p_compare <- function(a, b) {
z <- tryCatch(stats::anova(a, b), error = function(e) NULL)
if (is.null(z) || nrow(z) < 2L) return(NA_real_)
pcol <- grep("p-value|Pr\\(", names(z), value = TRUE, ignore.case = TRUE)
if (!length(pcol)) return(NA_real_)
.r4vn_long_safe_scalar(z[[pcol[1L]]][2L])
}
if (has_by) {
p_time <- extract_p_compare(fits$time_reduced, fits$add)
p_group <- extract_p_compare(fits$by_reduced, fits$add)
p_interaction <- extract_p_compare(fits$add, fits$full)
} else {
p_time <- extract_p_compare(fits$null, fits$full)
p_group <- NA_real_
p_interaction <- NA_real_
}
beta_vcov <- .r4vn_long_extract_beta_vcov(fits$full, engine)
singular <- FALSE
df_contrast <- Inf
} else {
# Base-R marginal models. For repeated subjects, cluster-robust sandwich
# variance is calculated internally by R4VN, so sandwich/geepack/lme4 are
# not required for routine binary, RR/PR, count, or Gaussian GEE analyses.
repeated_robust <- isTRUE(repeated)
robust <- repeated_robust || (identical(outcome_type, "binary") && effect_type %in% c("RR", "PR"))
engine <- if (repeated_robust) "cluster_robust" else "independent"
correlation <- if (repeated_robust) "working independence" else "independent"
if (isTRUE(repeated) && identical(outcome_type, "continuous") && !isTRUE(gee) &&
!requireNamespace("nlme", quietly = TRUE)) {
warning(
"Recommended package `nlme` is unavailable; using a marginal linear model with subject-clustered robust standard errors.",
call. = FALSE
)
package_note <- "nlme was unavailable; a marginal linear model with subject-clustered robust standard errors was used."
}
if (isTRUE(slope) && identical(outcome_type, "continuous") && !identical(engine, "mixed")) {
warning("`slope = TRUE` requires the recommended package `nlme`; the random slope was not fitted.", call. = FALSE)
}
fit_fun <- function(fixed) {
if (identical(outcome_type, "continuous")) {
stats::lm(fixed, data = data, na.action = stats::na.omit)
} else {
stats::glm(fixed, data = data, family = family, na.action = stats::na.omit)
}
}
fits$full <- fit_fun(fixed_full)
if (has_by) fits$add <- fit_fun(fixed_add)
cluster <- if (repeated_robust) data$.id else NULL
full_V <- if (robust) .r4vn_long_glm_robust_vcov(fits$full, cluster = cluster) else as.matrix(stats::vcov(fits$full))
if (is.null(full_V)) full_V <- as.matrix(stats::vcov(fits$full))
beta_vcov <- list(beta = stats::coef(fits$full), V = full_V)
if (has_by) {
add_V <- if (robust) .r4vn_long_glm_robust_vcov(fits$add, cluster = cluster) else as.matrix(stats::vcov(fits$add))
if (is.null(add_V)) add_V <- as.matrix(stats::vcov(fits$add))
add_beta <- stats::coef(fits$add)
time_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".time")
group_cols <- .r4vn_long_term_columns(fits$add, fixed_add, ".by")
interaction_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time:.by")
p_time <- .r4vn_long_wald(add_beta, add_V, time_cols)
p_group <- .r4vn_long_wald(add_beta, add_V, group_cols)
p_interaction <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, interaction_cols)
} else {
time_cols <- .r4vn_long_term_columns(fits$full, fixed_full, ".time")
p_time <- .r4vn_long_wald(beta_vcov$beta, beta_vcov$V, time_cols)
p_group <- NA_real_
p_interaction <- NA_real_
}
singular <- FALSE
if (repeated_robust) {
G <- length(unique(data$.id[!is.na(data$.id)]))
df_contrast <- if (G > 1L) G - 1L else Inf
} else {
df_contrast <- if (identical(outcome_type, "continuous")) stats::df.residual(fits$full) else Inf
}
}
list(
fit = fits$full,
fits = fits,
engine = engine,
family = family,
link = link,
fixed_formula = fixed_full,
beta = beta_vcov$beta,
V = beta_vcov$V,
p_time = p_time,
p_group = p_group,
p_interaction = p_interaction,
singular = singular,
df_contrast = df_contrast,
correlation = correlation,
package_note = package_note
)
}
.r4vn_long_grid_row <- function(data, fixed_formula, time_value, by_value = NULL) {
nd <- data.frame(.outcome = 0)
if (is.factor(data$.time)) {
nd$.time <- factor(as.character(time_value), levels = levels(data$.time))
} else {
nd$.time <- as.numeric(time_value)
}
if (".by" %in% names(data)) {
nd$.by <- factor(as.character(by_value), levels = levels(data$.by))
}
z_names <- grep("^\\.z[0-9]+$", names(data), value = TRUE)
for (z in z_names) {
if (is.factor(data[[z]])) {
nd[[z]] <- factor(levels(data[[z]])[1L], levels = levels(data[[z]]))
} else {
value <- mean(data[[z]], na.rm = TRUE)
if (!is.finite(value)) value <- 0
nd[[z]] <- value
}
}
if (".exposure" %in% names(data)) nd$.exposure <- 1
mm <- stats::model.matrix(fixed_formula, data = nd)
mm
}
.r4vn_long_contrast <- function(engine_fit, L, level = 0.95) {
beta <- engine_fit$beta
V <- engine_fit$V
common <- intersect(names(beta), colnames(L))
if (!length(common)) {
return(c(estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
}
L2 <- as.numeric(L[, common, drop = FALSE])
names(L2) <- common
b <- beta[common]
VV <- V[common, common, drop = FALSE]
estimate_link <- sum(L2 * b)
variance <- as.numeric(t(L2) %*% VV %*% L2)
if (!is.finite(variance) || variance < 0) {
return(c(estimate = NA_real_, lower = NA_real_, upper = NA_real_, p = NA_real_))
}
se <- sqrt(variance)
alpha <- 1 - level
df <- engine_fit$df_contrast
critical <- if (is.finite(df)) stats::qt(1 - alpha / 2, df = df) else stats::qnorm(1 - alpha / 2)
lower_link <- estimate_link - critical * se
upper_link <- estimate_link + critical * se
statistic <- if (se > 0) estimate_link / se else NA_real_
p <- if (!is.finite(statistic)) NA_real_ else {
if (is.finite(df)) 2 * stats::pt(-abs(statistic), df = df)
else 2 * stats::pnorm(-abs(statistic))
}
if (engine_fit$link %in% c("log", "logit")) {
c(
estimate = exp(estimate_link),
lower = exp(lower_link),
upper = exp(upper_link),
p = p
)
} else {
c(
estimate = estimate_link,
lower = lower_link,
upper = upper_link,
p = p
)
}
}
.r4vn_long_contrast_rows <- function(data, fit, time_display, time_reference,
by_display, time_continuous, change,
pairwise, level, adjust) {
has_by <- ".by" %in% names(data)
out <- list()
index <- 0L
add <- function(type, time, group1 = "", group2 = "", comparison = "", result) {
index <<- index + 1L
out[[index]] <<- data.frame(
type = type,
time = as.character(time),
group1 = as.character(group1),
group2 = as.character(group2),
comparison = as.character(comparison),
estimate = unname(result["estimate"]),
lower = unname(result["lower"]),
upper = unname(result["upper"]),
p = unname(result["p"]),
stringsAsFactors = FALSE
)
}
xrow <- function(time, group = NULL) {
.r4vn_long_grid_row(data, fit$fixed_formula, time, group)
}
if (isTRUE(time_continuous)) {
groups <- if (has_by) by_display else ""
base_time <- if (length(time_display)) min(as.numeric(time_display)) else 0
if (has_by && length(by_display) == 2L) {
for (tt in time_display) {
X1 <- xrow(tt, by_display[1L])
X2 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, X2 - X1, level)
add("between", tt, by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
}
for (g in groups) {
X0 <- xrow(base_time, if (has_by) g else NULL)
X1 <- xrow(base_time + 1, if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, X1 - X0, level)
add("slope", "Per 1 time unit", if (has_by) g else "", "", "Per 1 time unit", result)
}
if (has_by && length(by_display) == 2L) {
X00 <- xrow(base_time, by_display[1L])
X01 <- xrow(base_time + 1, by_display[1L])
X10 <- xrow(base_time, by_display[2L])
X11 <- xrow(base_time + 1, by_display[2L])
result <- .r4vn_long_contrast(fit, (X11 - X10) - (X01 - X00), level)
add("slope_difference", "Per 1 time unit", by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
} else {
if (has_by && length(by_display) == 2L) {
for (tt in time_display) {
X1 <- xrow(tt, by_display[1L])
X2 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, X2 - X1, level)
add("between", tt, by_display[1L], by_display[2L],
paste(by_display[2L], "vs", by_display[1L]), result)
}
}
if (isTRUE(change)) {
follow <- setdiff(time_display, as.character(time_reference))
groups <- if (has_by) by_display else ""
for (tt in follow) {
for (g in groups) {
X0 <- xrow(time_reference, if (has_by) g else NULL)
X1 <- xrow(tt, if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, X1 - X0, level)
add("change", tt, if (has_by) g else "", "",
paste(tt, "vs", time_reference), result)
}
if (has_by && length(by_display) == 2L) {
X00 <- xrow(time_reference, by_display[1L])
X01 <- xrow(tt, by_display[1L])
X10 <- xrow(time_reference, by_display[2L])
X11 <- xrow(tt, by_display[2L])
result <- .r4vn_long_contrast(fit, (X11 - X10) - (X01 - X00), level)
add("change_difference", tt, by_display[1L], by_display[2L],
paste0("(", by_display[2L], " change) vs (", by_display[1L], " change)"),
result)
}
}
}
if (isTRUE(pairwise)) {
if (length(time_display) > 1L) {
tp <- utils::combn(time_display, 2L, simplify = FALSE)
groups <- if (has_by) by_display else ""
for (pair in tp) {
for (g in groups) {
Xa <- xrow(pair[1L], if (has_by) g else NULL)
Xb <- xrow(pair[2L], if (has_by) g else NULL)
result <- .r4vn_long_contrast(fit, Xb - Xa, level)
add("pairwise_time", pair[2L], if (has_by) g else "", "",
paste(pair[2L], "vs", pair[1L]), result)
}
}
}
if (has_by && length(by_display) > 1L) {
gp <- utils::combn(by_display, 2L, simplify = FALSE)
for (tt in time_display) {
for (pair in gp) {
Xa <- xrow(tt, pair[1L])
Xb <- xrow(tt, pair[2L])
result <- .r4vn_long_contrast(fit, Xb - Xa, level)
add("pairwise_group", tt, pair[1L], pair[2L],
paste(pair[2L], "vs", pair[1L]), result)
}
}
}
}
}
if (!length(out)) {
return(data.frame(
type = character(), time = character(), group1 = character(),
group2 = character(), comparison = character(),
estimate = numeric(), lower = numeric(), upper = numeric(), p = numeric(),
p_adjusted = numeric(), stringsAsFactors = FALSE
))
}
result <- do.call(rbind, out)
if (!identical(adjust, "none")) {
result$p_adjusted <- stats::p.adjust(result$p, method = adjust)
} else {
result$p_adjusted <- result$p
}
result
}
.r4vn_long_summary_continuous <- function(x, type, digits, show_n) {
y <- suppressWarnings(as.numeric(x))
y <- y[is.finite(y)]
n <- length(y)
if (!n) return("")
if (identical(type, "median")) {
q <- stats::quantile(y, c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 7)
text <- paste0(
.r4vn_long_fmt(q[2L], digits), " (",
.r4vn_long_fmt(q[1L], digits), ", ",
.r4vn_long_fmt(q[3L], digits), ")"
)
} else if (identical(type, "full")) {
q <- stats::quantile(y, c(.25, .5, .75), na.rm = TRUE, names = FALSE, type = 7)
text <- paste0(
.r4vn_long_fmt(mean(y), digits), " (", .r4vn_long_fmt(stats::sd(y), digits), "); ",
.r4vn_long_fmt(q[2L], digits), " [", .r4vn_long_fmt(q[1L], digits), ", ",
.r4vn_long_fmt(q[3L], digits), "]; ",
.r4vn_long_fmt(min(y), digits), "-", .r4vn_long_fmt(max(y), digits)
)
} else {
text <- paste0(
.r4vn_long_fmt(mean(y), digits),
" (",
.r4vn_long_fmt(stats::sd(y), digits),
")"
)
}
if (isTRUE(show_n)) text <- paste0(text, " [n=", n, "]")
text
}
.r4vn_long_summary_binary <- function(x, event, digits) {
keep <- !is.na(x)
n <- sum(keep)
if (!n) return("")
cases <- sum(as.character(x[keep]) == as.character(event))
pct <- 100 * cases / n
paste0(cases, "/", n, " (", .r4vn_long_fmt(pct, digits), "%)")
}
.r4vn_long_summary_count <- function(x, digits, show_n) {
y <- suppressWarnings(as.numeric(x))
y <- y[is.finite(y)]
n <- length(y)
if (!n) return("")
text <- paste0(.r4vn_long_fmt(mean(y), digits), " (", .r4vn_long_fmt(stats::sd(y), digits), ")")
if (isTRUE(show_n)) text <- paste0(text, " [n=", n, "]")
text
}
.r4vn_long_observed_cell <- function(data, time_value, group_value, outcome_type,
summary_type, event, digits, show_n,
time_continuous) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(time_value)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(time_value)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(group_value)
}
x <- if (".outcome_display" %in% names(data)) {
data$.outcome_display[mask]
} else {
data$.outcome[mask]
}
if (identical(outcome_type, "continuous")) {
return(.r4vn_long_summary_continuous(x, summary_type, digits, show_n))
}
if (identical(outcome_type, "binary")) {
return(.r4vn_long_summary_binary(x, event, digits))
}
.r4vn_long_summary_count(x, digits, show_n)
}
.r4vn_long_effect_name <- function(outcome_type, effect_type) {
if (identical(outcome_type, "continuous")) return("Difference (95% CI)")
if (identical(outcome_type, "count")) return("IRR (95% CI)")
paste0(effect_type, " (95% CI)")
}
.r4vn_long_change_name <- function(outcome_type, effect_type) {
if (identical(outcome_type, "continuous")) return("Change (95% CI)")
if (identical(outcome_type, "count")) return("IRR (95% CI)")
paste0(effect_type, " (95% CI)")
}
.r4vn_long_build_rows <- function(data, outcome_label, outcome_type, summary_type,
event, effect_type, time_display, time_reference,
by_display, time_continuous, contrasts, tests,
change, digits, effect_digits, p_digits,
show_n) {
has_by <- ".by" %in% names(data)
two_groups <- has_by && length(by_display) == 2L
effect_header <- .r4vn_long_effect_name(outcome_type, effect_type)
rows <- list()
rid <- 0L
add_row <- function(values) {
rid <<- rid + 1L
rows[[rid]] <<- values
}
if (!has_by) {
for (tt in time_display) {
cell <- .r4vn_long_observed_cell(
data, tt, NULL, outcome_type, summary_type, event,
digits, show_n, time_continuous
)
if (isTRUE(time_continuous)) {
effect_text <- ""
p_text <- ""
} else {
effect_text <- if (as.character(tt) == as.character(time_reference)) {
"Ref"
} else {
hit <- contrasts[
contrasts$type == "change" &
contrasts$time == as.character(tt),
, drop = FALSE
]
if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
}
hitp <- contrasts[
contrasts$type == "change" &
contrasts$time == as.character(tt),
, drop = FALSE
]
p_text <- if (nrow(hitp)) .r4vn_long_fmt_p(hitp$p_adjusted[1L], p_digits) else ""
}
add_row(c(
Outcome = outcome_label,
Time = as.character(tt),
Summary = cell,
Effect = effect_text,
`p-value` = p_text
))
}
if (isTRUE(time_continuous)) {
hit <- contrasts[contrasts$type == "slope", , drop = FALSE]
add_row(c(
Outcome = "",
Time = "Change per 1 time unit",
Summary = "",
Effect = if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else "",
`p-value` = if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
))
}
add_row(c(
Outcome = "",
Time = "Overall time",
Summary = "",
Effect = "",
`p-value` = .r4vn_long_fmt_p(tests$p_time, p_digits)
))
} else {
for (tt in time_display) {
values <- c(Outcome = outcome_label, Time = as.character(tt))
for (g in by_display) {
values <- c(
values,
stats::setNames(
.r4vn_long_observed_cell(
data, tt, g, outcome_type, summary_type, event,
digits, show_n, time_continuous
),
g
)
)
}
if (two_groups) {
hit <- contrasts[contrasts$type == "between" & contrasts$time == as.character(tt), , drop = FALSE]
if (nrow(hit)) {
values <- c(
values,
stats::setNames(
.r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits),
effect_header
),
`p-value` = .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits)
)
} else {
values <- c(values, stats::setNames("", effect_header), `p-value` = "")
}
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
if (isTRUE(time_continuous)) {
slope_rows <- contrasts[contrasts$type == "slope", , drop = FALSE]
if (nrow(slope_rows)) {
values <- c(Outcome = "", Time = "Change per 1 time unit")
for (g in by_display) {
hit <- slope_rows[slope_rows$group1 == g, , drop = FALSE]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
values <- c(values, stats::setNames(txt, g))
}
if (two_groups) {
hit <- contrasts[contrasts$type == "slope_difference", , drop = FALSE]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
pv <- if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
values <- c(values, stats::setNames(txt, effect_header), `p-value` = pv)
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
} else if (isTRUE(change)) {
follow <- setdiff(as.character(time_display), as.character(time_reference))
for (tt in follow) {
values <- c(Outcome = "", Time = paste0("Change: ", tt, " vs ", time_reference))
for (g in by_display) {
hit <- contrasts[
contrasts$type == "change" &
contrasts$time == tt &
contrasts$group1 == g,
, drop = FALSE
]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
values <- c(values, stats::setNames(txt, g))
}
if (two_groups) {
hit <- contrasts[
contrasts$type == "change_difference" &
contrasts$time == tt,
, drop = FALSE
]
txt <- if (nrow(hit)) .r4vn_long_fmt_ci(hit$estimate[1L], hit$lower[1L], hit$upper[1L], effect_digits) else ""
pv <- if (nrow(hit)) .r4vn_long_fmt_p(hit$p_adjusted[1L], p_digits) else ""
values <- c(values, stats::setNames(txt, effect_header), `p-value` = pv)
} else {
values <- c(values, `p-value` = "")
}
add_row(values)
}
}
test_rows <- list(
c("Overall time", tests$p_time),
c("Overall group", tests$p_group),
c("Time x group", tests$p_interaction)
)
for (z in test_rows) {
values <- c(Outcome = "", Time = z[1L])
for (g in by_display) values <- c(values, stats::setNames("", g))
if (two_groups) values <- c(values, stats::setNames("", effect_header))
values <- c(values, `p-value` = .r4vn_long_fmt_p(as.numeric(z[2L]), p_digits))
add_row(values)
}
}
out <- as.data.frame(do.call(rbind, rows), stringsAsFactors = FALSE, check.names = FALSE)
rownames(out) <- NULL
if ("Effect" %in% names(out)) names(out)[names(out) == "Effect"] <- .r4vn_long_change_name(outcome_type, effect_type)
out
}
.r4vn_long_html_table <- function(data, title, notes, bold_p, p_bold) {
cols <- names(data)
header <- paste0("<th>", .r4vn_long_escape(cols), "</th>", collapse = "")
body <- character(nrow(data))
last_outcome <- NULL
for (i in seq_len(nrow(data))) {
cells <- as.character(data[i, , drop = TRUE])
outcome <- cells[1L]
if (nzchar(outcome) && identical(outcome, last_outcome)) cells[1L] <- ""
if (nzchar(outcome)) last_outcome <- outcome
td <- character(length(cells))
for (j in seq_along(cells)) {
value <- .r4vn_long_escape(cells[j])
class_attr <- ""
if (j == 1L && nzchar(cells[j])) {
value <- paste0('<span class="variable-name">', value, "</span>")
}
if (identical(cols[j], "p-value") && isTRUE(bold_p)) {
raw <- suppressWarnings(as.numeric(sub("^<", "", cells[j])))
if (is.finite(raw) && raw < p_bold) value <- paste0("<strong>", value, "</strong>")
}
td[j] <- paste0("<td", class_attr, ">", value, "</td>")
}
row_class <- if (grepl("^(Overall|Time x group)", cells[2L])) ' class="diagnostic-row"' else ""
body[i] <- paste0("<tr", row_class, ">", paste0(td, collapse = ""), "</tr>")
}
note_html <- if (length(notes)) {
paste0('<div class="model-note">', .r4vn_long_escape(notes), "</div>", collapse = "")
} else ""
paste0(
'<section class="r4vn-table">',
if (!is.null(title) && nzchar(title)) paste0('<div class="table-title">', .r4vn_long_escape(title), "</div>") else "",
"<table><thead><tr>", header, "</tr></thead><tbody>",
paste(body, collapse = ""),
"</tbody></table>",
note_html,
"</section>"
)
}
.r4vn_long_html_document <- function(table_html) {
css <- paste0(
"body{font-family:Arial,'Helvetica Neue',sans-serif;background:#fff;color:#111;margin:18px}",
".table-title{font-size:18px;font-weight:700;margin:0 0 8px}.section-title{font-size:15px;font-weight:700;margin:18px 0 6px}",
".r4vn-table{margin-bottom:28px;overflow-x:auto}",
".r4vn-table table{border-collapse:collapse;width:auto;min-width:760px;border-top:2px solid #111;border-bottom:2px solid #111}",
".r4vn-table th{padding:5px 9px;border-bottom:1.5px solid #111;font-weight:700;white-space:nowrap;text-align:center;background:#fff}",
".r4vn-table th:first-child{text-align:left;min-width:220px}",
".r4vn-table td{padding:4px 9px;vertical-align:top;border:0;white-space:nowrap;text-align:center}",
".r4vn-table td:first-child{text-align:left}",
".variable-name{font-weight:700}",
".diagnostic-row td{border-top:1px solid #aaa;font-style:italic}",
".model-note{font-size:12px;color:#333;margin-top:6px;line-height:1.35}",
".r4vn-plot{margin:18px 0 28px}.plot-wrap{max-width:1100px;overflow-x:auto}.plot-wrap svg{display:block;max-width:100%;height:auto;background:#fff}"
)
paste0(
'<!DOCTYPE html><html><head><meta charset="UTF-8">',
'<meta name="viewport" content="width=device-width,initial-scale=1">',
"<style>", css, "</style></head><body>",
table_html,
"</body></html>"
)
}
.r4vn_long_show_html <- function(file) {
if (!interactive()) return(invisible(file))
viewer <- getOption("viewer")
normalized <- normalizePath(file, winslash = "/", mustWork = FALSE)
if (is.function(viewer)) viewer(normalized) else utils::browseURL(normalized)
invisible(file)
}
.r4vn_long_descriptive_rows <- function(data, outcome_label, outcome_type,
summary_type, event, time_display,
by_display, time_continuous,
digits, show_n) {
groups <- if (".by" %in% names(data)) by_display else "Overall"
rows <- list()
k <- 0L
for (tt in time_display) {
for (gg in groups) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(tt)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(tt)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(gg)
}
# Work with row indices inside the current time x group cell. `observed`
# is cell-sized, whereas `mask` is data-sized; combining them directly
# would recycle vectors and can produce incorrect subject counts/warnings.
idx <- which(mask)
x <- if (".outcome_display" %in% names(data)) data$.outcome_display[idx] else data$.outcome[idx]
observed <- !is.na(x)
if (is.numeric(x)) observed <- observed & is.finite(x)
observed_idx <- idx[observed]
k <- k + 1L
rows[[k]] <- data.frame(
Outcome = outcome_label,
Time = as.character(tt),
Group = as.character(gg),
N = sum(observed),
Missing = sum(!observed),
Subjects = if (length(observed_idx)) {
length(unique(data$.id[observed_idx][!is.na(data$.id[observed_idx])]))
} else {
0L
},
Summary = .r4vn_long_observed_cell(
data, tt, if (".by" %in% names(data)) gg else NULL,
outcome_type, summary_type, event, digits, show_n, time_continuous
),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_tests_table <- function(results, digits = 3L) {
if (!length(results)) return(data.frame())
rows <- lapply(results, function(z) {
data.frame(
Outcome = z$label,
Test = c("Overall time", "Overall group", "Time x group"),
p = c(z$tests$p_time, z$tests$p_group, z$tests$p_interaction),
stringsAsFactors = FALSE,
check.names = FALSE
)
})
out <- do.call(rbind, rows)
out <- out[is.finite(out$p), , drop = FALSE]
rownames(out) <- NULL
if (!nrow(out)) return(out)
out$`p-value` <- vapply(out$p, .r4vn_long_fmt_p, character(1), digits = digits)
out$p <- NULL
out
}
.r4vn_long_contrasts_table <- function(results, effect_digits = 2L, p_digits = 3L) {
if (!length(results)) return(data.frame())
rows <- list()
k <- 0L
for (z in results) {
cc <- z$contrasts
if (is.null(cc) || !nrow(cc)) next
for (i in seq_len(nrow(cc))) {
k <- k + 1L
type_label <- switch(
as.character(cc$type[i]),
between = "Between groups",
change = "Change from reference time",
change_difference = "Difference in change",
pairwise_time = "Pairwise time comparison",
pairwise_group = "Pairwise group comparison",
slope = "Slope",
slope_difference = "Difference in slopes",
as.character(cc$type[i])
)
desc <- as.character(cc$comparison[i])
if (!nzchar(desc)) desc <- as.character(cc$type[i])
if (as.character(cc$type[i]) %in% c("between", "pairwise_group") && nzchar(as.character(cc$time[i]))) {
desc <- paste0(desc, " at ", cc$time[i])
}
if (as.character(cc$type[i]) %in% c("change", "pairwise_time") && nzchar(as.character(cc$group1[i]))) {
desc <- paste0(desc, " in ", cc$group1[i])
}
rows[[k]] <- data.frame(
Outcome = z$label,
Type = type_label,
Contrast = desc,
Effect = z$effect,
`Estimate (95% CI)` = .r4vn_long_fmt_ci(cc$estimate[i], cc$lower[i], cc$upper[i], effect_digits),
`p-value` = .r4vn_long_fmt_p(cc$p_adjusted[i], p_digits),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_diagnostics_one <- function(fit, model_data, outcome_label,
outcome_type, time_display,
time_continuous, repeated) {
model <- fit$fit
n_obs <- .r4vn_long_safe_scalar(tryCatch(stats::nobs(model), error = function(e) nrow(model_data)))
if (!is.finite(n_obs)) n_obs <- nrow(model_data)
subject_ok <- !is.na(model_data$.id)
n_subjects <- length(unique(model_data$.id[subject_ok]))
complete_subjects <- NA_integer_
if (isTRUE(repeated) && !isTRUE(time_continuous) && length(time_display)) {
z <- split(as.character(model_data$.time[subject_ok]), as.character(model_data$.id[subject_ok]))
complete_subjects <- sum(vapply(z, function(x) length(unique(x)) == length(time_display), logical(1)))
}
convergence <- "OK"
if (identical(fit$engine, "mixed")) {
if (inherits(model, "lme")) convergence <- "OK"
} else if (identical(fit$engine, "gee")) {
err <- .r4vn_long_safe_scalar(tryCatch(model$geese$error, error = function(e) NA_real_))
if (is.finite(err) && err != 0) convergence <- paste0("GEE error code ", err)
} else if (inherits(model, "glm") && isFALSE(model$converged)) {
convergence <- "Model did not converge"
}
# AIC/BIC are not reported for GEE or cluster-robust marginal inference.
# In particular, geeglm may return numeric(0), which previously triggered
# `if (is.finite(aic))` with "argument is of length zero".
if (fit$engine %in% c("gee", "cluster_robust")) {
aic <- NA_real_
bic <- NA_real_
} else {
aic <- .r4vn_long_safe_scalar(tryCatch(stats::AIC(model), error = function(e) NA_real_))
bic <- .r4vn_long_safe_scalar(tryCatch(stats::BIC(model), error = function(e) NA_real_))
}
icc <- NA_real_
if (identical(fit$engine, "mixed") && identical(outcome_type, "continuous")) {
if (inherits(model, "lme")) {
vc <- tryCatch(nlme::VarCorr(model), error = function(e) NULL)
if (!is.null(vc)) {
vals <- suppressWarnings(as.numeric(vc[, "Variance"]))
rn <- rownames(vc)
vi <- vals[grepl("Intercept", rn, fixed = TRUE)]
vr <- vals[grepl("Residual", rn, fixed = TRUE)]
if (length(vi) && length(vr) && is.finite(vi[1L] + vr[length(vr)]) && vi[1L] + vr[length(vr)] > 0) {
icc <- vi[1L] / (vi[1L] + vr[length(vr)])
}
}
}
}
engine_label <- switch(
fit$engine,
mixed = "Linear mixed model",
gee = "GEE",
cluster_robust = "Marginal robust regression",
independent = "Independent regression",
fit$engine
)
data.frame(
Outcome = outcome_label,
Type = outcome_type,
Engine = engine_label,
Correlation = if (!is.null(fit$correlation)) fit$correlation else "",
Observations = as.integer(n_obs),
Subjects = as.integer(n_subjects),
`Complete subjects` = complete_subjects,
Singular = if (identical(fit$engine, "mixed")) isTRUE(fit$singular) else NA,
Convergence = convergence,
AIC = if (is.finite(aic)) round(aic, 2) else NA_real_,
BIC = if (is.finite(bic)) round(bic, 2) else NA_real_,
ICC = if (is.finite(icc)) round(icc, 3) else NA_real_,
stringsAsFactors = FALSE,
check.names = FALSE
)
}
.r4vn_long_interpretation <- function(results, alpha = 0.05, p_digits = 3L) {
if (!length(results)) return(data.frame())
rows <- list()
k <- 0L
add <- function(outcome, section, text) {
k <<- k + 1L
rows[[k]] <<- data.frame(
Outcome = outcome, Section = section, Interpretation = text,
stringsAsFactors = FALSE, check.names = FALSE
)
}
for (z in results) {
pt <- z$tests$p_time
pg <- z$tests$p_group
pi <- z$tests$p_interaction
if (is.finite(pi)) {
add(
z$label, "Time x group",
if (pi < alpha) {
paste0("The change over time differs between groups (p=", .r4vn_long_fmt_p(pi, p_digits),
"). Interpret the time-specific and change contrasts rather than the main effects alone.")
} else {
paste0("There is no statistical evidence that the temporal pattern differs between groups (p=",
.r4vn_long_fmt_p(pi, p_digits), ").")
}
)
}
if (is.finite(pt)) {
add(
z$label, "Time",
paste0(if (pt < alpha) "There is evidence of an overall time effect" else "There is no statistical evidence of an overall time effect",
" (p=", .r4vn_long_fmt_p(pt, p_digits), ").")
)
}
if (is.finite(pg)) {
add(
z$label, "Group",
paste0(if (pg < alpha) "There is evidence of an overall group effect" else "There is no statistical evidence of an overall group effect",
" (p=", .r4vn_long_fmt_p(pg, p_digits), ").")
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_plot_rows <- function(data, outcome_label, outcome_type, event,
time_display, by_display, time_continuous,
level = 0.95) {
groups <- if (".by" %in% names(data)) by_display else "Overall"
alpha <- 1 - level
rows <- list()
k <- 0L
for (tt in time_display) {
for (gg in groups) {
if (isTRUE(time_continuous)) {
mask <- is.finite(data$.time) & data$.time == as.numeric(tt)
} else {
mask <- !is.na(data$.time) & as.character(data$.time) == as.character(tt)
}
if (".by" %in% names(data)) {
mask <- mask & !is.na(data$.by) & as.character(data$.by) == as.character(gg)
}
x_display <- if (".outcome_display" %in% names(data)) data$.outcome_display[mask] else data$.outcome[mask]
estimate <- lower <- upper <- NA_real_
n <- sum(!is.na(x_display))
measure <- "Observed value"
if (identical(outcome_type, "continuous")) {
y <- suppressWarnings(as.numeric(x_display))
y <- y[is.finite(y)]
n <- length(y)
if (n) {
estimate <- mean(y)
se <- if (n > 1L) stats::sd(y) / sqrt(n) else NA_real_
crit <- if (n > 1L) stats::qt(1 - alpha / 2, df = n - 1L) else NA_real_
if (is.finite(se) && is.finite(crit)) {
lower <- estimate - crit * se
upper <- estimate + crit * se
}
}
measure <- "Mean"
} else if (identical(outcome_type, "binary")) {
keep <- !is.na(x_display)
n <- sum(keep)
cases <- sum(as.character(x_display[keep]) == as.character(event))
if (n) {
estimate <- 100 * cases / n
ci <- tryCatch(suppressWarnings(stats::prop.test(cases, n, conf.level = level, correct = FALSE)$conf.int),
error = function(e) c(NA_real_, NA_real_))
lower <- 100 * ci[1L]
upper <- 100 * ci[2L]
}
measure <- paste0("Event percentage (", event, ")")
} else {
y <- suppressWarnings(as.numeric(x_display))
keep <- is.finite(y)
n <- sum(keep)
if (".exposure" %in% names(data)) {
e <- suppressWarnings(as.numeric(data$.exposure[mask]))
good <- keep & is.finite(e) & e > 0
events <- sum(y[good])
pt <- sum(e[good])
if (pt > 0) {
estimate <- events / pt
lower <- if (events > 0) stats::qchisq(alpha / 2, 2 * events) / (2 * pt) else 0
upper <- stats::qchisq(1 - alpha / 2, 2 * (events + 1)) / (2 * pt)
}
measure <- "Incidence rate per 1 person-time"
} else {
yy <- y[keep]
if (length(yy)) {
estimate <- mean(yy)
se <- if (length(yy) > 1L) stats::sd(yy) / sqrt(length(yy)) else NA_real_
crit <- if (length(yy) > 1L) stats::qt(1 - alpha / 2, df = length(yy) - 1L) else NA_real_
if (is.finite(se) && is.finite(crit)) {
lower <- max(0, estimate - crit * se)
upper <- estimate + crit * se
}
}
measure <- "Mean count"
}
}
k <- k + 1L
rows[[k]] <- data.frame(
outcome = outcome_label,
time = if (isTRUE(time_continuous)) as.numeric(tt) else as.character(tt),
time_order = match(as.character(tt), as.character(time_display)),
group = as.character(gg),
estimate = estimate,
lower = lower,
upper = upper,
n = n,
measure = measure,
time_continuous = isTRUE(time_continuous),
stringsAsFactors = FALSE,
check.names = FALSE
)
}
}
if (!length(rows)) return(data.frame())
out <- do.call(rbind, rows)
rownames(out) <- NULL
out
}
.r4vn_long_build_plot <- function(plot_data, ci = TRUE, plot_args = list()) {
if (is.null(plot_data) || !nrow(plot_data)) return(NULL)
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
if (length(plot_args) && is.null(names(plot_args))) stop("`plot_args` must be a named list.", call. = FALSE)
use_ci <- if (!is.null(plot_args$ci)) isTRUE(plot_args$ci) else isTRUE(ci)
structure(
list(data = plot_data, ci = use_ci, args = plot_args),
class = "r4vn_long_plot"
)
}
.r4vn_long_draw_plot <- function(plot_data, ci = TRUE, plot_args = list()) {
if (is.null(plot_data) || !nrow(plot_data)) return(invisible(NULL))
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
title <- if (!is.null(plot_args$title)) as.character(plot_args$title)[1L] else "Longitudinal profile"
xlab <- if (!is.null(plot_args$xlab)) as.character(plot_args$xlab)[1L] else "Time"
ylab <- if (!is.null(plot_args$ylab)) as.character(plot_args$ylab)[1L] else "Observed estimate (95% CI)"
line_width <- if (!is.null(plot_args$line_width)) as.numeric(plot_args$line_width)[1L] else 1.8
point_size <- if (!is.null(plot_args$point_size)) as.numeric(plot_args$point_size)[1L] else 1.05
base_size <- if (!is.null(plot_args$base_size)) as.numeric(plot_args$base_size)[1L] else 11
legend_position <- if (!is.null(plot_args$legend_position)) as.character(plot_args$legend_position)[1L] else "bottom"
font_family <- if (!is.null(plot_args$font_family)) as.character(plot_args$font_family)[1L] else "sans"
if (!nzchar(font_family)) font_family <- "sans"
use_ci <- if (!is.null(plot_args$ci)) isTRUE(plot_args$ci) else isTRUE(ci)
grid <- if (!is.null(plot_args$grid)) isTRUE(plot_args$grid) else TRUE
outcomes <- unique(as.character(plot_data$outcome))
n_out <- length(outcomes)
if (!n_out) return(invisible(NULL))
ncol_panels <- if (n_out <= 1L) 1L else if (n_out <= 4L) 2L else 3L
nrow_panels <- ceiling(n_out / ncol_panels)
old <- graphics::par(no.readonly = TRUE)
on.exit(graphics::par(old), add = TRUE)
graphics::par(
mfrow = c(nrow_panels, ncol_panels),
mar = c(4.3, 4.4, 3.2, 1.2),
oma = c(0, 0, if (n_out > 1L) 1.7 else 0, 0),
family = font_family,
cex = max(0.65, base_size / 11)
)
for (outcome in outcomes) {
pd <- plot_data[as.character(plot_data$outcome) == outcome, , drop = FALSE]
continuous_x <- all(pd$time_continuous)
if (continuous_x) {
x_all <- suppressWarnings(as.numeric(pd$time))
x_levels <- sort(unique(x_all[is.finite(x_all)]))
x_at <- x_levels
x_labels <- format(x_levels, trim = TRUE)
} else {
ord <- order(pd$time_order)
x_levels <- unique(as.character(pd$time[ord]))
x_all <- match(as.character(pd$time), x_levels)
x_at <- seq_along(x_levels)
x_labels <- x_levels
}
finite_y <- is.finite(pd$estimate)
y_values <- pd$estimate[finite_y]
if (use_ci) {
y_values <- c(y_values, pd$lower[is.finite(pd$lower)], pd$upper[is.finite(pd$upper)])
}
if (!length(y_values)) y_values <- c(0, 1)
yr <- range(y_values, finite = TRUE)
if (!all(is.finite(yr))) yr <- c(0, 1)
if (diff(yr) == 0) yr <- yr + c(-0.5, 0.5)
pad <- 0.07 * diff(yr)
ylim <- yr + c(-pad, pad)
if (continuous_x) {
xr <- range(x_all, finite = TRUE)
if (!all(is.finite(xr))) xr <- c(0, 1)
if (diff(xr) == 0) xr <- xr + c(-0.5, 0.5)
xlim <- xr
} else {
xlim <- c(0.6, max(1.4, length(x_levels) + 0.4))
}
panel_title <- if (n_out > 1L) outcome else title
graphics::plot(
NA_real_, NA_real_, type = "n", xlim = xlim, ylim = ylim,
xaxt = "n", xlab = xlab, ylab = ylab, main = panel_title,
bty = "l", las = 1
)
graphics::axis(1, at = x_at, labels = x_labels)
if (isTRUE(grid)) {
graphics::abline(h = graphics::axTicks(2), col = "grey90", lty = 1, lwd = 0.7)
}
groups <- unique(as.character(pd$group))
if (!length(groups)) groups <- "Overall"
colors <- plot_args$colors
if (is.null(colors) || !length(colors)) {
colors <- if (length(groups) == 1L) "black" else grDevices::hcl.colors(length(groups), palette = "Dark 3")
}
colors <- rep(colors, length.out = length(groups))
pch <- plot_args$point_shapes
if (is.null(pch) || !length(pch)) pch <- c(16, 17, 15, 18, 3, 4, 8)
pch <- rep(pch, length.out = length(groups))
lty <- plot_args$line_types
if (is.null(lty) || !length(lty)) lty <- seq_along(groups)
lty <- rep(lty, length.out = length(groups))
for (j in seq_along(groups)) {
gg <- groups[j]
z <- pd[as.character(pd$group) == gg, , drop = FALSE]
z <- z[order(z$time_order), , drop = FALSE]
xx <- if (continuous_x) suppressWarnings(as.numeric(z$time)) else match(as.character(z$time), x_levels)
good <- is.finite(xx) & is.finite(z$estimate)
if (sum(good) >= 1L) {
graphics::lines(xx[good], z$estimate[good], col = colors[j], lty = lty[j], lwd = line_width)
graphics::points(xx[good], z$estimate[good], col = colors[j], pch = pch[j], cex = point_size)
}
if (use_ci) {
ci_good <- is.finite(xx) & is.finite(z$lower) & is.finite(z$upper)
if (any(ci_good)) {
graphics::arrows(
xx[ci_good], z$lower[ci_good], xx[ci_good], z$upper[ci_good],
angle = 90, code = 3, length = 0.035,
col = colors[j], lwd = max(0.8, line_width * 0.65)
)
}
}
}
if (!(length(groups) == 1L && identical(groups, "Overall")) && !identical(legend_position, "none")) {
pos <- switch(
tolower(legend_position),
right = "topright", left = "topleft", top = "top", bottom = "bottom",
topright = "topright", topleft = "topleft", bottomright = "bottomright",
bottomleft = "bottomleft", "bottom"
)
graphics::legend(
pos, legend = groups, col = colors, lty = lty, pch = pch,
bty = "n", horiz = identical(pos, "bottom"), cex = 0.85
)
}
}
if (n_out > 1L && nzchar(title)) graphics::mtext(title, outer = TRUE, side = 3, line = 0.3, font = 2)
invisible(NULL)
}
.r4vn_long_base64 <- function(x) {
if (!length(x)) return("")
alphabet <- strsplit("ABCDEFGHIJKLMNOPQRSTUVWXYZabcdefghijklmnopqrstuvwxyz0123456789+/", "", fixed = TRUE)[[1L]]
n <- length(x)
pad <- (3L - n %% 3L) %% 3L
vals <- as.integer(x)
if (pad) vals <- c(vals, rep(0L, pad))
m <- matrix(vals, ncol = 3L, byrow = TRUE)
i1 <- bitwShiftR(m[, 1L], 2L)
i2 <- bitwOr(bitwShiftL(bitwAnd(m[, 1L], 3L), 4L), bitwShiftR(m[, 2L], 4L))
i3 <- bitwOr(bitwShiftL(bitwAnd(m[, 2L], 15L), 2L), bitwShiftR(m[, 3L], 6L))
i4 <- bitwAnd(m[, 3L], 63L)
out <- as.vector(rbind(alphabet[i1 + 1L], alphabet[i2 + 1L], alphabet[i3 + 1L], alphabet[i4 + 1L]))
if (pad >= 1L) out[length(out)] <- "="
if (pad == 2L) out[length(out) - 1L] <- "="
paste0(out, collapse = "")
}
.r4vn_long_plot_html <- function(graph) {
if (is.null(graph) || !inherits(graph, "r4vn_long_plot")) return("")
pd <- graph$data
if (is.null(pd) || !nrow(pd)) return("")
n_out <- length(unique(as.character(pd$outcome)))
width <- if (!is.null(graph$args$viewer_width)) as.numeric(graph$args$viewer_width)[1L] else 9.2
height <- if (!is.null(graph$args$viewer_height)) as.numeric(graph$args$viewer_height)[1L] else {
if (n_out <= 1L) 5.2 else 4.4 * ceiling(n_out / if (n_out <= 4L) 2 else 3)
}
if (!is.finite(width) || width <= 0) width <- 9.2
if (!is.finite(height) || height <= 0) height <- 5.2
# Prefer vector SVG in Viewer. If the local R build cannot create SVG,
# fall back to an inline PNG data URI; both paths use only base/recommended R.
tf <- tempfile(pattern = "r4vn-tablong-plot-", fileext = ".svg")
opened <- FALSE
ok <- tryCatch({
svg_family <- if (!is.null(graph$args$font_family)) as.character(graph$args$font_family)[1L] else "sans"
if (!nzchar(svg_family)) svg_family <- "sans"
grDevices::svg(tf, width = width, height = height, pointsize = 11, family = svg_family)
opened <- TRUE
.r4vn_long_draw_plot(pd, ci = graph$ci, plot_args = graph$args)
grDevices::dev.off()
opened <- FALSE
TRUE
}, error = function(e) FALSE)
if (opened) try(grDevices::dev.off(), silent = TRUE)
plot_markup <- ""
if (ok && file.exists(tf)) {
lines <- readLines(tf, warn = FALSE, encoding = "UTF-8")
start <- grep("<svg", lines, fixed = TRUE)[1L]
if (length(start) && !is.na(start)) {
plot_markup <- paste(lines[start:length(lines)], collapse = "\n")
}
}
unlink(tf)
if (!nzchar(plot_markup)) {
png_file <- tempfile(pattern = "r4vn-tablong-plot-", fileext = ".png")
png_open <- FALSE
png_ok <- tryCatch({
grDevices::png(
png_file,
width = max(900L, as.integer(width * 120)),
height = max(560L, as.integer(height * 120)),
res = 120,
type = if (.Platform$OS.type == "windows") "windows" else "cairo"
)
png_open <- TRUE
.r4vn_long_draw_plot(pd, ci = graph$ci, plot_args = graph$args)
grDevices::dev.off()
png_open <- FALSE
TRUE
}, error = function(e) FALSE)
if (png_open) try(grDevices::dev.off(), silent = TRUE)
if (png_ok && file.exists(png_file)) {
raw <- readBin(png_file, what = "raw", n = file.info(png_file)$size)
plot_markup <- paste0(
'<img alt="Longitudinal plot" src="data:image/png;base64,',
.r4vn_long_base64(raw), '" />'
)
}
unlink(png_file)
}
if (!nzchar(plot_markup)) return("")
paste0(
'<section class="r4vn-plot">',
'<div class="section-title">Longitudinal plot</div>',
'<div class="plot-wrap">', plot_markup, '</div>',
'</section>'
)
}
.r4vn_long_html_simple_table <- function(data, title) {
if (is.null(data) || !is.data.frame(data) || !nrow(data)) return("")
x <- data
# Viewer readability: a one/few-row table with many columns is easier to
# inspect after transposition. Keep the returned R object unchanged; this is
# presentation-only and follows the same compact Viewer rule used elsewhere
# in R4VN.
if (ncol(x) >= 8L && nrow(x) <= 4L) {
labels <- if ("Outcome" %in% names(x)) as.character(x$Outcome) else paste0("Result ", seq_len(nrow(x)))
labels <- make.unique(ifelse(is.na(labels) | !nzchar(labels), paste0("Result ", seq_along(labels)), labels))
payload <- if ("Outcome" %in% names(x)) x[, setdiff(names(x), "Outcome"), drop = FALSE] else x
tx <- data.frame(Metric = names(payload), stringsAsFactors = FALSE, check.names = FALSE)
for (i in seq_len(nrow(payload))) {
tx[[labels[i]]] <- vapply(payload, function(v) as.character(v[i]), character(1))
}
x <- tx
}
for (nm in names(x)) {
if (is.numeric(x[[nm]])) {
x[[nm]] <- ifelse(is.na(x[[nm]]), "", format(x[[nm]], trim = TRUE, scientific = FALSE))
} else {
x[[nm]][is.na(x[[nm]])] <- ""
}
}
header <- paste0("<th>", .r4vn_long_escape(names(x)), "</th>", collapse = "")
body <- vapply(seq_len(nrow(x)), function(i) {
vals <- vapply(as.list(x[i, , drop = FALSE]), function(z) as.character(z[1L]), character(1))
paste0("<tr>", paste0("<td>", .r4vn_long_escape(vals), "</td>", collapse = ""), "</tr>")
}, character(1))
paste0(
'<section class="r4vn-table secondary-table">',
'<div class="section-title">', .r4vn_long_escape(title), '</div>',
'<table><thead><tr>', header, '</tr></thead><tbody>',
paste0(body, collapse = ""), '</tbody></table></section>'
)
}
#' Longitudinal and Repeated-Measures Analysis
#'
#' Performs publication-ready longitudinal or repeated-measures analysis from
#' either long or wide data. `tablong()` is designed for the usual biomedical
#' workflow: describe each time point, test overall time and group effects,
#' test the time-by-group interaction, estimate clinically interpretable
#' contrasts, retain fitted models for advanced use, and optionally create a
#' longitudinal profile plot and a cautious interpretation table.
#'
#' The interface follows the R4VN principle of keeping routine analysis simple.
#' In most studies the essential call is only `vars()`, `time`, `id`, and
#' optionally `by`. Repeated continuous outcomes use a random-intercept model
#' when R's recommended `nlme` package is available. Repeated binary/count
#' outcomes use marginal regression with subject-clustered robust standard errors
#' calculated internally by R4VN, so routine analyses need no extra package.
#'
#' @param data Optional data frame. If `NULL`, the active data selected by
#' `usedf()` are used.
#' @param vars Outcome specification created by `vars()`. In long data, several
#' outcomes can be analyzed in one call. Continuous outcomes use the `c.`
#' prefix, `q.` requests median (IQR) descriptive display, and `f.` requests a
#' fuller continuous summary. An unprefixed two-level variable is treated as
#' binary. In wide data, the variables in `vars()` are repeated measurements
#' of the same outcome.
#' @param time In long data, an unquoted time variable. Use `c.timevar` for a
#' continuous linear time effect or `b2.timevar`, `b3.timevar`, etc. to select
#' the categorical reference level. In wide data, provide display labels such
#' as `c("Baseline", "Month 3", "Month 6")`; when omitted, the repeated
#' variable names are used as time labels.
#' @param id Optional subject identifier. In wide data it is optional because
#' each source row represents one subject. In long data, repeated IDs trigger
#' longitudinal analysis with within-subject correlation. If no ID is supplied,
#' or IDs do not repeat across time, observations are analyzed as repeated
#' cross-sectional samples.
#' @param by Optional grouping variable, for example treatment group. Prefix
#' `b2.`, `b3.`, etc. selects the reference group.
#' @param ref Optional categorical reference time label. This overrides a `bN.`
#' prefix supplied in `time`.
#' @param event Event level for binary outcomes. A single value applies to every
#' binary outcome; a named vector can specify a different event for each
#' outcome, for example `event = c(controlled = "Yes", admitted = "Yes")`.
#' @param adjusted Optional adjustment variables created by `vars()`. Use the
#' same R4VN prefixes as elsewhere, for example `vars(c.age, sex, b2.site)`.
#' @param gee Logical. Request a population-average marginal model for repeated
#' data. By default R4VN uses working-independence regression with subject-
#' clustered robust sandwich standard errors calculated internally. No extra
#' package is required. Set `ar1 = TRUE` to request AR(1) GEE through optional
#' package `geepack` when it is installed.
#' @param ar1 Logical. Request an AR(1) working correlation for a marginal GEE.
#' This advanced option uses optional package `geepack`. If it is unavailable,
#' `tablong()` warns and falls back to working-independence cluster-robust
#' inference instead of stopping the analysis.
#' @param slope Logical. With a mixed model and continuous time
#' (`time = c.month`), include a subject-specific random linear time slope in
#' addition to the random intercept.
#' @param or,rr,pr Logical effect switches for binary outcomes. Odds ratio is
#' the default when none is selected. `rr = TRUE` reports risk ratios and
#' `pr = TRUE` reports prevalence ratios using modified Poisson regression with
#' robust variance; repeated subjects use subject-clustered robust variance.
#' No external sandwich/GEE package is needed unless `ar1 = TRUE` is requested.
#' Only one switch may be `TRUE`.
#' @param count Logical. Treat numeric outcomes as non-negative counts and fit
#' Poisson models. Count outcomes report incidence-rate ratios (IRR).
#' @param exposure Optional positive exposure/person-time variable for count
#' models. In wide data it may also be `vars(exp0, exp1, ...)`, with one
#' exposure variable per repeated count variable.
#' @param change Logical. For categorical time, report change from the reference
#' time. With exactly two groups, also report the difference in change, i.e.
#' the usual difference-in-differences contrast. Default `TRUE`.
#' @param pairwise Logical. Calculate all available time and group pairwise
#' contrasts and retain them in `$contrasts` and `$contrasts_table`. The main
#' publication table remains compact. Default `FALSE`.
#' @param adjust Multiplicity adjustment applied to contrast p-values. Any method
#' accepted by `p.adjust()` may be used, including `"none"`, `"holm"`,
#' `"bonferroni"`, and `"BH"`.
#' @param missing Logical. Append cell-specific `n` to continuous/count summary
#' cells. Binary cells always show event/total. Detailed observed and missing
#' counts are always available in `$descriptive`.
#' @param level Confidence level, default 0.95.
#' @param digit Decimal places for descriptive summaries.
#' @param p_digit Decimal places for p-values.
#' @param effect_digit Decimal places for model effects and confidence intervals.
#' @param bold_p Logical. Bold p-values below `p_bold` in the HTML Viewer.
#' @param p_bold Threshold used by `bold_p`.
#' @param diagnostics Logical. Include the compact model-diagnostics table in the
#' HTML Viewer. Diagnostics are always retained in `$diagnostics`; default
#' `FALSE` keeps the primary Viewer concise.
#' @param diagnosis Logical singular alias for `diagnostics`, provided for consistency with other R4VN regression commands. Default `FALSE`. When explicitly supplied, it overrides `diagnostics`; omitting it preserves backward-compatible use of `diagnostics`.
#' @param interpretation Logical. Add a cautious deterministic interpretation
#' table. The default is `FALSE`. The interpretation emphasizes the
#' time-by-group interaction when present and does not replace substantive or
#' clinical interpretation by the researcher.
#' @param plot Logical. Create an observed longitudinal profile plot with 95%
#' confidence intervals, include the same plot directly in the HTML Viewer,
#' display it in the Plot pane, and store its specification in `$graph` and
#' `$plots$trajectory`. Plotting uses base R graphics; `ggplot2` is not required.
#' Default `FALSE`.
#' @param plot_args Named list controlling the profile plot. Supported entries
#' include `title`, `xlab`, `ylab`, `ci`, `line_width`, `point_size`,
#' `base_size`, `legend_position`, `font_family`, `colors`, `point_shapes`,
#' `line_types`, and `grid`. Generic `sans` is the default font for reliable
#' display in RStudio Viewer, browsers, Windows, macOS, and Linux.
#' @param name Logical. Display the original variable name after its variable
#' label in the main table.
#' @param title Optional table title.
#' @param file Optional HTML file path. When omitted, a temporary HTML file is
#' created. This file is the formatted Viewer report, not a replacement for
#' `tabexport()`.
#' @param raw Logical retained for backward compatibility. Raw models, tests,
#' contrasts, standardized long data, and reporting tables are always retained
#' in the returned object.
#' @param show Logical. Open the formatted HTML report in the Viewer/browser.
#' Default `TRUE`.
#'
#' @details
#' **Data format.** If `time` names a column in `data`, input is treated as long.
#' If `vars()` contains multiple repeated variables and `time` is a vector of
#' labels (or omitted), input is treated as wide and is reshaped internally.
#' The original data frame is never modified.
#'
#' **Continuous outcomes.** Repeated subjects use a random-intercept linear
#' mixed model through R's recommended `nlme` package when available.
#' `slope = TRUE` adds a random linear time slope when time is continuous. If
#' `nlme` is unavailable, R4VN falls back to a marginal linear model with
#' subject-clustered robust standard errors. Repeated cross-sectional data use
#' ordinary linear models. `gee = TRUE` explicitly requests the marginal model.
#'
#' **Binary outcomes.** The default effect is an odds ratio from logistic
#' regression. For repeated subjects, R4VN calculates subject-clustered robust
#' standard errors internally. `rr = TRUE` and `pr = TRUE` use modified Poisson
#' regression with robust variance and report RR or PR. This avoids requiring
#' `lme4`, `sandwich`, or `geepack` for routine binary longitudinal analysis.
#'
#' **Count outcomes.** `count = TRUE` fits a Poisson model and reports IRR.
#' Repeated subjects use subject-clustered robust standard errors calculated
#' internally. Supplying `exposure` adds `offset(log(exposure))` and the observed
#' plot is an incidence rate per one person-time unit.
#'
#' **Categorical time.** The model includes time, group when supplied, and the
#' time-by-group interaction. The main table reports observed summaries at each
#' time. With two groups it also reports the between-group effect at each time,
#' within-group change from the reference time, and the difference in change.
#' The omnibus `Time x group` p-value is the formal test that temporal changes
#' differ between groups.
#'
#' **Continuous time.** `time = c.month` estimates change per one time unit. With
#' a group variable, group-specific slopes and their difference are returned.
#'
#' **More than two groups.** The main table remains intentionally compact and
#' shows omnibus tests. Set `pairwise = TRUE` to obtain all model-based pairwise
#' comparisons in `$contrasts_table`; use `adjust` to control multiplicity.
#'
#' **Descriptive prefixes.** `q.` and `f.` affect the observed descriptive
#' summary only. Inferential effects remain based on the selected mean model;
#' they do not fit median regression.
#'
#' **Missing values.** Each model uses observations complete for that outcome,
#' time, ID/group, exposure if required, and adjustment variables. Descriptive
#' counts are retained separately so that missingness and attrition can be
#' reviewed before publication.
#'
#' **Package dependencies.** Routine `tablong()` analyses and plots are designed
#' to run with base/recommended R only. `nlme` is used for continuous random-
#' effects models and is bundled with standard R installations. `geepack` is
#' optional and used only when an AR(1) GEE is explicitly requested. `lme4`,
#' `sandwich`, `broom`, and `ggplot2` are not required by `tablong()`.
#'
#' **Returned reporting contract.** For programmatic reuse, the object includes
#' a flat main table plus descriptive, omnibus-test, contrast, diagnostics,
#' interpretation, plot, model, and metadata components. The object also
#' inherits from `r4vn_tab`, so existing `tabexport()` workflows continue to
#' work.
#'
#' @return Invisibly returns an object of class
#' `c("r4vn_tablong", "r4vn_tab")`. Important components are:
#' `$data` (main publication table), `$descriptive`, `$tests` and
#' `$tests_table`, `$contrasts` and `$contrasts_table`, `$diagnostics`,
#' `$interpretation`, `$tables`, `$graph`, `$plots`, `$models`, `$long_data`,
#' `$metadata`, `$html`, `$file`, `$results`, and `$call`.
#'
#' @seealso `vars()`, `tab()`, `tabexport()`, `usedf()`
#' @family R4VN tables
#' @export
#'
#' @examples
#' # -------------------------------------------------------------------------
#' # 1. Repeated cross-sectional continuous outcome: no optional package needed
#' # -------------------------------------------------------------------------
#' set.seed(11)
#' d <- data.frame(
#' period = factor(rep(c("Before", "After"), each = 80),
#' levels = c("Before", "After")),
#' group = factor(rep(rep(c("Control", "Intervention"), each = 40), 2)),
#' age = rnorm(160, 45, 10)
#' )
#' d$score <- 50 + 2 * (d$period == "After") +
#' 5 * (d$group == "Intervention") +
#' 4 * (d$period == "After" & d$group == "Intervention") +
#' 0.15 * d$age + rnorm(160, 0, 8)
#'
#' z <- tablong(
#' d, vars = vars(c.score), time = period, by = group,
#' adjusted = vars(c.age), show = FALSE
#' )
#' z$data
#' z$tests_table
#' z$contrasts_table
#'
#' \donttest{
#' # -----------------------------------------------------------------------
#' # 2. The same analysis with interpretation and a publication profile plot
#' # -----------------------------------------------------------------------
#' z2 <- tablong(
#' d, vars = vars(c.score), time = period, by = group,
#' adjusted = vars(c.age), interpretation = TRUE,
#' plot = TRUE, show = FALSE
#' )
#' z2$interpretation
#' z2$diagnostics
#' plot(z2)
#'
#' # -----------------------------------------------------------------------
#' # 3. No comparison group: change over time only
#' # -----------------------------------------------------------------------
#' z_time <- tablong(
#' d, vars = vars(c.score), time = period,
#' adjusted = vars(c.age), show = FALSE
#' )
#' z_time$tests_table
#'
#' # -----------------------------------------------------------------------
#' # 4. Change the reference time by value or by bN. prefix
#' # -----------------------------------------------------------------------
#' z_ref1 <- tablong(d, vars = vars(c.score), time = period,
#' by = group, ref = "After", show = FALSE)
#' z_ref2 <- tablong(d, vars = vars(c.score), time = b2.period,
#' by = group, show = FALSE)
#'
#' # -----------------------------------------------------------------------
#' # 5. Median/IQR or full descriptive display, while inference remains a
#' # mean model
#' # -----------------------------------------------------------------------
#' z_median <- tablong(d, vars = vars(q.score), time = period,
#' by = group, show = FALSE)
#' z_full <- tablong(d, vars = vars(f.score), time = period,
#' by = group, missing = TRUE, show = FALSE)
#'
#' # -----------------------------------------------------------------------
#' # 6. Long repeated data: linear mixed model
#' # -----------------------------------------------------------------------
#' set.seed(12)
#' n_subject <- 60
#' dl <- expand.grid(
#' id = seq_len(n_subject),
#' visit = factor(c("Baseline", "Month 3", "Month 6"),
#' levels = c("Baseline", "Month 3", "Month 6"))
#' )
#' dl <- dl[order(dl$id, dl$visit), ]
#' trt <- factor(sample(c("Control", "Intervention"), n_subject, TRUE),
#' levels = c("Control", "Intervention"))
#' dl$treatment <- rep(trt, each = 3)
#' u <- rnorm(n_subject, 0, 6)
#' dl$sbp <- 140 + u[dl$id] - 3 * (dl$visit == "Month 3") -
#' 5 * (dl$visit == "Month 6") -
#' 4 * (dl$treatment == "Intervention" & dl$visit == "Month 6") +
#' rnorm(nrow(dl), 0, 5)
#'
#' mixed <- tablong(
#' dl, vars = vars(c.sbp), time = visit, id = id,
#' by = treatment, diagnostics = TRUE, show = FALSE
#' )
#' mixed$models[[1]]
#' mixed$diagnostics
#'
#' # ---------------------------------------------------------------------
#' # 7. Wide repeated data: internally converted to long format
#' # ---------------------------------------------------------------------
#' dw <- data.frame(
#' id = seq_len(n_subject), treatment = trt,
#' sbp0 = rnorm(n_subject, 140, 10)
#' )
#' dw$sbp3 <- dw$sbp0 - 3 + rnorm(n_subject, 0, 4)
#' dw$sbp6 <- dw$sbp0 - 5 - 4 * (dw$treatment == "Intervention") +
#' rnorm(n_subject, 0, 4)
#' wide <- tablong(
#' dw, vars = vars(c.sbp0, c.sbp3, c.sbp6),
#' time = c("Baseline", "Month 3", "Month 6"),
#' id = id, by = treatment, show = FALSE
#' )
#' wide$input_format
#' head(wide$long_data)
#'
#' # ---------------------------------------------------------------------
#' # 8. Binary repeated outcome: marginal logistic model with clustered SE and OR
#' # ---------------------------------------------------------------------
#' p <- plogis(-1 + 0.4 * (dl$visit == "Month 6") +
#' 0.6 * (dl$treatment == "Intervention"))
#' dl$controlled <- factor(rbinom(nrow(dl), 1, p),
#' levels = c(0, 1), labels = c("No", "Yes"))
#' binary_or <- tablong(
#' dl, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", show = FALSE
#' )
#' binary_or$contrasts_table
#'
#' # ---------------------------------------------------------------------
#' # 9. Continuous time and random slope
#' # ---------------------------------------------------------------------
#' ds <- expand.grid(id = seq_len(50), month = c(0, 3, 6, 12))
#' ds <- ds[order(ds$id, ds$month), ]
#' ds$group <- factor(rep(sample(c("Control", "Intervention"), 50, TRUE), each = 4))
#' b0 <- rnorm(50, 0, 5)
#' b1 <- rnorm(50, 0, 0.15)
#' ds$score <- 50 + b0[ds$id] + (-0.3 + b1[ds$id]) * ds$month -
#' 0.25 * ds$month * (ds$group == "Intervention") + rnorm(nrow(ds), 0, 3)
#' slope_fit <- tablong(
#' ds, vars = vars(c.score), time = c.month, id = id,
#' by = group, slope = TRUE, show = FALSE
#' )
#' slope_fit$contrasts_table
#'
#' # ---------------------------------------------------------------------
#' # 10. Repeated count outcome with person-time offset -> IRR
#' # ---------------------------------------------------------------------
#' dc <- dl
#' dc$person_time <- runif(nrow(dc), 0.8, 1.2)
#' rate <- exp(0.2 + 0.2 * (dc$visit == "Month 6") -
#' 0.3 * (dc$treatment == "Intervention" & dc$visit == "Month 6"))
#' dc$events <- rpois(nrow(dc), rate * dc$person_time)
#' count_fit <- tablong(
#' dc, vars = vars(c.events), time = visit, id = id,
#' by = treatment, count = TRUE, exposure = person_time,
#' show = FALSE
#' )
#' count_fit$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 11. RR/PR without extra packages; optional AR(1) GEE
#' # -----------------------------------------------------------------------
#' set.seed(13)
#' n_subject <- 70
#' dg <- expand.grid(
#' id = seq_len(n_subject),
#' visit = factor(c("Baseline", "Month 6"),
#' levels = c("Baseline", "Month 6"))
#' )
#' dg <- dg[order(dg$id, dg$visit), ]
#' dg$treatment <- factor(
#' rep(sample(c("Control", "Intervention"), n_subject, TRUE), each = 2),
#' levels = c("Control", "Intervention")
#' )
#' p <- plogis(-1 + 0.3 * (dg$visit == "Month 6") +
#' 0.4 * (dg$treatment == "Intervention"))
#' dg$controlled <- factor(rbinom(nrow(dg), 1, p),
#' levels = c(0, 1), labels = c("No", "Yes"))
#'
#' fit_rr <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", rr = TRUE, show = FALSE
#' )
#' fit_pr <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", pr = TRUE,
#' show = FALSE
#' )
#' fit_rr$contrasts_table
#' fit_pr$contrasts_table
#'
#' # AR(1) is advanced and uses geepack only when explicitly requested.
#' if (requireNamespace("geepack", quietly = TRUE)) {
#' fit_pr_ar1 <- tablong(
#' dg, vars = vars(controlled), time = visit, id = id,
#' by = treatment, event = "Yes", pr = TRUE, gee = TRUE, ar1 = TRUE,
#' show = FALSE
#' )
#' }
#'
#' # -----------------------------------------------------------------------
#' # 12. Three or more groups: keep main table compact, request pairwise tests
#' # -----------------------------------------------------------------------
#' set.seed(14)
#' dm <- data.frame(
#' period = factor(rep(c("Baseline", "Follow-up"), each = 90),
#' levels = c("Baseline", "Follow-up")),
#' arm = factor(rep(rep(c("A", "B", "C"), each = 30), 2))
#' )
#' dm$score <- rnorm(nrow(dm), 50 + 2 * (dm$period == "Follow-up") +
#' 2 * (dm$arm == "B") + 4 * (dm$arm == "C"), 7)
#' multi_arm <- tablong(
#' dm, vars = vars(c.score), time = period, by = arm,
#' pairwise = TRUE, adjust = "holm", show = FALSE
#' )
#' multi_arm$tests_table
#' multi_arm$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 13. Several outcomes in one long-data analysis
#' # -----------------------------------------------------------------------
#' d$positive <- factor(
#' rbinom(nrow(d), 1, plogis(-1 + 0.5 * (d$period == "After"))),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#' multi_outcome <- tablong(
#' d, vars = vars(c.score, positive), time = period,
#' by = group, event = "Yes", show = FALSE
#' )
#' multi_outcome$data
#' multi_outcome$descriptive
#'
#' # -----------------------------------------------------------------------
#' # 14. Export remains compatible with ordinary R4VN table workflows
#' # -----------------------------------------------------------------------
#' export_data <- tabexport(z)
#' head(export_data)
#'
#' # -----------------------------------------------------------------------
#' # 15. Missing outcomes/attrition: inspect counts before publication
#' # -----------------------------------------------------------------------
#' d_missing <- d
#' d_missing$score[c(2, 7, 21, 100)] <- NA
#' miss_fit <- tablong(
#' d_missing, vars = vars(c.score), time = period, by = group,
#' missing = TRUE, diagnostics = TRUE, show = FALSE
#' )
#' miss_fit$descriptive
#' miss_fit$diagnostics
#'
#' # -----------------------------------------------------------------------
#' # 16. Several binary outcomes can use a named event vector
#' # -----------------------------------------------------------------------
#' d$admitted <- factor(
#' rbinom(nrow(d), 1, plogis(-1.4 + 0.4 * (d$period == "After"))),
#' levels = c(0, 1), labels = c("No", "Yes")
#' )
#' binary_set <- tablong(
#' d, vars = vars(positive, admitted), time = period, by = group,
#' event = c(positive = "Yes", admitted = "Yes"), show = FALSE
#' )
#' binary_set$tests_table
#'
#' # -----------------------------------------------------------------------
#' # 17. Continuous Gaussian GEE with AR(1) working correlation
#' # -----------------------------------------------------------------------
#' if (requireNamespace("geepack", quietly = TRUE)) {
#' set.seed(17)
#' dg2 <- expand.grid(id = seq_len(60), month = c(0, 3, 6, 12))
#' dg2 <- dg2[order(dg2$id, dg2$month), ]
#' dg2$group <- factor(rep(sample(c("Control", "Intervention"), 60, TRUE), each = 4))
#' dg2$score <- 55 - 0.2 * dg2$month -
#' 0.15 * dg2$month * (dg2$group == "Intervention") + rnorm(nrow(dg2), 0, 5)
#' gee_cont <- tablong(
#' dg2, vars = vars(c.score), time = c.month, id = id, by = group,
#' gee = TRUE, ar1 = TRUE, show = FALSE
#' )
#' gee_cont$contrasts_table
#' }
#'
#' # -----------------------------------------------------------------------
#' # 18. Wide count data can supply one person-time variable per time point
#' # -----------------------------------------------------------------------
#' set.seed(18)
#' nw <- 50
#' wc <- data.frame(
#' id = seq_len(nw),
#' group = factor(sample(c("Control", "Intervention"), nw, TRUE)),
#' pt0 = runif(nw, 0.8, 1.2),
#' pt6 = runif(nw, 0.8, 1.2)
#' )
#' wc$event0 <- rpois(nw, 1.2 * wc$pt0)
#' wc$event6 <- rpois(nw,
#' exp(log(1.2) - 0.25 * (wc$group == "Intervention")) * wc$pt6)
#' wide_count <- tablong(
#' wc, vars = vars(c.event0, c.event6),
#' time = c("Baseline", "Month 6"), id = id, by = group,
#' count = TRUE, exposure = vars(pt0, pt6), show = FALSE
#' )
#' wide_count$contrasts_table
#'
#' # -----------------------------------------------------------------------
#' # 19. Replot an existing result without refitting the statistical model
#' # -----------------------------------------------------------------------
#' plot(z, ci = FALSE, title = "Observed longitudinal profile",
#' base_size = 12, legend_position = "right")
#' }
tablong <- function(data = NULL, vars,
time = NULL, id = NULL, by = NULL,
ref = NULL, event = NULL, adjusted = NULL,
gee = FALSE, ar1 = FALSE, slope = FALSE,
or = FALSE, rr = FALSE, pr = FALSE,
count = FALSE, exposure = NULL,
change = TRUE, pairwise = FALSE,
adjust = "none", missing = FALSE,
level = 0.95, digit = 1, p_digit = 3, effect_digit = 2,
bold_p = TRUE, p_bold = 0.05,
diagnostics = FALSE, diagnosis = FALSE, interpretation = FALSE,
plot = FALSE, plot_args = list(),
name = FALSE, title = NULL, file = NULL,
raw = FALSE, show = TRUE) {
caller <- parent.frame()
data <- .r4vn_resolve_analysis_data(data)
meta <- .r4vn_long_meta(vars, "vars", data = data)
flags <- c(or = isTRUE(or), rr = isTRUE(rr), pr = isTRUE(pr))
if (sum(flags) > 1L) stop("Use only one of `or = TRUE`, `rr = TRUE`, or `pr = TRUE`.", call. = FALSE)
if (!is.logical(count) || length(count) != 1L || is.na(count)) stop("`count` must be TRUE or FALSE.", call. = FALSE)
if (isTRUE(count) && any(flags)) {
stop("`or`, `rr`, and `pr` are binary-outcome options and cannot be combined with `count = TRUE`.", call. = FALSE)
}
if (!is.logical(gee) || length(gee) != 1L || is.na(gee)) stop("`gee` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(ar1) || length(ar1) != 1L || is.na(ar1)) stop("`ar1` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(slope) || length(slope) != 1L || is.na(slope)) stop("`slope` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(change) || length(change) != 1L || is.na(change)) stop("`change` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(pairwise) || length(pairwise) != 1L || is.na(pairwise)) stop("`pairwise` must be TRUE or FALSE.", call. = FALSE)
if (!missing(diagnosis)) {
if (!is.logical(diagnosis) || length(diagnosis) != 1L || is.na(diagnosis)) stop("`diagnosis` must be TRUE or FALSE.", call. = FALSE)
diagnostics <- diagnosis
}
if (!is.logical(diagnostics) || length(diagnostics) != 1L || is.na(diagnostics)) stop("`diagnostics` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(interpretation) || length(interpretation) != 1L || is.na(interpretation)) stop("`interpretation` must be TRUE or FALSE.", call. = FALSE)
if (!is.logical(plot) || length(plot) != 1L || is.na(plot)) stop("`plot` must be TRUE or FALSE.", call. = FALSE)
if (!is.list(plot_args)) stop("`plot_args` must be a named list.", call. = FALSE)
if (length(plot_args) && is.null(names(plot_args))) stop("`plot_args` must be a named list.", call. = FALSE)
if (!is.numeric(level) || length(level) != 1L || !is.finite(level) || level <= 0 || level >= 1) stop("`level` must be between 0 and 1.", call. = FALSE)
if (!adjust %in% stats::p.adjust.methods) stop("`adjust` must be one of: ", paste(stats::p.adjust.methods, collapse = ", "), ".", call. = FALSE)
adjusted_meta <- .r4vn_long_adjusted_meta(
substitute(adjusted),
missing(adjusted),
caller,
data = data
)
time_source <- .r4vn_long_time_source(
substitute(time),
data,
nrow(meta),
caller
)
exposure_expr <- substitute(exposure)
if (identical(time_source$mode, "wide")) {
prep <- .r4vn_long_prepare_wide(
data = data,
meta = meta,
time_labels = time_source$labels,
id_expr = substitute(id),
id_missing = missing(id),
by_expr = substitute(by),
by_missing = missing(by),
adjusted_meta = adjusted_meta,
ref = ref,
exposure_expr = exposure_expr,
exposure_missing = missing(exposure),
env = caller
)
input_format <- "wide"
} else {
prep <- .r4vn_long_prepare_long_common(
data = data,
time_spec = time_source$spec,
id_expr = substitute(id),
id_missing = missing(id),
by_expr = substitute(by),
by_missing = missing(by),
adjusted_meta = adjusted_meta,
ref = ref,
env = caller
)
input_format <- "long"
if (!missing(exposure) && !identical(exposure_expr, quote(NULL))) {
if (!is.symbol(exposure_expr)) stop("In long data, `exposure` must be one unquoted variable.", call. = FALSE)
exposure_name <- as.character(exposure_expr)
if (!exposure_name %in% names(prep$data)) stop("`exposure` variable was not found.", call. = FALSE)
prep$data$.exposure <- prep$data[[exposure_name]]
}
}
if (isTRUE(slope) && !isTRUE(prep$time_continuous)) {
stop("`slope = TRUE` is available only with continuous time, for example `time = c.month`.", call. = FALSE)
}
if (isTRUE(ar1) && !isTRUE(gee) && !isTRUE(rr) && !isTRUE(pr)) {
warning("`ar1 = TRUE` is used only with GEE. Set `gee = TRUE`, `rr = TRUE`, or `pr = TRUE`.", call. = FALSE)
}
design <- if (isTRUE(prep$repeated)) "longitudinal" else "repeated cross-sectional"
design_note <- if (identical(input_format, "wide")) {
"Wide data were converted to long format internally; the original data were not modified."
} else if (identical(design, "repeated cross-sectional")) {
"No subject ID was observed at more than one time point; observations across periods were treated as independent repeated cross-sectional samples."
} else {
"Repeated observations from the same subject were modeled using within-subject correlation."
}
if (identical(input_format, "long") && !isTRUE(prep$repeated) && is.null(prep$id_name)) {
design_note <- "No `id` was supplied; observations across time were treated as independent repeated cross-sectional samples."
}
results <- list()
all_rows <- list()
all_models <- list()
all_tests <- list()
all_contrasts <- list()
all_descriptive <- list()
all_diagnostics <- list()
all_plot_data <- list()
long_storage <- list()
notes <- character()
if (identical(input_format, "wide")) {
outcome_iterations <- 1L
} else {
outcome_iterations <- seq_len(nrow(meta))
}
for (ii in outcome_iterations) {
if (identical(input_format, "wide")) {
od <- prep$data
outcome_name <- prep$outcome_name
summary_type <- meta$type[1L]
original_variable <- paste(meta$variable, collapse = ", ")
} else {
outcome_name <- meta$variable[ii]
if (!outcome_name %in% names(prep$data)) stop("Outcome `", outcome_name, "` was not found in `data`.", call. = FALSE)
od <- prep$data
summary_type <- meta$type[ii]
original_variable <- outcome_name
if (summary_type %in% c("mean", "median", "full")) {
od$.outcome <- .r4vn_long_numeric(od[[outcome_name]], outcome_name)
} else {
od$.outcome <- od[[outcome_name]]
}
}
if (identical(input_format, "wide")) {
outcome_label <- prep$outcome_name
} else {
outcome_label <- .r4vn_long_label(data[[outcome_name]], outcome_name)
}
if (isTRUE(name) && !identical(outcome_label, original_variable)) {
outcome_label <- paste0(outcome_label, " (", original_variable, ")")
}
# Preserve the observed outcome exactly as supplied for descriptive
# summaries. The model copy below may be recoded (for example Yes/No -> 0/1).
od$.outcome_display <- od$.outcome
outcome_type <- NULL
event_level <- NULL
if (isTRUE(count)) {
numeric_outcome <- suppressWarnings(as.numeric(od$.outcome))
nonmissing <- numeric_outcome[is.finite(numeric_outcome)]
if (length(nonmissing) && (any(nonmissing < 0) || any(abs(nonmissing - round(nonmissing)) > 1e-8))) {
stop("`count = TRUE` requires non-negative integer outcomes.", call. = FALSE)
}
od$.outcome <- numeric_outcome
outcome_type <- "count"
effect_type <- "IRR"
if (".exposure" %in% names(od)) {
od$.exposure <- suppressWarnings(as.numeric(od$.exposure))
if (any(!is.na(od$.exposure) & (!is.finite(od$.exposure) | od$.exposure <= 0))) {
stop("`exposure` must be positive wherever it is observed.", call. = FALSE)
}
}
} else if (summary_type %in% c("mean", "median", "full")) {
od$.outcome <- suppressWarnings(as.numeric(od$.outcome))
outcome_type <- "continuous"
effect_type <- "Difference"
} else {
observed <- .r4vn_long_levels(od$.outcome)
if (length(observed) != 2L) {
stop(
"Outcome `", outcome_name, "` is categorical with ", length(observed),
" observed levels. `tablong()` currently supports continuous, binary, and count outcomes.",
call. = FALSE
)
}
event_level <- .r4vn_long_event(od$.outcome, event, outcome_name)
od$.outcome <- as.integer(as.character(od$.outcome) == as.character(event_level))
attr(od$.outcome, "event") <- event_level
outcome_type <- "binary"
effect_type <- if (isTRUE(rr)) "RR" else if (isTRUE(pr)) "PR" else "OR"
}
model_vars <- c(".outcome", ".time", ".id", prep$covariates)
if (!is.null(prep$by_name)) model_vars <- c(model_vars, ".by")
if (".exposure" %in% names(od)) model_vars <- c(model_vars, ".exposure")
model_vars <- unique(model_vars)
complete_base <- stats::complete.cases(od[, model_vars, drop = FALSE])
model_data <- od[complete_base, , drop = FALSE]
if (nrow(model_data) < 5L) stop("Too few complete observations for outcome `", outcome_name, "`.", call. = FALSE)
if (identical(outcome_type, "binary") && length(unique(model_data$.outcome)) < 2L) {
stop("Binary outcome `", outcome_name, "` has no variation after removing missing model data.", call. = FALSE)
}
fit <- .r4vn_long_fit_one(
data = model_data,
outcome_type = outcome_type,
effect_type = effect_type,
repeated = prep$repeated,
gee = gee,
ar1 = ar1,
slope = slope,
time_continuous = prep$time_continuous,
covariates = prep$covariates,
exposure = ".exposure" %in% names(model_data)
)
if (isTRUE(prep$time_continuous)) {
display_times <- prep$time_display
} else {
display_times <- prep$time_display
}
contrasts <- .r4vn_long_contrast_rows(
data = model_data,
fit = fit,
time_display = display_times,
time_reference = prep$time_reference,
by_display = prep$by_display,
time_continuous = prep$time_continuous,
change = change,
pairwise = pairwise,
level = level,
adjust = adjust
)
tests <- list(
p_time = fit$p_time,
p_group = fit$p_group,
p_interaction = fit$p_interaction
)
rows <- .r4vn_long_build_rows(
data = od,
outcome_label = outcome_label,
outcome_type = outcome_type,
summary_type = summary_type,
event = event_level,
effect_type = effect_type,
time_display = display_times,
time_reference = prep$time_reference,
by_display = prep$by_display,
time_continuous = prep$time_continuous,
contrasts = contrasts,
tests = tests,
change = change,
digits = digit,
effect_digits = effect_digit,
p_digits = p_digit,
show_n = missing
)
descriptive_rows <- .r4vn_long_descriptive_rows(
data = od, outcome_label = outcome_label, outcome_type = outcome_type,
summary_type = summary_type, event = event_level,
time_display = display_times, by_display = prep$by_display,
time_continuous = prep$time_continuous, digits = digit, show_n = missing
)
diagnostic_row <- .r4vn_long_diagnostics_one(
fit = fit, model_data = model_data, outcome_label = outcome_label,
outcome_type = outcome_type, time_display = display_times,
time_continuous = prep$time_continuous, repeated = prep$repeated
)
plot_rows <- .r4vn_long_plot_rows(
data = od, outcome_label = outcome_label, outcome_type = outcome_type,
event = event_level, time_display = display_times,
by_display = prep$by_display, time_continuous = prep$time_continuous,
level = level
)
if (isTRUE(fit$singular)) {
notes <- c(notes, paste0(outcome_label, ": the mixed model produced a singular random-effects fit; inspect `$models` before publication."))
}
if (!is.null(fit$package_note) && length(fit$package_note) && nzchar(fit$package_note[1L])) {
notes <- c(notes, paste0(outcome_label, ": ", fit$package_note[1L]))
}
if (summary_type %in% c("median", "full") && identical(outcome_type, "continuous")) {
notes <- c(
notes,
paste0(
outcome_label,
": `", if (identical(summary_type, "median")) "q." else "f.",
"` changes the observed descriptive summary only; inferential effects are mean-model effects."
)
)
}
key <- make.unique(c(names(all_models), outcome_name))[length(all_models) + 1L]
all_rows[[key]] <- rows
all_models[[key]] <- fit$fit
all_tests[[key]] <- tests
all_contrasts[[key]] <- contrasts
all_descriptive[[key]] <- descriptive_rows
all_diagnostics[[key]] <- diagnostic_row
all_plot_data[[key]] <- plot_rows
long_storage[[key]] <- od
results[[key]] <- list(
outcome = outcome_name,
label = outcome_label,
type = outcome_type,
event = event_level,
effect = effect_type,
model_engine = fit$engine,
model = fit$fit,
tests = tests,
contrasts = contrasts,
data = rows
)
}
# Multiple outcomes may use different effect scales (for example a mean
# difference for SBP and an OR for a binary outcome). Standardize only the
# effect-column heading when needed so all outcome blocks can share one table.
effect_columns <- unique(unlist(lapply(all_rows, function(z) {
grep("\\(95% CI\\)$", names(z), value = TRUE)
}), use.names = FALSE))
if (length(effect_columns) > 1L) {
all_rows <- lapply(all_rows, function(z) {
hit <- grep("\\(95% CI\\)$", names(z), value = TRUE)
if (length(hit) == 1L) names(z)[names(z) == hit] <- "Model effect (95% CI)"
z
})
scale_notes <- vapply(results, function(z) {
scale <- if (identical(z$type, "continuous")) "mean difference" else z$effect
paste0(z$label, ": model effect is ", scale, ".")
}, character(1))
notes <- c(notes, scale_notes)
}
flat <- do.call(rbind, all_rows)
rownames(flat) <- NULL
descriptive_table <- if (length(all_descriptive)) do.call(rbind, all_descriptive) else data.frame()
if (nrow(descriptive_table)) rownames(descriptive_table) <- NULL
diagnostics_table <- if (length(all_diagnostics)) do.call(rbind, all_diagnostics) else data.frame()
if (nrow(diagnostics_table)) rownames(diagnostics_table) <- NULL
plot_data <- if (length(all_plot_data)) do.call(rbind, all_plot_data) else data.frame()
if (nrow(plot_data)) rownames(plot_data) <- NULL
tests_table <- .r4vn_long_tests_table(results, digits = p_digit)
contrasts_table <- .r4vn_long_contrasts_table(
results, effect_digits = effect_digit, p_digits = p_digit
)
interpretation_table <- if (isTRUE(interpretation)) {
.r4vn_long_interpretation(results, alpha = 0.05, p_digits = p_digit)
} else NULL
graph <- if (isTRUE(plot)) {
.r4vn_long_build_plot(plot_data, ci = TRUE, plot_args = plot_args)
} else NULL
base_notes <- c(
design_note,
if (identical(design, "longitudinal") && any(vapply(results, function(z) identical(z$model_engine, "mixed"), logical(1)))) {
"Continuous repeated outcomes used a random-intercept linear mixed model from R's recommended `nlme` package; `slope = TRUE` adds a random linear time slope when time is continuous."
} else NULL,
if (any(vapply(results, function(z) identical(z$model_engine, "cluster_robust"), logical(1)))) {
"Repeated binary/count outcomes and marginal models used working-independence regression with subject-clustered robust sandwich standard errors calculated internally by R4VN."
} else NULL,
if (any(vapply(results, function(z) identical(z$model_engine, "gee"), logical(1)))) {
"AR(1) GEE used robust sandwich standard errors through the optional `geepack` package."
} else NULL,
if (any(vapply(results, function(z) z$effect %in% c("RR", "PR") && identical(z$model_engine, "independent"), logical(1)))) {
"RR/PR in independent samples used modified Poisson regression with robust sandwich standard errors calculated internally by R4VN."
} else NULL,
if (!is.null(prep$by_name)) {
"Time and group p-values are omnibus main-effect tests from the additive model; the time x group p-value tests the interaction in the full model."
} else {
"The overall time p-value tests whether the time effect is zero."
},
if (isTRUE(change) && !isTRUE(prep$time_continuous)) {
paste0("Change contrasts use ", prep$time_reference, " as the reference time.")
} else NULL,
if (!identical(adjust, "none")) {
paste0("Displayed contrast p-values use ", adjust, " multiplicity adjustment; confidence intervals are unadjusted.")
} else NULL
)
notes <- unique(c(base_notes, notes))
notes <- notes[nzchar(notes)]
if (is.null(title) || !length(title) || is.na(title[1L]) || !nzchar(as.character(title[1L]))) {
title <- "Longitudinal / repeated-measures analysis"
} else {
title <- as.character(title[1L])
}
table_html <- .r4vn_long_html_table(flat, title, notes, bold_p, p_bold)
plot_html <- if (isTRUE(plot) && !is.null(graph)) .r4vn_long_plot_html(graph) else ""
secondary_html <- paste0(
plot_html,
.r4vn_long_html_simple_table(tests_table, "Omnibus tests"),
if (isTRUE(diagnostics)) .r4vn_long_html_simple_table(diagnostics_table, "Model diagnostics") else "",
if (isTRUE(interpretation)) .r4vn_long_html_simple_table(interpretation_table, "Interpretation") else ""
)
html <- .r4vn_long_html_document(paste0(table_html, secondary_html))
if (is.null(file) || !length(file) || is.na(file[1L]) || !nzchar(as.character(file[1L]))) {
file <- tempfile(pattern = "r4vn-tablong-", fileext = ".html")
} else {
file <- path.expand(as.character(file[1L]))
if (!grepl("\\.html?$", file, ignore.case = TRUE)) file <- paste0(file, ".html")
dir.create(dirname(file), recursive = TRUE, showWarnings = FALSE)
}
writeLines(enc2utf8(html), file, useBytes = TRUE)
file <- normalizePath(file, winslash = "/", mustWork = TRUE)
if (isTRUE(show)) .r4vn_long_show_html(file)
if (isTRUE(plot) && !is.null(graph) && interactive()) {
.r4vn_long_draw_plot(graph$data, ci = graph$ci, plot_args = graph$args)
}
combined_long <- if (length(long_storage) == 1L) {
long_storage[[1L]]
} else {
do.call(rbind, lapply(names(long_storage), function(nm) {
z <- long_storage[[nm]]
z$.r4vn_outcome_name <- nm
z
}))
}
rownames(combined_long) <- NULL
tables <- list(
`Main table` = flat,
Descriptive = descriptive_table,
`Omnibus tests` = tests_table,
Contrasts = contrasts_table,
Diagnostics = diagnostics_table
)
if (isTRUE(interpretation)) tables$Interpretation <- interpretation_table
result <- list(
data = flat,
descriptive = descriptive_table,
estimates = list(contrasts = contrasts_table),
table_html = table_html,
html = html,
file = file,
models = all_models,
tests = all_tests,
tests_table = tests_table,
contrasts = all_contrasts,
contrasts_table = contrasts_table,
diagnostics = diagnostics_table,
interpretation = interpretation_table,
tables = tables,
plot_data = plot_data,
graph = graph,
plots = list(trajectory = graph),
long_data = combined_long,
results = results,
input_format = input_format,
design = design,
time = list(
variable = prep$time_name,
continuous = prep$time_continuous,
levels = prep$time_display,
reference = prep$time_reference
),
by = prep$by_name,
id = prep$id_name,
adjusted = adjusted_meta,
metadata = list(
input_format = input_format, design = design,
repeated = isTRUE(prep$repeated), time_variable = prep$time_name,
time_continuous = isTRUE(prep$time_continuous),
time_levels = prep$time_display, time_reference = prep$time_reference,
group_variable = prep$by_name, id_variable = prep$id_name,
change = isTRUE(change), pairwise = isTRUE(pairwise), adjust = adjust,
confidence_level = level, interpretation = isTRUE(interpretation),
plot = isTRUE(plot)
),
call = match.call()
)
class(result) <- c("r4vn_tablong", "r4vn_tab", "list")
invisible(result)
}
#' Print an R4VN Longitudinal Table
#'
#' @param x Object returned by `tablong()`.
#' @param ... Additional arguments currently ignored.
#' @return The input object, invisibly.
#' @keywords internal
#' @method print r4vn_tablong
#' @export
print.r4vn_tablong <- function(x, ...) {
print(x$data, row.names = FALSE)
invisible(x)
}
#' Plot an R4VN Longitudinal Analysis
#'
#' Recreates the observed longitudinal profile stored by `tablong()`. This is
#' useful when the original analysis used `plot = FALSE` or when different plot
#' labels/sizing are wanted without refitting the statistical model.
#'
#' @param x Object returned by `tablong()`.
#' @param ci Show 95% confidence intervals. Default `TRUE`.
#' @param ... Named plot options accepted through `plot_args`, including
#' `title`, `xlab`, `ylab`, `line_width`, `point_size`, `base_size`, and
#' `legend_position`.
#' @return An R4VN plot specification, invisibly; the plot is drawn with base R graphics.
#' @method plot r4vn_tablong
#' @export
plot.r4vn_tablong <- function(x, ci = TRUE, ...) {
args <- list(...)
g <- .r4vn_long_build_plot(x$plot_data, ci = ci, plot_args = args)
if (!is.null(g)) .r4vn_long_draw_plot(g$data, ci = g$ci, plot_args = g$args)
invisible(g)
}
attr(tablong, "r4vn_version") <- "tablong-1.2.0-2026-09-02"
Try the R4VN package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.