Introduction

Percutaneous vertebroplasty (VP) and balloon kyphoplasty (BKP) are elective procedures used for the treatment of osteoporotic vertebral fractures. VP was first described in 1987 for the treatment of vertebral angioma [1]. Later, BKP, a modified VP technique that uses a balloon to erect the fractured vertebral body prior to the injection of bone cement, was introduced in 1998 [2]. After small-sized, uncontrolled trials showed substantial pain relief and few adverse events [3,4,5], the procedures rapidly gained ground in clinical practice for the treatment of osteoporotic vertebral fractures.

The evidence and recommendations for use of these procedures have evolved over time. Two high-quality randomized controlled trials (RCT) in 2009 showed no beneficial effect of VP compared to a sham procedure [6, 7]. Thereafter, the American Academy of Orthopaedic Surgeons (AAOS) recommended against the use of VP in the treatment of symptomatic osteoporotic spinal compression fractures [8], and the procedure rates declined by 62% in the USA until 2015 [9, 10]. Other professional societies, however, still recommended VP in patients with painful osteoporotic vertebral fractures refractory to medical therapy [11,12,13,14]. Two subsequent RCTs in patients with vertebral fractures of less than 6 weeks duration demonstrated more effective pain reduction after VP compared to a sham procedure [15, 16]. However, the VERTOS IV trial in 2018 showed no benefit of VP over a sham intervention [17], and a large meta-analysis in 2018 found no clinically relevant benefit of VP over sham procedures, with no difference based on duration of pain [18].

For BKP, recommendations were based on RCTs that showed greater pain reduction for BKP compared to conservative treatment [19,20,21,22], and improved correction of fracture-related kyphosis and a lower rate of cement leakage in BKP compared to VP [23, 24]. Despite the lack of sham procedure controlled trials, US and European guidelines both recommend the use of BKP for highly selected patients [8, 25]. As a consequence, BKP rates gradually increased from 2010 to 2015 [9, 26, 27].

It is ultimately unclear who will benefit from VP or BKP and the use of both procedures remain controversial [18]. While VP or BKP may lead to faster pain reduction and shorter hospital stay in selected patients, potential adverse events need to be considered. Leakage of bone cement is frequent (69% of VPs) and may result in myelopathy with need of surgical decompression (0.3% of patients with VP) [28]. Other rare but potentially life-threatening complications include cardiac or pulmonary cement embolism (3.9% and 3.5%, respectively [29]) as well as surgical site infections (0.4% [30]). There is an ongoing controversy on whether VP or BKP increases the risk of fractures of adjacent vertebra [31, 32].

In controversial procedures, such as VP/BKP, geographical variation may signal differing physician opinions more than differing medical need or patient preferences [33,34,35]. In Switzerland, guidelines recommend against the routine use of VP/BKP [11], and very wide variations in procedure rates across hospital service areas (HSAs) were observed between 2012 and 2013. [36]. Little is known about other potential determinants of the observed variation and whether these variations remain today. The aim of this study was to assess trends and regional variation in VP/BKP between 2013 and 2018 and examine potential determinants of the observed variation.

Methods

Data sources

We conducted a population-based, small area analysis based on routinely collected patient discharge data from all Swiss acute care hospitals for calendar years 2013 to 2018. The methods have been described previously in detail [36, 37]. Briefly, all Swiss hospitals are legally obligated to provide an anonymized, standardized dataset for each discharge to the Swiss Federal Statistical Office (SFSO). These data are centrally stored in the Swiss Hospital Discharge Masterfile. Recorded variables include age, sex, nationality, insurance status, type of admission (emergency vs. elective), up to 100 procedure codes based on the Swiss Classification of Operations (CHOP, an adaptation of the U.S. ICD-9-CM volume 3 procedure classification), and diagnostic codes based on the International Classification of Diseases, 10th revision, German Modification (ICD-10-GM) [38]. Furthermore, the patient’s area of residence and the location of the hospital within one of 705 Swiss administrative (MedStat) regions based on aggregated ZIP codes are recorded [39]. The SFSO reviews data integrity and completeness by means of a specifically designed software [40]. Since 2012, the Swiss Hospital Discharge Masterfile covers 100% of discharges, and data completeness for CHOP codes used in this analysis was high [41].

We used Swiss National Cohort data [42] from 2014 to determine the language region (Swiss German, Romance [Italian or French]), and data from the SFSO to determine the population density for each MedStat region. We used the average Swiss Socioeconomic Position (SSEP, version 2) as a measure of socioeconomic status. The SSEP version 2 was derived using 2012–2015 Swiss structural surveys data to rank Swiss neighborhoods based on four domains (median rent/m2, proportion of households with a person with no/low education, proportion with a person in manual/unskilled occupation, and mean crowding, all on the neighborhood level). The SSEP varies between zero (lowest) and 100 (highest) and correlates well with mortality [43]. We further obtained the density of orthopedic surgeons and radiologists per MedStat region for calendar year 2014 from the registry of physicians affiliated with the Swiss Medical Association (FMH) (orthopedic surgeons do the majority of BKP and radiologists do the majority of VP procedures in Switzerland). According to the Swiss Human Research Act, this study was exempted from ethics committee approval, because it was based on anonymized administrative data only.

Derivation of VP/BKP-specific hospital service areas

Although Swiss hospital care is primarily organized based on 26 geographic regions (the cantons), patients may utilize hospital services outside their canton of residence. Hence, the cantons are not suitable units to compare regional variation of VP/BKP utilization, because procedure rates may be skewed. We therefore generated reproducible HSAs by a fully automated method using all patients’ discharge data from calendar years 2013–2016, which has been previously described in detail [44, 45]. Briefly, in a first step, we identified 4,105,885 discharges of patients aged ≥18 years living in Switzerland from 155 Swiss acute care hospitals for calendar years 2013–2016. Across the 705 Swiss administrative (MedStat) regions, we identified regions that contain one or several hospitals as potential HSAs. Starting from these potential HSAs, in a centrifugal stepwise approach, we identified the geographically neighboring MedStat regions and merged them with the HSA if the majority of its residents were discharged from hospitals located in the specific HSA (plurality rule) [46]. HSAs with <50% of the patients being treated within the same HSA or <10 discharges overall were merged with the neighboring HSA which received most discharges until >50% and ≥10 discharges occurred within each HSA. This process yielded 63 general HSAs.

In a second step, we identified patient discharges with specific CHOP codes for VP (codes 81.65.x) and BKP (codes 81.66.x) from all Swiss acute care hospitals for calendar years 2013–2018 using the Swiss Hospital Discharge Masterfile. As VP and BKP are not performed in every hospital, we analyzed patient flows for VP/BKP. Using the procedure described above, HSAs were further aggregated into 31 VP/BKP-specific HSAs. We then drew visual maps of the 31 final HSAs using Geographical Information System (GIS)–compatible vector files.

Study population

We identified 18,522 discharges with specific CHOP codes for VP or BKP during calendar years 2013–2018 (Supplementary Fig 2). After delineating the HSAs, we excluded all discharges related to emergencies (n = 7169), and discharges with concomitant ICD-10-GM codes for accidents, intentional self-harm, or assault (ICD-10 GM code groups X, Y, Z; n = 2125). Because pathologic fractures due to vertebral metastases or multiple myeloma may be appropriate indications for VP/BKP [47, 48], we excluded 1373 discharges with ICD-10-GM codes for metastatic cancer or multiple myeloma (ICD-10 codes C77.x–C80.x, C90.0-, M49.5- and M82.0), leaving a final study population of 7855 patient discharges with elective VP or BKP. Discharges with both procedures were assigned to the BKP group.

Measures of variation

We calculated unadjusted and age- and sex-standardized VP/BKP procedure rates per 100,000 persons for each HSA using procedure counts and 2013–2018 census data for the adult Swiss population [49]. We used direct standardization with age bands of 18–49, 50–54, 55–59, 60–64, (...), 75–79, and ≥80 years. To examine the variation in procedure rates across the procedure-specific HSAs, we calculated three different measures of variation: We determined the extremal quotient (EQ), which is the highest, divided by the lowest procedure rate. While this is an intuitive measure of variation, it is prone to distortion by extreme values [34]. We further calculated the coefficient of variation (CV), which is the ratio of the standard deviation of the procedure rates to the mean rate, and the systematic component of variation (SCV) [34, 50, 51]. Although less intuitive, the SCV represents the non-random part of the variation in procedure rates while reducing the effect of extreme values [34, 50, 52]. It has been suggested that SCVs >5 are largely due to differences in practice styles or medical discretion and an SCV of >10 is considered indicative of a very high variation [34, 52].

Determinants of variation

Differences in illness incidences and socioeconomic factors are possible causes of variation [34]. Therefore, we explored the following domains, which potentially influence procedure rates: demographics (patients’ age), cultural and socioeconomic factors (language region, population density, socioeconomic status, insurance class, citizenship), population health (burden of disease), and supply factors (physician density). The language spoken by the majority of people living in a HSA was used to classify each HSA as either Swiss German or Romance (French/Italian) language region as a proxy for culture. We used population density as a proxy for the type of neighborhood a resident lives in (e.g., countryside areas have low density, cities generally have higher density). The socioeconomic status of each HSA was measured using the mean value of the SSEP of persons within an HSA. The percentage of discharges with (semi-)private insurance status and Swiss citizenship was used as an additional measure of the socioeconomic status of each HSA. As a proxy for the population burden of disease, we calculated age-standardized incidence rates of hip fractures (ICD 10 codes S720 - 22), colon (ICD 10 codes C18/19 and CHOP codes 457- 58 or 446) or lung cancer (ICD 10 codes C34 and CHOP codes 323 -26 or 329) treated surgically, acute myocardial infarctions (ICD 10 codes I21), or strokes (ICD 10 codes I63/64) for each HSA, as differences in these disease rates are likely to reflect true regional differences in burden of disease rather than differences in coding intensity or supply factors [53, 54]. The density of orthopedic surgeons and radiologists was used as a supply measure.

Statistical analyses

We used incrementally adjusted multilevel Poisson regression with a log link to model the procedure rates in each HSA using age bands of 18–49, 50–54, 55–59, 60–64, (...), 75–79, and ≥80 years. HSA was included as a random intercept in all models. In a stepwise approach, model 1 adjusted only for the calendar year of the procedure, model 2 was additionally adjusted for demographics, model 3 added socioeconomic and cultural factors, model 4 regional health, and model 5 supply factors. Variables included in the model were chosen a priori as we expected them to influence the rates. We used the models in three ways: (1) to assess to what extent VP/BKP rates in Switzerland can be explained by explanatory factors, (2) to assess the variance explained by the domains defined previously, and (3) to calculate intervention rates per 100,000 persons per HSA. For the first, we expressed the effect of explanatory factors on VP/BKP rates as incidence rate ratios (IRRs), defined as the VP/BKP rate in the defined category (e.g., Romance language region) relative to the estimated VP/BKP rate in the reference category (e.g., Swiss German language region). For the second use, we determined the percentage reduction in procedure variation across the 31 HSAs by examining the variance of the random intercept relative to model 1. We considered the residual, unexplained variation of the fully adjusted model a proxy for unwarranted variation. For the third use, we used the models to predict rates in each HSA. Sets of models were created for overall rates as well as for VP and BKP separately. We further assessed remaining variation in procedure rates across HSAs after full adjustment (model 5) using funnel plots. We plotted procedure rates against population size for each HSA. The mean procedure rates and the control limits of 2 and 3 standard deviations above and below the mean (95% and 99.8% confidence intervals), respectively, were calculated for all possible values for population size and used to create the funnel plot based on exact Poisson confidence intervals [55]. Statistical modeling was performed using Stata version 15.1 (StataCorp, College Station, TX, USA). HSAs were delineated and maps drawn using the R statistical software version 3.6.1 [56].

Results

Characteristics of VP/BKP-specific HSAs and the study population

For the 31 VP/BKP-specific HSAs delineated, the median population size per HSA was 785,112 persons (interquartile range [IQR] 355,009–2,456,318), with a median population density of 1580 persons per km2 (IQR 369–2712), a mean SSEP of 56 points (standard deviation [SD] 5.4), and a median density of orthopedic surgeons and radiologists of 21 per 10,000 persons (IQR 13–27). Overall, 23 HSAs were located in the Swiss German and 8 in Romance (4 French and 4 Italian) language regions.

Of the 7855 patients discharged with VP/BKP, 6110 (78%) underwent VP and 2413 (31%) BKP procedures (668 (9%) had both). The majority of patients were female (74%), ≥70 years old (72%), and Swiss citizens (92%, Table 1).

Table 1 Characteristics of patients receiving vertebroplasty (VP) and balloon kyphoplasty (BKP) in Switzerland 2013–2018
Full size table

Variation in procedure rates across HSAs

Across 2013–2018, the mean age/sex-standardized VP/BKP rate was 18 (range 2–35) per 100,000 persons, ranging across years from 16 per 100,000 persons per year in 2014 to 20 per 100,000 persons per year in 2018. After full adjustment for demographics, cultural, socioeconomic, health, and supply factors, and combining all years, the predicted VP/BKP rates varied between 4 and 40 per 100,000 persons per year across HSAs (Supplementary Fig. 1). For individual procedures, the mean age/sex-standardized VP rate was 12 (range 0–30 across HSAs) per 100,000 persons and the mean age/sex-standardized BKP rate was 6 (range 1–18 across HSAs) per 100,000 persons.

The overall EQ was 17.4, the CV 0.5, and the SCV 17.7, indicating very high variation across HSAs (Table 2). Both VP (EQ 73.1, CV 0.7, SCV 38.6) and BKP (EQ 28.9, CV 0.7, SCV 45.8) also showed very high regional variation. Between 2013 and 2018, the variation across HSAs decreased for both VP and BKP (by all measures, seen in Table 2), while the median procedure rates increased from 7 to 11 per 100,000 persons for VP and from 5 to 6 per 100,000 persons for BKP (Fig. 1). Despite the decrease in variation, it remained high in 2018 for both VP and BKP (e.g., SCV >6 for combined procedures and SCV >20 for individual procedures (Table 2).

Table 2 Measures of variation in age/sex-standardized vertebroplasty and balloon kyphoplasty procedure rates across Swiss HSAs through the years 2013–2018
Full size table
Fig. 1
figure 1

Age- and sex-standardized HSA vertebroplasty (A) and balloon kyphoplasty (B) rates (per 100,000 persons per year) through time. Each violin illustrates the procedure rates in a given year. Wider sections of the violin represent higher amounts of HSAs with a given procedure rate. The white dots inside the violins indicate median procedure rates, the black bars the interquartile ranges, and the black lines the lower/upper adjacent values. HSAs with fewer than 10 interventions in a given year were removed (i.e., the number of HSAs in each year changes. Consult Table 2 for the exact numbers)

Full size image

Determinants of variation in procedure rates

Between 2013 and 2018, overall VP/BKP rates increased by 7% per year (IRR 1.07, 95% confidence interval [CI] 1.05–1.10) in the fully adjusted model (Table 3). VP/BKP rates were highest in patients aged 75–79 years (IRR 23.37, CI 19.96–27.36 compared to patients aged 50–54 years). The rate in females was more than twice the rate in male patients (IRR 2.25, CI 2.14–2.36). Compared to the year-adjusted model (model 1), additional adjustment for demographic, cultural, and socioeconomic factors (model 3) explained 44.5% of the variation in procedure rates. After additional adjustment for health and supply factors (model 5), 51.6% of the variation in procedure rates was explained (i.e., 48.4% of the variation remained unexplained) (Table 3).

Table 3 Incidence rate ratios (IRR) of VP/BKP and remaining variance after incremental adjustment for determinants of variation across 31 Swiss HSAs 2013–2018
Full size table

The variance remaining in fully adjusted models was larger for the individual procedures (Supplementary Table 1). For VP and BKP, the remaining variance was 82.4% and 80.9%, respectively. The relationships between variables of interest and procedure rates were qualitatively similar for the individual procedures, with year, age, female sex associated with higher procedure rates (Supplementary Table 1)

When plotting the procedure rates against population size, several HSAs had procedure rates above or below the outer control limits of 3 standard deviations, indicating unusually high (HSAs 4, 6, and 8, around Bern) or low procedure rates (HSA 1 (Geneva) and 7 (Basel)) (Supplementary Fig. 3). HSAs with high procedure rates all lay in the greater Bern area, and tended to have a lower population and specialty-physician density compared to the HSAs with low procedure rates (Supplementary Table 2).

Discussion

In this analysis of all hospital discharges in Switzerland from 2013 to 2018, we found substantial regional variation in VP/BKP procedure rates, although the variation decreased over time. During the observed time period, while regional variation decreased, the procedure rates increased by 57% for VP and by 20% for BKP, respectively. We identified age, sex, and health and supply factors as statistically significant determinants of variation in VP/BKP rates. However, a large amount of variation remained unexplained by our models.

The variation in VP/BKP rates across HSAs steadily declined over time, but remained high by 2018. While previous studies also demonstrated substantial regional variation in VP/BKP rates across Switzerland [36] and across the USA [57-59], to our knowledge, our report is the first to specifically evaluate the time trend of the variation. We can only speculate on explanations for the observed trend. Between 2010 and 2016, many un-blinded trials that compared VP or BKP to conservative treatment demonstrated a statistically significant faster and/or greater pain reduction following an intervention [19, 22, 60-63]. Furthermore, the VAPOUR and VOPE trials showed better outcomes after VP in patients with acute osteoporotic vertebral fractures in 2016 [15, 16]. The results of these studies, together with the availability of an increasing number of international guidelines [12-14, 25, 64], might have convinced some physicians to use VP/BKP in patients with painful vertebral fractures refractory to medical therapy, even though the topic remained controversial. A large meta-analysis eventually showed no clinically relevant benefits of VP in 2018 [18] and, in early 2019, the American Society for Bone and Mineral Research (ASBMR) recommended against routine use of vertebral augmentation for osteoporotic vertebral fractures [65].

Large amounts of variation in procedure rates across Switzerland is not an isolated phenomenon for VP/BKP and has also been reported for other orthopedic interventions (e.g., SCV >33 for shoulder arthroscopy [66]) or hysterectomy (e.g., SCV >16 for abdominal hysterectomy [37]). The similar levels of regional variation may signal unequal access to surgery and widely differing physician’s opinions regarding surgical interventions across the country. It has previously been demonstrated that variation across geographically close small areas is intimately related to individual physicians’ practice style rather than to patient need or preferences [34].

While the regional variation decreased, we found a steady increase in both VP (+57%) and BKP rates (+20%) between 2013 and 2018. This increase is despite an unchanged recommendation from the Swiss medical board since 2011 against use of VK/BKP beyond trials. In comparison, in the USA, there was a substantial decrease in both VP and (less distinct) BKP rates in 2009 [58] after the publication of the first negative sham procedure controlled trials [6, 7], and the VP rates continued to decline between 2010 and 2018 while the BKP rates rose again during the mentioned time period [10]. In Germany and South Korea, VP/BKP rates also increased between 2005–2014 and 2012–2016, respectively [67, 68].

Several measured variables were associated with rates of VP and BKP in our analysis. Higher age and female sex are well-known risk factors for osteoporosis and the association with higher VP/BKP rates in this population group is evident. Less intuitively, we found an association between a higher density of specialty physicians and lower procedure rates. Previous studies demonstrated that higher rates of surgical interventions were not explained by differences in the number of surgeons within a specific area, but rather by enthusiasm for that procedure among a small number of high-volume surgeons in that region [69, 70]. Due to our study design, we were not able to assess surgeons’ enthusiasm for VP/BKP. Furthermore, we used the density of all registered radiologists and orthopedic surgeons (not only physicians performing VP/BKP) as a proxy for supply factors. Further studies to investigate physicians’ enthusiasm for VP/BKP might offer valuable clues to explaining regional variation in the procedure rates.

In HSAs with a lower burden of disease, we observed higher VP/BKP rates. One explanation might be that healthier patients usually are physically more active and thus are more likely to demand a surgical intervention for treatment of a vertebral fracture, hoping quickly to regain their active lifestyle. In previous studies, varying associations between comorbidities and VP/BKP rates were found [71-73]. It has to be noted that our study assessed burden of disease at a regional level, not on patient level.

In our analysis, nearly half of the regional variation in VP/BKP rates remained unexplained after adjustment for potential explanatory factors, highlighting clinical uncertainty about their use. VP and BKP remain controversial and it is not possible to determine the appropriate rate of such procedures, but the high regional variation indicates unequal access to the procedures across Switzerland. In practice, this may suggest overuse and unnecessarily putting patients at risk for adverse events in some regions, while patients in other regions might be left at risk for prolonged or more severe pain and disability. Better understanding is needed of how and why physicians and patients make these decisions about procedure use.

Our study has several limitations. We used administrative discharge data, and there might be errors due to inaccurate coding of main or comorbid diagnoses and procedures. Furthermore, data about the indications for procedures or about severity and age of the vertebral fractures were not available. Neither could we explore other potential determinants of regional variation, such as prevalence of osteoporosis, differences in physicians’ and patients’ preferences, or training of local physicians. Behavior/procedure rates of individual service providers, be they hospitals or physicians, could not be analyzed due to data protection issues, but will be important to assess in future work. We found associations between VP/BKP rates and several determinants on a regional level, but cannot infer causality. Furthermore, for the determinants of variation of interest, we used data of single calendar year (e.g., density of specialty physicians in 2014) and could not account for potential annual variation in these variables. Such covariates are also aggregated from finer scale measurements, which averages away some variation. In addition, our HSAs are of different sizes, and the regions they contain are not of uniform size (e.g., a relatively small HSA may be composed of a few large MedStat regions while a large HSA may be composed of many small and large MedStat regions). This effectively results in different sub-regions receiving different weights when considering the effects of the covariates if there are particularly high or low values in some sub-regions. As some MedStat regions (the smallest geographical unit of our analysis) contain several hospitals, we could not assess the procedure rates of and variation between individual hospitals in individual regions. Finally, because of lacking national data on outpatient procedures, we examined inpatient VP and BKP only. However, according to a major Swiss health insurer, the proportion of outpatient VP and BKP is very small (<10%).

In conclusion, we found a declining, but still high, level of variation in VP/BKP rates across Swiss HSAs between 2013 and 2018. Procedure rates during this time increased steadily. A substantial part of the regional variation remained unexplained by differences in demographic, cultural, socioeconomic, health, and supply factors.