inst/doc/simPreg.R

In simPreg: Simulate Pregnancy Data with Time-Varying Exposure

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

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

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

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

## -----------------------------------------------------------------------------
# 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)
)

## -----------------------------------------------------------------------------
# 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)
)

## -----------------------------------------------------------------------------
# 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)
)