Piecewise Linear Mixed Models With a Random Change Point

I’m also interested in fitting a linear mixed model with a random change point. However, my dataset has an additional categorical predictor, “group”, with 3 levels. I’m trying to incorporate a fixed effect of group into @paul.buerkner’s example code above, but am currently at a bit of a loss.

For reference, if we ignore the change point for a moment, my linear mixed model formula would look as follows:
y ~ age * group + (1 + age | person)

However, I have reason to believe that my data would be better fit by introducing a change point. Ideally, I’d like to estimate the effect of my categorical predictor “group” on the change point as well. As a start, I first tried adapting the code above by just including a main effect of “group” as follows:

bform <- bf(
    # intercept and main effect of group
    y ~ b0 + b1 * group +
        # pre-change slope
        b2 * (age - omega) * step(omega - age) +
        # post-change slope
        b3 * (age - omega) * step(age - omega),
    # intercept, slopes and change point varying by person
    b0 + b2 + b3 + omega ~ 1 + (1 | person),
    # fixed effect of group
    b1 ~ 1,
    nl = TRUE
)

However, this gave me the following error message:

Error: Factors with more than two levels are not allowed as covariates.

I saw that this issue came up previously, but I couldn’t really follow the solution in that thread. It seems I should incorporate the non-linear parameters in a linear formula, but I’m not sure how to go about this.

Thanks in advance for any advice!