APRV Use in PICU: Descriptive Analysis

Author

James Weitz

Published

December 16, 2025

1. Introduction

APRV (Airway Pressure Release Ventilation) is used in our PICU as an alternative mode of mechanical ventilation in selected patients with respiratory failure. This report describes how APRV has been used in our unit since 2020 and explores patient characteristics, ventilatory and gas-exchange parameters over time, and key clinical outcomes.

We use routinely collected data from the PICU information system (2020–present), including all mechanically ventilated admissions, both with and without APRV.

1.1 Objectives

The analysis aims to:

Describe demographic and clinical characteristics of patients receiving APRV.
Evaluate changes in oxygenation and respiratory parameters around APRV use.
Assess the association between APRV and neuromuscular blockade (paralysis), as a proxy for sedation intensity.
Report key clinical outcomes, including:
- Total duration of APRV
- Total duration of mechanical ventilation
- PICU length of stay
- 28-day mortality
Compare APRV vs non-APRV ventilated patients on high-level severity and outcome measures (e.g. PIM3 risk, mortality, LOS)

2. Methods

2.1 Data source

We used a time-series data-set exported from our unit electronic medical record (ICIP Carvue - Philips date) joined with data from a propriety database used to export data for the national PICU audit (PICANET, DB = mela - mediucs ….) . Data extraction was done using Microsoft SQL Server (T-SQL) and cleaned and analysed using R (ref and date))

Key variables included in the final dataset:

Ventilator and respiratory parameters: ventMode, FiO2, MAP, PEEP, PIP, PHigh, PLow, THigh, TLow, Pressure Support, PS above PEEP, PC above PEEP.
Physiology / blood gases: SpO2, Base Excess, HCO3, Lactate, pCO2, pH, pO2.
Other clinical variables: nmb_on (neuromuscular blockade on/off), pim3_death_risk, diagnoses, weight, admission/discharge times, etc.
Derived flags: ever_aprv (indicator of APRV use), ventilated(indicator of any mechanical ventilation).

2.2 Study population

We will include:

All PICU admissions between 2020-01-01 and the latest available date.
For ventilated cohort analyses, we restrict to encounters with ventilated == TRUE.
The APRV group is defined as ventilated patients with ever_aprv == TRUE.
The comparison group is ventilated patients with ever_aprv == FALSE.

Code

# Define ventilated cohort and APRV vs non-APRV groups

# vector of ventilated modes
vent_modes <- c("APRV", "APV", "ASV", "CMV", "HFO", "ManualVent",
                "NAVA", "PC", "PRVC", "SIMV", "VC", "VS")

# 1) Keep only time-points where the patient is actually on a ventilated mode
vent_ts <- aprv |>
  mutate(
    is_vent = ventMode %in% vent_modes
  ) |>
  filter(is_vent)


# 2) Create an APRV vs non-APRV grouping variable
#    (adjust the `ever_aprv_flag` line if your variable is coded differently)
vent_ts <- vent_ts |>
  mutate(
    ever_aprv_flag = ever_aprv %in% c(TRUE, 1, "Yes", "Y"),
    aprv_group = if_else(
      ever_aprv_flag,
      "APRV",
      "No APRV"
    ),
    aprv_group = factor(aprv_group, levels = c("No APRV", "APRV"))
  )

# 3) Create an encounter-level dataset (one row per PICU encounter)
#    We just take the first row per encounterId and keep relevant columns.
vent_encounters <- vent_ts |>
  arrange(encounterId, hour) |>
  distinct(encounterId, .keep_all = TRUE) |>
  select(
    encounterId,
    aprv_group,
    ever_aprv_flag,
    ventilated,
    decAge,
    weight_on_adm,
    pim3_death_risk,
    los_days,
    admitUTC,
    dischargeUTC,
    admission_source_,
    discharge_status_,
    PrimDiagnosis,
    OtherDiagnosis,
    Operations,
    diagCat
  )

vent_encounters %>% count(aprv_group)

2.3 Variables and definitions

2.3.1 Exposure: APRV

APRV use: defined by ever_aprv == TRUE.
APRV period: for APRV patients, we will later define when APRV starts and stops, e.g. based on ventMode or other mode indicators over time

Code

# Identify APRV episodes and derive APRV duration per encounter

# ---- 1) Define which vent modes count as APRV ----
# Update this list as needed (run: vent_ts |> count(ventMode) to inspect modes)
aprv_modes <- c("APRV")   # add other labels if needed, e.g. "BiLevel", etc.

# ---- 2) Create APRV time-series object ----
# This works on vent_ts (which already contains only truly ventilated rows)
aprv_ts <- vent_ts |>
  arrange(encounterId, hour) |>
  mutate(
    # TRUE if current row is APRV mode
    is_aprv = ventMode %in% aprv_modes,

    # TRUE when APRV starts at this timepoint
    aprv_start_flag = is_aprv & !lag(is_aprv, default = FALSE),

    # Number APRV episodes: 1, 2, 3, …
    aprv_episode = if_else(
      is_aprv,
      cumsum(aprv_start_flag),
      NA_integer_
    )
  )

# ---- 3) Summarise each APRV episode (start, end, duration) ----
aprv_episodes <- aprv_ts |>
  filter(is_aprv) |>
  group_by(encounterId, aprv_episode) |>
  summarise(
    aprv_start_time      = min(hour, na.rm = TRUE),
    aprv_end_time        = max(hour, na.rm = TRUE),
    aprv_duration        = difftime(aprv_end_time, aprv_start_time, units = "hours"),
    aprv_duration_hours  = as.numeric(aprv_duration),
    .groups              = "drop"
  )

# ---- 4) Summarise APRV per encounter ----
aprv_by_encounter <- aprv_episodes |>
  group_by(encounterId) |>
  summarise(
    aprv_total_hours      = sum(aprv_duration_hours, na.rm = TRUE),
    aprv_episodes_n       = n(),
    aprv_first_start_time = min(aprv_start_time, na.rm = TRUE),
    n_aprv_rows = n(),
    .groups               = "drop"
  )

Code

# Define time windows relative to APRV start and aggregate key variables
# Produces:
#   - aprv_windows_ts: row-level data with time_from_aprv_start_hours + period
#   - aprv_windows_enc_ITT: encounter × period medians (all rows after APRV start)
#   - aprv_windows_enc_PP: encounter × period medians (only when actually on APRV)

# ---- 1) Set window definitions (hours) ----
pre_window_hours   <- 6   # look-back period before APRV start
early_window_hours <- 6   # first 0–6 hours after APRV start
late_window_hours  <- 24  # 6–24 hours after APRV start
# anything >= late_window_hours will be "post"

# ---- 2) Add time-from-APRV-start and period labels to the time series ----
aprv_windows_ts <- aprv_ts |>
  left_join(
    aprv_by_encounter |>
      select(encounterId, aprv_first_start_time),
    by = "encounterId"
  ) |>
  mutate(
    # Time difference in HOURS relative to first APRV start
    time_from_aprv_start_hours = as.numeric(
      difftime(hour, aprv_first_start_time, units = "hours")
    ),

    period = case_when(
      # Never had APRV → no window period
      is.na(aprv_first_start_time) ~ NA_character_,

      # PRE window: within pre_window_hours hours before APRV start
      time_from_aprv_start_hours >= -pre_window_hours &
        time_from_aprv_start_hours < 0 ~ "pre",

      # EARLY: 0 to early_window_hours
      time_from_aprv_start_hours >= 0 &
        time_from_aprv_start_hours < early_window_hours ~ "early",

      # LATE: early_window_hours to late_window_hours
      time_from_aprv_start_hours >= early_window_hours &
        time_from_aprv_start_hours < late_window_hours ~ "late",

      # POST: anything after late_window_hours
      time_from_aprv_start_hours >= late_window_hours ~ "post",

      TRUE ~ NA_character_
    ),
    period = factor(period, levels = c("pre", "early", "late", "post"))
  )

# ---- 3) ITT-style summary: all data after APRV start, regardless of mode ----
aprv_windows_enc_ITT <- aprv_windows_ts |>
  filter(
    ever_aprv_flag,
    !is.na(period)
  ) |>
  group_by(encounterId, period) |>
  summarise(
    n_points      = n(),
    FiO2_median   = median(FiO2, na.rm = TRUE),
    SpO2_median   = median(SpO2, na.rm = TRUE),
    MAP_median    = median(MAP, na.rm = TRUE),
    pCO2_median   = median(pCO2, na.rm = TRUE),
    pH_median     = median(pH, na.rm = TRUE),
    .groups = "drop"
  )

# ---- 4) Per-protocol summary: only timepoints actually on APRV ----
aprv_windows_enc_PP <- aprv_windows_ts |>
  filter(
    ever_aprv_flag,
    !is.na(period),
    is_aprv == TRUE          # this is the key difference
  ) |>
  group_by(encounterId, period) |>
  summarise(
    n_points      = n(),
    FiO2_median   = median(FiO2, na.rm = TRUE),
    SpO2_median   = median(SpO2, na.rm = TRUE),
    MAP_median    = median(MAP, na.rm = TRUE),
    pCO2_median   = median(pCO2, na.rm = TRUE),
    pH_median     = median(pH, na.rm = TRUE),
    .groups = "drop"
  )

2.3.2 Time windows for respiratory analyses

For patients with APRV, we plan to summarise respiratory parameters over:

A pre-APRV window (e.g. X hours before APRV starts)
The early APRV period (e.g. first Y hours on APRV)
The later APRV period (e.g. after Y hours on APRV)
(Optional) A post-APRV window if patients switch back to conventional modes
The exact window definitions will be finalised once we look at how densely the time series is sampled.

2.3.3 Outcomes

We will consider:

APRV duration (hours): time between APRV start and stop.
Total duration of mechanical ventilation (hours): from first to last ventilated time.
PICU length of stay (days): los_days.
28-day mortality: death before or on day 28 from unit admission.
(Potentially) Use of neuromuscular blockade during ventilation (nmb_on).

Code

# Derive encounter-level outcomes:
#   - Ventilation duration (hours)
#   - PICU LOS (already present as los_days)
#   - 28-day PICU mortality (approximate)

# ---- 1) Ventilation duration from time-series data ----
# We assume `vent_ts` contains one row per timepoint while ventilated,
# with a datetime column `hour`.

vent_durations <- vent_ts |>
  group_by(encounterId) |>
  summarise(
    vent_start_time = min(hour, na.rm = TRUE),
    vent_end_time   = max(hour, na.rm = TRUE),
    vent_duration   = difftime(vent_end_time, vent_start_time, units = "hours"),
    vent_duration_hours = as.numeric(vent_duration),
    .groups = "drop"
  )


# ---- 3) Attach outcomes to encounter-level data ----
vent_encounters <- vent_encounters |>
  left_join(vent_durations, by = "encounterId") |>
  mutate(
    # Ensure LOS is numeric (if it's not already)
    los_days = as.numeric(los_days),

    # Time from PICU admission to PICU discharge (days)
    days_from_adm_to_discharge = as.numeric(
      difftime(dischargeUTC, admitUTC, units = "days")
    ),

    # Did this patient die in PICU?
    died_in_picu = discharge_status_ == "Dead",

    # 28-day PICU mortality (approximation):
    #   death in PICU AND discharge (death) within 28 days of unit admission.
    death_28d_picu = died_in_picu & days_from_adm_to_discharge <= 28
  )

2.4 Statistical approach

This is primarily a descriptive analysis:

Continuous variables: summarised using mean ± SD or median (IQR) as appropriate.
Categorical variables: summarised using counts and percentages.
Respiratory parameters: described over time and across defined windows.
APRV vs non-APRV ventilated patients:
- Compare baseline severity (e.g. PIM3)
- Compare key outcomes (vent duration, LOS, mortality) descriptively
- We may include simple unadjusted comparisons (e.g. boxplots, tables).

#| label: analysis-plan-sketch
#| echo: false
# You might later store your choices (e.g. transform variables, etc.)
# but no need to do anything here yet.

3. Results

3.1 Cohort overview

Describe:

Total number of PICU admissions
Number and proportion ventilated
Number and proportion with APRV

Code

# --- Basic counts ---

total_admissions <- aprv |>
  distinct(encounterId) |>
  tally() |>
  pull(n)

ventilated_admissions <- vent_encounters |>
  tally() |>
  pull(n)

aprv_admissions <- vent_encounters |>
  filter(ever_aprv_flag) |>
  tally() |>
  pull(n)

non_aprv_admissions <- ventilated_admissions - aprv_admissions

# --- Create a tidy summary table ---

cohort_overview_tbl <- tibble(
  Metric = c(
    "Total admissions",
    "Total Ventilated admissions",
    "Total Venilated with APRV"
  ),
  Count = c(
    total_admissions,
    ventilated_admissions,
    aprv_admissions
  )
)

cohort_overview_tbl |>
  gt() |>
  tab_header(title = "Cohort overview") |>
  cols_label(
    Metric = "Metric",
    Count  = "N"
  ) |>
  fmt_number(columns = Count, decimals = 0)

Table 1

Cohort overview
Metric	N
Total admissions	3,880
Total Ventilated admissions	1,486
Total Venilated with APRV	105

We identified X PICU admissions between 1 January 2020 and [end date].

Of these, Y (…%) were mechanically ventilated, and Z (…%) of ventilated

patients received APRV at some point during their PICU stay.

3.2 Demographic and clinical characteristics

Summarise baseline demographics and severity for:

APRV patients
Non-APRV ventilated patients

Key variables:

Age (decAge), weight_on_adm
PIM3 predicted mortality (pim3_death_risk)
Admission source, diagnoses, etc.

Code

# For now we load it here; you can move this to the setup chunk later if you like.
library(gtsummary)


baseline_data <- vent_encounters |>
  mutate(
    diagCat = forcats::fct_infreq(diagCat)
  )

tbl_baseline <- baseline_data |>
  select(
    aprv_group,
    decAge,
    weight_on_adm,
    pim3_death_risk,
    admission_source_,
    diagCat
  ) |>
  tbl_summary(
    by = aprv_group,
    label = list(
      decAge           ~ "Age (years)",
      weight_on_adm    ~ "Weight on admission (kg)",
      pim3_death_risk  ~ "PIM3 predicted mortality (%)",
      admission_source_~ "Admission source",
      diagCat          ~ "Primary diagnostic category"
    ),
    statistic = list(
      all_continuous() ~ "{median} ({p25}, {p75})",
      all_categorical() ~ "{n} ({p}%)"
    ),
    missing = "ifany"
  ) |>
  bold_labels()

tbl_baseline

Table 2

Characteristic	No APRV N = 1,381¹	APRV N = 105¹
Age (years)	1.7 (0.3, 7.5)	1.8 (0.8, 6.3)
Weight on admission (kg)	12 (5, 25)	12 (8, 21)
PIM3 predicted mortality (%)	2 (1, 5)	3 (1, 11)
Admission source
Other hospital	732 (53%)	64 (61%)
Same hospital	649 (47%)	41 (39%)
Primary diagnostic category
Respiratory	598 (43%)	87 (83%)
Neurology	276 (20%)	3 (2.9%)
Infection/Immune/Inflammatory	96 (7.0%)	4 (3.8%)
Neurosurgery	74 (5.4%)	0 (0%)
ENT	66 (4.8%)	2 (1.9%)
General Surgery	63 (4.6%)	1 (1.0%)
Cardiac	41 (3.0%)	2 (1.9%)
Trauma/Poisoning	39 (2.8%)	0 (0%)
Out of hospital Cardiac Arrest	34 (2.5%)	2 (1.9%)
Haematology/Oncology	26 (1.9%)	2 (1.9%)
Spinal	18 (1.3%)	0 (0%)
Endocrine/Metabolic/Electrolyte abnormalities	15 (1.1%)	1 (1.0%)
Genetic condition	13 (0.9%)	0 (0%)
Multi-organ failure	8 (0.6%)	0 (0%)
Other/Not Specified	5 (0.4%)	1 (1.0%)
Gastro	5 (0.4%)	0 (0%)
Unknown	4	0
¹ Median (Q1, Q3); n (%)

Text here briefly describing any notable differences.

3.3 Ventilation and respiratory parameters

3.3.1 Ventilator settings over time on APRV

Here you will show, for APRV patients:

Trends in PHigh, PLow, THigh, TLow, MAP, PEEP over time.
Possibly stratified by pre/early/late APRV windows.

#| label: fig-aprv-vent-params
#| echo: false
# Target: 1–2 plots:
#   - e.g. line plots or summary boxplots of PHigh, PLow, FiO2, etc.
#   - time could be "hours since APRV start"

3.3.2 Oxygenation and gas exchange

Focus on:

FiO2 and SpO2
pO2 / paO2 if available
pCO2 and pH
Lactate (if clinically relevant)

Code

# Ventilator settings (on APRV) across pre / early / late / post windows

# 1) Summarise APRV settings per encounter × period
aprv_vent_windows <- aprv_windows_ts |>
  filter(
    ever_aprv_flag,
    !is.na(period),
    is_aprv == TRUE              # only while actually on APRV
  ) |>
  group_by(encounterId, period) |>
  summarise(
    PHigh_median  = median(PHigh, na.rm = TRUE),
    PLow_median   = median(PLow,  na.rm = TRUE),
    THigh_median  = median(THigh, na.rm = TRUE),
    TLow_median   = median(TLow,  na.rm = TRUE),
    MAP_median    = median(MAP,   na.rm = TRUE),
    # PEEP_median   = median(PEEP,  na.rm = TRUE),
    .groups = "drop"
  )

# 2) Pivot longer for plotting
aprv_vent_long <- aprv_vent_windows |>
  tidyr::pivot_longer(
    cols = ends_with("_median"),
    names_to = "parameter",
    values_to = "value"
  ) |>
  mutate(
    parameter = dplyr::recode(
      parameter,
      PHigh_median = "PHigh",
      PLow_median  = "PLow",
      THigh_median = "THigh",
      TLow_median  = "TLow",
      MAP_median   = "MAP",
      #PEEP_median  = "PEEP"
    )
  )

# 3) Boxplots of settings by period, facetted by parameter
ggplot(aprv_vent_long, aes(x = period, y = value)) +
  geom_boxplot() +
  facet_wrap(~ parameter, scales = "free_y") +
  labs(
    title = "Ventilator settings on APRV by time window",
    x = "Window relative to APRV start",
    y = "Median value per encounter"
  )

Code

# Physiological changes across APRV windows (ITT-style: all data after APRV start)

# 1) Prepare ITT-style physiology data with SF ratio and OSI
aprv_phys_ITT <- aprv_windows_enc_ITT |>
  mutate(
    # Your FiO2 appears to be % (e.g. 40–80), so convert to fraction
    FiO2_frac = FiO2_median / 100,

    # SF ratio (SpO2 / FiO2 fraction)
    SF_ratio = SpO2_median / FiO2_frac,

    # OSI = (FiO2 * MAP * 100) / SpO2, where FiO2 is fraction and SpO2 is %
    OSI = (FiO2_frac * MAP_median * 100) / SpO2_median
  )

# ---- Oxygenation panel: FiO2, SpO2, SF ratio, OSI ----
ox_long <- aprv_phys_ITT |>
  select(
    encounterId,
    period,
    FiO2_median,
    SpO2_median,
    SF_ratio,
    OSI
  ) |>
  tidyr::pivot_longer(
    cols = c(FiO2_median, SpO2_median, SF_ratio, OSI),
    names_to = "measure",
    values_to = "value"
  ) |>
  mutate(
    measure = dplyr::recode(
      measure,
      FiO2_median = "FiO2 (%)",
      SpO2_median = "SpO2 (%)",
      SF_ratio    = "SF Ratio (SpO2/FiO2)",
      OSI         = "OSI"
    ),
    measure = factor(
      measure,
      levels = c("FiO2 (%)", "SpO2 (%)", "SF Ratio (SpO2/FiO2)", "OSI")
    )
  )

p_oxygenation <- ggplot(ox_long, aes(x = period, y = value)) +
  geom_boxplot() +
  facet_wrap(~ measure, scales = "free_y") +
  labs(
    title = "Oxygenation metrics by window relative to APRV start (ITT)",
    x = "Window (pre / early / late / post)",
    y = "Median value per encounter"
  )

# ---- Ventilation panel: pCO2 and pH ----
vent_long <- aprv_phys_ITT |>
  select(
    encounterId,
    period,
    pCO2_median,
    pH_median
  ) |>
  tidyr::pivot_longer(
    cols = c(pCO2_median, pH_median),
    names_to = "measure",
    values_to = "value"
  ) |>
  mutate(
    measure = dplyr::recode(
      measure,
      pCO2_median = "pCO2",
      pH_median   = "pH"
    )
  )

p_ventilation <- ggplot(vent_long, aes(x = period, y = value)) +
  geom_boxplot() +
  facet_wrap(~ measure, scales = "free_y") +
  labs(
    title = "Ventilation metrics by window relative to APRV start (ITT)",
    x = "Window (pre / early / late / post)",
    y = "Median value per encounter"
  )

# Print both plots
p_oxygenation
p_ventilation

Text narrative here: e.g. “Median FiO₂ decreased from … to … over the early APRV period, while SpO₂ remained stable/increased.”

3.4 Neuromuscular blockade (paralysis)

Although sedation depth is not available, we can describe neuromuscular blockade use:

Proportion of APRV vs non-APRV ventilated patients with any nmb_on == TRUE
Duration/proportion of ventilated time with paralysis (if feasible)

Code

# Neuromuscular blockade summary for APRV vs non-APRV patients

# APRV time-series only, with NMB flag
nmb_aprv_ts <- aprv_ts |>
  filter(is_aprv) |>
  mutate(nmb_flag = nmb_on %in% c(TRUE, 1, "Yes", "Y"))

# One row per encounter: any NMB on APRV and proportion of APRV time paralysed
nmb_aprv_by_encounter <- nmb_aprv_ts |>
  group_by(encounterId) |>
  summarise(
    any_nmb_on_aprv  = any(nmb_flag, na.rm = TRUE),
    prop_nmb_on_aprv = mean(nmb_flag, na.rm = TRUE),
    .groups = "drop"
  )

# Restrict to APRV encounters
aprv_encounters_nmb <- vent_encounters |>
  filter(ever_aprv_flag) |>
  left_join(nmb_aprv_by_encounter, by = "encounterId")

# Summary numbers
n_total <- nrow(aprv_encounters_nmb)

n_any <- aprv_encounters_nmb |>
  summarise(n = sum(any_nmb_on_aprv, na.rm = TRUE)) |>
  pull(n)

exposed <- aprv_encounters_nmb |>
  filter(any_nmb_on_aprv, !is.na(prop_nmb_on_aprv))

med <- median(exposed$prop_nmb_on_aprv, na.rm = TRUE)
q1  <- quantile(exposed$prop_nmb_on_aprv, 0.25, na.rm = TRUE)
q3  <- quantile(exposed$prop_nmb_on_aprv, 0.75, na.rm = TRUE)
mx  <- max(exposed$prop_nmb_on_aprv, na.rm = TRUE)

# Build a clean presentation table
nmb_aprv_tbl <- tibble::tibble(
  Measure = c(
    "Any neuromuscular blockade during APRV",
    "Proportion of APRV time paralysed (among those receiving NMB)",
    "Maximum proportion of APRV time paralysed (among those receiving NMB)"
  ),
  Value = c(
    paste0(n_any, " / ", n_total, " (", round(100 * n_any / n_total), "%)"),
    paste0(round(100 * med), "% (", round(100 * q1), "–", round(100 * q3), "%)"),
    paste0(round(100 * mx), "%")
  )
)

nmb_aprv_tbl |>
  gt() |>
  tab_header(title = "Neuromuscular blockade during APRV") |>
  cols_label(Measure = "Measure", Value = "Summary")

Table 3

Neuromuscular blockade during APRV
Measure	Summary
Any neuromuscular blockade during APRV	24 / 105 (23%)
Proportion of APRV time paralysed (among those receiving NMB)	35% (25–63%)
Maximum proportion of APRV time paralysed (among those receiving NMB)	100%

Text narrative interpreting whether APRV is associated with more/less paralysis.

3.5 Clinical outcomes

3.5.1 APRV duration and mechanical ventilation duration

Code

# Summaries of APRV duration and ventilation duration
# APRV duration applies only to APRV patients
# Ventilation duration applies to both APRV and non-APRV groups

library(gtsummary)

duration_data <- vent_encounters |>
  left_join(vent_durations, by = "encounterId") |>
  left_join(aprv_by_encounter, by = "encounterId") |>
  transmute(
    aprv_group,
    vent_duration_hours = vent_duration_hours.x,  # or whichever is correct
    aprv_total_hours
  )

duration_table_vent <- duration_data |>
  select(aprv_group, vent_duration_hours) |>
  gtsummary::tbl_summary(
    by = aprv_group,
    statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"),
    missing = "no",
    label = list(
      vent_duration_hours ~ "Duration of mechanical ventilation (hours)"
    )
  ) |>
  gtsummary::modify_header(label ~ "**Duration metric**") |>
  gtsummary::bold_labels()

duration_table_aprv <- duration_data |>
  filter(aprv_group == "APRV") |>
  select(aprv_total_hours) |>
  gtsummary::tbl_summary(
    statistic = list(all_continuous() ~ "{median} ({p25}, {p75})"),
    missing = "no",
    label = list(
      aprv_total_hours ~ "Duration of APRV (hours)"
    )
  ) |>
  gtsummary::modify_header(label ~ "**APRV duration (APRV patients only)**") |>
  gtsummary::bold_labels()

duration_table_aprv
duration_table_vent

Table 4

APRV duration (APRV patients only)	N = 105¹
Duration of APRV (hours)	33 (13, 75)
¹ Median (Q1, Q3)

Duration metric	No APRV N = 1,381¹	APRV N = 105¹
Duration of mechanical ventilation (hours)	65 (18, 125)	229 (134, 346)
¹ Median (Q1, Q3)

3.5.2 PICU length of stay and 28-day mortality

Code

# Summaries of outcomes for APRV vs non-APRV groups:
#   - PICU length of stay (days)
#   - 28-day PICU mortality
#   - Optionally: PIM3-risk–stratified descriptive summaries

library(gtsummary)

# ---- 1) Select relevant outcome variables ----
outcomes_data <- vent_encounters |>
  select(
    aprv_group,
    los_days,
    death_28d_picu,
    pim3_death_risk
  )

# ---- 2) Create standard table of outcomes ----
tbl_outcomes <- outcomes_data |>
  gtsummary::tbl_summary(
    by = aprv_group,
    missing_stat = "no",
    missing_text = "no",
    statistic = list(
      los_days ~ "{median} ({p25}, {p75})",
      pim3_death_risk ~ "{median} ({p25}, {p75})",
      death_28d_picu ~ "{n} ({p}%)"
    ),
    label = list(
      los_days = "PICU Length of Stay (days)",
      death_28d_picu = "28 day mortality",
      pim3_death_risk = "PIM3 predicted mortality (%)"
    ),
    missing = "no"
  ) |>
  gtsummary::modify_spanning_header(c("stat_1", "stat_2") ~ "**Group**") |>
  gtsummary::modify_header(label ~ "**Outcome**") |>
  gtsummary::bold_labels()

tbl_outcomes

Table 5

Outcome	Group
Outcome	No APRV N = 1,381¹	APRV N = 105¹
PICU Length of Stay (days)	6 (3, 10)	15 (8, 24)
28 day mortality	66 (4.8%)	9 (8.6%)
PIM3 predicted mortality (%)	2 (1, 5)	3 (1, 11)
¹ Median (Q1, Q3); n (%)

Interpret patterns: e.g. “APRV patients had longer ventilation and LOS, consistent with higher baseline severity (higher PIM3). Raw mortality was X% vs Y%.”

4. Discussion

Summarise:

How APRV is being used in your unit (who gets it, when, and for how long)
What happens to oxygenation/ventilation parameters when APRV is started
Apparent relationships with paralysis use
How outcomes compare with non-APRV ventilated patients (with caution about confounding)

Include:

Strengths: whole-population, high-resolution time-series data.
Limitations: observational, selection bias (sicker patients may get APRV), incomplete sedation data, etc.
Clinical implications and future work.

5. Appendices

5.1 Variable dictionary (selected)

#| label: tbl-variable-dictionary
#| echo: false
# Target: small table documenting:
#   - key variables used
#   - their definitions / units
#   - any transformations

5.2 Sensitivity analyses (optional)

Placeholder for any extra analyses you might add later, e.g.:

Limiting to first PICU admission per patient
Restricting to patients with certain diagnoses (e.g. ARDS vs other)
Alternative definitions of time windows

#| label: sensitivity-analyses
#| echo: false
# Placeholder for any future extended analyses.