1 Install the packages and upload the 3 datasets and join them

# Install Bioconductor package manager
if (!requireNamespace("BiocManager", quietly = TRUE)) {
  install.packages("BiocManager")
}

# Install Icens
BiocManager::install("Icens", ask = FALSE, update = FALSE)

library(tidyverse)
library(survival)
library(survminer)
library(lubridate)
library(icenReg)
library(interval)
library(Icens)


pregmpox<-read.csv("C:/Users/msaleeb/Downloads/DM-reports-PREGMPOX/Project isaac/Survival-New/pregmpox.csv", stringsAsFactors = T, na.strings = "",sep=",")
dates<-read.csv("C:/Users/msaleeb/Downloads/DM-reports-PREGMPOX/Project isaac/Survival-New/dates.csv", stringsAsFactors = T, na.strings = "",sep=",")
PCR<-read.csv("C:/Users/msaleeb/Downloads/DM-reports-PREGMPOX/Project isaac/Survival-New/PCR.csv",stringsAsFactors = T, na.strings = "",sep=",", 
              fileEncoding = "UTF-8-BOM",
              check.names = FALSE)



Trans_data<- dates|> left_join(PCR, by="id")
Final_data <- Trans_data|> left_join(pregmpox, by="id")

2 Clean the final dataset

Final_data$datvisite1<- dmy(Final_data$datvisite1)
Final_data$Start_date<-dmy(Final_data$Start_date)
Final_data$End_date<- dmy(Final_data$End_date)
Final_data$Start_date_episode2<- dmy(Final_data$Start_date_episode2)
Final_data$End_date_episode2<-dmy(Final_data$End_date_episode2)
Final_data$Interval_censoring_date1<- dmy(Final_data$Interval_censoring_date1)
Final_data$Interval_censoring_date2<- dmy(Final_data$Interval_censoring_date2)
Final_data$Last_seen_date<- dmy(Final_data$Last_seen_date)


Final_data$ndays_to_healing_first<- Final_data$End_date-Final_data$Start_date
Final_data$ndays_healing_second<- Final_data$End_date_episode2-Final_data$Start_date_episode2


Final_data$ndays_to_healing_first <- as.numeric(gsub(" days", "",Final_data$ndays_to_healing_first ))
Final_data$ndays_healing_second <- as.numeric(gsub(" days", "",Final_data$ndays_healing_second))



errors1<-Final_data|> filter(ndays_to_healing_first<=15)|>select(id, ndays_to_healing_first,Start_date, End_date)
errors2<-Final_data|> filter(ndays_to_healing_first>67)|>select(id, ndays_to_healing_first,Start_date, End_date)

Time to complete disappearance of mpox lesions was calculated from the date of lesion onset.

We have 3 groups of patients here in this dataset:


1-For patients with a documented exact healing date, the lower and upper bounds were identical and corresponded to the number of days between lesion onset and complete healing.
2-For interval-censored observations, the lower bound was defined as the time from lesion onset to the last assessment at which lesions were present, while the upper bound was defined as the time from lesion onset to the first assessment confirming complete healing.


3-For right-censored observations, the lower bound represented the time from lesion onset to the last assessment at which complete healing had not occurred (I extracted from redcap the last date when we saw the woman and it was reported that she is still having the lesions), and the upper bound was set to infinity.


These bounds will be used to estimate the survival function using the Turnbull nonparametric maximum likelihood estimator.

Final_data <- Final_data %>%
  mutate(
    censoring_type = case_when(
      Healed == 1 & !is.na(End_date) ~
        "Exact event",
      
      Healed == 1 &
        !is.na(Interval_censoring_date1) &
        !is.na(Interval_censoring_date2) ~
        "Interval censored",
      
      Healed == 0 & !is.na(Last_seen_date) ~
        "Right censored",
      
      TRUE ~ "Missing information"
    ),
    
    Lower = case_when(
      censoring_type == "Exact event" ~
        as.numeric(End_date - Start_date),
      
      censoring_type == "Interval censored" ~
        as.numeric(Interval_censoring_date1 - Start_date),
      
      censoring_type == "Right censored" ~
        as.numeric(Last_seen_date - Start_date),
      
      TRUE ~ NA_real_
    ),
    
    Upper = case_when(
      censoring_type == "Exact event" ~
        as.numeric(End_date - Start_date),
      
      censoring_type == "Interval censored" ~
        as.numeric(Interval_censoring_date2 - Start_date),
      
      censoring_type == "Right censored" ~ Inf,
      
      TRUE ~ NA_real_
    )
  )


table(Final_data$censoring_type)
## 
##       Exact event Interval censored    Right censored 
##               172                 5                11

We have only 5 cases we are not sure about their time of healing (interval censoring ) and 11 who are lost to follow-up

Turnbull is a nonparametric estimator that extends Kaplan–Meier to accommodate interval-censored data, where the event is known to have occurred between two time points but its exact timing is unknown

errors2
##               id ndays_to_healing_first Start_date   End_date
## 1 AM/HNYA/06/339                    411 2025-05-29 2026-07-14
## 2      BT/KV/013                    376 2025-06-23 2026-07-04
## 3      JK/01/013                    104 2025-01-11 2025-04-25
## 4      KM/MT/023                     73 2025-05-17 2025-07-29
## 5    MMG/CUB/005                     89 2024-11-17 2025-02-14
## 6    MMG/CUB/009                    380 2025-03-05 2026-03-20
## 7    NK/HNYA/002                     98 2025-04-25 2025-08-01
## 8    NK/HNYA/003                     71 2025-04-07 2025-06-17

Before doing the survival analysis we have to remove the 8 cases shown above (women with very high number of days to healing)

  Final_data2 <- Final_data %>%
  anti_join(errors2 %>% select(id), by = "id")

3 Survival curve by using Turbull estimator

set.seed(123)

fit_turnbull <- interval::icfit(
  L = Final_data2$Lower,
  R = Final_data2$Upper,
  conf.int = TRUE,
  control = interval::icfitControl(
    conf.level = 0.95,
    B = 500
  )
)


# 1b. Fit survfit object only to extract the median
fit_median <- survfit(
  Surv(Lower, Upper, type = "interval2") ~ 1,
  data = Final_data2
)

# Extract median time
median_obj <- quantile(fit_median, probs = 0.5, conf.int = TRUE)
median_time <- as.numeric(median_obj$quantile)

# Optional: median CI
median_lower <- as.numeric(median_obj$lower)
median_upper <- as.numeric(median_obj$upper)


# 2. Identify exact events and right censoring
exact_times <- Final_data2 %>%
  filter(Lower == Upper) %>%
  pull(Lower)

censored_times <- Final_data2 %>%
  filter(is.infinite(Upper)) %>%
  pull(Lower)

# 3. Turnbull curve and confidence intervals
plot(
  fit_turnbull,
  dtype = "survival",
  conf.int = TRUE,
  shade = TRUE,
  XLAB = "Days since lesion onset",
  YLAB = "Probability of remaining unhealed",
  main = "Time to complete disappearance of mpox lesions",
  estpar = list(
    col = "#176B8A",
    lwd = 3,
    lty = 1
  ),
  cipar = list(
    col = "#A5B7C5",
    lwd = 1.5,
    lty = 2
  )
)

# 4. Add marks for observed event and censoring times
rug(
  exact_times,
  side = 1,
  col = adjustcolor("#176B8A", alpha.f = 0.3),
  ticksize = 0.025
)

rug(
  censored_times,
  side = 3,
  col = "#B85C38",
  ticksize = 0.025,
  lwd = 1.5
)

# 5. Add median line and label
if (!is.na(median_time) && is.finite(median_time)) {
  
  # horizontal line at S(t)=0.5
  segments(
    x0 = 0, y0 = 0.5,
    x1 = median_time, y1 = 0.5,
    col = "#D55E00",
    lty = 2,
    lwd = 2
  )
  
  # vertical line at median time
  segments(
    x0 = median_time, y0 = 0,
    x1 = median_time, y1 = 0.5,
    col = "#D55E00",
    lty = 2,
    lwd = 2
  )
  
  # point at the median
  points(
    median_time, 0.5,
    pch = 19,
    col = "#D55E00",
    cex = 1.3
  )
  
  # text label
  text(
    x = median_time,
    y = 0.56,
    labels = paste0("Median = ", round(median_time, 1), " days"),
    col = "#D55E00",
    font = 2,
    pos = 4,
    cex = 0.95
  )
}

# 6. Add legend
legend(
  "topright",
  legend = c(
    "Turnbull estimate",
    "95% confidence limits",
    "Right-censored observations",
    paste0("Median = ", round(median_time, 1), " days")
  ),
  col = c("#176B8A", "#A5B7C5", "#B85C38", "#D55E00"),
  lty = c(1, 2, NA, 2),
  lwd = c(3, 1.5, NA, 2),
  pch = c(NA, NA, 124, 19),
  bty = "n",
  cex = 0.95
)

The figure illustrates the probability of remaining unhealed over time, estimated using the Turnbull estimator with 95% confidence intervals, accounting for exact, interval-censored, and right-censored observations. The survival probability decreased progressively throughout follow-up, reflecting the increasing occurrence of complete lesion healing. The median time to complete lesion healing was 27 days, corresponding to the time at which the estimated probability of remaining unhealed reached 50%. Approximately 50% of patients remained unhealed at Day 27, decreasing to around 20% by Day 45 and approximately 10% by Day 55. By approximately Day 65–70, the estimated probability of remaining unhealed approached zero, indicating that nearly all patients were estimated to have achieved complete lesion healing. The dashed lines represent the pointwise 95% confidence limits, while the orange marks at the top indicate right-censored observations.

4 Cumulative risk curve by using Turnbull estimator

# Turnbull cumulative healing curve with 95% CI
plot(
  fit_turnbull,
  dtype = "cdf",
  conf.int = TRUE,
  shade = TRUE,
  XLAB = "Days since lesion onset",
  YLAB = "Cumulative probability of complete lesion healing",
  main = "Time to complete disappearance of mpox lesions",
  estpar = list(
    col = "#176B8A",
    lwd = 3,
    lty = 1
  ),
  cipar = list(
    col = "#A5B7C5",
    lwd = 1.5,
    lty = 2
  )
)

# Exact healing times
rug(
  exact_times,
  side = 1,
  col = adjustcolor("#176B8A", alpha.f = 0.3),
  ticksize = 0.025
)

# Right-censored observations
rug(
  censored_times,
  side = 3,
  col = "#B85C38",
  ticksize = 0.025,
  lwd = 1.5
)

# Median healing time
if (length(median_time) == 1 &&
    !is.na(median_time) &&
    is.finite(median_time)) {

  segments(
    0, 0.5, median_time, 0.5,
    col = "#D55E00",
    lty = 2,
    lwd = 2
  )

  segments(
    median_time, 0, median_time, 0.5,
    col = "#D55E00",
    lty = 2,
    lwd = 2
  )

  points(
    median_time, 0.5,
    pch = 19,
    col = "#D55E00",
    cex = 1.3
  )

  text(
    median_time, 0.56,
    labels = paste0(
      "Median = ", round(median_time, 1), " days"
    ),
    col = "#D55E00",
    font = 2,
    pos = 4,
    cex = 0.95
  )
}

# Legend
has_median <- length(median_time) == 1 &&
  !is.na(median_time) && is.finite(median_time)

legend(
  "bottomright",
  legend = c(
    "Cumulative healing probability",
    "95% confidence limits",
    "Right-censored observations",
    if (has_median)
      paste0("Median = ", round(median_time, 1), " days")
  ),
  col = c(
    "#176B8A",
    "#A5B7C5",
    "#B85C38",
    if (has_median) "#D55E00"
  ),
  lty = c(1, 2, NA, if (has_median) 2),
  lwd = c(3, 1.5, NA, if (has_median) 2),
  pch = c(NA, NA, 124, if (has_median) 19),
  bty = "n",
  cex = 0.95
)

The figure presents the cumulative probability of complete disappearance of mpox lesions over time, estimated using the Turnbull nonparametric maximum likelihood estimator with 95% confidence intervals. The cumulative probability of healing increased progressively with time, indicating that an increasing proportion of patients achieved complete lesion healing during follow-up. The median time to complete lesion healing was 27 days, meaning that approximately 50% of patients were estimated to have completely healed by Day 27 after lesion onset. The cumulative probability of healing reached approximately 80% by Day 45, 90% by Day 55, and approached 100% by approximately Day 65–70. The dashed lines represent the pointwise 95% confidence limits around the estimated cumulative healing probability, while the orange marks at the top indicate right-censored observations. Overall, the figure suggests that most patients achieved complete lesion healing within approximately two months of lesion onset.

Important note: Both figures represent the same underlying Turnbull analysis and median healing time, but from complementary perspectives. Figure 1 shows the probability of remaining unhealed whereas Figure 2 shows cumulative probability of healing. I recommend using Figure 2 (cumulative healing probability) as the main figure because it directly reflects your outcome of interest: complete disappearance of mpox lesions.

5 If you need to compare between the levels of any risk factor (I selected vaccination status)

# Prepare data
data_vacc <- Final_data2 %>%
  filter(
    !is.na(vacc_mpox),
    !is.na(Lower),
    !is.na(Upper)
  ) %>%
  mutate(
    vacc_mpox = factor(vacc_mpox)
  )



# Fit stratified Turnbull estimator
fit_vacc <- interval::icfit(
  Surv(Lower, Upper, type = "interval2") ~ vacc_mpox,
  data = data_vacc,
  conf.int = TRUE,
  control = interval::icfitControl(
    conf.level = 0.95,
    B = 500,
    seed = 123
  )
)
cols <- c("#176B8A", "#D55E00")

plot(
  fit_vacc,
  dtype = "cdf",
  conf.int = TRUE,
  shade = FALSE,

  XLAB = "Days since lesion onset",
  YLAB = "Cumulative probability of complete lesion healing",
  main = "Time to complete lesion healing by vaccination status",

  LEGEND = FALSE,

  estpar = list(
    col = cols,
    lwd = 3,
    lty = 1
  ),

  cipar = list(
    col = c("#9AB7C5", "#E8AD8C"),
    lwd = 1.2,
    lty = 2
  )
)

legend(
  "bottomright",
  legend = levels(data_vacc$vacc_mpox),
  col = cols,
  lty = 1,
  lwd = 3,
  bty = "n",
  title = "Vaccination status"
)

6 Logrank test to compare between the 2 groups

test_vacc <- interval::ictest(
  Surv(Lower, Upper, type = "interval2") ~ vacc_mpox,
  data = data_vacc,
  scores = "logrank1"
)

print(test_vacc)
## 
##  Asymptotic Logrank two-sample test (permutation form), Sun's scores
## 
## data:  Surv(Lower, Upper, type = "interval2") by vacc_mpox 
## Z = 0.2769, p-value = 0.7818
## alternative hypothesis: survival distributions not equal
## 
##                 n Score Statistic*
## vacc_mpox=Non 169        0.8243716
## vacc_mpox=Oui  11       -0.8243716
## * like Obs-Exp, positive implies earlier failures than expected

According to this test, no significant difference in the time to healing between vaccinated and not vaccinated people (p-values=0.7818)

7 Cox proportional hazard model with 2 covariates

set.seed(123)
fit_cox <- ic_sp(
  Surv(Lower, Upper, type = "interval2") ~
    age+ vacc_mpox,
  data = Final_data2,
  model = "ph",
  bs_samples = 200
)

summary(fit_cox)
## 
## Model:  Cox PH
## Dependency structure assumed: Independence
## Baseline:  semi-parametric 
## Call: ic_sp(formula = Surv(Lower, Upper, type = "interval2") ~ age + 
##     vacc_mpox, data = Final_data2, model = "ph", bs_samples = 200)
## 
##              Estimate Exp(Est) Std.Error z-value       p
## age          -0.03204   0.9685   0.01317 -2.4330 0.01498
## vacc_mpoxOui -0.06642   0.9357   0.65790 -0.1009 0.91960
## 
## final llk =  -629.8693 
## Iterations =  50 
## Bootstrap Samples =  200
# Hazard ratios
exp(coef(fit_cox))
##          age vacc_mpoxOui 
##    0.9684692    0.9357407

Age (HR = 0.9685, p = 0.024) For each additional year of age, the instantaneous risk of complete lesion healing decreases by approximately 3.15%, after adjusting for mpox vaccination status. This suggests that older participants tend to heal more slowly than younger participants. The association is statistically significant at the 5% level.

Mpox vaccination status (HR = 0.9357, p = 0.926) However, the association is not statistically significant (p = 0.926). Therefore, these results do not provide evidence that vaccination status is associated with faster or slower lesion healing. Importantly, this does not mean vaccination has no effect. It means this model did not detect a statistically significant association.

8 How to check the proportionality assumption???

You should do that for every risk factor put in the model separately, we need to see that the lines are parrallel to each other

diag_covar(
  fit_cox,
  varName = "vacc_mpox",
  factorSplit = TRUE,
  main = "Proportional hazards assumption: Drinking",
  lgdLocation = "topright"
)

Here in this figure we can clearly see that the proportionality of hazard doesn’t hold across the vaccination status.

For the numerical variables as age, you should categorize it to 3 categories and see whether the lines will be parrallel or they will cross each other.

# Calculate the median age
median_age <- median(
  Final_data2$age,
  na.rm = TRUE
)

print(median_age)
## [1] 24
# Examine vaccination distribution
table(
  cut(Final_data2$age,
      breaks = c(-Inf, median_age, Inf)),
  Final_data2$vacc_mpox
)
##            
##             Non Oui
##   (-Inf,24]  87   4
##   (24, Inf]  82   7
# Attempt adjusted PH diagnostic
icenReg::diag_covar(
  fit_cox,
  varName = "age",
  factorSplit = FALSE,
  numericCuts = median_age,
  main = "PH assumption: Age",
  lgdLocation = "topright",
  col = c("blue", "red")
)

The proportionality of hazards works here, but we cannot use the above model, we have to switch to the AFT

9 AFT Model

fit_aft <- ic_par(
  Surv(Lower, Upper, type = "interval2") ~
    age+vacc_mpox,
  data = Final_data2,
  model = "aft",
  dist = "weibull"
)

summary(fit_aft)
## 
## Model:  Accelerated Failure Time
## Dependency structure assumed: Independence
## Baseline:  weibull 
## Call: ic_par(formula = Surv(Lower, Upper, type = "interval2") ~ age + 
##     vacc_mpox, data = Final_data2, model = "aft", dist = "weibull")
## 
##              Estimate Exp(Est) Std.Error z-value        p
## log_shape     0.64070    1.898   0.06074  10.550 0.000000
## log_scale     3.56000   35.180   0.04269  83.400 0.000000
## age           0.01514    1.015   0.00736   2.057 0.039690
## vacc_mpoxOui  0.54720    1.728   0.18110   3.021 0.002519
## 
## final llk =  -692.5018 
## Iterations =  7
# Time ratios for covariates
exp(coef(fit_aft))
##    log_shape    log_scale          age vacc_mpoxOui 
##     1.897843    35.175103     1.015255     1.728428

Age (Time Ratio = 1.015, p = 0.0397) For every additional year of age, the time required for complete lesion healing increases by approximately 1.5%, after adjusting for vaccination status. This suggests that older participants take longer to heal than younger participants. The association is statistically significant, although the estimated effect per year is relatively small. Mpox vaccination (Time Ratio = 1.728, p = 0.0025) Vaccinated participants have an estimated 72.8% longer time to complete lesion healing compared with unvaccinated participants of the same age. (I know that this is counter intuitive!!)

10 Sensitivity analysis

10.1 Removing the 5 cases with interval censoring

Sensitivity_data <- Final_data2 %>%
  filter(
    censoring_type %in% c(
      "Exact event",
      "Right censored"
    )
  ) %>%
  mutate(
    # Time until event or right censoring
    time = Lower,
    
    # Event indicator
    # 1 = Complete healing
    # 0 = Right censored
    event = ifelse(
      censoring_type == "Exact event", 1, 0
    )
  ) %>%
  filter(
    !is.na(time),
    time >= 0
  )

# Verify the sample


table(Sensitivity_data$censoring_type)
## 
##    Exact event Right censored 
##            164             11

10.2 Kaplan-Meier curve

fit <- survfit(Surv(time, event) ~ vacc_mpox, data = Sensitivity_data)


ggsurvplot(
  fit, data = Sensitivity_data,
  pval = TRUE,                 # overall log-rank p-value on the plot
  conf.int = TRUE,
  risk.table = TRUE,
  xlab = "Time (days)",
  ylab = "probability of not healing",
  legend.title = "Vaccination status"
)

10.3 Cumulative risk curve

ggsurvplot(
  fit,
  data = Sensitivity_data,
  fun = "event",
  pval = TRUE,
  conf.int = TRUE,
  risk.table = TRUE,
  surv.scale = "percent",

  # Add median healing time
  surv.median.line = "hv",

  palette = c("blue", "red"),
  xlab = "Time (days)",
  ylab = "Cumulative probability of complete healing (%)",
  legend.title = "Vaccination status",
  ggtheme = theme_minimal()
)

# Extract median healing time for each group
surv_median(fit)
##          strata median lower upper
## 1 vacc_mpox=Non   29.0    24    35
## 2 vacc_mpox=Oui   17.5    14    NA

11 Cox proportional hazard model

m1 <- coxph(Surv(time, event) ~ age  + vacc_mpox ,data = Sensitivity_data)

summary(m1)
## Call:
## coxph(formula = Surv(time, event) ~ age + vacc_mpox, data = Sensitivity_data)
## 
##   n= 175, number of events= 164 
## 
##                  coef exp(coef) se(coef)      z Pr(>|z|)  
## age          -0.03473   0.96587  0.01432 -2.425   0.0153 *
## vacc_mpoxOui -0.01524   0.98487  0.35716 -0.043   0.9660  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
##              exp(coef) exp(-coef) lower .95 upper .95
## age             0.9659      1.035    0.9391    0.9934
## vacc_mpoxOui    0.9849      1.015    0.4891    1.9833
## 
## Concordance= 0.571  (se = 0.028 )
## Likelihood ratio test= 6.28  on 2 df,   p=0.04
## Wald test            = 6  on 2 df,   p=0.05
## Score (logrank) test = 6.03  on 2 df,   p=0.05

Age (HR = 0.965, p = 0.015) For each additional year of age, the instantaneous risk of complete lesion healing decreases by approximately 3.15%, after adjusting for mpox vaccination status. This suggests that older participants tend to heal more slowly than younger participants. The association is statistically significant at the 5% level.

Mpox vaccination status (HR = 0.984, p = 0.966) However, the association is not statistically significant (p = 0.926). Therefore, these results do not provide evidence that vaccination status is associated with faster or slower lesion healing.

I don’t see big difference between this model and the one with the interval censoring

11.1 Check the proportionality assumption

cox.zph(m1)
##           chisq df      p
## age        0.55  1 0.4582
## vacc_mpox  8.49  1 0.0036
## GLOBAL     9.58  2 0.0083

It doesn’t hold for the vacc_mpox as we noticed in the model with interval censoring.

We could not ue this model and must shift to th AFT

11.2 AFT model with Weibul distribution

 aft1 <- survreg(Surv(time, event) ~ age+vacc_mpox, data=Sensitivity_data,
                     dist = "weibul") 

summary(aft1)
## 
## Call:
## survreg(formula = Surv(time, event) ~ age + vacc_mpox, data = Sensitivity_data, 
##     dist = "weibul")
##                 Value Std. Error      z      p
## (Intercept)   3.11191    0.19051  16.33 <2e-16
## age           0.01669    0.00751   2.22 0.0263
## vacc_mpoxOui  0.55267    0.18711   2.95 0.0031
## Log(scale)   -0.63739    0.06108 -10.44 <2e-16
## 
## Scale= 0.529 
## 
## Weibull distribution
## Loglik(model)= -690   Loglik(intercept only)= -698.5
##  Chisq= 16.98 on 2 degrees of freedom, p= 0.00021 
## Number of Newton-Raphson Iterations: 7 
## n= 175

Age (TR = 1.017, p = 0.026) For every additional year of age, the estimated time required for complete lesion healing increases by approximately 1.7%, after adjusting for vaccination status. This suggests that older participants tend to take longer to heal than younger participants.

Mpox vaccination (TR = 1.738, p = 0.003) Vaccinated participants have an estimated 73.8% longer time to complete lesion healing compared with unvaccinated participants of the same age.

The same interpretation as in the AFT model done above with the 5 cases with interval censoring

You can do the same and add the risk factors that you need and compare between the 2 models (the one with interval censoring and the one after removing the interval censoring)