Models and methods for analysing clustered recurrent hospitalisations in the presence of COVID-19 effects

Abstract Recurrent events such as hospitalisations are outcomes that can be used to monitor dialysis facilities’ quality of care. However, current methods are not adequate to analyse data from many facilities with multiple hospitalisations, especially when adjustments are needed for multiple time scales. It is also controversial whether direct or indirect standardisation should be used in comparing facilities. This study is motivated by the need of the Centers for Medicare and Medicaid Services to evaluate US dialysis facilities using Medicare claims, which involve almost 8,000 facilities and over 500,000 dialysis patients. This scope is challenging for current statistical software’s computational power. We propose a method that has a flexible baseline rate function and is computationally efficient. Additionally, the proposed method shares advantages of both indirect and direct standardisation. The method is evaluated under a range of simulation settings and demonstrates substantially improved computational efficiency over the existing R package survival. Finally, we illustrate the method with an important application to monitoring dialysis facilities in the U.S., while making time-dependent adjustments for the effects of COVID-19.


Introduction
In order to identify poor performance and intervene as necessary, indices of a healthcare facility's performance are computed and compared to a benchmark or norm by the Centers for Medicare and Medicaid Services (CMS).
For facilities with poor outcomes relative to the national norm, surveillance may be reinforced, and financial penalties imposed, in rare cases leading to loss of accreditation.Because these evaluations are high stakes, and the quality assessment depends on the accuracy of the evaluations, it is crucial to use appropriate statistical methods to analyse facility-level outcomes while accounting for differences in the characteristics of treated patients.
This work is motivated by the development of quality measures for dialysis facilities.The Medicare data for 2020 involves 7,979 US facilities and 509,609 patients who were diagnosed with kidney failure prior to 30 September 2020.The quality of dialysis treatment greatly affects dialysis patients' quality of life.
Hospitalisations are an important indicator of patient morbidity and quality of life.They are also very costly, accounting for approximately 32% of total Medicare expenditures for dialysis patients in the U.S. (United States Renal Data System, 2020).On average, dialysis patients had 1.58 hospital admissions per person-year in 2018 (United States Renal Data System, 2021).The frequency of hospitalisations is thus an important measurement in monitoring dialysis facilities.To accurately assess facility performance and reduce patient hospital admissions, a risk-adjusted standardised measure, the standardised hospitalisation ratio (SHR), has been used by CMS on the dialysis facility reports site (DialysisData, 2022).Numerically, the SHR = O/E, where O is the observed number of hospital admissions of a facility and E is the expected number of hospitalisations, given patient characteristics, that would be observed if the facility had the same hospitalisation rate as the national norm.An SHR greater than 1.0 indicates that the facility's hospitalisation rate is worse than expected based on overall national rates with adjustment for the measured characteristics of patients of this facility.In particular, the SHR has been implemented by CMS in the End-Stage Renal Disease Quality Incentive Program (ESRD QIP), a value-based purchasing program that links a dialysis facility's payment to the quality of the care it provides.
A facility's 'expected' count (the denominator of the SHR), typically stems from a proportional rates model (Kalbfleisch & Prentice, 2002;Lawless & Nadeau, 1995;Lin et al., 2000), which is a recurrent event analogue of the well-known Cox relative risk model (Cox, 1972).While successful for moderate sample sizes and small-dimensional data, existing statistical software for these methods does not scale to our motivating data set because of the large sample size and high-dimensional facility effects.For example, when implementing the proportional rates model, existing software typically requires dividing each patient's follow-up time into a set of records with at most one recurrent event in each record.Duplicating the covariates in each record is time-consuming and results in a large amount of redundant information.In our motivating example, the data set had 1,007,094 observations and 257 variables, which can become unduly large if further split by events.
Estimation of the effects of many covariates, especially time-dependent ones, and estimation of high-dimensional facility effects further increase the computational burden.
Another important aspect of this example is that patient hospital admissions may vary along two time scales, and modelling these time scales simultaneously adds another level of complexity.First, the frequency of hospitalisation depends on time since the initiation of kidney failure, as shown by Liu et al. (2012), by the United States Renal Data System (2021), and by our estimated model discussed below.Second, the rate of hospital admissions varies with calendar time, with more hospital admissions in winter and fewer in the summer.In addition, in 2020, as an effect of the coronavirus disease 2019 (COVID-19) pandemic, the hospitalisation rate dropped sharply in early spring and gradually picked up in summer.Both the seasonal effect and pandemic impact are evident in Figure 1.Such features require a model to account for the variation over calendar time.
To accommodate high-dimensional clustered recurrent events data, Liu et al. (2012) proposed a two-stage model with a piecewise constant baseline.In Stage I, this method fits a model stratified by facilities, and in Stage II, fits an unstratified model to estimate the overall baseline rate function, using an offset defined with the estimated β obtained from Stage I.This model solves the high-dimensional problem by not estimating facility effects directly but computing the expected number E (defined above) using the baseline rate function from Stage II and β from Stage I, and adjusts for the calendar-time effect with some time interval indicators.The two-stage model could allow for a time-varying rate over calendar time by introducing more intervals, but an already very complex data structure would become even more complicated.
The SHR = O/E is an example of indirect standardisation (Breslow & Day, 1985;He & Schaubel, 2014;Inskip, 2005;Roessler et al., 2021), where both O and E are based on the characteristics of the patients in a facility.Although indirect standardisation is a useful measure for facilities to evaluate themselves or for a governing body to evaluate a facility's events compared to the national norm, it has been argued that indirect standardised measures for different facilities should not be compared, since each facility-specific estimator is essentially adjusted to a different covariate distribution (Breslow & Day, 1985;He & Schaubel, 2014;Inskip, 2005;Roessler et al., 2021).To overcome this disadvantage, direct standardisation can be used, which gives SHR * = E * /O * , where O * is the total observed number of hospitalisations in all facilities and E * is the expected number of hospitalisations in the whole population if all facilities had hospitalisation rates as in the given facility.Because the same standard population is applied to all facilities, direct standardised measures are directly comparable (He & Schaubel, 2014;Inskip, 2005).However, direct standardisation is more difficult for users to understand and may seem less relevant, since the comparison is based on a different patient mix than they treat.In contrast, indirect standardisation is focused on the given patient mix and is more intuitive.Since both ways of standardisation have advantages and drawbacks, it would be useful to have an approach that admits both direct and indirect interpretations.
To address these issues, we propose a computationally efficient modelling procedure for clustered recurrent events data, which considers both the calendar-time effect and time since initiation of kidney failure.It makes the following contributions: First, the proposed method fits in large data settings.It updates the high-dimensional facility effects using a fixed-point algorithm, which works well even when the number of facilities is large, and estimates the model parameters directly without substantial manipulation or replication of the data.This saves computation time and memory.Second, by using calendar time as the baseline, the proposed method captures various trends of hospitalisations in multiple time scales.The baseline rate is a non-parametric function of calendar time that can reflect seasonal or other trends.The time since the start of dialysis is then modelled as a risk adjustment.Third, the proposed model obtains the same facility effect estimates (or SHRs) from indirect and direct standardisation, and thus shares the classic indirect standardisation's ease of computation and interpretation, but also provides estimates that are comparable across facilities.
The rest of the article is organised as follows: Section 2 presents the model and estimation algorithm.In Section 3, we conduct simulations to evaluate the time-and-memory usage and accuracy of the proposed method.In Section 4, we apply the proposed model and method to a real-data application of hospitalisations of US Medicare dialysis patients in the presence of COVID-19.Some concluding remarks are given in Section 5. Additional simulation results and real-data summary and estimates are relegated to the online supplementary material.Weekly hospitalization rate in 2018 to 2020 2 Model specification and estimation

Set-up and notation
Let j = 1, . . ., F index the facility, with facility size n 1 , . . ., n F , where F is the total number of facilities.Let i = 1, . . ., n index the patient, where n =  F j=1 n j is the total number of patients.Let B i , C i , and D i denote the calendar time of beginning to be at risk, right censoring, and death for patient i, respectively.Patients are said to be 'left truncated', if they enter the study after the start of the observation period, which is slightly different from the usual definition.Define the follow-up time (a, b).Define the at-risk process by Y i (t) = 1(B i ≤ t ≤ X i ), with 1(•) being the indicator function and t being a time point in a year.Let N * i (t) denote the cumulative number of events for patient i in time interval (0, t], and let N i (t)= ∫ t 0 Y i (s) dN * i (s) denote the counting process for the observed number of events, where dN * i (t) is the number of events of the patient i at time t.Because the recurrent event of hospitalisation is stopped by the terminal event of death, Assume one patient can have at most one event at a time; i.e. dN * i (t) ≤ 1.Let Z i (t) = (Z i1 (t), Z i2 (t), . . ., Z ip (t)) T be the possibly time-dependent p-dimensional covariate vector.Let G i be the facility index associated with patient i, and, finally, define Let N j (t) be the counting process for the number of events in facility j.Let dN.(t) =  F j=1 dN j (t) = n i=1 dN i (t) be the total observed number of events at time t.
With reference to Lin et al. (2000), we assume a semi-parametric proportional rates model.Conditional on patient i from facility j being alive, where dμ 0 (t) is the baseline rate function, α = (α 1 , . . ., α F ) is the facility effect parameter, and β = (β 1 , . . ., β p ) T is the coefficient for covariate effects.Under the assumption of independent left truncation and censoring, we have (2)

Estimation
Fitting risk adjustment models to recurrent event data with very large samples is computationally challenging, especially when the numbers of healthcare facilities and the candidate adjustment covariates are large and time dependent.Most existing statistical methods fail in large-scale settings because of lack of computational power.As there are almost 8,000 facilities in our motivating data set, treating facilities as categorical variables in the model comes with computation problems, and it requires special methods to calculate the inverse of the observed information matrix.The method we propose circumvents including thousands of indicator variables for facilities by using a fixed-point algorithm.
We first estimate covariate effects β with a stratified model, and then update α by an explicit formula in each iteration.
Let t 1 < t 2 < • • • < t M be the unique event time, t 1 ≥ 0, and τ ≥ t M be the end of the study period.Let F j,m be the set of patients in facility j who have an event at t m , m = 1, 2, . . ., M, and dN j (t m ) be the number of events in facility j at t m ; i.e. dN j (t m ) is the cardinality of F j,m .Let dN.(t m ) =  F j=1 dN j (t m ) be the total number of events at t m .
To estimate β, we account for tied event times using the Breslow (1974) estimator for a stratified model.Any of the standard approximation would work here because the number of events compared to the number of at-risk patients is small.
The estimator is the solution to the following estimating equation: where Next, given α and β, an estimating equation for dμ 0 (t) is dN.(t) = dμ 0 (t)S (0) (t), where , which gives rise to the Nelson-Aalen estimate.
and dμ 0 (t) = 0 elsewhere.For given dμ 0 (t) and β, one can also estimate α from the unbiased estimating equation Plugging in dμ 0 (t) and β, one can obtain where O j is the observed number of hospital admissions of facility j and E j is the (estimated) expected number of hospital admissions if the facility j has event rates as in the national norm: However, the baseline rate function dμ 0 (t) and facility effect α cannot be identified, because for any solution dμ 0 (t) and α, dμ 0 (t) = xdμ 0 (t) and α = α/x is also a solution for any positive constant x.Therefore, we need an appropriate constraint.The constraint might be, for example,  F j=1 n j α j = 0, or median(α j ) = 0. We, however, will consider the constraint,  F j=1 O j =  F j=1 E j , for which the SHR for the whole population equals to 1.Under this constraint, the baseline rate function dμ 0 (t), the limit of dμ (s)  0 (t) defined below, can be interpreted as the national norm.Let where the superscript (s) is for estimation from the sth iteration.We then estimate α(s) j by a fixedpoint algorithm which satisfies the constraint j in each iteration.In this algorithm, C (s) converges to 1, and exp (α j ) = O j /E j = SHR j .The proposed estimation procedure is summarised in Algorithm 1, where our pre-specified convergence criterion is that the absolute difference in αj between the current and previous iterations is smaller than or equal to 10 −6 , for any j.
The log likelihood function of α j is a continuous and concave function, which guarantees the convergence of this fixed-point algorithm (Armijo, 1966;Goldstein, 1967;Wu et al., 2022).We are optimising the log likelihood function, and only use the scaling factor C to recenter the estimated facility effects.

Standard error of β
Since β is the maximum likelihood estimator (MLE) of the stratified Cox model, which is well studied in the literature, β is consistent, and the asymptotic distribution is provided by Kalbfleisch and Prentice (2002).Thus, under certain regularity conditions, where Σ(β) is the asymptotic variance-covariance matrix.
If we consider a complete intensity model, the asymptotic variance-covariance matrix of β can be estimated by I ( β) −1 , where I ( β) (provided in the online supplementary material) is the observed information matrix for a stratified model evaluated at β = β.The standard error of each βk , k = 1, . . ., p, is estimated by the square root of the corresponding diagonal entry of I ( β) −1 .This, however, is based on the very strong assumption that given the full history up to time t, including the events and the paths of covariates, the rates depend only on the covariates in Z i (t) and therefore often referred to as the naive standard error.
In this analysis, we consider the marginal rate model (Lin et al., 2000).The marginal rate model only models the marginal effect of the covariates but does not condition on the prior events.In this model, we use the sandwich variance estimator of the asymptotic variance-covariance matrix of β where is the weighted average of the covariates of patients in facility G i , the facility associated with patient i, at t m , and Ŝ(0) Gi (t m ) is a norming factor for facility G i at t m , evaluated at β = β.The sandwich estimator is a robust version of the estimated variance-covariance matrix for β.
Algorithm 1 The fixed-point algorithm estimation procedure.
1: Estimate β with a stratified model.
Since exp (α j ) = O j /E j and . A robust estimator of the standard error of αj is the sandwich estimator: In this algorithm, we use a constraint so that the proposed method and the standard approach are fitting the same model and optimising the same likelihood function.In our approach, β and dμ 0 (t) are estimated from  F j=1 n j = n patients, which converge at the rate of n −1/2 .In contrast, α j is estimated by a pseudo-MLE (Gong & Samaniego, 1981) based on only n j patients, which converges at the rate of n −1/2 j ; if F is large, then β and dμ 0 (t) can be considered as fixed and their variation can be ignored, compared to that of αj .We note that, in most profiling applications, the number of facilities and patients is very large compared to n j , so that β and dμ 0 (t) can be precisely estimated.This assumption is common in provider profiling applications, which usually involve data sets with very large sample sizes (He et al., 2019;Jones & Spiegelhalter, 2011;Kalbfleisch & Wolfe, 2013).If there are only a small number of facilities, the methods can be adapted to account for the variability in the estimates of β and dμ 0 (t).

The equivalence between indirect and direct standardisation
As discussed above, it is important to adjust for the differences in patient characteristics before assessing and comparing facilities.An important advantage of the proposed method is that it can be interpreted as either indirect or direct standardisation.This equivalence can be verified as follows: Let E * (j) be the expected number of hospitalisations in the whole population if all patients were treated in the given facility j.Direct standardisation of facility-specific hospitalisation outcomes is defined as where we have imposed the constraint Thus, the direct standardised measure reduces to SHR under the proposed model.

Simulation study
To compare the performance of the proposed method with the current R package survival (version 3.2-7), we implemented the proposed algorithm in R, with functions written in Rcpp language, and assessed its performance using different simulation settings.Because the standard method is not capable of estimating the model parameters in our example data, we carried out simulations on smaller data sets.All scenarios were repeated 1,000 times.While the estimation results from both algorithms were almost the same, the computation time and memory used by our proposed method were much smaller than the survival package.
Five continuous variables were generated from independent normal distributions with mean 0 and variance 0.09, and five binary predictors were generated from independent Bernoulli distributions which take the value 1 with probabilities equal to 0.2, 0.28, 0.36, 0.44, and 0.52.The facility effect α 1 = 0 for the first facility, and α j ∼ Normal(0, 0.2 2 ) for the rest.We generated the data to mimic the motivating data example, that patients arrived on various days over the year.The number of days in the study was τ = 365.Let B i denote the calendar time of beginning to be at risk.For each i, P(B i = 0) = 0.8, and otherwise, B i was chosen randomly from {0, . . ., τ − 1}.Independently of B i , The time at risk for the ith patient was t i = min (τ − B i , t * i ) where P(t * i = τ) = 0.75 and otherwise t * i was chosen randomly from {1, . . ., τ}, so that time to exit was no larger than the number of days in the study.Time to exit X i = B i + t i .We assumed entrances happened at the beginning of a day and events and exits happened at the end of a day, with events preceding exits.
For each patient i in facility j, we generated hospital admissions from a Poisson process with rate ρ i by generating the gap time T l i , l = 1, 2, . . ., as an exponential distribution with rate ρ i = ρ 0 × exp(α j + Z i β), with ρ 0 = 0.003.This model resulted in an event rate similar to our real data application.T l i , l = 1, 2, . . ., was rounded up to the next integer.The process of generating T l i , l = 1, 2, . . ., was repeated until  L l=1 T l i ≥ t i , for some L, and T L i was discarded if  L l=1 T l i > t i .We consider the following six scenarios: Scenario 1: Eight different numbers of facilities were tested: 30, 50, 100, 200, 300, 500, 1,000, and 2,000.In each facility, the number of patients followed a Poisson distribution with mean 100.Scenario 2: The number of facilities was 100, and the number of patients in each facility followed a Poisson distribution with mean from 50 to 5,000.Scenario 3: With other configurations being the same as in Scenario 2, we fixed the size of the first facility and the mean facility size as 50.Scenario 4: We explored the accuracy of β with correlation within patients.With other configurations being the same as in Scenario 3, the events generating model became: gap time T i ∼ exp (ρ ij ) for patient i in facility j, where and W i ∼ Gamma(1, 1) was a patient level random effect.Scenario 5: We fixed the facility effect of the last facility as 0, changing the size of the last facility from 50 to 5,000, as in Scenario 2. Other configurations were the same as in Scenario 3. Scenario 6: We explored the accuracy of α with correlation within patients.The true data generating system was again the equation ( 11).Other configurations were the same as in Scenario 5.
As an illustration, in Scenario 1, the percentage of patients with at risk time t i smaller than the number of days (365) ranged from 39.74 to 40.06, and the percentage of patients who dropped before the end of the study (X i < τ) ranged from 22.32 to 22.71.A patient on average had 1.25 to 1.38 events, and the average number of events per 365 days was 1.58-1.73.Results in other scenarios were similar and are omitted here.For the simulation and the real data application, we set α j = −10 when estimating the effect of facilities with no events, to avoid taking the log of zero.

Evaluation of computation time and memory
We compared the computation time and memory of the proposed method with the survival package in Scenarios 1 and 2. Experiments were conducted on an Intel Ⓡ Xeon Ⓡ Gold 6254 quadprocessor with max frequency 4 GHz and RAM 576 GB.The computation time and memory used were recorded on R (version 3.6.3)using R base function proc.time and R function profmem in R package profmem (version 0.6.0),respectively.The computation time and memory used are summarised in Figure 2. In both scenarios, our proposed method was substantially faster and required much less memory, compared to the survival package; furthermore, as sample size increased, the time used by our method increased more slowly than the survival package.When the number of facilities was 2,000, our method was over 3,000 times faster than the survival package.In addition, the survival package used around 35 GB of memory when the number of facilities was 2,000, which is larger than most computers' memory.In contrast, we had almost 8,000 facilities in our application, and without revision, the survival package could not fit a model on our real data, while our proposed method was able to.Our proposed method was also computationally more efficient when the average facility size increased.

Accuracy of β estimation
We studied the accuracy of estimation of covariate effects in Scenarios 3 and 4. In each repetition, we fitted the model with our proposed method and the survival package, and recorded the bias and standard error.A summary of these statistics is listed in Table 1, where the rows of 'Independent' are the results of Scenario 3, and the rows of 'Correlated' are the results of Scenario 4. Our proposed method had results similar to the survival package, with the sandwich asymptotic standard error (sASE) close to the empirical standard deviation (ESD) and a small bias in estimation; The top two panels are from Scenario 1, in which the mean facility size was fixed at 100, and the number of facilities varied between 30 and 2,000; and the bottom two panels are from Scenario 2, in which the number of facilities was fixed at 100, and the mean facility size varied between 50 and 5,000.
however, our proposed method saved computation time and memory and thus could handle a much larger data set.As is typically the case, the robust estimator gave better results when there was correlation.More results can be found in Tables S1 and S2 of the online supplementary material.

Accuracy of α estimation
We were most interested in the accuracy of the estimates of the facility effects.Therefore, in Table 2, we compared the performance of the proposed method and the survival package in estimating the effect of the last facility, α F , for Scenarios 5 (Independent) and 6 (Correlated).As mentioned above, the facility effects were not completely identifiable since the baseline event rate was flexible.Comparison of the two approaches requires a common choice of reference.In particular, when estimating the facility effects, the 'factor' argument in the survival package will treat facilities as covariates, and let the first facility have 0 effect (as reference), so we deducted the first facility's estimated effect from all estimated facility effects.Therefore, both methods satisfied the constraint that the first facility has 0 effect.Once this was done, the two methods gave similar results in the estimation of the facility effects and their standard errors, except that the variance estimates were better with the robust estimates when correlation was present.Surprisingly, in Tables 1 and 2, the coverage probabilities of the robust confidence interval tends to be slightly smaller than the nominal level, for which we do not have a reason; nevertheless, as Tables 1 and 2 show, the robust method is much better than the naive approach in the correlated scenario.

Application
We applied the proposed method to the study of hospital admissions among Medicare patients on chronic dialysis for kidney failure in US Medicare-certified dialysis facilities, derived from the CMS claims and clinical and administrative databases.Our analyses were based on data from the calendar year 2020, with a non-parametric baseline rate function based on the calendar year.This is especially useful for analysing data from 2020, as there was a drop of hospital admissions in early spring of 2020 due to the effect of COVID-19, as elective medical services were delayed, patients and physicians deferred medical care out of caution, and the US government required hospitals to redistribute resources to COVID-19 cases and people to stay at home (Birkmeyer et al., 2020;Kendzerska et al., 2021;The Council of State Governments, 2021).The change over calendar time can be captured by the proposed non-parametric baseline rate function.
The primary outcome of interest was the hospitalisation sequence.Each patient could have multiple hospital admissions (i.e.recurrent events).Hospital admissions were considered to happen at the end of each day, with arrivals and changes of COVID-19 stages preceding any admissions if on the same day.Patients were subject to left truncation if they started to be at risk after the start of the observation period, 1 January 2020.Subjects were followed until the earliest of the following: death, loss to follow-up, 3 days prior to transplant, or 31 December 2020.
The impact of COVID-19 was allowed to vary with time.After exploration, we arrived at this model: We divided the post diagnosis period into four stages and assumed a separate constant effect for each stage; after diagnosis, the patient remained in the 'COVID' state until there were no further reported COVID-19 diagnoses for a consecutive 21-day period.We divided this 'COVID' state into the first 10 days (COVID1) and the rest (COVID2), with the hypothesis that the impact was most severe during the first few days after infection.After the 21-day period, the patient entered the post-COVID state and remained in that state unless there was another hospitalisation with a COVID-19 diagnosis, in which case he/she became a late-COVID patient.Patients who had an observed previous COVID-19 diagnosis were referred to as in the 'any-COVID' group, i.e. the any-COVID group was the combination of the four stages.We assumed that the COVID-19's interactions with other covariates were the same in all four stages.All patients started in the no-COVID group, i.e. the patients without an observed previous COVID-19 diagnosis.As mentioned above, the Medicare data involve 509,609 patients in 7,979 facilities.Among these patients, 63,521 (12.5%) entered the 'COVID' state during the observation period, 34,375 later became post-COVID patients and 1,900 entered the late-COVID state on or before 31 December 2020.More details can be found in Figure S1 of the online supplementary material.
Other variables included in the model were patient age, gender, race (White, Black, Asian/Pacific Islander, and other), ethnicity (Hispanic, non-Hispanic, and unknown), body mass index (BMI) as a categorical variable, time since the onset of dialysis as a categorical variable, proportion of time with Medicare Advantage (MA) coverage, nursing home status, 13 incident comorbidities, 41 conditions for prevalent comorbidities, an indicator of whether a patient had less than 6 months of Medicare coverage in 2019, and some interaction terms and fixed effects for facilities.MA is supplementary insurance to Medicare, and the proportion of time with MA coverage is a time dependent variable that specifies the percentage of time a patient was on Medicare Advantage in the previous 12 months.Nursing home status had three categories: no nursing home care (0 days), short-term nursing home care (1 day to 89 days), and long-term nursing home care (90 days or more).We defined incident comorbidities as those present when the patient first started dialysis and prevalent comorbidities as those observed in inpatient claims within the year 2019.For any patient who had at least six months of claims in 2019, the comorbidities were identified through Medicare claims; for patients with less than 6 months of Medicare coverage in 2019, no prevalent comorbidities were recorded, but an indicator variable accounted for this situation.
While allowing the effect of COVID-19 to be time-dependent, we examined the proportional rate assumption for other covariates graphically by plotting the Schoenfeld residuals for tied data, (Park & Hendry, 2015;Schoenfeld, 1982), over calendar time for each covariate.Taking age as an example, the formula means that for every time point, we sum the age of all patients who have events at that time point, and minus the number of events times the average age of all at-risk patients at that time point.No other important covariates violate greatly the proportional rate assumption.Example Schoenfeld residual plots can be found in the online supplementary material (Figure S2).The variance of the estimated covariate effects was calculated by the sandwich estimator specified in equation ( 9).
The patients' characteristics are summarised in Table S5 of the online supplementary material, and the model results are summarised in Table 3 and more completely in Table S6 of the online supplementary material.Table 4 provides a comparison of the main characteristics' effect in the no-COVID and COVID-19 group, together with the 95% confidence intervals.As expected, the effect of COVID-19 varied significantly in the four stages.A patient in COVID1 had an extremely high relative rate of hospitalisations of 20.18 (p < 0.001) compared to patients in the no-COVID group.Patients in COVID2 had a still quite large relative rate of 2.10 (p < 0.001).The relative rates for post-COVID and late-COVID patients were respectively 1.37 and 2.04 (p < 0.001 for both), still substantial increases in the rate of hospitalisations, as compared to a no-COVID patient.
Among other adjusted variables, age had a U-shaped effect on hospitalisation, where the 65-to-80-year-old patients had the lowest risk.The estimated age effect curve and baseline rate function are provided in Figure 3. Females on average were more likely to be hospitalised (RR = 1.07, p< 0.001 for all comparisons in this paragraph) than males, if they had no COVID-19 history, but their interactions with previous COVID-19 diagnosis had lower risks (RR = 0.89).In the COVID-19 group, the relative rate of females vs. males was smaller than 1 (RR=1.07× 0.89 = 0.95).Similarly, compared to White patients, Black and Asian/Pacific Islander had lower risks of hospitalisations (RR = 0.93 and 0.77, respectively) in the no-COVID group, but similar risks (RR = 0.93 × 1.09 = 1.01 and RR = 0.77 × 1.40 = 1.08, respectively) in the COVID-19 group; compared to non-Hispanics, Hispanics had a lower relative rate of hospitalisation (RR = 0.88) in the no-COVID group, and a relative rate of about 1.10 = 0.88 × 1.25 in the COVID-19 group.Table 4 provides a comparison of the main characteristics' effect in the no-COVID and COVID-19 group, together with the 95% confidence intervals.While most prevalent comorbidities increased the risks of hospitalisation, the following led to at least a 25% increase in relative rate and had at least 0.5% prevalence in the no-COVID group (RR, prevalence in the no-COVID group): opioid dependence (1.44, 0.7%), pancreatitis (1.44, 0.5%), bipolar disorder (1.30, 0.8%), other liver disease (1.26, 1.0%), diabetes with complications (1.25, 29.4%), and major depressive affective disorder (1.25, 8.1%).Patients with less than 6 months of Medicare coverage in the prior calendar year, 2019, also had a higher risk of hospitalisation (RR: 1.54, 17.4% prevalence in the no-COVID group) due to unknown comorbidities.Patients with 100% Medicare Advantage coverage on average had a lower risk of hospital admissions (RR = 0.93), compared to patients with no Medicare Advantage coverage.
A weekly averaged O and E plot over time for all patients is given in Figure 3, where O is the observed number of hospital admissions, and E is the expected number obtained from the fitted model, assuming that all facilities have the same performance as the national average (i.e. the baseline).A similar O and E plot is given in Figure 4 for New York City and the rest of New York State.The O and E matched well for the rest of New York State whereas the prediction for New York City had a short lag.New York City is especially of interest as it was greatly impacted at the beginning of the pandemic and had a sharp increase of hospitalisations in the spring of 2020.Estimation results for other states were also explored.Plots for Illinois and Texas are provided in Figure 5.The fit was very good.If the model assumption is correct, (O − E)/ �� E √ would  S3 in the online supplementary material, in which we exclude facilities with expected number of hospitalisations under the null (E) greater than 800 or smaller than 5, or the observed number of hospitalisations (O) equal to 0.
Computational efficiency is also an important contribution of our analysis.It took less than 17 hours for the proposed method to fit the model.Parallel computing would further shorten the time used, depending on the number of threads used to compute the results.Table S4 in the online supplementary material shows the time, speedup, and efficiency estimated from simulation data, which provides an idea of how fast the algorithm can be when using parallel computation.We also developed code for models without time-dependent covariates; it took about 3 hours to fit the model for the real data.
During the observation period, many patients died, and it is instructive to see what effect COVID-19 had on mortality and to compare the estimates between the hospitalisation and mortality models.A model with the same covariates was fitted with mortality events being the outcome.Most coefficients had a similar pattern as in the hospitalisation model, but there were some differences.The model estimates are shown in Table 3 and in the online supplementary material (Table S6).On one hand, similar to the hospitalisation model, compared to the no-COVID group, COVID1 patients had the highest hazard ratio (11.77, p< 0.001), followed by COVID2 (6.25, p< 0.001), late-COVID (4.14, p< 0.001), and post-COVID (1.20, p < 0.001).On the other hand, the baseline hazard of mortality increased in early spring (see Figure S4 of the online supplementary material), while the baseline hospitalisation rate decreased.

Discussion
In this paper, we proposed a new modelling procedure for clustered recurrent events data.We developed a computationally efficient algorithm for estimation of high-dimensional facility effects and model parameters.The proposed fixed-point algorithm is very fast and provides essentially the same results as the slower partial likelihood approach.The proposed method used calendar time as the baseline time scale, which has only a fixed number of distinct values (365 or 366) and is easier to estimate, and therefore is more suitable for modelling multiple time scales; other time scales, for example, time since the onset of kidney failure, can be readily adjusted for as a categorical covariate.In addition, the proposed model allows a flexible change of baseline rate over time.This analysis also shows the equivalence between the indirect and direct standardisation under a recommended constraint.The model can be extended to the discrete case (Prentice & Kalbfleisch, 2003).Let 0 < t 1 < t 2 < • • • < t M be the possible event times and t 0 = 0 be the start of the study period.Under the assumption of independent left truncation and censoring, the semi-parametric proportional rates
model is where γ 0k is the instant event probability at time t k for a patient with covariates equal to zero and treated by a facility at the national norm.As mentioned above, we assume the number of events for a patient at a time point is at most 1.The discrete model requires some restrictions, usually not of practical importance, in order to assure the probabilities in equation ( 12) are at most 1.The discrete model ( 12) would provide results nearly identical to those in the continuous model.Note that the Breslow approximation for ties is exactly correct in the discrete model (Prentice & Kalbfleisch, 2003), and the variance adjustment makes little difference.The equivalence between the indirect and direct standardisation also presents in the discrete model.We have focused on frequency of hospitalisations, using the number of hospital admissions as our outcome, although the model also works well on analysing other recurrent events, e.g.duration of hospitalisations, transfusion, and emergency department visits.
Our model assumes that the facility effect is constant over the year but in fact, it is possible that the effect varies somewhat; if so, our estimate would be an 'average' effect.If this were deemed to be important, it would be possible to allow separate effects for various periods during the year, a possible direction for future work.

Figure 1 .
Figure 1.Unadjusted weeklyhospitalisation rates among Medicare dialysis patients at risk on 1 January of each year from 2018 to 2020.The decrease at the end of the year for 2020 is potentially due to reporting delays.

Figure 2 .
Figure 2. Comparison of the computation time and memory used in the survival package and our proposed method.The top two panels are from Scenario 1, in which the mean facility size was fixed at 100, and the number of facilities varied between 30 and 2,000; and the bottom two panels are from Scenario 2, in which the number of facilities was fixed at 100, and the mean facility size varied between 50 and 5,000.

Figure
Figure (a) The estimated baseline hospitalisation rate function averaged by week.(b) The estimated age effect curve.(c) Weekly average number of hospital admissions, observed (O) and expected (E), for all patients.(d) (O − E)/ �� E √ plot for all patients.If the model assumption is correct, (O − E)/ �� E √ would approximately follow the standard normal distribution.

Figure 4 .
Figure 4. Weekly average number of hospital admissions, observed (O) and expected (E), for New York City (NYC) and the rest of New York (NY) State.If the model assumption is correct, (O − E)/ �� E √ would approximately follow the standard normal distribution.Both NYC and the rest of NY State have this value within ±2, so the variation can be considered due to randomness.(a) O and E for NYC facilities in 2020, weekly average.(b) (O − E)/ �� E √ for NYC facilities in 2020.(c) O and E for facilities in the rest of NY state in 2020.(d) (O − E)/ �� E √ for facilities in the rest of NY.

Table 1 .
Accuracy of β estimation of the proposed method in Scenario 3 (Independent) and 4 (Correlated) Note.Results were from 1,000 repetitions.ESD = empirical standard deviation; MSE = mean square root error; ASE = the naive estimate of the asymptotic standard error; CP = coverage probability; sASE = sandwich estimator of the asymptotic standard error; sCP = coverage probability based on the sandwich estimator.

Table 2 .
Results of α F estimation using the proposed method in Scenario 5 (Independent) and 6 (Correlated) Note. Results were from 1,000 repetitions.ESD = empirical standard deviation; ASE = the naive estimate of the asymptotic standard error; CP = coverage probability; sASE = sandwich estimator of the asymptotic standard error; sCP = coverage probability based on the sandwich estimator.

Table 3 .
Model fitting results for the main covariates The facility effects are important for surveillance.As an illustration of different facility effects, we provide a histogram and a funnel plot in Figure

Table 3 .
Continued Other coefficients are listed in Table S5 of the online supplementary material.a b c means RR (or HR, the hazard ratio) is b with (a, c) as the 95% CI.BMI = body mass index; ESRD = end-stage renal disease.

Table 4 .
Estimated relative rate (RR) and 95% confidence interval (CI) for patients with and without a previously observed COVID-19 diagnosis, assuming COVID-19's interaction effect is the same in all four stages