การวัดซ้ำ: เลือกวิธีตามชนิดของผลลัพธ์ และตามว่าต้องการคำตอบระดับประชากรหรือระดับผู้ป่วย

* Simulated longitudinal trial: a symptom score (lower is better) at weeks 0, 4, 8, 12 and a clinician-rated
* responder (0/1) at weeks 4, 8, 12, in 300 participants randomised 1:1 to treatment or control.
* Simulated data: not evidence about any real drug or patient.
* The simulated trial file (not published) is read from a folder two levels above this script.
version 18
clear all
set more off
set linesize 120

* Record whether each command used below is available (all are built into Stata).
foreach c in mixed contrast lincom testparm xtgee melogit anova {
    capture which `c'
    if _rc == 0 {
        display "VERIFY `c' available"
    }
    else {
        display "VERIFY `c' missing"
    }
}

* Load the simulated trial file (long format: one row per participant-visit).
import delimited using ../../datasets/W4/W4.csv, clear asdouble
label define armlab 0 "control" 1 "treatment"
label values arm armlab

* Sample sizes.
quietly count
display "CANON w4.n_obs " r(N)
egen byte first = tag(id)
quietly count if first
display "CANON w4.n " r(N)
quietly count if first & arm == 0
display "CANON w4.n.ctrl " r(N)
quietly count if first & arm == 1
display "CANON w4.n.trt " r(N)
quietly count if !missing(responder)
display "CANON w4.n_obs_resp " r(N)

* ---------------------------------------------------------------- 2 x T means table
* Observed mean score in each arm at each visit (the cells a saturated mean model reproduces).
table arm week, statistic(mean score) statistic(sd score) nformat(%7.2f)
forvalues a = 0/1 {
    local g = cond(`a' == 1, "trt", "ctrl")
    forvalues v = 0/3 {
        quietly summarize score if arm == `a' & visit == `v'
        display "CANON w4.mean.`g'.wk" 4*`v' " " strtrim(string(r(mean), "%12.4f"))
    }
}

* ---------------------------------------------------------------- mixed model, arm x visit
* Categorical visit, arm-by-visit interaction, random intercept and random slope on week (REML).
mixed score i.arm##i.visit || id: week, covariance(unstructured) reml
* beta1 = 1.arm: the arm difference at week 0 (zero in truth, because arms are randomised).
lincom 1.arm
display "CANON w4.lmm.b1 " strtrim(string(r(estimate), "%12.4f"))
display "CANON w4.lmm.b1.se " strtrim(string(r(se), "%12.4f"))
display "CANON w4.lmm.b1.lo " strtrim(string(r(lb), "%12.4f"))
display "CANON w4.lmm.b1.hi " strtrim(string(r(ub), "%12.4f"))
* beta3 at each visit = 1.arm#v.visit: the difference in differences from week 0.
display "CANON w4.lmm.b3.wk4 " strtrim(string(_b[1.arm#1.visit], "%12.4f"))
display "CANON w4.lmm.b3.wk8 " strtrim(string(_b[1.arm#2.visit], "%12.4f"))
lincom 1.arm#3.visit
display "CANON w4.lmm.b3.wk12 " strtrim(string(r(estimate), "%12.4f"))
display "CANON w4.lmm.b3.wk12.se " strtrim(string(r(se), "%12.4f"))
display "CANON w4.lmm.b3.wk12.lo " strtrim(string(r(lb), "%12.4f"))
display "CANON w4.lmm.b3.wk12.hi " strtrim(string(r(ub), "%12.4f"))
* Arm contrast at each visit (beta1 + beta3 at that visit): what a reader usually wants.
contrast r.arm@visit, effects
forvalues v = 0/3 {
    if `v' == 0 {
        quietly lincom 1.arm
    }
    else {
        quietly lincom 1.arm + 1.arm#`v'.visit
    }
    display "CANON w4.lmm.diff.wk" 4*`v' " " strtrim(string(r(estimate), "%12.4f"))
    display "CANON w4.lmm.diff.wk" 4*`v' ".se " strtrim(string(r(se), "%12.4f"))
    display "CANON w4.lmm.diff.wk" 4*`v' ".lo " strtrim(string(r(lb), "%12.4f"))
    display "CANON w4.lmm.diff.wk" 4*`v' ".hi " strtrim(string(r(ub), "%12.4f"))
}
* Joint Wald test that the arm difference is the same at every visit (3 interaction terms).
testparm 1.arm#i.visit
display "CANON w4.lmm.int.chi2 " strtrim(string(r(chi2), "%12.4f"))
display "CANON w4.lmm.int.df " r(df)
display "CANON w4.lmm.int.p " strtrim(string(r(p), "%12.4f"))
* Variance components on the SD scale (slope SD is per week).
display "CANON w4.lmm.sd_int " strtrim(string(exp(_b[lns1_1_2:_cons]), "%12.4f"))
display "CANON w4.lmm.sd_slope " strtrim(string(exp(_b[lns1_1_1:_cons]), "%12.4f"))
display "CANON w4.lmm.corr " strtrim(string(tanh(_b[atr1_1_1_2:_cons]), "%12.4f"))
display "CANON w4.lmm.sd_resid " strtrim(string(exp(_b[lnsig_e:_cons]), "%12.4f"))

* ---------------------------------------------------------------- constrained baseline (cLDA)
* Same random-effects structure, but one common baseline mean for both arms.
generate byte trt_wk4 = arm * (visit == 1)
generate byte trt_wk8 = arm * (visit == 2)
generate byte trt_wk12 = arm * (visit == 3)
mixed score i.visit trt_wk4 trt_wk8 trt_wk12 || id: week, covariance(unstructured) reml
lincom trt_wk12
display "CANON w4.clda.diff.wk12 " strtrim(string(r(estimate), "%12.4f"))
display "CANON w4.clda.diff.wk12.se " strtrim(string(r(se), "%12.4f"))
display "CANON w4.clda.diff.wk12.lo " strtrim(string(r(lb), "%12.4f"))
display "CANON w4.clda.diff.wk12.hi " strtrim(string(r(ub), "%12.4f"))

* ---------------------------------------------------------------- week-12 ANCOVA and change score
* Baseline score copied to every row of the participant.
bysort id (visit): generate double score0 = score[1]
generate double change = score - score0
* ANCOVA: week-12 score on arm, adjusted for the baseline score (OLS).
regress score i.arm score0 if visit == 3
lincom 1.arm
display "CANON w4.ancova.wk12 " strtrim(string(r(estimate), "%12.4f"))
display "CANON w4.ancova.wk12.se " strtrim(string(r(se), "%12.4f"))
display "CANON w4.ancova.wk12.lo " strtrim(string(r(lb), "%12.4f"))
display "CANON w4.ancova.wk12.hi " strtrim(string(r(ub), "%12.4f"))
display "CANON w4.ancova.baseline_coef " strtrim(string(_b[score0], "%12.4f"))
* Change score: (week-12 minus week-0) on arm, no baseline adjustment (OLS).
regress change i.arm if visit == 3
lincom 1.arm
display "CANON w4.change.wk12 " strtrim(string(r(estimate), "%12.4f"))
display "CANON w4.change.wk12.se " strtrim(string(r(se), "%12.4f"))
display "CANON w4.change.wk12.lo " strtrim(string(r(lb), "%12.4f"))
display "CANON w4.change.wk12.hi " strtrim(string(r(ub), "%12.4f"))

* ---------------------------------------------------------------- repeated-measures ANOVA
* Classical RM-ANOVA (sphericity assumed), with the Greenhouse-Geisser correction.
anova score arm / id|arm visit arm#visit, repeated(visit)
* Term 4 is arm#visit; e(gg1) and e(hf1) are the epsilons for the repeated variable.
display "CANON w4.rmanova.f_int " strtrim(string(e(F_4), "%12.4f"))
display "CANON w4.rmanova.df1 " e(df_4)
display "CANON w4.rmanova.df2 " e(df_r)
display "CANON w4.rmanova.p_int " strtrim(string(Ftail(e(df_4), e(df_r), e(F_4)), "%12.4f"))
display "CANON w4.rmanova.gg_eps " strtrim(string(e(gg1), "%12.4f"))
display "CANON w4.rmanova.hf_eps " strtrim(string(e(hf1), "%12.4f"))
display "CANON w4.rmanova.p_int_gg " strtrim(string(Ftail(e(gg1)*e(df_4), e(gg1)*e(df_r), e(F_4)), "%12.4f"))

* ---------------------------------------------------------------- binary responder
* Observed responder proportion by arm and visit (weeks 4, 8, 12).
table arm week if visit >= 1, statistic(mean responder) nformat(%6.3f)
forvalues a = 0/1 {
    local g = cond(`a' == 1, "trt", "ctrl")
    forvalues v = 1/3 {
        quietly summarize responder if arm == `a' & visit == `v'
        display "CANON w4.resp.prop.`g'.wk" 4*`v' " " strtrim(string(r(mean), "%12.4f"))
    }
}
xtset id visit
* GEE: population-averaged OR, exchangeable working correlation, robust (sandwich) SE.
xtgee responder i.arm i.visit if visit >= 1, family(binomial) link(logit) corr(exchangeable) vce(robust)
matrix R = e(R)
scalar gee_b = _b[1.arm]
scalar gee_se = _se[1.arm]
display "CANON w4.gee.logor " strtrim(string(gee_b, "%12.4f"))
display "CANON w4.gee.logor.se " strtrim(string(gee_se, "%12.4f"))
display "CANON w4.gee.or " strtrim(string(exp(gee_b), "%12.4f"))
display "CANON w4.gee.or.lo " strtrim(string(exp(gee_b - invnormal(0.975)*gee_se), "%12.4f"))
display "CANON w4.gee.or.hi " strtrim(string(exp(gee_b + invnormal(0.975)*gee_se), "%12.4f"))
display "CANON w4.gee.alpha " strtrim(string(R[1,2], "%12.4f"))
* GLMM: subject-specific OR from a random-intercept logistic model (12-point adaptive quadrature).
melogit responder i.arm i.visit if visit >= 1 || id:, intpoints(12)
scalar glmm_b = _b[1.arm]
scalar glmm_se = _se[1.arm]
scalar var_u = _b[/var(_cons[id])]
display "CANON w4.glmm.logor " strtrim(string(glmm_b, "%12.4f"))
display "CANON w4.glmm.logor.se " strtrim(string(glmm_se, "%12.4f"))
display "CANON w4.glmm.or " strtrim(string(exp(glmm_b), "%12.4f"))
display "CANON w4.glmm.or.lo " strtrim(string(exp(glmm_b - invnormal(0.975)*glmm_se), "%12.4f"))
display "CANON w4.glmm.or.hi " strtrim(string(exp(glmm_b + invnormal(0.975)*glmm_se), "%12.4f"))
display "CANON w4.glmm.var_u " strtrim(string(var_u, "%12.4f"))
* Attenuation: beta_PA is approximately beta_SS / sqrt(1 + 0.346 sigma^2).
scalar atten = sqrt(1 + 0.346*var_u)
display "CANON w4.atten.factor " strtrim(string(atten, "%12.4f"))
display "CANON w4.atten.pred_or_pa " strtrim(string(exp(glmm_b/atten), "%12.4f"))
display "CANON w4.atten.obs_ratio " strtrim(string(glmm_b/gee_b, "%12.4f"))

* ---------------------------------------------------------------- simulated truth
* The values the trial was simulated with, and the targets they imply for the estimates above.
* Score: mean 50 - 0.5 x week in control and 50 - 1 x week in treatment; random intercept SD 6,
* random slope SD 0.5 per week, intercept-slope correlation -0.2, residual SD 4.
scalar t_mu0 = 50
scalar t_slope_ctrl = -0.5
scalar t_slope_trt = -1
scalar t_sd_int = 6
scalar t_sd_slope = 0.5
scalar t_corr = -0.2
scalar t_sd_resid = 4
foreach g in ctrl trt {
    forvalues v = 0/3 {
        display "CANON w4.truth.score.mean.`g'.wk" 4*`v' " " strtrim(string(t_mu0 + t_slope_`g'*4*`v', "%12.4f"))
    }
}
* true arm difference at each visit; beta1 (week 0) and beta3 (difference in differences from week 0)
forvalues v = 0/3 {
    display "CANON w4.truth.score.diff.wk" 4*`v' " " strtrim(string((t_slope_trt - t_slope_ctrl)*4*`v' + 0, "%12.4f"))
}
display "CANON w4.truth.score.b1 " strtrim(string(0, "%12.4f"))
forvalues v = 1/3 {
    display "CANON w4.truth.score.b3.wk" 4*`v' " " strtrim(string((t_slope_trt - t_slope_ctrl)*4*`v' + 0, "%12.4f"))
}
* ANCOVA, change score and constrained baseline all target the week-12 difference in a randomised trial
display "CANON w4.truth.score.ancova_wk12 " strtrim(string((t_slope_trt - t_slope_ctrl)*12, "%12.4f"))
display "CANON w4.truth.score.change_wk12 " strtrim(string((t_slope_trt - t_slope_ctrl)*12, "%12.4f"))
display "CANON w4.truth.score.clda_wk12 " strtrim(string((t_slope_trt - t_slope_ctrl)*12, "%12.4f"))
display "CANON w4.truth.score.sd_int " strtrim(string(t_sd_int, "%12.4f"))
display "CANON w4.truth.score.sd_slope " strtrim(string(t_sd_slope, "%12.4f"))
display "CANON w4.truth.score.corr " strtrim(string(t_corr, "%12.4f"))
display "CANON w4.truth.score.sd_resid " strtrim(string(t_sd_resid, "%12.4f"))
* implied covariance of the four visits, its SD at each visit and the true Greenhouse-Geisser epsilon,
* tr(SP)^2 / (3 tr(SPSP)) with P the centring matrix
mata:
G = (st_numscalar("t_sd_int")^2, st_numscalar("t_corr")*st_numscalar("t_sd_int")*st_numscalar("t_sd_slope") \ st_numscalar("t_corr")*st_numscalar("t_sd_int")*st_numscalar("t_sd_slope"), st_numscalar("t_sd_slope")^2)
Z = (J(4, 1, 1), (0 \ 4 \ 8 \ 12))
S = Z*G*Z' + st_numscalar("t_sd_resid")^2*I(4)
P = I(4) - J(4, 4, 1/4)
st_matrix("t_sdwk", sqrt(diagonal(S))')
st_numscalar("t_gg", trace(S*P)^2 / (3*trace(S*P*S*P)))
end
forvalues v = 0/3 {
    display "CANON w4.truth.score.sd.wk" 4*`v' " " strtrim(string(t_sdwk[1, `v' + 1], "%12.4f"))
}
display "CANON w4.truth.score.gg_eps " strtrim(string(t_gg, "%12.4f"))
* Responder: logit P(response) = alpha(week) + log(2.5) x treatment + u, u normal with variance 4;
* alpha = -1.5, -1, -0.5 at weeks 4, 8, 12 (control log odds for a participant whose u is 0)
scalar t_or_ss = 2.5
scalar t_var_u = 4
matrix t_alpha = (-1.5, -1, -0.5)
display "CANON w4.truth.responder.or_ss " strtrim(string(t_or_ss, "%12.4f"))
display "CANON w4.truth.responder.logor_ss " strtrim(string(ln(t_or_ss), "%12.4f"))
display "CANON w4.truth.responder.var_u " strtrim(string(t_var_u, "%12.4f"))
display "CANON w4.truth.responder.sd_u " strtrim(string(sqrt(t_var_u), "%12.4f"))
forvalues v = 1/3 {
    display "CANON w4.truth.responder.alpha.wk" 4*`v' " " strtrim(string(t_alpha[1, `v'], "%12.4f"))
}
* the attenuation approximation at the true variance
scalar t_atten = sqrt(1 + 0.346*t_var_u)
display "CANON w4.truth.responder.atten_factor " strtrim(string(t_atten, "%12.4f"))
display "CANON w4.truth.responder.approx_logor_pa " strtrim(string(ln(t_or_ss)/t_atten, "%12.4f"))
display "CANON w4.truth.responder.approx_or_pa " strtrim(string(exp(ln(t_or_ss)/t_atten), "%12.4f"))
* exact population-averaged response probabilities and odds ratio by week: the logistic curve averaged over
* the normal random intercept (trapezoid rule on a fine grid)
mata:
real scalar pa_prob(real scalar a, real scalar v)
{
    real colvector u, f
    real scalar h
    h = 0.0005
    u = range(-40, 40, h)
    f = invlogit(a :+ u) :* normalden(u, 0, sqrt(v))
    return((sum(f) - (f[1] + f[rows(f)])/2)*h)
}
A = st_matrix("t_alpha")
b = ln(st_numscalar("t_or_ss"))
v = st_numscalar("t_var_u")
R = J(3, 3, .)
for (j = 1; j <= 3; j++) {
    R[j, 1] = pa_prob(A[j], v)
    R[j, 2] = pa_prob(A[j] + b, v)
    R[j, 3] = (R[j, 2]/(1 - R[j, 2])) / (R[j, 1]/(1 - R[j, 1]))
}
st_matrix("t_pa", R)
end
forvalues v = 1/3 {
    display "CANON w4.truth.responder.exact_prob.ctrl.wk" 4*`v' " " strtrim(string(t_pa[`v', 1], "%12.4f"))
    display "CANON w4.truth.responder.exact_prob.trt.wk" 4*`v' " " strtrim(string(t_pa[`v', 2], "%12.4f"))
    display "CANON w4.truth.responder.exact_or_pa.wk" 4*`v' " " strtrim(string(t_pa[`v', 3], "%12.4f"))
}
* one summary over the three weeks: the geometric mean of the three exact odds ratios
display "CANON w4.truth.responder.exact_or_pa_mean " strtrim(string(exp((ln(t_pa[1, 3]) + ln(t_pa[2, 3]) + ln(t_pa[3, 3]))/3), "%12.4f"))

. * GEE: population-averaged OR, exchangeable working correlation, robust (sandwich) SE.
. xtgee responder i.arm i.visit if visit >= 1, family(binomial) link(logit) corr(exchangeable) vce(robust)

Iteration 1:  Tolerance = .0005793
Iteration 2:  Tolerance = 2.742e-06
Iteration 3:  Tolerance = 5.583e-09

GEE population-averaged model                        Number of obs    =    900
Group variable: id                                   Number of groups =    300
Family: Binomial                                     Obs per group:
Link:   Logit                                                     min =      3
Correlation: exchangeable                                         avg =    3.0
                                                                  max =      3
                                                     Wald chi2(3)     =  21.23
Scale parameter = 1                                  Prob > chi2      = 0.0001

                                     (Std. err. adjusted for clustering on id)
------------------------------------------------------------------------------
             |               Robust
   responder | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         arm |
  treatment  |   .3932917   .1782105     2.21   0.027     .0440055    .7425778
             |
       visit |
          2  |   .0980098   .1337247     0.73   0.464    -.1640858    .3601053
          3  |   .5054341     .12858     3.93   0.000      .253422    .7574462
             |
       _cons |  -.6476142   .1508334    -4.29   0.000    -.9432422   -.3519862
------------------------------------------------------------------------------

# Simulated longitudinal trial: a symptom score (lower is better) at weeks 0, 4, 8, 12 and a clinician-rated
# responder (0/1) at weeks 4, 8, 12, in 300 participants randomised 1:1 to treatment or control.
# Simulated data: not evidence about any real drug or patient.
# The simulated trial file (not published) is read from a folder two levels above this script.

# Packages: lme4 (mixed models), geepack (GEE); install into the user library when missing.
for (p in c("lme4", "geepack")) {
  if (!requireNamespace(p, quietly = TRUE)) {
    try(install.packages(p, repos = "https://cloud.r-project.org"))
  }
  cat(sprintf("VERIFY %s %s\n", p, if (requireNamespace(p, quietly = TRUE)) "available" else "missing"))
}
suppressPackageStartupMessages({
  library(lme4)
  library(geepack)
})

# One CANON line per quantity: 4 decimals for estimates, integers for counts.
canon <- function(key, x, int = FALSE) {
  cat("CANON w4.", key, " ", if (int) format(x) else formatC(x, format = "f", digits = 4), "\n", sep = "")
}
z975 <- qnorm(0.975)

# Load the simulated trial file (long format: one row per participant-visit).
w4 <- read.csv(file.path("..", "..", "datasets", "W4", "W4.csv"))
w4 <- w4[order(w4$id, w4$visit), ]
w4$visitf <- factor(w4$visit)

# Sample sizes.
first <- !duplicated(w4$id)
canon("n_obs", nrow(w4), int = TRUE)
canon("n", sum(first), int = TRUE)
canon("n.ctrl", sum(first & w4$arm == 0), int = TRUE)
canon("n.trt", sum(first & w4$arm == 1), int = TRUE)
canon("n_obs_resp", sum(!is.na(w4$responder)), int = TRUE)

# ---------------------------------------------------------------- 2 x T means table
# Observed mean score in each arm at each visit (the cells a saturated mean model reproduces).
means <- tapply(w4$score, list(arm = w4$arm, week = w4$week), mean)
print(round(means, 2))
for (a in 0:1) for (wk in c(0, 4, 8, 12)) {
  canon(sprintf("mean.%s.wk%d", if (a == 1) "trt" else "ctrl", wk), means[as.character(a), as.character(wk)])
}

# ---------------------------------------------------------------- mixed model, arm x visit
# Categorical visit, arm-by-visit interaction, random intercept and random slope on week (REML).
lmm <- lmer(score ~ arm * visitf + (week | id), data = w4, REML = TRUE)
print(summary(lmm), correlation = FALSE)
b <- fixef(lmm)
V <- as.matrix(vcov(lmm))
# Wald estimate, SE and 95% CI (normal) of a linear combination L'b.
lc <- function(L) {
  est <- sum(L * b)
  se <- sqrt(drop(t(L) %*% V %*% L))
  c(est = est, se = se, lo = est - z975 * se, hi = est + z975 * se)
}
L_of <- function(names) {
  L <- setNames(numeric(length(b)), names(b))
  L[names] <- 1
  L
}
# beta1 = arm: the arm difference at week 0 (zero in truth, because arms are randomised).
r <- lc(L_of("arm"))
canon("lmm.b1", r["est"]); canon("lmm.b1.se", r["se"]); canon("lmm.b1.lo", r["lo"]); canon("lmm.b1.hi", r["hi"])
# beta3 at each visit = arm:visit: the difference in differences from week 0.
canon("lmm.b3.wk4", b["arm:visitf1"])
canon("lmm.b3.wk8", b["arm:visitf2"])
r <- lc(L_of("arm:visitf3"))
canon("lmm.b3.wk12", r["est"]); canon("lmm.b3.wk12.se", r["se"])
canon("lmm.b3.wk12.lo", r["lo"]); canon("lmm.b3.wk12.hi", r["hi"])
# Arm contrast at each visit (beta1 + beta3 at that visit): what a reader usually wants.
for (v in 0:3) {
  r <- lc(L_of(if (v == 0) "arm" else c("arm", paste0("arm:visitf", v))))
  wk <- 4 * v
  canon(sprintf("lmm.diff.wk%d", wk), r["est"])
  canon(sprintf("lmm.diff.wk%d.se", wk), r["se"])
  canon(sprintf("lmm.diff.wk%d.lo", wk), r["lo"])
  canon(sprintf("lmm.diff.wk%d.hi", wk), r["hi"])
}
# Joint Wald test that the arm difference is the same at every visit (3 interaction terms).
idx <- paste0("arm:visitf", 1:3)
chi2 <- drop(t(b[idx]) %*% solve(V[idx, idx]) %*% b[idx])
canon("lmm.int.chi2", chi2)
canon("lmm.int.df", length(idx), int = TRUE)
canon("lmm.int.p", pchisq(chi2, length(idx), lower.tail = FALSE))
# Variance components on the SD scale (slope SD is per week).
vc <- VarCorr(lmm)$id
canon("lmm.sd_int", attr(vc, "stddev")[["(Intercept)"]])
canon("lmm.sd_slope", attr(vc, "stddev")[["week"]])
canon("lmm.corr", attr(vc, "correlation")[1, 2])
canon("lmm.sd_resid", sigma(lmm))

# ---------------------------------------------------------------- constrained baseline (cLDA)
# Same random-effects structure, but one common baseline mean for both arms.
w4$trt_wk4 <- w4$arm * (w4$visit == 1)
w4$trt_wk8 <- w4$arm * (w4$visit == 2)
w4$trt_wk12 <- w4$arm * (w4$visit == 3)
clda <- lmer(score ~ visitf + trt_wk4 + trt_wk8 + trt_wk12 + (week | id), data = w4, REML = TRUE)
est <- fixef(clda)[["trt_wk12"]]
se <- sqrt(as.matrix(vcov(clda))["trt_wk12", "trt_wk12"])
canon("clda.diff.wk12", est); canon("clda.diff.wk12.se", se)
canon("clda.diff.wk12.lo", est - z975 * se); canon("clda.diff.wk12.hi", est + z975 * se)

# ---------------------------------------------------------------- week-12 ANCOVA and change score
# One row per participant: baseline and week-12 scores side by side.
wide <- data.frame(id = w4$id[w4$visit == 0], arm = w4$arm[w4$visit == 0],
                   score0 = w4$score[w4$visit == 0], score12 = w4$score[w4$visit == 3])
# ANCOVA: week-12 score on arm, adjusted for the baseline score (OLS).
anc <- lm(score12 ~ arm + score0, data = wide)
print(summary(anc))
ci <- confint(anc)["arm", ]
canon("ancova.wk12", coef(anc)[["arm"]])
canon("ancova.wk12.se", coef(summary(anc))["arm", "Std. Error"])
canon("ancova.wk12.lo", ci[[1]]); canon("ancova.wk12.hi", ci[[2]])
canon("ancova.baseline_coef", coef(anc)[["score0"]])
# Change score: (week-12 minus week-0) on arm, no baseline adjustment (OLS).
chg <- lm(I(score12 - score0) ~ arm, data = wide)
ci <- confint(chg)["arm", ]
canon("change.wk12", coef(chg)[["arm"]])
canon("change.wk12.se", coef(summary(chg))["arm", "Std. Error"])
canon("change.wk12.lo", ci[[1]]); canon("change.wk12.hi", ci[[2]])

# ---------------------------------------------------------------- repeated-measures ANOVA
# Classical RM-ANOVA (sphericity assumed): arm between participants, visit within.
rma <- aov(score ~ factor(arm) * visitf + Error(factor(id) / visitf), data = w4)
tab <- summary(rma)[["Error: factor(id):visitf"]][[1]]
print(tab)
row <- grep("factor(arm):visitf", rownames(tab), fixed = TRUE)
f_int <- tab[row, "F value"]; df1 <- tab[row, "Df"]; df2 <- tab[nrow(tab), "Df"]
canon("rmanova.f_int", f_int)
canon("rmanova.df1", df1, int = TRUE)
canon("rmanova.df2", df2, int = TRUE)
canon("rmanova.p_int", pf(f_int, df1, df2, lower.tail = FALSE))
# Greenhouse-Geisser and Huynh-Feldt epsilons from the within-arm pooled covariance of the 4 visits.
Y <- matrix(w4$score, ncol = 4, byrow = TRUE)
grp <- w4$arm[w4$visit == 0]
S <- ((sum(grp == 0) - 1) * cov(Y[grp == 0, ]) + (sum(grp == 1) - 1) * cov(Y[grp == 1, ])) / (nrow(Y) - 2)
C <- qr.Q(qr(cbind(1, poly(1:4, 3))))[, -1]
M <- t(C) %*% S %*% C
gg <- sum(diag(M))^2 / (3 * sum(diag(M %*% M)))
hf <- (nrow(Y) * 3 * gg - 2) / (3 * (nrow(Y) - 2 - 3 * gg))
canon("rmanova.gg_eps", gg)
canon("rmanova.hf_eps", min(hf, 1))
canon("rmanova.p_int_gg", pf(f_int, gg * df1, gg * df2, lower.tail = FALSE))

# ---------------------------------------------------------------- binary responder
d1 <- w4[w4$visit >= 1, ]
d1$visitf <- factor(d1$visit)
# Observed responder proportion by arm and visit (weeks 4, 8, 12).
props <- tapply(d1$responder, list(arm = d1$arm, week = d1$week), mean)
print(round(props, 3))
for (a in 0:1) for (wk in c(4, 8, 12)) {
  canon(sprintf("resp.prop.%s.wk%d", if (a == 1) "trt" else "ctrl", wk), props[as.character(a), as.character(wk)])
}
# GEE: population-averaged OR, exchangeable working correlation, robust (sandwich) SE.
gee <- geeglm(responder ~ arm + visitf, id = id, data = d1, family = binomial, corstr = "exchangeable")
print(summary(gee))
gee_b <- coef(gee)[["arm"]]
gee_se <- summary(gee)$coefficients["arm", "Std.err"]
canon("gee.logor", gee_b); canon("gee.logor.se", gee_se)
canon("gee.or", exp(gee_b))
canon("gee.or.lo", exp(gee_b - z975 * gee_se)); canon("gee.or.hi", exp(gee_b + z975 * gee_se))
canon("gee.alpha", gee$geese$alpha[[1]])
# GLMM: subject-specific OR from a random-intercept logistic model (12-point adaptive quadrature).
glmm <- glmer(responder ~ arm + visitf + (1 | id), data = d1, family = binomial, nAGQ = 12)
print(summary(glmm), correlation = FALSE)
glmm_b <- fixef(glmm)[["arm"]]
glmm_se <- sqrt(as.matrix(vcov(glmm))["arm", "arm"])
var_u <- VarCorr(glmm)$id[1, 1]
canon("glmm.logor", glmm_b); canon("glmm.logor.se", glmm_se)
canon("glmm.or", exp(glmm_b))
canon("glmm.or.lo", exp(glmm_b - z975 * glmm_se)); canon("glmm.or.hi", exp(glmm_b + z975 * glmm_se))
canon("glmm.var_u", var_u)
# Attenuation: beta_PA is approximately beta_SS / sqrt(1 + 0.346 sigma^2).
atten <- sqrt(1 + 0.346 * var_u)
canon("atten.factor", atten)
canon("atten.pred_or_pa", exp(glmm_b / atten))
canon("atten.obs_ratio", glmm_b / gee_b)

# ---------------------------------------------------------------- simulated truth
# The values the trial was simulated with, read from the simulation's settings file (not published), and
# the targets they imply for the estimates above.
tj <- jsonlite::read_json(file.path("..", "..", "datasets", "W4", "truth.json"))
sc <- tj$score
for (g in c("ctrl", "trt")) for (wk in c(0, 4, 8, 12)) {
  canon(sprintf("truth.score.mean.%s.wk%d", g, wk), sc$mu0 + sc[[paste0("slope_", g, "_per_week")]] * wk)
}
# true arm difference at each visit; beta1 (week 0) and beta3 (difference in differences from week 0)
for (wk in c(0, 4, 8, 12)) {
  canon(sprintf("truth.score.diff.wk%d", wk), (sc$slope_trt_per_week - sc$slope_ctrl_per_week) * wk + 0)
}
canon("truth.score.b1", sc$beta1_baseline_difference)
for (wk in c(4, 8, 12)) canon(sprintf("truth.score.b3.wk%d", wk), sc$beta3_did[[paste0("wk", wk)]])
# ANCOVA, change score and constrained baseline all target the week-12 difference in a randomised trial
canon("truth.score.ancova_wk12", sc$ancova_wk12_effect)
canon("truth.score.change_wk12", sc$change_score_wk12_effect)
canon("truth.score.clda_wk12", sc$constrained_baseline_wk12_effect)
canon("truth.score.sd_int", sc$sd_int); canon("truth.score.sd_slope", sc$sd_slope)
canon("truth.score.corr", sc$corr_int_slope); canon("truth.score.sd_resid", sc$sd_resid)
# implied covariance of the four visits, its SD at each visit and the true Greenhouse-Geisser epsilon,
# tr(SP)^2 / (3 tr(SPSP)) with P the centring matrix
cv <- sc$corr_int_slope * sc$sd_int * sc$sd_slope
G <- matrix(c(sc$sd_int^2, cv, cv, sc$sd_slope^2), 2)
Z <- cbind(1, c(0, 4, 8, 12))
S_true <- Z %*% G %*% t(Z) + diag(sc$sd_resid^2, 4)
P <- diag(4) - matrix(1 / 4, 4, 4)
sd_wk <- sqrt(diag(S_true))
for (k in 1:4) canon(sprintf("truth.score.sd.wk%d", 4 * (k - 1)), sd_wk[k])
canon("truth.score.gg_eps", sum(diag(S_true %*% P))^2 / (3 * sum(diag(S_true %*% P %*% S_true %*% P))))
# Responder: logit P(response) = alpha(week) + log(OR_SS) x treatment + u, u normal with variance var_u
rs <- tj$responder
alpha <- unlist(rs$alpha_ctrl_wk4_wk8_wk12)
canon("truth.responder.or_ss", rs$or_ss)
canon("truth.responder.logor_ss", log(rs$or_ss))
canon("truth.responder.var_u", rs$var_u)
canon("truth.responder.sd_u", sqrt(rs$var_u))
for (k in 1:3) canon(sprintf("truth.responder.alpha.wk%d", 4 * k), alpha[k])
# the attenuation approximation at the true variance
atten_t <- sqrt(1 + 0.346 * rs$var_u)
canon("truth.responder.atten_factor", atten_t)
canon("truth.responder.approx_logor_pa", log(rs$or_ss) / atten_t)
canon("truth.responder.approx_or_pa", exp(log(rs$or_ss) / atten_t))
# exact population-averaged response probabilities and odds ratio by week: the logistic curve averaged over
# the normal random intercept (numerical integration)
pa_prob <- function(a, v) integrate(function(u) plogis(a + u) * dnorm(u, 0, sqrt(v)), -Inf, Inf,
                                    rel.tol = 1e-10)$value
pa_or <- numeric(3)
for (k in 1:3) {
  p0 <- pa_prob(alpha[k], rs$var_u); p1 <- pa_prob(alpha[k] + log(rs$or_ss), rs$var_u)
  pa_or[k] <- (p1 / (1 - p1)) / (p0 / (1 - p0))
  canon(sprintf("truth.responder.exact_prob.ctrl.wk%d", 4 * k), p0)
  canon(sprintf("truth.responder.exact_prob.trt.wk%d", 4 * k), p1)
  canon(sprintf("truth.responder.exact_or_pa.wk%d", 4 * k), pa_or[k])
}
# one summary over the three weeks: the geometric mean of the three exact odds ratios
canon("truth.responder.exact_or_pa_mean", exp(mean(log(pa_or))))

> glmm <- glmer(responder ~ arm + visitf + (1 | id),
+     data = d1, family = binomial, nAGQ = 12)

> print(summary(glmm), correlation = FALSE)
Generalized linear mixed model fit by maximum likelihood (Adaptive
  Gauss-Hermite Quadrature, nAGQ = 12) [glmerMod]
 Family: binomial  ( logit )
Formula: responder ~ arm + visitf + (1 | id)
   Data: d1

      AIC       BIC    logLik -2*log(L)  df.resid
   1126.3    1150.3    -558.1    1116.3       895

Scaled residuals:
   Min     1Q Median     3Q    Max
-1.609 -0.539 -0.341  0.582  1.669

Random effects:
 Groups Name        Variance Std.Dev.
 id     (Intercept) 3.43     1.85
Number of obs: 900, groups:  id, 300

Fixed effects:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)   -1.016      0.240   -4.24  2.2e-05 ***
arm            0.615      0.281    2.19  0.02859 *
visitf2        0.154      0.210    0.73  0.46401
visitf3        0.795      0.213    3.74  0.00019 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1