Notes on the paper by Guo, Xu, and Bradley P. Carlin, “Separate and Joint Modeling of Longitudinal and Event Time Data Using Standard Computer Packages” (2004).

Overview of Joint Models

What is a Joint Model?

A joint model simultaneously analyzes:

Longitudinal outcomes (repeated measurements over time, e.g., CD4 counts)
Survival outcomes (time-to-event data, e.g., time to death)

The key insight is that these two processes are interconnected:

The longitudinal trajectory can predict or influence the risk of an event
The event (death) causes dropout from the longitudinal study
Shared subject-specific characteristics affect both outcomes

Why Not Analyze Separately?

Analysis Type	Problem
Longitudinal only	Ignores informative dropout; assumes missing data is random
Survival only	Ignores time-varying biomarkers
Two-stage approach	Underestimates uncertainty; introduces bias

The Guo-Carlin (2004) Model

Background

Data: AIDS clinical trial comparing ddI vs ddC
Longitudinal: CD4 cell counts over time
Survival: Time to death
Joint modeling approach: Henderson et al. (2000) framework with bivariate Gaussian random effects

Model Specification

Longitudinal Submodel

\[ y_{ij} = \mathbf{x}_{1i}^T(s_{ij})\boldsymbol{\beta}_1 + U_{1i} + U_{2i}s_{ij} + \epsilon_{ij} \]

\[ \epsilon_{ij} \sim N(0, \sigma_\epsilon^2) \]

where: - \(y_{ij}\) = CD4 count for subject \(i\) at visit \(j\) - \(\mathbf{x}_{1i}^T(s_{ij})\boldsymbol{\beta}_1\) = fixed effects (includes drug, time, drug×time) - \(U_{1i}\) = random intercept (baseline deviation) - \(U_{2i}\) = random slope (rate of change) - \(\epsilon_{ij}\) = measurement error

Understanding \(s_{ij}\) in the Longitudinal Model

What is \(s_{ij}\)?

\(s_{ij}\) represents the time of the \(j\)-th measurement for subject \(i\). The notation uses \(s\) (instead of \(t\)) to denote measurement times (longitudinal observations), while \(t\) is reserved for event times (survival outcomes). \(s_{ij}\) is the actual time point (e.g., in months) when the \(j\)-th CD4 measurement was taken for subject \(i\).

Why use \(s\) instead of \(t\)?

Distinguishes measurement times from event times:
- \(s\) = time of longitudinal measurement (e.g., clinic visit)
- \(t\) = time of event (death or censoring)
Different scales: Longitudinal measurements (\(s_{ij}\)) can occur at any time, but are typically measured at scheduled visits (e.g., months 0, 3, 6, 12). Event times (\(t_i\)) can occur at any time, not necessarily at visit times.
Flexible evaluation: The longitudinal trajectory can be evaluated at any time \(s\), while the survival model uses time \(t\) to model the hazard. The association structure links them through functions like \(W_{2i}(s)\) or \(W_{2i}(t)\).

Example:

If subject \(i\) has CD4 measurements at months 0, 3, 6, and 12: - \(s_{i1} = 0\) (baseline visit) - \(s_{i2} = 3\) (3-month visit) - \(s_{i3} = 6\) (6-month visit) - \(s_{i4} = 12\) (12-month visit)

Then the longitudinal submodel for this subject would be: - \(y_{i1} = \mathbf{x}_{1i}^T(0)\boldsymbol{\beta}_1 + U_{1i} + U_{2i}(0) + \epsilon_{i1}\) (Baseline CD4) - \(y_{i2} = \mathbf{x}_{1i}^T(3)\boldsymbol{\beta}_1 + U_{1i} + U_{2i}(3) + \epsilon_{i2}\) (3-month CD4) - \(y_{i3} = \mathbf{x}_{1i}^T(6)\boldsymbol{\beta}_1 + U_{1i} + U_{2i}(6) + \epsilon_{i3}\) (6-month CD4) - \(y_{i4} = \mathbf{x}_{1i}^T(12)\boldsymbol{\beta}_1 + U_{1i} + U_{2i}(12) + \epsilon_{i4}\) (12-month CD4)

The Role of \(s_{ij}\) in the Association Structure:

The true trajectory \(U_{1i} + U_{2i}s_{ij}\) at measurement time \(s_{ij}\) enters the survival model through the association structure: \[ W_{2i}(t) = \gamma_1 U_{1i} + \gamma_2 U_{2i} + \gamma_3 (U_{1i} + U_{2i}t) \]

This links the longitudinal measurements (taken at times \(s_{ij}\)) to the survival outcome (at time \(t\)).

Key Takeaway: The observed CD4 value \(y_{ij}\) is used to estimate the random effects (\(U_{1i}, U_{2i}\)), but the survival model uses the true trajectory \(U_{1i} + U_{2i}t\) evaluated at any time \(t\), not the noisy observed values \(y_{ij}\).

Survival Submodel

\[ t_i \sim \text{Weibull}(p, \mu_i(t)) \]

\[ \log(\mu_i(t)) = \mathbf{x}_{2i}^T(t)\boldsymbol{\beta}_2 + W_{2i}(t) \]

where: - \(t_i\) = time to death - \(p > 0\) = Weibull shape parameter - \(\mathbf{x}_{2i}^T(t)\boldsymbol{\beta}_2\) = fixed effects (includes drug, age, etc.) - \(W_{2i}(t)\) = association structure

Association Structure

\[ W_{2i}(t) = \gamma_1 U_{1i} + \gamma_2 U_{2i} + \gamma_3 (U_{1i} + U_{2i}t) \]

Random Effects

\[ \begin{pmatrix} U_{1i} \\ U_{2i} \end{pmatrix} \sim N\left( \begin{pmatrix} 0 \\ 0 \end{pmatrix}, \begin{pmatrix} \sigma_1^2 & \sigma_{12} \\ \sigma_{12} & \sigma_2^2 \end{pmatrix} \right) \]

How Are the Longitudinal and Survival Models Linked?

The Role of \(y_{ij}\) in the Survival Model

Important distinction: \(y_{ij}\) (the observed CD4 measurement) does NOT directly enter the survival model!

Instead: - The survival model uses the underlying true trajectory (not the observed values with measurement error) - The true trajectory is: \(U_{1i} + U_{2i}s_{ij}\) (the random intercept + random slope component) - This is then weighted by the association parameters \(\gamma_1, \gamma_2, \gamma_3\)

The Two-Stage Approach for Model XI

Since brms doesn’t easily allow excluding \(\gamma_3\) while keeping \(\gamma_1\) and \(\gamma_2\) in a joint model, we use a two-stage approach:

Stage 1: Fit the longitudinal model to estimate the random effects (\(U_{1i}, U_{2i}\)) Stage 2: Use the estimated random effects as covariates in the survival model

This gives us Model XI: \(W_{2i}(t) = \gamma_1 U_{1i} + \gamma_2 U_{2i}\) (no \(\gamma_3\))

Summary: The Flow of Information

Observed data: \(y_{ij}\) (CD4 measurements at times \(s_{ij}\))
Estimates: Random effects \(U_{1i}\) and \(U_{2i}\) from the longitudinal model
Prediction: True trajectory \(U_{1i} + U_{2i}t\) at any time \(t\)
Link: The true trajectory enters the survival model through \(W_{2i}(t)\)
Survival: \(W_{2i}(t)\) affects the hazard of death

Model Comparison: Full vs Model XI

The “Model XI” Naming Convention

Why is it called “Model XI”?

In the original Guo and Carlin (2004) paper, they systematically compared 11 different joint models (Models I through XI) with varying association structures.

Model XI (Roman numeral for 11) was the eleventh and final model they considered, which had the most parsimonious association structure.

The Evolution of Models I-XI

In their paper, Guo and Carlin systematically tested different ways to link the longitudinal and survival processes:

Model	Association Structure	Description
Models I-V	Various combinations	Testing different functional forms
Models VI-X	More complex structures	Including current value, slope, etc.
Model XI	\(\gamma_1 U_{1i} + \gamma_2 U_{2i}\)	Only random effects (intercept + slope)

Why Model XI Was Special

Model XI was the most restrictive model because: 1. It excluded \(\gamma_3\) (the current CD4 value effect) 2. It assumed that only the underlying trajectory (baseline deviation + rate of change) matters for survival 3. The current CD4 measurement at event time provides no additional prognostic information

The Key Finding

Despite being the simplest model, Model XI had the lowest DIC (Deviance Information Criterion), meaning:

“The extra parameter \(\gamma_3\) does not improve fit enough to justify the complexity”

This was a clinically important finding because it suggested that: - Doctors don’t need to know the current CD4 value to predict survival - The patient’s trajectory (intercept and slope) contains all the prognostic information - This simplifies clinical monitoring

Full Model (with \(U_{3i} = 0\))

\[ W_{2i}(t) = \gamma_1 U_{1i} + \gamma_2 U_{2i} + \gamma_3 (U_{1i} + U_{2i}t) \]

Association Parameters:

\(\gamma_1\): Effect of random intercept (baseline deviation) on survival
\(\gamma_2\): Effect of random slope (rate of change) on survival
\(\gamma_3\): Effect of fitted longitudinal value at event time on survival

Model XI (\(\gamma_3 = 0\))

\[ W_{2i}(t) = \gamma_1 U_{1i} + \gamma_2 U_{2i} \]

Key Difference: - Model XI assumes the current CD4 value does NOT directly affect survival beyond the patient’s underlying trajectory - Only baseline deviation and rate of change matter

Which Model Fits Better?

According to Guo and Carlin (2004), Model XI has lower DIC, indicating:

The extra parameter \(\gamma_3\) does not improve fit enough to justify the complexity
Current CD4 value provides no additional prognostic information beyond the trajectory

Complete R Code with Simulated Data using brms

# ============================================================
# COMPLETE JOINT MODELING WITH BRMS
# Guo-Carlin (2004) Model Implementation
# ============================================================

# Install required packages (if not already installed)
# install.packages(c("brms", "dplyr", "tidyr", "ggplot2", "survival", "MASS", "loo", "bayesplot"))

library(brms)
library(dplyr)
library(tidyr)
library(ggplot2)
library(survival)
library(MASS)
library(loo)        # For model comparison
library(bayesplot)  # For diagnostic plots

# Set seed for reproducibility
set.seed(2026)

# Set brms options
options(mc.cores = parallel::detectCores())

# ============================================================
# STEP 1: SIMULATE DATA
# ============================================================

simulate_joint_data <- function(n = 100, n_visits = 5, max_time = 36) {
  
  # ---------- Parameters ----------
  # Longitudinal model
  beta10 <- 200   # Baseline CD4 (ddC group)
  beta11 <- -15   # Time effect (monthly decline)
  beta12 <- 30    # Drug effect (ddI vs ddC, baseline difference)
  beta13 <- 2     # Drug × Time interaction (ddI has slower decline)
  sigma_e <- 30   # Measurement error
  
  # Survival model
  beta20 <- -5    # Baseline log hazard
  beta21 <- 0.3   # Drug effect on survival
  beta22 <- 0.01  # Age effect
  p_shape <- 1.5  # Weibull shape parameter
  
  # Random effects parameters
  sigma1 <- 40    # SD of random intercept
  sigma2 <- 8     # SD of random slope
  rho <- -0.5     # Correlation between intercept and slope
  
  # Association parameters
  gamma1 <- -0.01 # Intercept effect on survival
  gamma2 <- 0.05  # Slope effect on survival
  gamma3 <- -0.02 # Current value effect on survival
  
  # Create covariance matrix for random effects
  Sigma <- matrix(c(sigma1^2, rho*sigma1*sigma2,
                    rho*sigma1*sigma2, sigma2^2), nrow = 2)
  
  # ---------- Generate Subject-Level Data ----------
  id <- 1:n
  drug <- rbinom(n, 1, 0.5)
  
  # Generate TRUE random effects
  U <- MASS::mvrnorm(n, mu = c(0, 0), Sigma = Sigma)
  U1 <- U[, 1]
  U2 <- U[, 2]
  
  age <- round(rnorm(n, 40, 10), 1)
  
  # ---------- Step 1: Generate Survival Times ----------
  surv_data <- data.frame()
  
  for(i in 1:n) {
    log_mu <- beta20 + beta21*drug[i] + beta22*age[i] +
              gamma1*U1[i] + gamma2*U2[i] + gamma3*(U1[i] + U2[i]*max_time)
    
    mu <- exp(log_mu)
    U_surv <- runif(1)
    death_time <- (-log(U_surv) / mu)^(1/p_shape)
    
    censored <- death_time > max_time
    event_time <- min(death_time, max_time)
    
    surv_data <- rbind(surv_data, data.frame(
      id = i,
      drug = drug[i],
      age = age[i],
      survtime = event_time,
      status = as.integer(!censored)
    ))
  }
  
  # ---------- Step 2: Generate Longitudinal Data ----------
  long_data <- data.frame()
  
  for(i in 1:n) {
    event_time <- surv_data$survtime[i]
    candidate_times <- sort(runif(n_visits * 3, 0, max_time))
    visit_times <- c(0, candidate_times[candidate_times < event_time & candidate_times > 0])
    
    if(length(visit_times) < 2) {
      visit_times <- c(0, min(event_time * 0.5, max_time * 0.1))
    }
    
    visit_times <- sort(unique(visit_times))[1:min(n_visits, length(visit_times))]
    
    for(j in 1:length(visit_times)) {
      time <- visit_times[j]
      mu_fixed <- beta10 + beta11*time + beta12*drug[i] + beta13*drug[i]*time
      U_part <- U1[i] + U2[i]*time
      epsilon <- rnorm(1, 0, sigma_e)
      y <- mu_fixed + U_part + epsilon
      
      long_data <- rbind(long_data, data.frame(
        id = i,
        drug = drug[i],
        age = age[i],
        time = time,
        CD4 = y
      ))
    }
  }
  
  # Store TRUE values for validation
  true_values <- data.frame(
    id = 1:n,
    U1_true = U1,
    U2_true = U2
  )
  
  return(list(
    long_data = long_data,
    surv_data = surv_data,
    true_values = true_values,
    true_params = list(
      beta10 = beta10, beta11 = beta11, beta12 = beta12, beta13 = beta13,
      beta20 = beta20, beta21 = beta21, beta22 = beta22,
      sigma1 = sigma1, sigma2 = sigma2, rho = rho,
      gamma1 = gamma1, gamma2 = gamma2, gamma3 = gamma3,
      sigma_e = sigma_e, p_shape = p_shape
    )
  ))
}

# Generate the data
sim_data <- simulate_joint_data(n = 150, n_visits = 6)
long_data <- sim_data$long_data
surv_data <- sim_data$surv_data
true_values <- sim_data$true_values
true_params <- sim_data$true_params

# longitudinal data structure
print(head(long_data))

##   id drug  age       time       CD4
## 1  1    1 38.2 0.00000000 291.91785
## 2  1    1 38.2 0.05082128 285.57667
## 3  2    1 44.9 0.00000000 171.22451
## 4  2    1 44.9 0.05760438 125.19121
## 5  2    1 44.9 2.40817487  94.73255
## 6  2    1 44.9 6.52086902 154.90172

# survival data structure 
print(head(surv_data))

##   id drug  age   survtime status
## 1  1    1 38.2  0.1016426      1
## 2  2    1 44.9 36.0000000      0
## 3  3    0 26.4 36.0000000      0
## 4  4    0 36.3  6.5477075      1
## 5  5    1 47.6  2.3593613      1
## 6  6    0 41.5 20.2602790      1

# Validation check
merged_check <- long_data %>%
  left_join(surv_data, by = c("id", "drug", "age")) %>%
  group_by(id) %>%
  summarise(
    max_time = max(time),
    event_time = first(survtime),
    valid = max_time <= event_time
  )

print(head(merged_check))

## # A tibble: 6 × 4
##      id max_time event_time valid
##   <int>    <dbl>      <dbl> <lgl>
## 1     1   0.0508      0.102 TRUE 
## 2     2   7.78       36     TRUE 
## 3     3   6.29       36     TRUE 
## 4     4   6.46        6.55  TRUE 
## 5     5   1.91        2.36  TRUE 
## 6     6   8.11       20.3   TRUE

# note, here max_time is the last measurement time,  event_time is death time. last measure time should always less than event time.

# ============================================================
# STEP 2: EXPLORATORY PLOTS
# ============================================================

# Longitudinal trajectories
ggplot(long_data, aes(x = time, y = CD4, color = factor(drug), group = id)) +
  geom_line(alpha = 0.3) +
  geom_smooth(aes(group = factor(drug)), method = "lm", se = FALSE, size = 1.5) +
  labs(x = "Time (months)", y = "CD4 Count",
       title = "Longitudinal CD4 Trajectories by Treatment",
       color = "Treatment (0=ddC, 1=ddI)") +
  theme_minimal()

# Kaplan-Meier survival curves
km_fit <- survfit(Surv(survtime, status) ~ drug, data = surv_data)
plot(km_fit, col = c("blue", "red"), lwd = 2,
     xlab = "Time (months)", ylab = "Survival Probability",
     main = "Kaplan-Meier Survival Curves by Treatment")
legend("topright", legend = c("ddC", "ddI"), col = c("blue", "red"), lwd = 2)

# ============================================================
# STEP 3: FIT MODELS USING TWO-STAGE APPROACH WITH BRMS
# ============================================================

# Fitting models using two-stage approach with brms#

# ============================================================
# STAGE 1: Fit Longitudinal Model
# ============================================================

#STAGE 1: Fitting Longitudinal Model

fit_long <- brm(
  CD4 ~ drug * time + (time | id),
  data = long_data,
  family = gaussian(),
  prior = c(
    prior(normal(0, 100), class = "b"),
    prior(cauchy(0, 10), class = "sd"),
    prior(cauchy(0, 10), class = "sigma")
  ),
  chains = 2,
  iter = 1000,
  warmup = 500,
  seed = 123,
  refresh = 0
)


# ============================================================
# STAGE 2: Extract Random Effects for Model XI
# ============================================================

# Extract the random effects (U1 and U2) from the longitudinal model
ranef_long <- ranef(fit_long)
ranef_df <- as.data.frame(ranef_long$id)

# Create subject-level dataset with random effects
subject_data <- data.frame(
  id = as.numeric(rownames(ranef_df)),
  U1_est = ranef_df$Estimate.Intercept,
  U2_est = ranef_df$Estimate.time
)

# Merge with survival data
surv_data_with_re <- surv_data %>%
  left_join(subject_data, by = "id")

# ============================================================
# STAGE 3: Fit Full Model (Joint) with brms
# ============================================================

# Full model using mvbrmsformula with correlated random effects
joint_data <- long_data %>%
  left_join(surv_data, by = c("id", "drug", "age"))

fit_full <- brm(
  mvbrmsformula(
    CD4 ~ drug * time + (time | id),  # Longitudinal model
    survtime | cens(status) ~ drug + age + (1 | id)  # Survival model
  ),
  data = joint_data,
  family = list(gaussian(), brmsfamily("weibull")),
  prior = c(
    prior(normal(0, 100), class = "b", resp = "CD4"),
    prior(normal(0, 100), class = "b", resp = "survtime"),
    prior(cauchy(0, 10), class = "sd", resp = "CD4"),
    prior(cauchy(0, 10), class = "sd", resp = "survtime"),
    prior(cauchy(0, 10), class = "sigma", resp = "CD4")
  ),
  chains = 2,
  iter = 1000,
  warmup = 500,
  seed = 123,
  refresh = 0
)

# ============================================================
# STAGE 4: Fit Model XI (gamma3 = 0) using Two-Stage Approach
# ============================================================

# Model XI: Survival with estimated random effects (U1 and U2)
# This gives us: W(t) = gamma1*U1 + gamma2*U2 (no gamma3)
fit_xi <- brm(
  survtime | cens(status) ~ drug + age + U1_est + U2_est,
  data = surv_data_with_re,
  family = brmsfamily("weibull"),
  prior = c(
    prior(normal(0, 100), class = "b"),
    prior(cauchy(0, 10), class = "shape")
  ),
  chains = 2,
  iter = 1000,
  warmup = 500,
  seed = 123,
  refresh = 0
)

# ============================================================
# STEP 4: MODEL COMPARISON
# ============================================================

# Extract WAIC for all models
waic_full <- waic(fit_full)
waic_xi <- waic(fit_xi)

# WAIC Values

#Full Model (Joint)

round(waic_full$estimates["waic", "Estimate"], 2)

## [1] 4756.51

# Model XI (U1+U2): 

round(waic_xi$estimates["waic", "Estimate"], 2)

## [1] -577.74

# Compare Full vs Model XI
diff_xi_full <- waic_xi$estimates["waic", "Estimate"] - waic_full$estimates["waic", "Estimate"]
cat("\n\nDifference (XI - Full):", round(diff_xi_full, 2))

## 
## 
## Difference (XI - Full): -5334.25

## 
## 
## --- CONCLUSION ---

## 
## ✅ Model XI fits better (lower WAIC).
##    The current CD4 value does NOT add prognostic information.
##    ✅ This matches Guo & Carlin's (2004) finding!

# ============================================================
# STEP 5: RESULTS SUMMARY
# ============================================================

full_summary <- summary(fit_full)
print(full_summary)

##  Family: MV(gaussian, weibull) 
##   Links: mu = identity
##          mu = log 
## Formula: CD4 ~ drug * time + (time | id) 
##          survtime | cens(status) ~ drug + age + (1 | id) 
##    Data: joint_data (Number of observations: 635) 
##   Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
##          total post-warmup draws = 1000
## 
## Multilevel Hyperparameters:
## ~id (Number of levels: 150) 
##                             Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS
## sd(CD4_Intercept)               0.81      0.49     0.32     1.30 2.65        2
## sd(CD4_time)                    2.93      1.10     1.82     4.04 2.74        2
## sd(survtime_Intercept)          0.00      0.00     0.00     0.00 2.24        3
## cor(CD4_Intercept,CD4_time)    -0.03      0.91    -0.95     0.88 2.55        2
##                             Tail_ESS
## sd(CD4_Intercept)                 14
## sd(CD4_time)                      13
## sd(survtime_Intercept)            12
## cor(CD4_Intercept,CD4_time)       19
## 
## Regression Coefficients:
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## CD4_Intercept          4.17      3.85     0.32     8.02 1.90        3       11
## survtime_Intercept     3.58      0.00     3.58     3.58 1.00     1052      639
## CD4_drug              -0.04      0.68    -0.72     0.65 2.20        3       30
## CD4_time              -0.67      0.85    -1.52     0.18 2.04        3       11
## CD4_drug:time         -0.23      0.08    -0.32    -0.14 2.12        3       12
## survtime_drug          0.00      0.00    -0.00     0.00 1.01      769      610
## survtime_age           0.00      0.00    -0.00     0.00 1.01      989      606
## 
## Further Distributional Parameters:
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_CD4        188.97      1.29   187.40   190.48 2.43        2       22
## shape_survtime 27186.22   1572.06 25329.93 30159.02 1.85        3       16
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

# Model XI summary
#Model XI (U1+U2) ---\n")
xi_summary <- summary(fit_xi)
print(xi_summary)

##  Family: weibull 
##   Links: mu = log 
## Formula: survtime | cens(status) ~ drug + age + U1_est + U2_est 
##    Data: surv_data_with_re (Number of observations: 150) 
##   Draws: 2 chains, each with iter = 1000; warmup = 500; thin = 1;
##          total post-warmup draws = 1000
## 
## Regression Coefficients:
##           Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept     3.58      0.00     3.58     3.58 1.02      351      124
## drug          0.00      0.00    -0.00     0.00 1.02      132      268
## age           0.00      0.00    -0.00     0.00 1.02      940      460
## U1_est       -0.00      0.00    -0.00     0.00 1.02      801      755
## U2_est        0.00      0.00    -0.00     0.00 1.01      644      362
## 
## Further Distributional Parameters:
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## shape 20263.07   3344.28 12288.56 24259.01 1.65        3       11
## 
## Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

# Extract Model XI parameters
xi_fixed <- posterior_summary(fit_xi)

print(xi_fixed)

##                  Estimate    Est.Error          Q2.5         Q97.5
## b_Intercept  3.583489e+00 3.064138e-05  3.583430e+00  3.583553e+00
## b_drug       6.436575e-07 1.403062e-05 -2.705986e-05  2.825057e-05
## b_age        2.910676e-08 6.451974e-07 -1.304955e-06  1.265246e-06
## b_U1_est    -1.564041e-11 2.247034e-07 -4.275241e-07  4.429984e-07
## b_U2_est     1.271203e-07 1.805156e-06 -3.343061e-06  3.847305e-06
## shape        2.026307e+04 3.344284e+03  1.228856e+04  2.425901e+04
## Intercept    3.583491e+00 1.001814e-05  3.583470e+00  3.583507e+00
## lprior      -4.218700e+01 3.623932e-01 -4.257766e+01 -4.121742e+01
## lp__         2.576683e+02 9.878451e+00  2.318615e+02  2.693245e+02

# Find specific parameters
param_names <- rownames(xi_fixed)
print(param_names)

## [1] "b_Intercept" "b_drug"      "b_age"       "b_U1_est"    "b_U2_est"   
## [6] "shape"       "Intercept"   "lprior"      "lp__"

# ============================================================
# STEP 6: DIAGNOSTIC PLOTS 
# ============================================================

# Extract R-hat for all parameters (fixed + random)
# Option A: Using the summary object directly
rhat_full <- full_summary$fixed[, "Rhat"]
rhat_full

## [1] 1.900209 1.001955 2.196299 2.041891 2.121594 1.005552 1.014508

rhat_xi <- xi_summary$fixed[, "Rhat"]
rhat_xi

## [1] 1.020319 1.023870 1.023475 1.015814 1.006959

# ============================================================
# TRACE PLOTS
# ============================================================

# For Model XI - get parameter names
xi_param_names <- rownames(posterior_summary(fit_xi))
drug_params <- grep("drug", xi_param_names, value = TRUE)
u1_params <- grep("U1_est", xi_param_names, value = TRUE)
u2_params <- grep("U2_est", xi_param_names, value = TRUE)

params_to_plot <- c(drug_params, u1_params, u2_params)
params_to_plot <- params_to_plot[!is.na(params_to_plot)]

if (length(params_to_plot) > 0) {
  cat("Plotting parameters:", paste(params_to_plot, collapse=", "), "\n")
  # Use plot method for brms objects
  plot(fit_xi, pars = params_to_plot[1:min(4, length(params_to_plot))])
}

## Plotting parameters: b_drug, b_U1_est, b_U2_est

# For Full Model - plot key parameters
full_param_names <- rownames(posterior_summary(fit_full))
full_drug_params <- grep("drug", full_param_names, value = TRUE)
# Get all drug parameter names
# View all of them

print(full_drug_params)

## [1] "b_CD4_drug"      "b_CD4_drug:time" "b_survtime_drug"

# ============================================================
# POSTERIOR PREDICTIVE CHECKS
# ============================================================
pp_check(fit_full, resp = "CD4")

pp_check(fit_xi)

INTERPRETATION OF DRUG EFFECTS

Longitudinal Model (\(\beta_{12}\), \(\beta_{13}\))

\(\beta_{12}\) = Baseline CD4 difference (ddI vs ddC)
\(\beta_{13}\) = CD4 slope difference (ddI vs ddC)
- If \(\beta_{13} > 0\): ddI has slower CD4 decline than ddC
- If \(\beta_{13} < 0\): ddI has faster CD4 decline than ddC

Survival Model (\(\beta_{21}\))

\(\beta_{21}\) = Log-hazard ratio (ddI vs ddC)
- If \(\beta_{21} > 0\): ddI has higher mortality hazard than ddC
- If \(\beta_{21} < 0\): ddI has lower mortality hazard than ddC
- Hazard ratio \(= \exp(\beta_{21})\)

Association Parameters

\[ \begin{aligned} \gamma_1 \text{ (intercept)} &= \text{Coefficient of } U1_{\text{est}} \\ \gamma_2 \text{ (slope)} &= \text{Coefficient of } U2_{\text{est}} \\ \gamma_3 \text{ (current value)} &= 0 \text{ in Model XI (excluded)} \end{aligned} \]

Key Takeaways

Joint Model Summary

Component	Description
Longitudinal model	Linear mixed effects model for CD4 trajectories
Survival model	Weibull model for time to death
Association	Shared random effects (\(U_{1i}, U_{2i}\)) with parameters \(\gamma_1, \gamma_2, \gamma_3\)
Drug effects	\(\beta_{12}, \beta_{13}\) (CD4), \(\beta_{21}\) (survival)

Model Comparison

Criterion	Full Model	Model XI (\(\gamma_3 = 0\))
Association	\(\gamma_1, \gamma_2, \gamma_3\)	\(\gamma_1, \gamma_2\)
Complexity	Higher	Lower
WAIC	Varies	Varies

The Two-Stage Approach for Model XI

Stage 1: Fit longitudinal model to estimate random effects

fit_long <- brm(CD4 ~ drug * time + (time | id), data = long_data)

Stage 2: Extract estimated random effects and use in survival model

# Extract U1 and U2 estimates
ranef_long <- ranef(fit_long)
U1_est <- ranef_long$id[, "Estimate.Intercept"]
U2_est <- ranef_long$id[, "Estimate.time"]

# Fit Model XI: No gamma3 (current value)
fit_xi <- brm(survtime | cens(status) ~ drug + age + U1_est + U2_est, 
              data = surv_data, family = brmsfamily("weibull"))

Why This Approach Works

True Model XI: Excludes \(\gamma_3\) (current value effect)
Only includes \(\gamma_1\) and \(\gamma_2\): Through U1 and U2
Can compare with Full Model: Using WAIC/LOOIC
Validates the hypothesis: If Model XI has lower WAIC, current CD4 value doesn’t add information

Joint Modeling of Longitudinal and Survival Data: Guo-Carlin (2004) Model using brms

2026-07-18

Overview of Joint Models

What is a Joint Model?

Why Not Analyze Separately?

The Guo-Carlin (2004) Model

Background

Model Specification

Longitudinal Submodel

Understanding \(s_{ij}\) in the Longitudinal Model

Survival Submodel

Association Structure

Random Effects

How Are the Longitudinal and Survival Models Linked?

The Role of \(y_{ij}\) in the Survival Model

The Two-Stage Approach for Model XI

Summary: The Flow of Information

Model Comparison: Full vs Model XI

The “Model XI” Naming Convention

The Evolution of Models I-XI

Why Model XI Was Special

The Key Finding

Full Model (with \(U_{3i} = 0\))

Model XI (\(\gamma_3 = 0\))

Which Model Fits Better?

Complete R Code with Simulated Data using brms

INTERPRETATION OF DRUG EFFECTS

Longitudinal Model (\(\beta_{12}\), \(\beta_{13}\))

Survival Model (\(\beta_{21}\))

Association Parameters

Key Takeaways

Joint Model Summary

Model Comparison

The Two-Stage Approach for Model XI

Why This Approach Works