simPreg
In simPreg: Simulate Pregnancy Data with Time-Varying Exposure

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

The R package simPreg simulates pregnancy data using hazard rates derived from Norwegian registry data as default inputs. Pregnancies can end in spontaneous or non-spontaneous live birth from week 22, or in late miscarriage/stillbirth from week 12. Because the distribution of early miscarriages (before week 12) is uncertain, these are omitted from the simulation. A time-varying exposure is also simulated, with default hazard rates derived from the distribution of COVID-19 vaccination timing during pregnancy in Norway in 2020–2023.

Simulating pregnancy outcome proportions

The function simPregProp() calculates pregnancy outcome proportions across gestational age and exposure timing within a time-to-event framework.

library(simPreg)
df_prop <- simPregProp()
head(df_prop)

Sampling pregnancy data

The function simPregSamp() generates a sample of user-specified size n from the proportions produced by simPregProp(). The logical argument expand determines whether the sampled frequencies are returned as grouped counts or expanded to one row per pregnancy.

Example using expand = FALSE:

set.seed(33)
df_samp_grp <- simPregSamp(
  df = df_prop,
  n = 100000,
  expand = FALSE
)
head(df_samp_grp)

Example using expand = TRUE:

set.seed(33)
df_samp_one <- simPregSamp(
  df = df_prop,
  n = 100000,
  expand = TRUE
)
head(df_samp_one)
table(df_samp_one$Outcome)

All hazard rate vectors can be replaced by user-specified vectors, provided that they have the same length. The default vectors have length 301 and are derived by survival analysis of real-world data.

Specifying hazard ratios

By default, no effect of exposure on any pregnancy outcome is used in the simulations. Users can specify hazard ratios that are constant throughout gestation (proportional hazards) or vary over gestation.

Specifying an increased (proportional) hazard of late miscarriage/stillbirth following exposure:

# HR = 5
df_prop_sb <- simPregProp(hr.late.miscarriage.stillbirth = rep(5, 301))
set.seed(34)
df_samp_sb <- simPregSamp(
  df = df_prop_sb,
  n = 100000,
  expand = TRUE
)

# Calculate proportions among exposed and unexposed pregnancies
prop_sb_exposed <- mean(
  df_samp_sb$Outcome[!is.na(df_samp_sb$ExpGA)] ==
    "late_miscarriage_stillbirth"
)
prop_sb_unexposed <- mean(
  df_samp_sb$Outcome[is.na(df_samp_sb$ExpGA)] ==
    "late_miscarriage_stillbirth"
)
c(
  exposed = round(prop_sb_exposed, 3),
  unexposed = round(prop_sb_unexposed, 3)
)

A non-null hazard ratio for preterm birth following exposure can be specified by setting the HR to a value different from 1 during the preterm period.

# HR = 2 from gestational day 154 through day 258, and 1 otherwise
hr_preterm <- c(rep(1, 153), rep(2, 105), rep(1, 43))
df_prop_pt <- simPregProp(
  hr.spont.livebirth = hr_preterm,
  hr.nonspont.livebirth = hr_preterm
)
set.seed(35)
df_samp_pt <- simPregSamp(
  df = df_prop_pt,
  n = 100000,
  expand = TRUE
)

# Calculate proportions among exposed and unexposed pregnancies
prop_pt_exposed <- mean(
  df_samp_pt$GA[!is.na(df_samp_pt$ExpGA)] >= 154 &
    df_samp_pt$GA[!is.na(df_samp_pt$ExpGA)] <= 258
)
prop_pt_unexposed <- mean(
  df_samp_pt$GA[is.na(df_samp_pt$ExpGA)] >= 154 &
    df_samp_pt$GA[is.na(df_samp_pt$ExpGA)] <= 258
)
c(
  exposed = round(prop_pt_exposed, 3),
  unexposed = round(prop_pt_unexposed, 3)
)

Analyzing simulated data

The simulated data can be analyzed with a Cox proportional hazards model from the survival package, treating exposure as time-varying. Pregnancies are considered unexposed before exposure and exposed thereafter, with follow-up represented in start–stop format.

# Exclude late miscarriage/stillbirth
df_samp_pt <- subset(
  df_samp_pt,
  Outcome != "late_miscarriage_stillbirth"
)

# Define preterm birth
df_samp_pt$preterm <- as.integer(df_samp_pt$GA <= 258)

# Start follow-up at day 153 to allow events from day 154 through day 258,
# and censor at day 258 otherwise
df_samp_pt$tstart <- 153
df_samp_pt$tstop <- pmin(df_samp_pt$GA, 258)

# Assign unique identifiers
df_samp_pt$id <- seq_len(nrow(df_samp_pt))

# Create start-stop data with exposure as a time-varying covariate
df_samp_pt_tv <- survival::tmerge(
  data1 = df_samp_pt,
  data2 = df_samp_pt,
  id = id,
  tstart = tstart,
  tstop = tstop,
  preterm = event(tstop, preterm),
  exposed = tdc(ExpGA)
)

# Fit Cox proportional hazards model
cox_pt <- survival::coxph(
  survival::Surv(tstart, tstop, preterm) ~ exposed,
  data = df_samp_pt_tv
)

# Display estimated HR and 95% CI
c(
  HR = round(exp(coef(cox_pt)), 3),
  lower_95 = round(exp(confint(cox_pt)[1]), 3),
  upper_95 = round(exp(confint(cox_pt)[2]), 3)
)

This example demonstrates how simPreg can be used to simulate pregnancy data with time-varying exposures and specified exposure effects. The simulated data can be analyzed using standard time-to-event methods.