*2.2. Data Analysis*

Data from 2010 to 2019 (baseline years) and population forecasts up to the year 2040 were used to project the annual incidence of shoulder arthroplasty in Germany. A linear (Poisson, "classic approach") regression analysis was performed to estimate the expected incidence with the calendar year, sex, and patient age as covariates. The incidence was calculated by dividing the estimated number of arthroplasties for the national total and for each age subgroup by the corresponding official population forecast. An offset variable for the size of the population was chosen to ensure that the procedural number did not exceed the total population number. To overcome overdispersion problems that could result in an underestimation of the variance, we used a robust sandwich covariance matrix estimator for variance calculation. To minimize the error of variance underestimation of the estimated parameter because of overdispersion, we applied a quasi-Poisson regression to our data in accordance with the theory of quasi-likelihood.

As regressions based on logarithm, or an exponent, like Poisson, will only fit optimally when that is the exact nature of the true relationship, they might not be economically plausible, as in principle the projected counts can rise to infinity. In contrast, it seems quite reasonable to imagine that there might rather be an asymptotic or curvilinear relationship. To overcome this issue, an alternative estimate of the TSA projections was also determined by fitting the incidence rates (counts per 100,000 persons and year) with a negative binomial regression model using a monotone B-spline approach ("new approach") for modeling time effects and accounting for respective gender and age groups. Splines are used in statistics in order to mathematically reproduce flexible shapes. Knots are placed at several places within the data range, to identify the points where adjacent functional pieces join each other. The advantage of using splines for yearly data compared to the traditional approach is the more accurate curve estimation for the nonlinear trend changes and the simple way of modeling interactions between the time variables.

To compare the prediction accuracy of each model, the dataset was then split into training (years 2010–2017) and testing subsets (years 2018–2019). Both models were analyzed regarding common forecast accuracy measurement instruments (mean squared error (MSE), root mean squared error (RMSE), the mean absolute percentage error (MAPE) and Theil's U inequality coefficients, of which the first (U1) is a measure of forecast accuracy and the second (U2) is a measure of forecast quality [21]), which showed lower and thus more accurate values for the negative binominal approach using monotone B-splines (Table 1).

**Table 1.** Two-year forecast accuracy using an out-of-sample training-test validation set (minor numbers indicating greater forecast quality).


MSE: mean squared error; RMSE: root mean squared error; MAPE: mean absolute percentage error; U1/2: Thiel's U inequality coefficients.

Because of the anonymization of the diagnosis-related group DRG data, arthroplasty patients who underwent a revision (replacement or explanation) could not be individually followed and therefore, actual revision rates could not be calculated. Instead, we estimated the revision burden (RB) by dividing the number of revisions in the form of replacements or explanations by the number of all primary and revision shoulder arthroplasties (HA and TSA), as described previously [22]. Future projections for revision arthroplasty were also calculated, as mentioned above. For better comparison, a "constant rate" approach was used, based on the average revision burden during the baseline years and the projections of shoulder replacements in 2040.

All statistical analyses were performed using R Version 3.4.0 (R Development Core Team, The R Foundation for Statistical Computing, Vienna, Austria).
