Modelling cumulative mean/median rehabilitation in time

I have long-format data including patients’ received rehabilitation within one-year time period. I am interested in modelling (showing) two things simultaneously: (1) total cumulative mean/median received rehabilitation at the end of the follow up, (2) duration of the rehabilitation provision. Would it be possible to model such thing with brms or rstanarm and how to do it?

My dream is to get something like this from it

Data structure


df = read_csv("")


Plotted cumulative hours for one patient

df %>% 
  filter(ID == 327) %>% 
  ggplot(aes(day_of_rehabilitation, cumsum(hours_received)))+
  ggtitle("Line for patient number 327")