Load Data

data <- read.csv("data/trab_clean.csv")

Prepare Long Format for Modeling

iop_long_model <- data %>%
  dplyr::select(Patient_ID, Age, surgery_clean, iop_preop, iop_1dpo, iop_6mo, iop_1yr) %>%
  tidyr::pivot_longer(
    cols = starts_with("iop_"),
    names_to = "timepoint",
    values_to = "IOP"
  ) %>%
  dplyr::filter(!is.na(IOP)) %>%
  dplyr::mutate(
    timepoint = factor(timepoint,
      levels = c("iop_preop", "iop_1dpo", "iop_6mo", "iop_1yr", "iop_2yr"),
      labels = c("Pre-op", "1DPO", "6MO", "1YR", "2YR")
    )
  )

Fit Linear Mixed Model

# LMM: IOP as outcome, timepoint as fixed effect, random intercept for each patient
model_iop_lmm <- lmer(IOP ~ timepoint + (1 | Patient_ID), data = iop_long_model)

# Summary of fixed effects
summary(model_iop_lmm)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: IOP ~ timepoint + (1 | Patient_ID)
##    Data: iop_long_model
## 
## REML criterion at convergence: 566.2
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.3078 -0.5242 -0.1992  0.2311  3.9955 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Patient_ID (Intercept) 30.68    5.539   
##  Residual               56.59    7.523   
## Number of obs: 81, groups:  Patient_ID, 30
## 
## Fixed effects:
##               Estimate Std. Error      df t value Pr(>|t|)    
## (Intercept)     24.489      1.758  65.356  13.927  < 2e-16 ***
## timepoint1DPO  -13.799      2.178  53.650  -6.336 5.03e-08 ***
## timepoint6MO   -10.693      2.208  56.591  -4.842 1.03e-05 ***
## timepoint1YR    -6.615      3.048  58.566  -2.170   0.0341 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tm1DPO tmp6MO
## timepnt1DPO -0.528              
## timepont6MO -0.539  0.424       
## timepont1YR -0.378  0.306  0.321

Plot Model Predictions

emm <- emmeans(model_iop_lmm, ~ timepoint)
plot(emm, comparisons = TRUE) +
  labs(
    title = "Estimated Mean IOP with 95% CI",
    x = "Timepoint",
    y = "Predicted IOP (mmHg)"
  ) +
  theme_minimal()

Add Covariates to the IOP Model

Let’s test if surgery type or age influences IOP trajectory:

# Linear Mixed Model with Surgery Type and Age
model_cov <- lmer(IOP ~ timepoint + surgery_clean + Age + (1 | Patient_ID), data = iop_long_model)

summary(model_cov)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: IOP ~ timepoint + surgery_clean + Age + (1 | Patient_ID)
##    Data: iop_long_model
## 
## REML criterion at convergence: 559.7
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -1.4471 -0.4695 -0.1506  0.3112  3.8888 
## 
## Random effects:
##  Groups     Name        Variance Std.Dev.
##  Patient_ID (Intercept) 25.61    5.061   
##  Residual               55.46    7.447   
## Number of obs: 81, groups:  Patient_ID, 30
## 
## Fixed effects:
##                              Estimate Std. Error        df t value Pr(>|t|)    
## (Intercept)                  24.37795    4.76457  32.29480   5.117 1.38e-05 ***
## timepoint1DPO               -13.94130    2.15395  55.09901  -6.472 2.73e-08 ***
## timepoint6MO                -10.59871    2.19573  57.57542  -4.827 1.06e-05 ***
## timepoint1YR                 -7.01505    3.01741  59.66518  -2.325   0.0235 *  
## surgery_cleanTrabeculectomy  -4.59665    3.08361  29.25474  -1.491   0.1468    
## Age                           0.07111    0.07094  30.12217   1.003   0.3241    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) tm1DPO tmp6MO tmp1YR srgr_T
## timepnt1DPO -0.202                            
## timepont6MO -0.086  0.420                     
## timepont1YR -0.081  0.306  0.324              
## srgry_clnTr -0.746  0.023 -0.105 -0.014       
## Age         -0.861 -0.001 -0.100 -0.079  0.501

KM Analysis

Define Failure and Time-to-Event Variables - Assume we define failure as: - IOP > 21 or - Meds restarted (not med-free).

You’ll need these columns:

# Create survival variables
 data <- data %>%
  mutate(
    # Use last known IOP and meds status for defining failure
    time_to_event = case_when(
      !is.na(iop_6mo) & iop_6mo > 21 ~ 6,
      !is.na(iop_1yr) & iop_1yr > 21 ~ 12,
      !is.na("Current.medication") & "Current.medication" == 1 ~ 18,
      TRUE ~ NA_real_
    ),
    status = if_else(!is.na(time_to_event), 1L, 0L),
    time_to_event = if_else(is.na(time_to_event), 18, time_to_event)
  )

Build and Plot Kaplan-Meier Curve

# Create the survival object
km_object <- Surv(time = data$time_to_event, event = data$status)

# Fit the Kaplan-Meier survival curve
km_fit <- survfit(Surv(time_to_event, status) ~ 1, data = data)

km_fit
## Call: survfit(formula = Surv(time_to_event, status) ~ 1, data = data)
## 
##       n events median 0.95LCL 0.95UCL
## [1,] 33      4     NA      NA      NA
# Plot
km_fit %>%
  ggsurvfit() +
  labs(
    title = "Kaplan-Meier Curve: Time to Failure",
    x = "Months Since Surgery",
    y = "Probability of Success"
  ) +
  theme_minimal()

Save Cleaned + Long Data

#Export cleaned files for reproducibility:
write.csv(data, "data/trab_survival.csv", row.names = FALSE)
write.csv(iop_long_model, "data/iop_long.csv", row.names = FALSE)