Arjas Plots in Survival Analysis I would like to construct an Arjas plot to test the proportional hazards assumption. The example data set I am using can be found here: https://stats.oarc.ucla.edu/sas/seminars/sas-survival. I have googled "Arjas plot", I have tried contacting one person who is familiar with them and my email has gone unanswered, and I tried contacting Klein and Moeschberger who authored Survival Analysis... (They both have died.)

Following Section 11.4 of the second edition of Klein and Moeschberger, an Arjas plot starts with a Cox model of survival as a function of covariates $Z^*$, then plots the model-predicted against the observed number of events for each level $g$ of an additional categorical predictor $Z_1$ as time increases. A continuous predictor $Z_1$ needs to be broken down into groups.
The idea is that if the predictor $Z_1$ isn't important, then the plots for all the groups should follow the identity line. If $Z_1$ is important and meets the proportional hazards assumption, then each plot should be close to a straight line but at a different angle. If it doesn't meet the proportional hazards assumption, the plots should be other than straight lines, indicating changes in hazard over event times.
The observed number of events through time $t_i$ for cases $j$ in group $g$ is simply:
$$N_g(t_i)=\sum_{Z_{1j}=g} \delta_j I(T_j \le t_i),$$
where $\delta_j$ is the case's 0/1 censored/event indicator and $T_j$ is the censoring/event time. That's just a running sum of events in that group through each event time.
The predicted total number of events through time $t_i$ for cases $j$ in group $g$ is based on the cumulative hazard $\hat H$, conditional on the covariates $Z_j^*$, at the lower of $t_i$ or $T_j$:
$$\text{TOT}_g(t_i) =  \sum_{Z_{1j}=g} \hat H(\min(t_i,T_j)|Z_j^*).$$
In other words, the sum at any time $t_i$ is over all cases in the group regardless of time, but the cumulative hazard for case $j$ after its censoring/event time $T_j$ is kept constant at its value at $T_j$. In R, you can get the full matrix of estimated cumulative hazards for all cases from a coxph model based on a dataset originalData as:
cumhazModel <- survfit(coxphModel,newdata=originalData)$cumhaz

which gives one row per event time and one column for each case. It is simplest if you put the originalData into the order of increasing event times to start.
Collect the columns corresponding a group $g$. Then, in each column $j$, replace all cumulative hazards after the case's censoring/event time $T_j$ with the value in place at its $T_j$. Row sums of the resulting matrix provide the predicted number of events for group $g$, $\text{TOT}_g$, for each event time. Going up in time, plot those values at each event time against the corresponding number of total observed events for the group.
This github page shows an implementation in R that reproduces the plots in Klein and Moeschberger.
