Additive hazard model: Estimating a martingale residual process Consider an additive hazard model with a single categorical covariate
$$ \alpha(t | x_i) = \beta_0(t) + \beta_1(t) x_i, $$
where e.g., $x_i \in \{0, 1\}$.
To assess the goodness-of-fit of this model, I would like to estimate the martingale residual process
$$M(t) = N(t) - \hat{\Lambda}(t),$$
where $N(t)$ and $\hat{\Lambda}(t)$ are $n$-dimensional vectors of the observed counting processes and the estimated cumulative intensity processes, respectively.
When I try to estimate the martingale residual process $M(t)$ in R using timereg::cum.residuals, I get an error message indicating that I should specify the additive model above with a continuous covariate.
Why is it required to include a continuous covariate in an additive hazard model before estimating a martingale residual process?
Here's a minimal example in R:
df <- data.frame(
  t = c(5, 10, 40, 80, 120, 400, 600),
  x = c(12, 10, 3, 5, 3, 4, 1),
  d = c(0, 1, 0, 0, 1, 1, 0)
)

df$x_cat <- cut(df$x, c(0, 5, 15))

fit <- timereg::aalen(survival::Surv(t, d) ~ x_cat, df, residuals = 1, n.sim = 0) # error

resid <- timereg::cum.residuals(fit, df, cum.resid = 1)

plot(resid, score = 2)

which produces the following error message
Error in timereg::cum.residuals(fit, df, cum.resid = 1) : 
   No continous covariates given to cumulate residuals 

Try running the code above using the continuous covariate x instead of the factor x_cat when specifying the regression model.
 A: These covariate-associated cumulative residuals are typically used for evaluating the functional form of a covariate in the model. For example, is the general single-predictor form that you wrote
$$ \alpha(t | x_1) = \beta_0(t) + \beta_1(t) x_1 $$
correct as it stands? Or does the predictor $x_1$ need to be transformed in some way to describe the data properly? For example, perhaps the true functional form is $f(x_1)$, for a model:
$$ \alpha(t | x_1) = \beta_0(t) + \beta_1(t) f(x_1). $$
That's clearly an issue with a continuous predictor. Plots of cumulative martingale residuals associated with a continuous predictor can indicate whether or not an untransformed $x_1$ is adequate.* See the Martinussen and Scheike text "Dynamic Regression Models for Survival Data," Springer (2006), in particular Section 5.7 on goodness-of-fit procedures.
If $x_1$ is binary, however, then any transformation  of it will simply be equivalent to a rescaling of $\beta_1(t)$. There is thus no transformation to evaluate for a binary predictor. With dummy/treatment coding of a multi-level categorical predictor, you just have a whole set of effectively binary predictors, with the same implication.
The warning from the code isn't exactly correct. Try breaking your $x$ predictor into a 3-level ordered factor:
df$x_cat2 <- cut(df$x, c(0, 3, 5, 15))
df$x_cat2 <- ordered(df$x_cat2)

and your code run with x_cat2 won't give the error found with x_cat . Alternatively, if a "continuous" predictor only happens to have 2 distinct values, the cum.residuals() function should return the same error (although I haven't tried that myself). So what's technically required to use cum.residuals() is to have an ordered predictor of some sort with more than 2 distinct values.
This isn't a limitation on estimating the martingale residual process, however. The aalen() function called with residuals = 1 returns an object with a "residuals" list containing a vector of event times and a matrix with the estimated martingale increments for each event time and case. You can process that information on residuals any way you wish.

*Similarly, in Cox models the smoothed martingale residuals plotted against values of a continuous predictor omitted from the model can indicate a helpful functional form for a transformation of that predictor.
