In this vignette we briefly illustrate how to fit PAMMs for left-truncated data. As an example we use the data already discussed in the left-truncation section of the data transformation vignette.

We first we load some libraries as well as the data.

library(eha)
library(pammtools)
library(ggplot2)
theme_set(theme_bw())
data(infants, package = "eha")
head(infants)
##   stratum enter exit event mother age  sex      parish   civst    ses year
## 1       1    55  365     0   dead  26  boy Nedertornea married farmer 1877
## 2       1    55  365     0  alive  26  boy Nedertornea married farmer 1870
## 3       1    55  365     0  alive  26  boy Nedertornea married farmer 1882
## 4       2    13   76     1   dead  23 girl Nedertornea married  other 1847
## 5       2    13  365     0  alive  23 girl Nedertornea married  other 1847
## 6       2    13  365     0  alive  23 girl Nedertornea married  other 1848

The header shows that we have information on entry (enter) and exit (exit) times into the riskset. Survival is indicated by the event variable. The variable of interest is mother, i.e., whether mother was alive or dead. For each child who’s mother died, 2 children of the same age and same values for the other covariates were included into the study who’s mother was still alive at inclusion age. The maximum follow up was 365 days.

As a proof of concept we show that we can recreate the effects found by the Cox model for left truncated data as implemented in eha::coxreg:

# fit cox model for left-truncated data
fit <- coxreg(
  formula = Surv(enter, exit, event) ~  mother,
  data = infants)
# fit pam to left-truncated data
infants_ped <- infants %>% as_ped(Surv(enter, exit, event)~.)
pam <- pamm(
  formula = ped_status ~ s(tend)+ mother,
  trafo_args = list(formula = Surv(enter, exit, event)~.),
  data = infants)
# compare coefficients
cbind(coef(pam)[2], coef(fit))
##               [,1]    [,2]
## motherdead 1.68071 1.69185

The coefficients of the mother variable are very similar and indicate that children with diseased mother’s are about 5.37 times as likely to die within one year after birth.

## compare (baseline) cumulative hazard probabilities
base <- survival::basehaz(fit)
ndf <- make_newdata(infants_ped, tend = unique(tend), mother = c("alive")) %>%
  add_cumu_hazard(pam)
ggplot(ndf, aes(x = tend, y = cumu_hazard)) +
  geom_line() +
  geom_ribbon(aes(ymin = cumu_lower, ymax = cumu_upper), alpha = .3) +
  geom_step(data = base, aes(x = time, y = hazard), col = 2)