Pre-apneic capnography waveform abnormalities during procedural sedation and analgesia

For reproducibility, this notebook contains the data and code used to perform statistical analysis for the journal article 'Pre-apneic capnography waveform abnormalities during procedural sedation and analgesia'.

Environment set-up

install.packages(c("coxme", "qwraps2", "grid", "gridExtra"))

78.0s

library(dplyr)
library(coxme)
library(qwraps2)
library(ggplot2)
library(grid)
library(gridExtra)

2.1s

Code to produce table of participant characteristics (Table 1 from the article)

1.3s

# Code to produce table of participant characteristics

survival.data <- read.csv(  data.csv
)

Sample <- survival.data %>% 
  distinct(ID, Sex, Age, BMI, OSA) %>% 
  select(Sex, Age, BMI, OSA) %>%
  mutate(Sex = as.factor(Sex), OSA = as.factor(OSA))

myCapnoData <- readr::read_csv(file=  data.csv
)

tribtable <- tribble(~Characteristic, ~"Mean (SD) or Frequency (%)",
                     "Age", mean_sd(myCapnoData$Age, denote_sd = "paren"),
                     "Female", n_perc0(myCapnoData$Sex==1),
                      "Body mass index", mean_sd(myCapnoData$BMI,  na_rm = TRUE, show_n = "never", denote_sd = "paren"),
                     "Diagnosis of obstructive sleep apnea", n_perc0(myCapnoData$OSA==1),
                     "Charlson comorbidity index", mean_sd(myCapnoData$CCI, na_rm = TRUE, show_n = "never", denote_sd = "paren"),
                     "Permanent pacemaker implant or generator change", n_perc0(myCapnoData$procedure=="PPM"),
                     "Implantable cardioverter defibrillator implant or generator change", n_perc0(myCapnoData$procedure=="ICD"),
                     "Cardiac resynchronisation therapy", n_perc0(myCapnoData$procedure=="CRT"),
                     "Atrial flutter ablation", n_perc0(myCapnoData$procedure=="flutter"),
                     "Other arrhythmia ablation", n_perc0(myCapnoData$procedure=="RFA"),
                     "Diagnostic electrophysiology study", n_perc0(myCapnoData$procedure=="EPS"),
                     "Loop recorder implant", n_perc0(myCapnoData$procedure=="Loop recorder implant"),
                     "ASA Class I", n_perc0(myCapnoData$ASA==1),
                     "ASA Class II", n_perc0(myCapnoData$ASA==2),
                     "ASA Class III", n_perc0(myCapnoData$ASA==3),
                     "ASA Class IV", n_perc0(myCapnoData$ASA==4),
                     "Midazolam total dose (mg)", mean_sd(myCapnoData$Midazolam,  na_rm = TRUE, show_n = "never", denote_sd = "paren"),
                     "Fetanyl total dose (mg)", mean_sd(myCapnoData$Fentanyl,  na_rm = TRUE, show_n = "never", denote_sd = "paren"))

tribtable

1.3s

0 items

Code to produce figure displaying the distribution of apneic episodes (Figure 1)

0.8s

# Format data to obtain survival times between apneic episodes
survival.times <- survival.data %>% 
  group_by(ID) %>%
  mutate(dur = lead(stop-start)) %>%
  filter(event==1 ) %>%
  mutate(dur = stop - lag(stop + dur)) %>%
  mutate(dur = ifelse(is.na(dur), stop, dur)) %>%
  mutate(episode = 1:n()) %>%
  select(episode, dur) %>% 
  ungroup()
  

survival.times  <- survival.times[order(survival.times$episode, -survival.times$dur),]

survival.times$orders <- 1:nrow(survival.times)

median.time <- survival.times %>%
  group_by(episode) %>%
  dplyr::summarize(med=median(dur))

median.epi <- survival.times %>% 
  group_by(episode) %>%
  dplyr::summarize(med=median(orders))

# Plotting
ggplot(data=survival.times, aes(x=dur, y=orders)) +
geom_rect(data=survival.times, mapping = aes(xmin=0,xmax=dur,ymin=orders,
                        ymax=orders+1, fill=factor(episode)))  +
geom_segment(x=median.time$med[1],xend=median.time$med[1],
          y=survival.times %>% filter(episode == 1) %>% pull(orders) %>% min() ,
          yend=survival.times %>% filter(episode == 1) %>% pull(orders) %>% max(),
          colour="black",linetype="dashed") +
  geom_segment(x=median.time$med[2],xend=median.time$med[2],
          y=survival.times %>% filter(episode == 2) %>% pull(orders) %>% min() ,
          yend=survival.times %>% filter(episode == 2) %>% pull(orders) %>% max(),
          colour="black",linetype="dashed") +
  geom_segment(x=median.time$med[3],xend=median.time$med[3],
          y=survival.times %>% filter(episode == 3) %>% pull(orders) %>% min() ,
          yend=survival.times %>% filter(episode == 3) %>% pull(orders) %>% max(),
          colour="black",linetype="dashed") +
  geom_segment(x=median.time$med[4],xend=median.time$med[4],
          y=survival.times %>% filter(episode == 4) %>% pull(orders) %>% min() ,
          yend=survival.times %>% filter(episode == 4) %>% pull(orders) %>% max(),
          colour="black",linetype="dashed") +
  geom_segment(x=median.time$med[5],xend=median.time$med[5],
          y=survival.times %>% filter(episode == 5) %>% pull(orders) %>% min() ,
          yend=survival.times %>% filter(episode == 5) %>% pull(orders) %>% max(),
          colour="black",linetype="dashed") +
scale_y_continuous(breaks=median.epi$med[c(1,2,3,4,5,10,15,20,25)], labels = c(1,2,3,4,5,10,15,20,25)) +
  coord_cartesian(xlim=c(0,1000)) +
xlab("Time (s)") + ylab("Episode") +
ggtitle("Survival time by episode") +
 theme_bw() +
    theme(panel.grid.major.y =element_blank(),
          panel.grid.minor.y =element_blank(),
          panel.background=element_blank(),
          legend.position="none",
          axis.line.x = element_line(color = "black"),
          axis.line.y = element_line(color = "black"))

0.8s

Code to produce results for the mixed effects Cox model (Table 2)

2.0s

# Load dataset for first model
survival.data <- read.csv(survival_recurrent_20.csv
)

survival.data <- survival.data %>% filter(state!=4)

survival.data$state <- survival.data$state %>% as.factor()
survival.data$Sex <- survival.data$Sex %>% as.factor()
survival.data$OSA <- survival.data$OSA %>% as.factor()

# Subtract 1 so that the baseline parameterization is at episodes = 0
survival.data <- survival.data %>% group_by(ID) %>% mutate(episode = episode - 1)  %>%
  ungroup()

# Load dataset for second model, where apneic episode is defined as at least 20 seconds or more of apneic breathing
survival.data.20 <- read.csv(file=  survival_recurrent_20.csv
)

survival.data.20 <- survival.data.20 %>% filter(!(state==4 & stop-start>=20))

survival.data.20$state <- survival.data.20$state %>% as.factor()
survival.data.20$Sex <- survival.data.20$Sex %>% as.factor()
survival.data.20$OSA <- survival.data.20$OSA %>% as.factor()

# Subtract 1 so that the baseline parameterization is at episodes = 0 and n.interventions=0
survival.data.20 <- survival.data.20 %>% group_by(ID) %>% mutate(episode = episode - 1) %>%
  ungroup()
# Fit models

model1 <- coxme(Surv(start, stop, event) ~ state + n.intervention + episode + (1|ID), data=survival.data)
model2 <- coxme(Surv(start, stop, event) ~ state + n.intervention + episode + (1|ID), data=survival.data.20)

# Extract information about models to present in a table

coef1 <- model1$coefficients
coef2 <- model2$coefficients

stderr1 <- vcov(model1) %>% diag() %>% sqrt()
stderr2 <- vcov(model2) %>% diag() %>% sqrt()

HR1 <- coef1  %>% exp()
HR2 <- coef2  %>% exp()

lower.HR1 <- (coef1 - 1.96*stderr1) %>% exp()
lower.HR2 <- (coef2 - 1.96*stderr2) %>% exp()

upper.HR1 <- (coef1 + 1.96*stderr1) %>% exp()
upper.HR2 <- (coef2 + 1.96*stderr2) %>% exp()
  
summary1 <- data.frame( predictor=c("hypopnea", "bradypnea", "sedation dose", "episode"), 
                        HR = exp(coef1),
                        `2.5 %` = lower.HR1,
                        `97.5 %` = upper.HR1)
summary2 <- data.frame( predictor=c("hypopnea", "bradypnea", "apnea", "sedation dose", "episode"), 
                        HR = exp(coef2),
                        `2.5 %` = lower.HR2,
                        `97.5 %` = upper.HR2)

rownames(summary1) <- NULL
rownames(summary2) <- NULL

model.table <- right_join(summary1, summary2, by="predictor" )
model.table$predictor <- c("Hypopnea", "Bradypnea", "Apnea", "Sedation dose", "Number of apneic episodes")
names(model.table) <- c(" ", "HR", "lower", "upper",  "HR", "lower", "upper")
model.table

2.0s

0 items

Code to produce figure demonstrating change in risk for apnea throughout procedures (Figure 2)

1.5s

# Define upper and lower confidence intervals for predictions from model:
predict.conf <- function(model, id = 27){
  patient <- survival.data %>% mutate(linear.predictor =  predict(model)) %>% filter(ID==id)
  
  se <- vcov(model) %>% diag() %>% sqrt()
  
  patient <- patient %>%
    mutate(upper = linear.predictor + 1.96*(se[1]*(state == "2") + se[2]*(state == "3") + se[3]*n.intervention + se[4]*episode)) %>%
    mutate(lower = linear.predictor - 1.96*(se[1]*(state == "2") + se[2]*(state == "3") + se[3]*n.intervention + se[4]*episode))
  
  return(patient)
}

patient.27 <- predict.conf(model1, id=27)
patient.27 <- patient.27 %>% mutate(lead_lin = lead(linear.predictor))

# Plot for linear predictor
model.pred <- ggplot(patient.27 , aes(Time, linear.predictor)) +
  geom_segment(aes(x = start, xend = stop, y = linear.predictor, yend = linear.predictor)) + 
  geom_segment(data = patient.27 %>% filter(event==0),  
                             aes(x = stop, xend = stop, y = linear.predictor, yend = lead_lin), linetype=2) +
  geom_hline(yintercept = 0, linetype=2) +
  geom_rect(aes(xmin=start, xmax=stop, ymin=lower, ymax=upper), alpha=0.15) +
  scale_x_continuous(expand = c(0, 0)) +
  ylab("log(HR)") + 
  ggtitle("Predicted log-Hazard Ratio for a patient ") +
  theme_bw() +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(margin=margin(t = 0, r = 9, b = 0, l = 0)))
  
patient.27$state <- patient.27$state %>% as.numeric()

col.palette <- c("1" = "#2b83ba","2" = "#fdae61", "3" = "#abdda4", "4" = "#d7191c")
  
inter.plot <- ggplot(patient.27, aes(Time, n.intervention)) +
  geom_step(aes(x = Time, y= n.intervention)) +
  geom_segment(aes(x=start,xend=stop, y= n.intervention, yend= n.intervention)) +
  scale_y_continuous(limits=c(1,3), breaks=c(1,2,3)) +
  scale_x_continuous(expand = c(0,0)) +
  theme_bw() +
  theme(axis.title.x = element_blank(),
        axis.title.y = element_text(margin=margin(t = 0, r = 15, b = 0, l = 0))) +
  ylab("Sedation\n Doses")


state.plot <- ggplot(patient.27, aes(Time, state)) +
  geom_rect(aes(xmin = start, xmax = stop, ymin = 0, ymax = 1, fill = factor(patient.27$state))) +
  scale_x_continuous(expand = c(0, 0)) +
  scale_y_continuous(expand = c(0,0), limits=c(0,1)) +
  theme_classic() +
  theme(axis.ticks.y = element_blank(),
        axis.text.y = element_blank(),  
        axis.title.y = element_text(margin=margin(t = 0, r = 15, b = 0, l = 0)),
        panel.background = element_rect(fill = col.palette[4]),
        legend.position = "bottom") +
  ylab("State") + 
  xlab("Time (seconds)") +
  scale_fill_manual(values = col.palette,
                    labels = c("Normal","Hypopnea", "Bradypnea", "Apnea"), 
                    limits = c("1", "2","3","4"),
                    name = "State legend:")

grid.newpage()
grid.arrange(ggplotGrob(model.pred), rbind(ggplotGrob(inter.plot), ggplotGrob(state.plot) , size = "last"))

1.5s

References

data.csv

survival_recurrent_25.csv

survival_recurrent_15.csv

survival_recurrent_10.csv

survival_recurrent_20.csv

survival_data.csv

survival_recurrent_0.csv

survival_recurrent.csv