RCTs with baseline and follow-up measurements of a continuous outcome

Marco Cattaneo

Department of Clinical Research, University of Basel

18 November 2025

pre-post design

statistical analysis of trt effect in RCT with bsl and fup measurements of a continuous outcome (y):

fup ~ trt: outcome without considering baseline
fup - bsl ~ trt: change from baseline
fup ~ trt + bsl: ANCOVA (with or without interaction)
y ~ trt * time + pat: longitudinal with random/fixed effects for pat (with or without trt effect at baseline)

all methods can be extended in various ways (covariate adjustment, heteroskedasticity, additional correlations due to centers/patients, …)

methods comparison

all methods estimate the same quantity in RCTs (under the corresponding model assumptions)
if we have an unbalanced baseline in an RCT, then completely ignoring the baseline can be problematic, and considering only the change from baseline can be problematic as well (ceiling/floor effect)
in observational studies, it is in general not clear what exactly we want to estimate and whether we should adjust for the baseline value (see e.g. Lord’s paradox)
there is a large, partially conflicting literature comparing the methods in RCTs and observational studies (authors often simulate from one model and find out that the corresponding method is best…)
I simulated some non-normal bivariate data that may resemble data often met in RCTs: all methods are based on incorrect distributional assumptions!

simulations

load("res.rda")
library(MASS)
library(ggplot2)
library(lme4)
library(lmerTest)
library(pbkrtest)
library(future.apply)
set.seed(0)

n = 1000000
mu0 = c(0,2)
mu1 = c(0,3)
Sigma = matrix(c(2,1,1,2),2)
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
round(colMeans(dat[1:(n/2),]),2)  # means bsl/fup for trt0

[1] 5.00 6.78

round(cov(dat[1:(n/2),]), 2)  # cov matrix bsl/fup for trt0

     [,1] [,2]
[1,] 1.74 0.76
[2,] 0.76 1.36

round(colMeans(dat[n/2+1:(n/2),]),2)  # means bsl/fup for trt1

[1] 5.00 7.54

round(cov(dat[n/2+1:(n/2),]),2)  # cov matrix bsl/fup for trt1

     [,1] [,2]
[1,] 1.74 0.65
[2,] 0.65 1.01

es = mean(dat[n/2+1:(n/2),2])-mean(dat[1:(n/2),2])
round(es, 2)  # effect size of trt (difference in mean fup)

[1] 0.76

n = 50
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
db.wide = data.frame(pat = 1:n, bsl = dat[,1], fup = dat[,2], trt = rep(0:1, each = n/2))
db.long = reshape(db.wide, c("bsl","fup"), "y", "time", times = 0:1, direction = "long")
db.long$trt.b = factor(db.long$trt + 2*(1-db.long$time), 0:3, c(0:1,"b","b"))
ggplot(db.wide, aes(bsl, fup, color = factor(trt))) +
  geom_point()

# outcome without considering baseline
summary(lm(fup ~ trt, db.wide))


Call:
lm(formula = fup ~ trt, data = db.wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.33255 -0.42923  0.00668  0.68734  1.45364 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   7.2616     0.1740  41.729   <2e-16 ***
trt           0.4122     0.2461   1.675      0.1    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.8701 on 48 degrees of freedom
Multiple R-squared:  0.05521,   Adjusted R-squared:  0.03553 
F-statistic: 2.805 on 1 and 48 DF,  p-value: 0.1005

# change from baseline
summary(lm((fup - bsl) ~ trt, db.wide))


Call:
lm(formula = (fup - bsl) ~ trt, data = db.wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-2.8555 -0.6115 -0.2386  0.6748  1.8775 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)   1.7384     0.2038   8.531 3.51e-11 ***
trt           0.9616     0.2882   3.337  0.00164 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.019 on 48 degrees of freedom
Multiple R-squared:  0.1883,    Adjusted R-squared:  0.1714 
F-statistic: 11.13 on 1 and 48 DF,  p-value: 0.001642

# ANCOVA with one slope
summary(lm(fup ~ trt + bsl, db.wide))


Call:
lm(formula = fup ~ trt + bsl, data = db.wide)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.6291 -0.3667  0.1306  0.4569  1.2626 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  4.99799    0.47595  10.501 6.46e-14 ***
trt          0.63735    0.20618   3.091  0.00335 ** 
bsl          0.40984    0.08223   4.984 8.88e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7112 on 47 degrees of freedom
Multiple R-squared:  0.3819,    Adjusted R-squared:  0.3556 
F-statistic: 14.52 on 2 and 47 DF,  p-value: 1.231e-05

# ANCOVA with two slopes
summary(lm(fup ~ trt * scale(bsl), db.wide))


Call:
lm(formula = fup ~ trt * scale(bsl), data = db.wide)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.64173 -0.37086  0.08824  0.46115  1.36494 

Coefficients:
               Estimate Std. Error t value Pr(>|t|)    
(Intercept)      7.1825     0.1450  49.542  < 2e-16 ***
trt              0.6309     0.2045   3.085  0.00344 ** 
scale(bsl)       0.3646     0.1546   2.359  0.02263 *  
trt:scale(bsl)   0.2788     0.2077   1.342  0.18614    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7052 on 46 degrees of freedom
Multiple R-squared:  0.4052,    Adjusted R-squared:  0.3664 
F-statistic: 10.44 on 3 and 46 DF,  p-value: 2.333e-05

# fixed-eff longitudinal with one baseline (corresponds to change from baseline)
summary(lm(y ~ trt.b + factor(pat), db.long))


Call:
lm(formula = y ~ trt.b + factor(pat), data = db.long)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4277 -0.3257  0.0000  0.3257  1.4277 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    8.35573    0.51952  16.084  < 2e-16 ***
trt.b1         0.96161    0.28818   3.337 0.001642 ** 
trt.bb        -1.73838    0.20377  -8.531 3.51e-11 ***
factor(pat)2   0.03067    0.72045   0.043 0.966225    
factor(pat)3  -1.60129    0.72045  -2.223 0.030985 *  
factor(pat)4  -0.75546    0.72045  -1.049 0.299616    
factor(pat)5  -0.68529    0.72045  -0.951 0.346268    
factor(pat)6  -1.26744    0.72045  -1.759 0.084909 .  
factor(pat)7  -1.77220    0.72045  -2.460 0.017552 *  
factor(pat)8  -1.77219    0.72045  -2.460 0.017552 *  
factor(pat)9  -2.48139    0.72045  -3.444 0.001198 ** 
factor(pat)10 -1.30084    0.72045  -1.806 0.077255 .  
factor(pat)11 -0.65782    0.72045  -0.913 0.365770    
factor(pat)12 -0.19312    0.72045  -0.268 0.789808    
factor(pat)13 -2.19266    0.72045  -3.043 0.003787 ** 
factor(pat)14 -0.65437    0.72045  -0.908 0.368269    
factor(pat)15 -0.44374    0.72045  -0.616 0.540857    
factor(pat)16  0.13527    0.72045   0.188 0.851859    
factor(pat)17 -2.00798    0.72045  -2.787 0.007596 ** 
factor(pat)18 -2.01766    0.72045  -2.801 0.007330 ** 
factor(pat)19 -2.69980    0.72045  -3.747 0.000480 ***
factor(pat)20  0.35317    0.72045   0.490 0.626219    
factor(pat)21 -1.63116    0.72045  -2.264 0.028124 *  
factor(pat)22 -1.07479    0.72045  -1.492 0.142283    
factor(pat)23 -1.08937    0.72045  -1.512 0.137070    
factor(pat)24 -0.56002    0.72045  -0.777 0.440784    
factor(pat)25 -1.01276    0.72045  -1.406 0.166240    
factor(pat)26 -0.40800    0.73471  -0.555 0.581258    
factor(pat)27 -1.84039    0.73471  -2.505 0.015697 *  
factor(pat)28 -0.80240    0.73471  -1.092 0.280228    
factor(pat)29 -4.05415    0.73471  -5.518 1.35e-06 ***
factor(pat)30 -0.49340    0.73471  -0.672 0.505083    
factor(pat)31 -0.65682    0.73471  -0.894 0.375791    
factor(pat)32 -2.04415    0.73471  -2.782 0.007695 ** 
factor(pat)33 -3.36597    0.73471  -4.581 3.30e-05 ***
factor(pat)34 -1.54324    0.73471  -2.100 0.040966 *  
factor(pat)35 -3.13185    0.73471  -4.263 9.40e-05 ***
factor(pat)36 -1.52205    0.73471  -2.072 0.043699 *  
factor(pat)37 -1.74678    0.73471  -2.377 0.021463 *  
factor(pat)38 -1.21384    0.73471  -1.652 0.105036    
factor(pat)39 -1.08659    0.73471  -1.479 0.145694    
factor(pat)40 -2.18232    0.73471  -2.970 0.004635 ** 
factor(pat)41 -1.00132    0.73471  -1.363 0.179282    
factor(pat)42 -0.40932    0.73471  -0.557 0.580039    
factor(pat)43 -2.49290    0.73471  -3.393 0.001393 ** 
factor(pat)44 -1.64385    0.73471  -2.237 0.029938 *  
factor(pat)45 -1.52581    0.73471  -2.077 0.043202 *  
factor(pat)46  0.04269    0.73471   0.058 0.953912    
factor(pat)47 -2.02597    0.73471  -2.757 0.008214 ** 
factor(pat)48 -2.56599    0.73471  -3.493 0.001038 ** 
factor(pat)49 -0.49123    0.73471  -0.669 0.506950    
factor(pat)50 -2.88272    0.73471  -3.924 0.000277 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7204 on 48 degrees of freedom
Multiple R-squared:  0.8963,    Adjusted R-squared:  0.786 
F-statistic: 8.131 on 51 and 48 DF,  p-value: 8.409e-12

# fixed-eff longitudinal with two baselines (corresponds to change from baseline)
summary(lm(y ~ trt * time + factor(pat), db.long))


Call:
lm(formula = y ~ trt * time + factor(pat), data = db.long)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.4277 -0.3257  0.0000  0.3257  1.4277 

Coefficients: (1 not defined because of singularities)
              Estimate Std. Error t value Pr(>|t|)    
(Intercept)    6.61735    0.51952  12.737  < 2e-16 ***
trt           -2.88272    0.73471  -3.924 0.000277 ***
time           1.73838    0.20377   8.531 3.51e-11 ***
factor(pat)2   0.03067    0.72045   0.043 0.966225    
factor(pat)3  -1.60129    0.72045  -2.223 0.030985 *  
factor(pat)4  -0.75546    0.72045  -1.049 0.299616    
factor(pat)5  -0.68529    0.72045  -0.951 0.346268    
factor(pat)6  -1.26744    0.72045  -1.759 0.084909 .  
factor(pat)7  -1.77220    0.72045  -2.460 0.017552 *  
factor(pat)8  -1.77219    0.72045  -2.460 0.017552 *  
factor(pat)9  -2.48139    0.72045  -3.444 0.001198 ** 
factor(pat)10 -1.30084    0.72045  -1.806 0.077255 .  
factor(pat)11 -0.65782    0.72045  -0.913 0.365770    
factor(pat)12 -0.19312    0.72045  -0.268 0.789808    
factor(pat)13 -2.19266    0.72045  -3.043 0.003787 ** 
factor(pat)14 -0.65437    0.72045  -0.908 0.368269    
factor(pat)15 -0.44374    0.72045  -0.616 0.540857    
factor(pat)16  0.13527    0.72045   0.188 0.851859    
factor(pat)17 -2.00798    0.72045  -2.787 0.007596 ** 
factor(pat)18 -2.01766    0.72045  -2.801 0.007330 ** 
factor(pat)19 -2.69980    0.72045  -3.747 0.000480 ***
factor(pat)20  0.35317    0.72045   0.490 0.626219    
factor(pat)21 -1.63116    0.72045  -2.264 0.028124 *  
factor(pat)22 -1.07479    0.72045  -1.492 0.142283    
factor(pat)23 -1.08937    0.72045  -1.512 0.137070    
factor(pat)24 -0.56002    0.72045  -0.777 0.440784    
factor(pat)25 -1.01276    0.72045  -1.406 0.166240    
factor(pat)26  2.47472    0.72045   3.435 0.001231 ** 
factor(pat)27  1.04232    0.72045   1.447 0.154460    
factor(pat)28  2.08031    0.72045   2.888 0.005807 ** 
factor(pat)29 -1.17143    0.72045  -1.626 0.110503    
factor(pat)30  2.38931    0.72045   3.316 0.001743 ** 
factor(pat)31  2.22589    0.72045   3.090 0.003330 ** 
factor(pat)32  0.83857    0.72045   1.164 0.250195    
factor(pat)33 -0.48326    0.72045  -0.671 0.505579    
factor(pat)34  1.33947    0.72045   1.859 0.069131 .  
factor(pat)35 -0.24914    0.72045  -0.346 0.730999    
factor(pat)36  1.36067    0.72045   1.889 0.064990 .  
factor(pat)37  1.13594    0.72045   1.577 0.121430    
factor(pat)38  1.66887    0.72045   2.316 0.024849 *  
factor(pat)39  1.79613    0.72045   2.493 0.016166 *  
factor(pat)40  0.70040    0.72045   0.972 0.335838    
factor(pat)41  1.88139    0.72045   2.611 0.011997 *  
factor(pat)42  2.47340    0.72045   3.433 0.001238 ** 
factor(pat)43  0.38981    0.72045   0.541 0.590960    
factor(pat)44  1.23886    0.72045   1.720 0.091953 .  
factor(pat)45  1.35691    0.72045   1.883 0.065709 .  
factor(pat)46  2.92540    0.72045   4.061 0.000180 ***
factor(pat)47  0.85675    0.72045   1.189 0.240215    
factor(pat)48  0.31672    0.72045   0.440 0.662185    
factor(pat)49  2.39148    0.72045   3.319 0.001728 ** 
factor(pat)50       NA         NA      NA       NA    
trt:time       0.96161    0.28818   3.337 0.001642 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7204 on 48 degrees of freedom
Multiple R-squared:  0.8963,    Adjusted R-squared:  0.786 
F-statistic: 8.131 on 51 and 48 DF,  p-value: 8.409e-12

# mixed-eff longitudinal with one baseline
summary(lmer(y ~ trt.b + (1|pat), db.long))

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ trt.b + (1 | pat)
   Data: db.long

REML criterion at convergence: 284.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.87718 -0.49717  0.00964  0.52858  2.39153 

Random effects:
 Groups   Name        Variance Std.Dev.
 pat      (Intercept) 0.6688   0.8178  
 Residual             0.5242   0.7240  
Number of obs: 100, groups:  pat, 50

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   7.1076     0.2006 96.0272  35.439  < 2e-16 ***
trt.b1        0.7202     0.2558 72.7491   2.815  0.00627 ** 
trt.bb       -1.8591     0.1932 57.4925  -9.622 1.39e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr) trt.b1
trt.b1 -0.638       
trt.bb -0.693  0.662

summary(lmer(y ~ trt.b + (1|pat), db.long), ddf = "Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: y ~ trt.b + (1 | pat)
   Data: db.long

REML criterion at convergence: 284.6

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.87718 -0.49717  0.00964  0.52858  2.39153 

Random effects:
 Groups   Name        Variance Std.Dev.
 pat      (Intercept) 0.6688   0.8178  
 Residual             0.5242   0.7240  
Number of obs: 100, groups:  pat, 50

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   7.1076     0.2017 96.0569  35.236  < 2e-16 ***
trt.b1        0.7202     0.2594 73.3136   2.776  0.00698 ** 
trt.bb       -1.8591     0.1944 58.2228  -9.562 1.53e-13 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
       (Intr) trt.b1
trt.b1 -0.638       
trt.bb -0.693  0.662

# mixed-eff longitudinal with two baselines (corresponds to change from baseline)
summary(lmer(y ~ trt * time + (1|pat), db.long))

Linear mixed model fit by REML. t-tests use Satterthwaite's method ['lmerModLmerTest']
Formula: y ~ trt * time + (1 | pat)
   Data: db.long

REML criterion at convergence: 282

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.73662 -0.49960  0.00951  0.58134  2.24256 

Random effects:
 Groups   Name        Variance Std.Dev.
 pat      (Intercept) 0.6387   0.7992  
 Residual             0.5190   0.7204  
Number of obs: 100, groups:  pat, 50

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   5.5233     0.2152 73.6000  25.666  < 2e-16 ***
trt          -0.5494     0.3043 73.6000  -1.805  0.07510 .  
time          1.7384     0.2038 48.0000   8.531 3.51e-11 ***
trt:time      0.9616     0.2882 48.0000   3.337  0.00164 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
         (Intr) trt    time  
trt      -0.707              
time     -0.473  0.335       
trt:time  0.335 -0.473 -0.707

summary(lmer(y ~ trt * time + (1|pat), db.long), ddf = "Kenward-Roger")

Linear mixed model fit by REML. t-tests use Kenward-Roger's method ['lmerModLmerTest']
Formula: y ~ trt * time + (1 | pat)
   Data: db.long

REML criterion at convergence: 282

Scaled residuals: 
     Min       1Q   Median       3Q      Max 
-1.73662 -0.49960  0.00951  0.58134  2.24256 

Random effects:
 Groups   Name        Variance Std.Dev.
 pat      (Intercept) 0.6387   0.7992  
 Residual             0.5190   0.7204  
Number of obs: 100, groups:  pat, 50

Fixed effects:
            Estimate Std. Error      df t value Pr(>|t|)    
(Intercept)   5.5233     0.2152 73.6000  25.666  < 2e-16 ***
trt          -0.5494     0.3043 73.6000  -1.805  0.07510 .  
time          1.7384     0.2038 48.0000   8.531 3.51e-11 ***
trt:time      0.9616     0.2882 48.0000   3.337  0.00164 ** 
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Correlation of Fixed Effects:
         (Intr) trt    time  
trt      -0.707              
time     -0.473  0.335       
trt:time  0.335 -0.473 -0.707

sim = function(n, mu0, mu1, Sigma) {
  dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
    (\(x) tanh(x/5)*5+5)()
  db.wide = data.frame(pat = 1:n, bsl = dat[,1], fup = dat[,2], trt = rep(0:1, each = n/2))
  db.long = reshape(db.wide, c("bsl","fup"), "y", "time", times = 0:1, direction = "long")
  db.long$trt.b = factor(db.long$trt + 2*(1-db.long$time), 0:3, c(0:1,"b","b"))
  c(summary(lm(fup ~ trt, db.wide))$coef["trt",c(1,4)],
    summary(lm((fup - bsl) ~ trt, db.wide))$coef["trt",c(1,4)],
    summary(lm(fup ~ trt + bsl, db.wide))$coef["trt",c(1,4)],
    summary(lm(fup ~ trt * scale(bsl), db.wide))$coef["trt",c(1,4)],
    summary(lmer(y ~ trt.b + (1|pat), db.long))$coef["trt.b1",c(1,5)],
    summary(lmer(y ~ trt.b + (1|pat), db.long), ddf = "Kenward-Roger")$coef["trt.b1",c(1,5)])}

summarize = function(res, es) {
  tab = round(t(sapply(1:(nrow(res)/2), \(x){
    c(mean(res[2*x-1,])-es,sd(res[2*x-1,]),sqrt(mean((res[2*x-1,]-es)^2)),mean(res[2*x,]<0.05))})), 2)
  colnames(tab) = c("bias", "SD estimate", "RMSE", "power")
  rownames(tab) = c("outcome without considering baseline", 
                    "change from baseline", 
                    "ANCOVA with one slope", 
                    "ANCOVA with two slopes", 
                    "mixed-eff longitudinal with one baseline and Satterthwaite approx", 
                    "mixed-eff longitudinal with one baseline and Kenward-Roger approx")
  print(tab)}

plan(multisession)
# res1.1 = future_replicate(10000, sim(n, mu0, mu1, Sigma), future.seed = 0)
summarize(res1.1, es)  # change bad, outcome not good, ANCOVA best

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.31 0.31  0.68
change from baseline                                                 0        0.34 0.34  0.56
ANCOVA with one slope                                                0        0.27 0.27  0.78
ANCOVA with two slopes                                               0        0.27 0.27  0.78
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.27 0.27  0.73
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.27 0.27  0.71

# res0.1 = future_replicate(10000, sim(n, mu0, mu0, Sigma), future.seed = 0)
summarize(res0.1, 0)

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.33 0.33  0.05
change from baseline                                                 0        0.35 0.35  0.05
ANCOVA with one slope                                                0        0.29 0.29  0.05
ANCOVA with two slopes                                               0        0.29 0.29  0.05
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.29 0.29  0.04
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.29 0.29  0.04

n = 1000000
mu0 = c(2,2)
mu1 = c(2,3)
Sigma = matrix(c(3,2,2,3),2)
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
round(colMeans(dat[1:(n/2),]),2)  # means bsl/fup for trt0

[1] 6.73 6.73

round(cov(dat[1:(n/2),]), 2)  # cov matrix bsl/fup for trt0

     [,1] [,2]
[1,] 1.96 1.29
[2,] 1.29 1.96

round(colMeans(dat[n/2+1:(n/2),]),2)  # means bsl/fup for trt1

[1] 6.73 7.48

round(cov(dat[n/2+1:(n/2),]),2)  # cov matrix bsl/fup for trt1

     [,1] [,2]
[1,] 1.97 1.14
[2,] 1.14 1.51

es = mean(dat[n/2+1:(n/2),2])-mean(dat[1:(n/2),2])
round(es, 2)  # effect size of trt (difference in mean fup)

[1] 0.75

n = 50
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
db.wide = data.frame(pat = 1:n, bsl = dat[,1], fup = dat[,2], trt = rep(0:1, each = n/2))
ggplot(db.wide, aes(bsl, fup, color = factor(trt))) +
  geom_point()

# res1.2 = future_replicate(10000, sim(n, mu0, mu1, Sigma), future.seed = 0)
summarize(res1.2, es)  # outcome bad, change not good

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.38 0.38  0.51
change from baseline                                                 0        0.31 0.31  0.63
ANCOVA with one slope                                                0        0.28 0.28  0.73
ANCOVA with two slopes                                               0        0.28 0.28  0.73
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.28 0.28  0.72
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.28 0.28  0.72

# res0.2 = future_replicate(10000, sim(n, mu0, mu0, Sigma), future.seed = 0)
summarize(res0.2, 0)

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.40 0.40  0.05
change from baseline                                                 0        0.32 0.32  0.05
ANCOVA with one slope                                                0        0.30 0.30  0.05
ANCOVA with two slopes                                               0        0.30 0.30  0.05
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.30 0.30  0.05
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.30 0.30  0.05

n = 1000000
mu0 = c(5,7)
mu1 = c(5,8)
Sigma = matrix(c(3,1,1,3),2)
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
round(colMeans(dat[1:(n/2),]),2)  # means bsl/fup for trt0

[1] 8.62 9.31

round(cov(dat[1:(n/2),]), 2)  # cov matrix bsl/fup for trt0

     [,1] [,2]
[1,] 0.66 0.12
[2,] 0.12 0.22

round(colMeans(dat[n/2+1:(n/2),]),2)  # means bsl/fup for trt1

[1] 8.62 9.52

round(cov(dat[n/2+1:(n/2),]),2)  # cov matrix bsl/fup for trt1

     [,1] [,2]
[1,] 0.67 0.08
[2,] 0.08 0.11

es = mean(dat[n/2+1:(n/2),2])-mean(dat[1:(n/2),2])
round(es, 2)  # effect size of trt (difference in mean fup)

[1] 0.21

n = 50
dat = rbind(mvrnorm(n/2, mu0, Sigma), mvrnorm(n/2, mu1, Sigma)) |>
  (\(x) tanh(x/5)*5+5)()
db.wide = data.frame(pat = 1:n, bsl = dat[,1], fup = dat[,2], trt = rep(0:1, each = n/2))
ggplot(db.wide, aes(bsl, fup, color = factor(trt))) +
  geom_point()

# res1.3 = future_replicate(10000, sim(n, mu0, mu1, Sigma), future.seed = 0)
summarize(res1.3, es)  # change and longitudinal bad (2% singular), ANCOVA best

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.12 0.12  0.46
change from baseline                                                 0        0.22 0.22  0.15
ANCOVA with one slope                                                0        0.11 0.11  0.49
ANCOVA with two slopes                                               0        0.11 0.11  0.50
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.11 0.11  0.12
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.11 0.11  0.11

# res0.3 = future_replicate(10000, sim(n, mu0, mu0, Sigma), future.seed = 0)
summarize(res0.3, 0)  # longitudinal conservative (2% singular)

                                                                  bias SD estimate RMSE power
outcome without considering baseline                                 0        0.13 0.13  0.05
change from baseline                                                 0        0.22 0.22  0.05
ANCOVA with one slope                                                0        0.13 0.13  0.05
ANCOVA with two slopes                                               0        0.13 0.13  0.05
mixed-eff longitudinal with one baseline and Satterthwaite approx    0        0.13 0.13  0.01
mixed-eff longitudinal with one baseline and Kenward-Roger approx    0        0.13 0.13  0.01

# save(res1.1, res0.1, res1.2, res0.2, res1.3, res0.3, file = "res.rda")

conclusions

fup ~ trt (outcome without considering baseline):
- loss of statistical power when high correlation between baseline and follow-up ❌
- issues with unbalanced baseline ❌
- compatibility with nonparametric analysis ✔️
fup - bsl ~ trt (change from baseline):
- loss of statistical power when low correlation between baseline and follow-up ❌
- issues with unbalanced baseline (ceiling/floor effect) ❌
- typically more interesting effect size ✔️
- compatibility with nonparametric analysis ✔️
fup ~ trt +/* scale(bsl) (ANCOVA with or without interaction):
- simple model, quite robust to deviations from assumptions ✔️
- typically more interesting effect size as adjusted analysis of change ✔️
y ~ trt.b + (1|pat) (mixed-effect longitudinal with one baseline):
- more complex model, not so robust to deviations from assumptions ❌
- partial observations possible (but limited practical utility) ✔️