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