# ============================================================================= # 01_simulate_coa_data.R # The Psychometrician's Handbook, Chapter 1 (Correlations) # Simulated COA validation dataset, per Ch1_Dataset_Spec_v0.3 # # Generates `adcorr`: one long CDISC ADaM BDS data frame, 18 parameters, # 52 rows per subject, 10,400 rows at the canonical n = 200. # # Every relationship below comes from a generative model with a stated # population or latent truth. Sample estimates land near those truths and # vary around them. No correlation is reverse-engineered to a target. # ============================================================================= pkgs <- c("dplyr","tidyr","ggplot2","psych","rmcorr","lme4","nlme", "corrr","corrplot","blandr","polycor","mirt","MASS") for (p in pkgs) if (!requireNamespace(p, quietly = TRUE)) install.packages(p, repos = "https://cloud.r-project.org") suppressPackageStartupMessages({ library(dplyr); library(tidyr); library(MASS); library(mirt) }) # --- Locked truths (spec 5, plus Addendum A2) -------------------------------- STUDYID_ <- "PHB-CH1" SEED <- 2026 N_CANON <- 200 # 5a. Population correlations for the four continuous scales POP_R <- matrix(c(1.00, 0.72, 0.65, 0.18, 0.72, 1.00, 0.58, 0.22, 0.65, 0.58, 1.00, 0.20, 0.18, 0.22, 0.20, 1.00), 4, 4, dimnames = list(c("TARGTOT","RELATEDA","RELATEDB","UNRELATD"), c("TARGTOT","RELATEDA","RELATEDB","UNRELATD"))) # 5b + A2. Ordinal items: mirt GRM on a moderately skewed latent trait. # Loadings give latent r(i,j) = lambda_i * lambda_j, averaging about 0.73. LAMBDA <- c(0.84, 0.87, 0.82, 0.86) THRESH <- list(c(-1.4,-0.4,0.5,1.5), # 4 thresholds -> 5 categories, c(-1.6,-0.5,0.4,1.4), # mildly asymmetric as symptom items run c(-1.2,-0.3,0.6,1.6), c(-1.5,-0.4,0.5,1.5)) SKEW_SHAPE <- 8 # gamma shape; skew = 2/sqrt(8) ~ 0.71 # mirt works in the logistic metric. Convert the factor-model parameters: # a_normal = lambda / sqrt(1 - lambda^2), a_mirt = 1.702 * a_normal # threshold tau on y* maps to IRT b = tau / lambda, so d = -a_mirt * b D_CONST <- 1.702 A_MIRT <- D_CONST * LAMBDA / sqrt(1 - LAMBDA^2) D_MIRT <- t(sapply(seq_along(LAMBDA), function(i) -A_MIRT[i] * (THRESH[[i]] / LAMBDA[i]))) # Calibrated generative constants (see spec 5c to 5g) TUNE <- list( sev_p = 0.35, # SEVGRP prevalence, genuine dichotomy sev_rpb = 0.40, # 5d target point-biserial with TARGTOT dich_rho = 0.50, # 5d latent severity correlation, recovered by biserial dich_p = 0.35, # DICHOTOM prevalence at the clinical cutoff fat_s1 = 1.0, # 5c FATVAS beta shapes -> right skew ~1.35 fat_s2 = 7.0, walk_rho = 0.955, # 5c latent fatigue to walk, before the skew transform pgis_rho = 0.635, # 5e -> Spearman ~0.60 with TARGTOT after a 5-cat cut pgic_rho = 0.470, # 5e -> r ~0.45 with true change after a 7-cat cut retest_r = 0.85, # 5e ICC target chg_mu = -0.35, # mean improvement at Week 12, in SD units chg_sd = 0.55, b_var = 0.70, # 5f/5g between-person variance share (ICC ~0.7) w_var = 0.30, fin_rb = 0.71, # 5f between-person fin_rw = 0.24, # 5f within-person dia_rb = 0.00, # 5g between-person, deliberately null dia_rw = 0.55, # 5g within-person, the hidden relationship dia_ar1 = 0.40 # 5g AR(1) on within-person residuals ) std <- function(x) (x - mean(x)) / stats::sd(x) cut_k <- function(x, k) as.integer(cut(x, breaks = c(-Inf, qnorm(seq_len(k - 1) / k), Inf), labels = FALSE)) # AR(1) series with correlated innovations: contemporaneous correlation of the # two series equals the innovation correlation when both share the same phi. ar1_pair <- function(n_t, phi, rho_w, sd_w) { s_innov <- sd_w * sqrt(1 - phi^2) E <- MASS::mvrnorm(n_t, c(0, 0), matrix(c(1, rho_w, rho_w, 1), 2) * s_innov^2) out <- matrix(0, n_t, 2) out[1, ] <- E[1, ] / sqrt(1 - phi^2) for (t in 2:n_t) out[t, ] <- phi * out[t - 1, ] + E[t, ] out } #' Simulate the Chapter 1 COA validation dataset #' #' @param n number of subjects #' @param seed optional seed; the canonical dataset uses 2026 #' @return a long ADaM BDS data frame simulate_coa_data <- function(n = N_CANON, seed = NULL, tune = TUNE) { if (!is.null(seed)) set.seed(seed) # -- 5a/5d. Four continuous scales, with SEVGRP built into the target ------ # SEVGRP is a genuine dichotomy, so it is drawn first and the target PRO is # constructed to carry a real group difference. The Cholesky factor keeps # every population correlation in POP_R exact. SEVGRP <- rbinom(n, 1, tune$sev_p) g <- std(SEVGRP) v1 <- tune$sev_rpb * g + sqrt(1 - tune$sev_rpb^2) * rnorm(n) v <- rbind(v1, matrix(rnorm(3 * n), nrow = 3)) z <- t(chol(POP_R)) %*% v TARGTOT_z <- as.numeric(z[1, ]); RELA_z <- as.numeric(z[2, ]) RELB_z <- as.numeric(z[3, ]); UNRE_z <- as.numeric(z[4, ]) # DICHOTOM: a latent severity score cut at a clinical cutoff, not a median # split. Biserial recovers dich_rho; point-biserial is attenuated. sev_lat <- tune$dich_rho * TARGTOT_z + sqrt(1 - tune$dich_rho^2) * rnorm(n) DICHOTOM <- as.integer(sev_lat > qnorm(1 - tune$dich_p)) # -- 5b. Ordinal items from a mirt GRM on a skewed latent ------------------ theta <- { gg <- rgamma(n, SKEW_SHAPE, 1); (gg - SKEW_SHAPE) / sqrt(SKEW_SHAPE) } items <- mirt::simdata(a = matrix(A_MIRT, ncol = 1), d = D_MIRT, itemtype = "graded", Theta = matrix(theta, ncol = 1)) + 1L colnames(items) <- c("ITEM01","ITEM02","ITEM03","ITEM04") # -- 5c. Fatigue and walk: strong, monotonic, curved ----------------------- # WALK6MIN is linear in latent fatigue, so the curvature the Pearson # coefficient trips over comes from the skewed FATVAS marginal. fat_lat <- rnorm(n) FATVAS <- 100 * qbeta(pnorm(fat_lat), tune$fat_s1, tune$fat_s2) WALK6MIN <- 400 - 110 * (tune$walk_rho * fat_lat + sqrt(1 - tune$walk_rho^2) * rnorm(n)) WALK6MIN <- pmax(WALK6MIN, 60) # -- 5e. Anchors, retest and change ---------------------------------------- PGIS <- cut_k(tune$pgis_rho * TARGTOT_z + sqrt(1 - tune$pgis_rho^2) * rnorm(n), 5) RETEST_z <- tune$retest_r * TARGTOT_z + sqrt(1 - tune$retest_r^2) * rnorm(n) chg_true <- tune$chg_mu + tune$chg_sd * rnorm(n) W12_z <- TARGTOT_z + chg_true PGIC <- cut_k(tune$pgic_rho * std(chg_true) + sqrt(1 - tune$pgic_rho^2) * rnorm(n), 7) # -- 5f. Monthly pair: the relationship lives between persons -------------- Sb <- matrix(c(1, tune$fin_rb, tune$fin_rb, 1), 2) * tune$b_var Sw <- matrix(c(1, tune$fin_rw, tune$fin_rw, 1), 2) * tune$w_var b_fin <- MASS::mvrnorm(n, c(0, 0), Sb) fin_long <- do.call(rbind, lapply(1:4, function(m) { w <- MASS::mvrnorm(n, c(0, 0), Sw) data.frame(SUBJID = 1:n, AVISITN = m, FINTOX = b_fin[, 1] + w[, 1], SYMBUR = b_fin[, 2] + w[, 2]) })) # -- 5g. Daily diary: the relationship lives within persons ---------------- Sb2 <- matrix(c(1, tune$dia_rb, tune$dia_rb, 1), 2) * tune$b_var b_dia <- MASS::mvrnorm(n, c(0, 0), Sb2) dia_long <- do.call(rbind, lapply(1:n, function(i) { w <- ar1_pair(14, tune$dia_ar1, tune$dia_rw, sqrt(tune$w_var)) data.frame(SUBJID = i, AVISITN = 1:14, DIARY1 = b_dia[i, 1] + w[, 1], DIARY2 = b_dia[i, 2] + w[, 2]) })) # -- Assemble long BDS ------------------------------------------------------ sid <- sprintf("%s-%03d", STUDYID_, 1:n) to_score <- function(z_, mean_, sd_) round(pmin(pmax(mean_ + sd_ * z_, 0), 100)) blk <- function(paramcd, avisit, avisitn, aval) data.frame(SUBJID = 1:n, PARAMCD = paramcd, AVISIT = avisit, AVISITN = avisitn, AVAL = aval) rows <- rbind( blk("TARGTOT", "Baseline", 0, to_score(TARGTOT_z, 50, 15)), blk("TARGTOT", "Retest", 1, to_score(RETEST_z, 50, 15)), blk("TARGTOT", "Week 12", 12, to_score(W12_z, 50, 15)), blk("RELATEDA", "Baseline", 0, to_score(RELA_z, 50, 12)), blk("RELATEDB", "Baseline", 0, to_score(RELB_z, 50, 12)), blk("UNRELATD", "Baseline", 0, to_score(UNRE_z, 50, 12)), blk("ITEM01", "Baseline", 0, items[, 1]), blk("ITEM02", "Baseline", 0, items[, 2]), blk("ITEM03", "Baseline", 0, items[, 3]), blk("ITEM04", "Baseline", 0, items[, 4]), blk("FATVAS", "Baseline", 0, round(FATVAS)), blk("WALK6MIN", "Baseline", 0, round(WALK6MIN)), blk("SEVGRP", "Baseline", 0, SEVGRP), blk("DICHOTOM", "Baseline", 0, DICHOTOM), blk("PGIS", "Baseline", 0, PGIS), blk("PGIC", "Week 12", 12, PGIC), do.call(rbind, lapply(1:4, function(m) { s <- fin_long[fin_long$AVISITN == m, ] rbind(data.frame(SUBJID = s$SUBJID, PARAMCD = "FINTOX", AVISIT = paste("Month", m), AVISITN = m, AVAL = round(pmin(pmax(45 + 14 * s$FINTOX, 0), 100), 1)), data.frame(SUBJID = s$SUBJID, PARAMCD = "SYMBUR", AVISIT = paste("Month", m), AVISITN = m, AVAL = round(pmin(pmax(52 + 14 * s$SYMBUR, 0), 100), 1))) })), do.call(rbind, lapply(1:14, function(dy) { s <- dia_long[dia_long$AVISITN == dy, ] rbind(data.frame(SUBJID = s$SUBJID, PARAMCD = "DIARY1", AVISIT = paste("Day", dy), AVISITN = dy, AVAL = round(pmin(pmax(4.8 + 1.6 * s$DIARY1, 0), 10), 1)), data.frame(SUBJID = s$SUBJID, PARAMCD = "DIARY2", AVISIT = paste("Day", dy), AVISITN = dy, AVAL = round(pmin(pmax(5.1 + 1.6 * s$DIARY2, 0), 10), 1))) })) ) meta <- data.frame( PARAMCD = c("TARGTOT","RELATEDA","RELATEDB","UNRELATD","ITEM01","ITEM02", "ITEM03","ITEM04","FATVAS","WALK6MIN","SEVGRP","DICHOTOM", "PGIS","PGIC","FINTOX","SYMBUR","DIARY1","DIARY2"), PARAMN = 1:18, PARAM = c( "Target PRO - total score (continuous)", "Related measure A - total score (convergent, continuous)", "Related measure B - total score (convergent, continuous)", "Unrelated measure - total score (divergent, continuous)", "Symptom item 1 - ordinal 5-point Likert", "Symptom item 2 - ordinal 5-point Likert", "Symptom item 3 - ordinal 5-point Likert", "Symptom item 4 - ordinal 5-point Likert", "Fatigue VAS - right-skewed continuous (0-100)", "Six-minute walk distance - approximately normal continuous (m)", "Clinical status - genuine dichotomy (0=absent, 1=present)", "Disease present by score cutoff - dichotomised continuum", "Patient Global Impression of Severity - ordinal anchor (5-point)", "Patient Global Impression of Change - ordinal anchor (7-point)", "Financial toxicity - repeated continuous (monthly)", "Symptom burden - repeated continuous (monthly)", "Daily symptom diary assessment 1 - repeated continuous (daily)", "Daily symptom diary assessment 2 - repeated continuous (daily)"), stringsAsFactors = FALSE) adcorr <- rows %>% left_join(meta, by = "PARAMCD") %>% mutate(STUDYID = STUDYID_, USUBJID = sid[SUBJID]) %>% group_by(PARAMCD, SUBJID) %>% arrange(AVISITN, .by_group = TRUE) %>% # BASE is the first visit of any parameter measured more than once mutate(BASE = if (n() > 1) first(AVAL) else NA_real_, CHG = if (n() > 1) AVAL - BASE else NA_real_) %>% ungroup() %>% dplyr::select(STUDYID, USUBJID, SUBJID, PARAMCD, PARAM, PARAMN, AVISIT, AVISITN, AVAL, BASE, CHG) %>% arrange(SUBJID, PARAMN, AVISITN) as.data.frame(adcorr) } #' Replication wrapper, for the sampling-variability figures only (spec 7) replicate_estimates <- function(n_grid = c(25, 50, 100, 200), R = 500, seed = SEED) { set.seed(seed) do.call(rbind, lapply(n_grid, function(nn) { do.call(rbind, lapply(seq_len(R), function(rep_i) { d <- simulate_coa_data(n = nn) w <- d %>% filter(AVISIT == "Baseline", PARAMCD %in% c("TARGTOT","RELATEDA")) %>% dplyr::select(USUBJID, PARAMCD, AVAL) %>% pivot_wider(names_from = PARAMCD, values_from = AVAL) ct <- cor.test(w$TARGTOT, w$RELATEDA) data.frame(n = nn, replicate = rep_i, r = unname(ct$estimate), lcl = ct$conf.int[1], ucl = ct$conf.int[2]) })) })) } # --- Canonical dataset ------------------------------------------------------- set.seed(SEED) adcorr <- simulate_coa_data(n = N_CANON) # ============================================================================= # EXPORT # ============================================================================= OUTDIR <- Sys.getenv("CH1_OUT", ".") write.csv(adcorr, file.path(OUTDIR, "adcorr.csv"), row.names = FALSE, na = "") # ============================================================================= # DECK CODE BLOCKS # Each header names the slide in Chapter 1 v0.204 that the block appears on. # On-slide snippets are tight excerpts of the code below. # ============================================================================= suppressPackageStartupMessages({ library(psych); library(polycor); library(rmcorr) library(nlme); library(lme4); library(corrr); library(blandr); library(ggplot2) }) say <- function(...) cat(sprintf(...), "\n") # CDISC BDS is long. Most blocks need one row per subject, so pivot first. wide <- adcorr %>% filter(AVISIT == "Baseline") %>% pivot_wider(id_cols = USUBJID, names_from = PARAMCD, values_from = AVAL) cat("\n########## SLIDE 23 cross-sectional example ##########\n") ct <- cor.test(wide$TARGTOT, wide$RELATEDA) say("r = %.3f 95%% CI [%.3f, %.3f] p = %.2g", ct$estimate, ct$conf.int[1], ct$conf.int[2], ct$p.value) cat("\n########## SLIDE 27 change-score example ##########\n") chg <- adcorr %>% filter(PARAMCD == "TARGTOT", AVISIT == "Week 12") %>% dplyr::select(USUBJID, CHG) anc <- adcorr %>% filter(PARAMCD == "PGIC") %>% dplyr::select(USUBJID, PGIC = AVAL) d27 <- inner_join(chg, anc, by = "USUBJID") ct27 <- cor.test(d27$CHG, d27$PGIC) say("r(change, anchor) = %.3f 95%% CI [%.3f, %.3f]", ct27$estimate, ct27$conf.int[1], ct27$conf.int[2]) cat("\n########## SLIDE 34 person-mean centering ##########\n") d34 <- adcorr %>% filter(PARAMCD %in% c("FINTOX","SYMBUR")) %>% pivot_wider(id_cols = c(USUBJID, AVISITN), names_from = PARAMCD, values_from = AVAL) %>% group_by(USUBJID) %>% mutate(fin_w = FINTOX - mean(FINTOX), sym_w = SYMBUR - mean(SYMBUR)) %>% ungroup() pm34 <- d34 %>% group_by(USUBJID) %>% summarise(f = mean(FINTOX), s = mean(SYMBUR), .groups = "drop") say("naive r = %.3f", cor(d34$FINTOX, d34$SYMBUR)) say("within r = %.3f", cor(d34$fin_w, d34$sym_w)) say("between r = %.3f (person means)", cor(pm34$f, pm34$s)) cat("\n########## SLIDE 35 rmcorr ##########\n") fit35 <- rmcorr(participant = USUBJID, measure1 = FINTOX, measure2 = SYMBUR, dataset = as.data.frame(d34)) say("rmcorr r = %.3f 95%% CI [%.3f, %.3f] df = %d", fit35$r, fit35$CI[1], fit35$CI[2], fit35$df) cat("\n########## SLIDE 36 mixed effects with AR(1) ##########\n") d36 <- adcorr %>% filter(PARAMCD %in% c("DIARY1","DIARY2")) %>% pivot_wider(id_cols = c(USUBJID, AVISITN), names_from = PARAMCD, values_from = AVAL) %>% as.data.frame() say("naive pooled r = %.3f", cor(d36$DIARY1, d36$DIARY2)) pm36 <- d36 %>% group_by(USUBJID) %>% summarise(a = mean(DIARY1), b = mean(DIARY2), .groups = "drop") say("person-mean r = %.3f (between persons)", cor(pm36$a, pm36$b)) m36 <- lme(DIARY2 ~ DIARY1, random = ~ 1 | USUBJID, correlation = corAR1(form = ~ AVISITN | USUBJID), data = d36) # within-person correlation from the variance components d36c <- d36 %>% group_by(USUBJID) %>% mutate(x = DIARY1 - mean(DIARY1), y = DIARY2 - mean(DIARY2)) %>% ungroup() say("within-person r = %.3f (centred)", cor(d36c$x, d36c$y)) say("AR(1) phi estimate = %.3f", coef(m36$modelStruct$corStruct, unconstrained = FALSE)) cat("\n########## SLIDE 44 Pearson ##########\n") say("Pearson r = %.3f", cor(wide$TARGTOT, wide$RELATEDA)) cat("\n########## SLIDE 45 Spearman and Kendall ##########\n") say("Spearman rho = %.3f", cor(wide$ITEM02, wide$ITEM03, method = "spearman")) say("Kendall tau = %.3f", cor(wide$ITEM02, wide$ITEM03, method = "kendall")) say("Pearson r = %.3f (for contrast)", cor(wide$ITEM02, wide$ITEM03)) cat("\n########## SLIDE 46 polyserial ##########\n") say("polyserial = %.3f (PGIS ordinal x TARGTOT continuous)", polyserial(wide$TARGTOT, wide$PGIS)) say("Pearson = %.3f (for contrast)", cor(wide$TARGTOT, wide$PGIS)) cat("\n########## SLIDE 47 polychoric ##########\n") say("polychoric = %.3f", polychor(wide$ITEM02, wide$ITEM03)) pc <- psych::polychoric(as.data.frame(wide[, c("ITEM02","ITEM03")]), progress = FALSE)$rho say("psych::polychoric = %.3f", pc[1,2]) say("latent truth = %.3f", LAMBDA[2] * LAMBDA[3]) cat("\n########## SLIDE 48 point-biserial and biserial ##########\n") say("point-biserial, SEVGRP (true dichotomy) = %.3f", cor(wide$TARGTOT, wide$SEVGRP)) say("point-biserial, DICHOTOM (cut continuum) = %.3f", cor(wide$TARGTOT, wide$DICHOTOM)) say("biserial, DICHOTOM = %.3f", polyserial(wide$TARGTOT, wide$DICHOTOM)) cat("\n########## SLIDE 49 ICC ##########\n") rt <- adcorr %>% filter(PARAMCD == "TARGTOT", AVISIT %in% c("Baseline","Retest")) %>% pivot_wider(id_cols = USUBJID, names_from = AVISIT, values_from = AVAL) icc <- psych::ICC(as.data.frame(rt[, c("Baseline","Retest")]))$results say("ICC(2,1) = %.3f [%.3f, %.3f]", icc$ICC[2], icc$`lower bound`[2], icc$`upper bound`[2]) say("ICC(3,1) = %.3f [%.3f, %.3f]", icc$ICC[3], icc$`lower bound`[3], icc$`upper bound`[3]) cat("\n########## SLIDE 57 scatter ##########\n") say("plotted TARGTOT vs RELATEDA, n = %d", nrow(wide)) cat("\n########## SLIDE 58 correlation matrix ##########\n") m58 <- wide %>% dplyr::select(TARGTOT, RELATEDA, RELATEDB, UNRELATD) print(round(cor(m58, use = "pairwise.complete.obs"), 2)) cat("\n########## SLIDE 59 Bland-Altman ##########\n") ba <- blandr.statistics(rt$Baseline, rt$Retest) say("bias = %.2f LoA [%.2f, %.2f]", ba$bias, ba$lowerLOA, ba$upperLOA) cat("\n########## SLIDE 60 spaghetti ##########\n") sp60 <- adcorr %>% filter(PARAMCD == "SYMBUR") say("SYMBUR trajectories: %d subjects x %d months", dplyr::n_distinct(sp60$USUBJID), dplyr::n_distinct(sp60$AVISITN)) say("mean by month: %s", paste(sprintf("%.1f", tapply(sp60$AVAL, sp60$AVISITN, mean)), collapse = ", "))