---
title: "simPreg"
subtitle: "Simulate Pregnancy Data with Time-Varying Exposure"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{simPreg}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
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.

```{r}
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`:

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

Example using `expand = TRUE`:

```{r}
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:

```{r}
# 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.

```{r}
# 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.

```{r}
# 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.