3  Part 3：Figure 1 — MMRM 跑 BCVA + CST

預計時間：40 分鐘（最長的一章）。 deliverable：兩張 paper 等級的曲線圖（BCVA、CST 隨時間變化）。

📚 MMRM 方法學的 robust citations

本章 MMRM pipeline 對齊以下 RCT 縱向資料分析的 methodological consensus： - MMRM > LOCF / ANOVA：Mallinckrodt 20081、Mallinckrodt 2004 covariance structure2、Gueorguieva & Krystal 2004 Move Over ANOVA3 - MMRM bias / type I error：Bell & Rabe 2020 Trials4 - Missing data 處理 / MAR 假設 / FDA 觀點：Little 2012 NRC report NEJM5、LaVange & Permutt 2016 regulatory perspective6 - ICH E9 (R1) estimands framework：Mehrotra et al. 20167

AI 對 MMRM 的細節（covariance 結構選擇、Kenward-Roger DoF 校正、estimand 對應 ITT vs treatment policy）容易給看似合理但偏離 ICH/FDA 推薦的答案。請以教材的「參考程式碼」為準，並對 MMRM 結構選擇查上方 ref 自行驗證。

我有兩個 R data frame：baseline（patient_id, arm, study, region, bcva_baseline, cst_baseline, bcva_strat, lld_strat…）和 followup（patient_id, week, bcva, cst, irf, srf）。請給我 R code：

  1. 把兩份資料 merge：在 followup 加上 arm、bcva_baseline、cst_baseline、study、bcva_strat、lld_strat、region 欄位
  2. 算 bcva_change = bcva - bcva_baseline、cst_change = cst - cst_baseline
  3. 把 week 變成 factor（levels = c(“4”,“8”,“12”)，當 visit 用）
  4. 用 mmrm::mmrm() 跑兩個模型：
    1. bcva_change ~ arm + visit + arm:visit + study + bcva_strat + lld_strat + region + us(visit | patient_id)
    2. 同模型但 outcome 換成 cst_change
  5. 用 emmeans::emmeans(model, ~ arm | visit) 取 adjusted means + 95% CI
  6. 用 ggplot 畫成 line plot：x = visit、y = adjusted mean change、color by arm、加 error bar
  7. 用 patchwork 把兩張圖合在一起，標題分別 (A) BCVA、(B) CST

請加中文註解。我用 Posit.Cloud。

3.1 任務 10：什麼是 MMRM？

📋 複製這段話，貼給 AI（🆕 開新對話）：

請用 100 字以下、不用數學公式，跟一個沒學過統計的眼科醫師解釋：什麼是 mixed model for repeated measures (MMRM)？為什麼眼科 longitudinal trial 都用 MMRM 而不是 paired t-test 或 ANOVA？「unstructured covariance」是什麼意思？

速記
  • paired t-test：只能比 2 個時間點
  • ANOVA：把每個時間點當獨立 group，沒利用「同一隻眼睛」的 within-subject correlation——這正是 Gueorguieva & Krystal 2004 Archives Gen Psychiatry 〈Move Over ANOVA〉所主張要被 MMRM 取代的原因3
  • MMRM：明確 model 「同一隻眼睛在不同 visit 之間相關」1,2，並且自動處理 missing（implicit imputation, MAR assumption）5
  • unstructured covariance：不假設「兩 visit 距離越遠相關越低」，每對 visit 各有自己的 covariance；最 flexible，paper 都用這個2,4

3.2 任務 11：把 baseline 與 followup 接起來

📋 複製這段話，貼給 AI（🆕 開新對話）：

我有兩個 R data frame：baseline（一筆病人一列，欄位 patient_id, arm, study, region, bcva_baseline, cst_baseline, bcva_strat, lld_strat 等）和 followup（long format，欄位 patient_id, week, bcva, cst, irf, srf）。請給我 R code：用 dplyr::left_join() 把 baseline 的 arm、bcva_baseline、cst_baseline、study、bcva_strat、lld_strat、region 加到 followup 上，並算 bcva_change = bcva - bcva_baseline、cst_change = cst - cst_baseline。

3.2.1 參考程式碼

# ── 前置（可單獨執行）：套件 + 資料 ──
library(readr); library(dplyr); library(tidyr)
library(mmrm); library(emmeans)
library(ggplot2); library(patchwork)
arm_colours <- c("faricimab" = "#1F5C8B", "aflibercept" = "#C75D38")
baseline <- read_csv("data/faricimab_baseline.csv", show_col_types = FALSE)
followup <- read_csv("data/faricimab_followup.csv", show_col_types = FALSE)
fu_long <- followup |>
  left_join(
    baseline |>
      select(patient_id, arm, study, region,
             bcva_baseline, cst_baseline,
             bcva_strat, lld_strat),
    by = "patient_id"
  ) |>
  mutate(
    bcva_change = bcva - bcva_baseline,
    cst_change  = cst - cst_baseline,
    visit       = factor(week, levels = c(4, 8, 12)),
    arm         = factor(arm, levels = c("aflibercept", "faricimab")),
    study       = factor(study),
    bcva_strat  = factor(bcva_strat),
    lld_strat   = factor(lld_strat),
    region      = factor(region)
  ) |>
  filter(!is.na(bcva_change) | !is.na(cst_change))

glimpse(fu_long)
Rows: 3,853
Columns: 16
$ patient_id    <chr> "P0001", "P0001", "P0001", "P0002", "P0002", "P0002", "P…
$ week          <dbl> 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, …
$ bcva          <dbl> 69, 74, 71, 86, 78, 83, 80, 73, 78, 79, 85, 85, 83, 85, …
$ cst           <dbl> 142, 110, 148, 368, 306, 286, 166, 143, 128, 187, 250, 2…
$ irf           <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0,…
$ srf           <dbl> 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1,…
$ arm           <fct> faricimab, faricimab, faricimab, aflibercept, aflibercep…
$ study         <fct> LUCERNE, LUCERNE, LUCERNE, LUCERNE, LUCERNE, LUCERNE, TE…
$ region        <fct> US-Canada, US-Canada, US-Canada, US-Canada, US-Canada, U…
$ bcva_baseline <dbl> 64, 64, 64, 79, 79, 79, 76, 76, 76, 73, 73, 73, 74, 74, …
$ cst_baseline  <dbl> 300, 300, 300, 381, 381, 381, 290, 290, 290, 382, 382, 3…
$ bcva_strat    <fct> 55-73, 55-73, 55-73, >=74, >=74, >=74, >=74, >=74, >=74,…
$ lld_strat     <fct> >=33, >=33, >=33, <33, <33, <33, <33, <33, <33, >=33, >=…
$ bcva_change   <dbl> 5, 10, 7, 7, -1, 4, 4, -3, 2, 6, 12, 12, 9, 11, 18, 6, 9…
$ cst_change    <dbl> -158, -190, -152, -13, -75, -95, -124, -147, -162, -195,…
$ visit         <fct> 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, 12, 4, 8, …
Tip

為什麼把 aflibercept 設為 reference level？ 這樣 model 跑出來的 coefficient 是「faricimab 比 aflibercept 多 X letters」，方向直觀。

📖 名詞｜factor / level / reference
  • factor = R 對「類別變項」的型別（例如 arm 只能是 faricimab 或 aflibercept）
  • level = 類別變項的「可能選項」
  • reference level = 第一個 level，model 把它當基準。其他 level 的 coefficient 是「跟 reference 比」
  • 想換基準：factor(arm, levels = c("aflibercept", "faricimab")) 把 aflibercept 排第一

3.3 任務 12：跑 MMRM — BCVA

📋 複製這段話，貼給 AI（↩︎️ 續前對話）：

請用 R 套件 mmrm 跑這個模型：

mmrm(
  formula = bcva_change ~ arm + visit + arm:visit + study + bcva_strat + lld_strat + region + us(visit | patient_id),
  data = fu_long
)

跑完印 summary。然後告訴我：（1）arm:visit interaction term 在說什麼？（2）us(visit | patient_id) 是什麼意思？

3.3.1 參考程式碼

# ── 前置（可單獨執行）：套件 + 資料 ──
library(readr); library(dplyr); library(tidyr)
library(mmrm); library(emmeans)
library(ggplot2); library(patchwork)
arm_colours <- c("faricimab" = "#1F5C8B", "aflibercept" = "#C75D38")
baseline <- read_csv("data/faricimab_baseline.csv", show_col_types = FALSE)
followup <- read_csv("data/faricimab_followup.csv", show_col_types = FALSE)
# ⚠️ 沿用前面任務的物件：fu_long（請先跑完任務11，本段不重算）
m_bcva <- mmrm(
  formula = bcva_change ~ arm + visit + arm:visit +
    study + bcva_strat + lld_strat + region +
    us(visit | patient_id),
  data = fu_long |> filter(!is.na(bcva_change))
)

summary(m_bcva)
mmrm fit

Formula:     bcva_change ~ arm + visit + arm:visit + study + bcva_strat +  
    lld_strat + region + us(visit | patient_id)
Data:        filter(fu_long, !is.na(bcva_change)) (used 3853 observations from 
1329 subjects with maximum 3 timepoints)
Covariance:  unstructured (6 variance parameters)
Method:      Satterthwaite
Vcov Method: Asymptotic
Inference:   REML

Model selection criteria:
     AIC      BIC   logLik deviance 
 24310.1  24341.3 -12149.1  24298.1 

Coefficients: 
                       Estimate Std. Error         df t value Pr(>|t|)    
(Intercept)           4.631e+00  4.871e-01  1.436e+03   9.508  < 2e-16 ***
armfaricimab          6.850e-01  3.975e-01  1.319e+03   1.723    0.085 .  
visit8                1.020e+00  2.263e-01  1.286e+03   4.510 7.09e-06 ***
visit12               2.159e+00  2.240e-01  1.276e+03   9.639  < 2e-16 ***
studyTENAYA          -2.616e-01  3.589e-01  1.320e+03  -0.729    0.466    
bcva_strat>=74       -5.397e-01  5.771e-01  1.320e+03  -0.935    0.350    
bcva_strat55-73      -1.079e-01  3.910e-01  1.320e+03  -0.276    0.783    
lld_strat>=33        -1.052e-02  3.810e-01  1.321e+03  -0.028    0.978    
regionUS-Canada      -2.577e-01  3.607e-01  1.320e+03  -0.714    0.475    
armfaricimab:visit8  -1.483e-03  3.201e-01  1.288e+03  -0.005    0.996    
armfaricimab:visit12  9.242e-02  3.162e-01  1.272e+03   0.292    0.770    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Covariance estimate:
         4       8      12
4  51.8242 36.5658 36.5880
8  36.5658 54.1432 38.2153
12 36.5880 38.2153 53.0027
模型 syntax 拆解
  • arm + visit + arm:visit — main effects 與 interaction（讓兩 arm 的時間趨勢可以不一樣）
  • study + bcva_strat + lld_strat + region — 分層因子（paper 也是這樣）；多重 prognostic factor 的 adjustment 原則見 Kahan & Morris 20138
  • us(visit | patient_id) — 「同一個 patient 跨 visit 的 covariance 是 unstructured」；UN 在 RCT longitudinal data 是 FDA 接受的預設選擇2,6
  • Estimand 對應：本書 demo 採「treatment policy」估計量（如 ICH E9 R1）7
📖 名詞｜interaction / unstructured covariance / emmeans
  • interaction (arm:visit) — 「兩 arm 的時間趨勢可以不一樣」。沒有 interaction 的話 model 強制兩條線平行，那就抓不到「藥效隨時間變化」
  • us(visit | patient_id) — 念作「visit given patient_id is unstructured」。意思是「同一個病人的 4 個 visit 之間互相相關，但相關有多強讓 model 自己估，不假設規則」
  • emmeans() — estimated marginal means，「校正過其他變項後的邊際平均」。比直接 mean() 更接近 paper 寫的「adjusted mean」

3.4 任務 13：取 adjusted means

📋 複製這段話，貼給 AI（↩︎️ 續前對話）：

用 emmeans::emmeans() 從上面的 mmrm model 取出每個 arm 在每個 visit 的 adjusted mean change from baseline，加上 95% CI。然後算每個 visit 兩 arm 的差異（faricimab − aflibercept）。

3.4.1 參考程式碼

# ── 前置（可單獨執行）：套件 + 資料 ──
library(readr); library(dplyr); library(tidyr)
library(mmrm); library(emmeans)
library(ggplot2); library(patchwork)
arm_colours <- c("faricimab" = "#1F5C8B", "aflibercept" = "#C75D38")
baseline <- read_csv("data/faricimab_baseline.csv", show_col_types = FALSE)
followup <- read_csv("data/faricimab_followup.csv", show_col_types = FALSE)
# ⚠️ 沿用前面任務的物件：m_bcva（請先跑完任務11、12，本段不重算）
emm_bcva <- emmeans(m_bcva, ~ arm | visit)
emm_bcva
visit = 4:
 arm         emmean    SE   df lower.CL upper.CL
 aflibercept   4.15 0.305 1385     3.55     4.75
 faricimab     4.84 0.312 1396     4.22     5.45

visit = 8:
 arm         emmean    SE   df lower.CL upper.CL
 aflibercept   5.17 0.312 1398     4.56     5.78
 faricimab     5.85 0.318 1406     5.23     6.48

visit = 12:
 arm         emmean    SE   df lower.CL upper.CL
 aflibercept   6.31 0.310 1382     5.70     6.92
 faricimab     7.09 0.317 1401     6.47     7.71

Results are averaged over the levels of: study, bcva_strat, lld_strat, region 
Confidence level used: 0.95 
contrast(emm_bcva, method = "revpairwise") |>
  summary(infer = TRUE)
visit = 4:
 contrast                estimate    SE   df lower.CL upper.CL t.ratio p.value
 faricimab - aflibercept    0.685 0.397 1319  -0.0947     1.46   1.723  0.0850

visit = 8:
 contrast                estimate    SE   df lower.CL upper.CL t.ratio p.value
 faricimab - aflibercept    0.684 0.407 1315  -0.1157     1.48   1.678  0.0936

visit = 12:
 contrast                estimate    SE   df lower.CL upper.CL t.ratio p.value
 faricimab - aflibercept    0.777 0.405 1308  -0.0161     1.57   1.922  0.0548

Results are averaged over the levels of: study, bcva_strat, lld_strat, region 
Confidence level used: 0.95

對照 paper Fig 1A：week 12 faricimab 應該 ~+6.7 letters、aflibercept ~+5.9，差 ~+0.8。

3.5 任務 14：畫 Figure 1A — BCVA

📋 複製這段話，貼給 AI（↩︎️ 續前對話）：

請用 ggplot2 畫一張 line plot：x 軸是 visit (4, 8, 12 週)、y 軸是 adjusted mean change from baseline、不同 arm 用不同顏色（faricimab 深藍 #1F5C8B、aflibercept 赭橙 #C75D38）、加 error bar (95% CI)、加 point。標題 “(A) BCVA change from baseline”，y 軸 label “Adjusted Mean BCVA Change (Letters)”。

3.5.1 參考程式碼

# ── 前置（可單獨執行）：套件 + 資料 ──
library(readr); library(dplyr); library(tidyr)
library(mmrm); library(emmeans)
library(ggplot2); library(patchwork)
arm_colours <- c("faricimab" = "#1F5C8B", "aflibercept" = "#C75D38")
baseline <- read_csv("data/faricimab_baseline.csv", show_col_types = FALSE)
followup <- read_csv("data/faricimab_followup.csv", show_col_types = FALSE)
# ⚠️ 沿用前面任務的物件：emm_bcva（請先跑完任務11~13，本段不重算）

emm_bcva_df <- as.data.frame(emm_bcva) |>
  mutate(week = as.numeric(as.character(visit)))

# 加一筆 week=0 的 baseline reference (change = 0)
emm_bcva_df <- bind_rows(
  data.frame(arm = c("aflibercept","faricimab"),
             week = 0, emmean = 0, SE = 0,
             df = NA, lower.CL = 0, upper.CL = 0,
             visit = factor("0", levels = c("0","4","8","12"))),
  emm_bcva_df
)

p_bcva <- ggplot(emm_bcva_df,
       aes(x = week, y = emmean, colour = arm, group = arm)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = 0.4, linewidth = 0.8) +
  scale_colour_manual(values = arm_colours, name = NULL,
                      labels = c("aflibercept" = "Aflibercept Q8W",
                                 "faricimab"   = "Faricimab up to Q16W")) +
  scale_x_continuous(breaks = c(0, 4, 8, 12)) +
  labs(title = "(A) BCVA change from baseline",
       x = "Time (Weeks)",
       y = "Adjusted Mean BCVA Change\n(ETDRS Letters)") +
  theme(legend.position = "bottom")

p_bcva

Figure 1A. Adjusted mean change from baseline in BCVA over 12-week head-to-head dosing phase.

3.6 任務 15：複製整段，跑 CST → Figure 1B

📋 複製這段話，貼給 AI（↩︎️ 續前對話）：

把上面 BCVA 那段 mmrm + emmeans + ggplot 完整 code 複製一份，但 outcome 換成 cst_change。最後用 patchwork 套件把 Figure 1A (BCVA) 與 Figure 1B (CST) 並排成 2 panel，模仿 paper Fig 1。

3.6.1 參考程式碼

# ── 前置（可單獨執行）：套件 + 資料 ──
library(readr); library(dplyr); library(tidyr)
library(mmrm); library(emmeans)
library(ggplot2); library(patchwork)
arm_colours <- c("faricimab" = "#1F5C8B", "aflibercept" = "#C75D38")
baseline <- read_csv("data/faricimab_baseline.csv", show_col_types = FALSE)
followup <- read_csv("data/faricimab_followup.csv", show_col_types = FALSE)
# ⚠️ 沿用前面任務的物件：fu_long（請先跑完任務11，本段不重算）
m_cst <- mmrm(
  formula = cst_change ~ arm + visit + arm:visit +
    study + bcva_strat + lld_strat + region +
    us(visit | patient_id),
  data = fu_long |> filter(!is.na(cst_change))
)

emm_cst <- emmeans(m_cst, ~ arm | visit)
emm_cst
visit = 4:
 arm         emmean   SE   df lower.CL upper.CL
 aflibercept   -116 1.63 1417     -119     -113
 faricimab     -127 1.66 1434     -130     -124

visit = 8:
 arm         emmean   SE   df lower.CL upper.CL
 aflibercept   -128 1.61 1406     -131     -125
 faricimab     -138 1.64 1416     -141     -135

visit = 12:
 arm         emmean   SE   df lower.CL upper.CL
 aflibercept   -134 1.63 1407     -137     -130
 faricimab     -145 1.66 1434     -148     -141

Results are averaged over the levels of: study, bcva_strat, lld_strat, region 
Confidence level used: 0.95 
emm_cst_df <- as.data.frame(emm_cst) |>
  mutate(week = as.numeric(as.character(visit))) |>
  bind_rows(
    data.frame(arm = c("aflibercept","faricimab"),
               week = 0, emmean = 0, SE = 0,
               df = NA, lower.CL = 0, upper.CL = 0,
               visit = factor("0", levels = c("0","4","8","12"))),
    x = _
  ) |>
  arrange(arm, week)

p_cst <- ggplot(emm_cst_df,
       aes(x = week, y = emmean, colour = arm, group = arm)) +
  geom_line(linewidth = 1) +
  geom_point(size = 3) +
  geom_errorbar(aes(ymin = lower.CL, ymax = upper.CL),
                width = 0.4, linewidth = 0.8) +
  scale_colour_manual(values = arm_colours, name = NULL,
                      labels = c("aflibercept" = "Aflibercept Q8W",
                                 "faricimab"   = "Faricimab up to Q16W")) +
  scale_x_continuous(breaks = c(0, 4, 8, 12)) +
  labs(title = "(B) CST change from baseline",
       x = "Time (Weeks)",
       y = "Adjusted Mean CST Change\n(μm)") +
  theme(legend.position = "bottom")

p_cst

Figure 1B. Adjusted mean change from baseline in CST over 12-week head-to-head dosing phase.

3.6.2 兩張合成 paper Fig 1

p_bcva + p_cst + plot_layout(guides = "collect") &
  theme(legend.position = "bottom")

Figure 1. Reproduced from Cheung et al. 2025 — Adjusted mean change from baseline in (A) BCVA and (B) CST.
🎯 你做出 paper 圖了

回去看一下 paper Figure 1（refs/tenaya-lucerne-paper.pdf 第 4 頁），你的圖跟 paper 走勢應該一致： - BCVA：兩 arm 一路上升，幾乎重疊，week 12 達 ~+6 letters - CST：兩 arm 都急速下降，faricimab 比 aflibercept 多降 ~12–14 μm

本章重點
  1. MMRM = 「同一隻眼被測多次」的標配模型
  2. unstructured covariance + MAR implicit imputation = paper 標準寫法
  3. mmrm::mmrm() 一行就能跑、emmeans 取 adjusted mean、ggplot 出圖
  4. 這個 pipeline 跟你後面換院內資料時是一模一樣的 code
1.
Mallinckrodt CH, Lane PW, Schnell D, Peng Y, Mancuso JP. Recommendations for the primary analysis of continuous endpoints in longitudinal clinical trials. Drug Information Journal. 2008;42(4):303-319.
2.
Mallinckrodt CH, Kaiser CJ, Watkin JG, Molenberghs G, Carroll RJ. The effect of correlation structure on treatment contrasts estimated from incomplete clinical trial data with likelihood-based repeated measures compared with last observation carried forward ANOVA. Clinical Trials. 2004;1(6):477-489. doi:10.1191/1740774504cn049oa
3.
Gueorguieva R, Krystal JH. Move over ANOVA: Progress in analyzing repeated-measures data and its reflection in papers published in the Archives of General Psychiatry. Archives of General Psychiatry. 2004;61(3):310-317. doi:10.1001/archpsyc.61.3.310
4.
Bell ML, Rabe BA. The mixed model for repeated measures for cluster randomized trials: A simulation study investigating bias and type I error with missing continuous data. Trials. 2020;21:148. doi:10.1186/s13063-020-4114-9
5.
Little RJ, D’Agostino R, Cohen ML, et al. The prevention and treatment of missing data in clinical trials. New England Journal of Medicine. 2012;367(14):1355-1360. doi:10.1056/NEJMsr1203730
6.
LaVange LM, Permutt T. A regulatory perspective on missing data in the aftermath of the NRC report. Statistics in Medicine. 2016;35(17):2853-2864. doi:10.1002/sim.6840
7.
Mehrotra DV, Hemmings RJ, Russek-Cohen E. Seeking harmony: Estimands and sensitivity analyses for confirmatory clinical trials. Clinical Trials. 2016;13(4):456-458. doi:10.1177/1740774516633115
8.
Kahan BC, Morris TP. Adjusting for multiple prognostic factors in the analysis of randomised trials. BMC Medical Research Methodology. 2013;13:99. doi:10.1186/1471-2288-13-99