How get survival function estimates from Aalen's additive regression?

Question

The Aalen model assumes that the cumulative hazard H(t) for a subject can be expressed as a(t) + X B(t), where a(t) is a time-dependent intercept term, X is the vector of covariates for the subject (possibly time-dependent), and B(t) is a time-dependent matrix of coefficients.

My intuition was that coefficient plots from plot.aareg produce cumulative hazards for each explanatory variable and intercept (or time-dependent baseline hazard). But should the cumulative hazard function be non-decreasing?

Anyhow, is there any way to get the survival function from this model?

library(survival)
library(tidyr)
library(dplyr)
lfit <- aareg(Surv(time, status) ~ sex , data=lung,
              nmin=1)
plot(lfit)
tibble(t=lfit$times, coef=lfit$coefficient[,'Intercept']) %>%
  group_by(t) %>%
  summarise(baseline_coef_t=sum(coef)) %>%
  mutate(baseline_cumhaz=cumsum(baseline_coef_t),
         baseline_survival=exp(-baseline_cumhaz))

Answer 1

Klein and Moeschberger (Second Edition) devote Chapter 10 to additive models. Practical Notes to Section 10.2 discuss this matter.

3. The estimates of the baseline hazard rate are not constrained to be nonnegative by this least-squares estimation procedure...

Similarly, with respect to cumulative hazards and associated survival curves either via Nelson-Aalen (your estimate) or product-limit approaches, they warn:

6. ... Some care is needed in interpreting either estimate because $\hat H(t | Z)$ [estimated cumulative hazard] need not be monotone over the time interval.

What you found is thus to be expected. You might consider using the slope values returned for the model as linear approximations to the cumulative hazard. From the manual page for summary.aareg():

The slope is based on a weighted linear regression to the cumulative coefficient plot, and may be a useful measure of the overall size of the effect... (Of course the plots [of coefficients over time] are often highly non-linear, so it is only a rough substitute).

You might consider flexible fits of the cumulative hazards that enforce monotonicity. This page has some suggestions. You probably need to weight the individual estimates over time similarly to how the slope values are estimated by summary.aareg().

My sense, as a non-expert in additive models, is that they are best used for evaluating associations of predictors with outcome in a highly flexible way, rather than for modeling survival functions.