The health effects of demand-side cost-sharing in European health insurance

We read in the data. We drop rows (combination of calendar year, NUTS 2 region, age and gender) where mortality measures (deaths, population or lagged_mortality) or oop (HF3_PC_CHE) and fraction of population that postponed treatment because it was too expensive (TOOEXP) are missing.

We select ages between age_min and age_max and calendar years between first_year and last_year. We drop rows where the number of deaths during a year exceeds the population at January 1st in that age-gender-year-region category.

age_min = 35
age_max = 85
age_range = np.arange(age_max-age_min+1)[:,np.newaxis]
plot_age = np.arange(age_min,age_max+1)
first_year = 2009
last_year = 2019

df = pd.read_csv('./data/data_deaths_by_age_nuts_2.csv')

df.rename(columns={'at risk of poverty':'poverty',\
                   'percentage_material_deprivation':'deprivation',\
                   'UNMET':'unmet'},inplace=True)

df.dropna(subset=['deaths','population', 'TOOEXP',\
                  'HF3_PC_CHE','lagged_mortality'],
                    axis=0, how ='any',inplace=True)
df = df[(df.population > df.deaths) & (df.age >= age_min) & \
        (df.age <= age_max) & (df.year <= last_year) &\
        (df.year >= first_year)]
df['mortality'] = df.deaths/df.population*100 \
    # mortality as a percentage

# lagged mortality as fraction of mean lagged mortality
# per age/gender group
df['lagged_mortality_s'] = (df['lagged_mortality'])/\
    df.groupby(['age','sex'])['lagged_mortality'].\
    transform('mean')

#len(df)

As we use a different python kernel here from the one used in the rest of the paper (due to conflicting supporting packages at the time of the analysis), we need to import some libraries and the data again. We use geopandas to plot the map of the NUTS 2 regions where the color indicates mortality per 100k population.

import numpy as np
import geopandas as gpd
import matplotlib.pyplot as plt
plt.style.use('Solarize_Light2')
import pandas as pd
# import altair as alt

# read the NUTS shapefile and extract
# the polygons for a individual countries
nuts=gpd.read_file('./SHP/NUTS_RG_60M_2021_4326_LEVL_2.shp')

age_min = 35
age_max = 85
plot_age = np.arange(age_min,age_max+1)
first_year = 2009
last_year = 2019
df = pd.read_csv('./data/data_deaths_by_age_nuts_2.csv')
df.dropna(subset=['deaths','population', 'TOOEXP',\
                  'HF3_PC_CHE','lagged_mortality'],
                    axis=0, how ='any',inplace=True)
df = df[(df.population > df.deaths) & (df.age >= age_min) & \
        (df.age <= age_max) & (df.year <= last_year) &\
        (df.year >= first_year)]
df['mortality'] = df.deaths/df.population*100000
df = df[(df.year==2018) & (df.age==40) & (df.sex=='F')]


nuts = nuts.to_crs(epsg=3035)
nuts['centroids'] = nuts.centroid
nuts = nuts.merge(df, how='inner',\
                  left_on = 'NUTS_ID', right_on = 'nuts2')

nuts[nuts.sex=='F'].plot(column='mortality',
                         legend=True,
                         figsize=(16,16),
                         # vmin = 71,
                         vmax = 0.002*100000,
                         missing_kwds={'color': 'lightgrey'},
                         legend_kwds={'label': "Mortality rate",
                                      'orientation': "vertical"})
# adjust plot domain to focus on EU region
plt.xlim(0.25e7, 0.6e7)
plt.ylim(1.3e6, 5.5e6)
plt.xticks([],[])
plt.yticks([],[])
plt.title(\
  'Mortality 40 year old females in 2018 (per 100,000 population)');
# plt.tight_layout()
# plt.legend('right');

Consider a population (of a certain age and gender in a particular year) in an EU region where a fraction \(\alpha \in \langle 0,1 \rangle\) has low income \(y^l\) and fraction \(1-\alpha\) high income \(y^h\). We think of \(\alpha\) as the poverty rate. Let \(\pi^j\) denote the probability that someone with income \(y^j, j=l,h\) falls ill. As is well known, low income people tend to have a lower health status (Cutler, Lleras-Muney, and Vogl 2011). We capture this by assuming \(\pi^l > \pi^h\). People on low income may have a less healthy diet, exercise less etc. due to either the cost or knowledge of healthy lifestyle choices. This makes it more likely that they fall ill. Thus we separate the direct health effect of income (\(\pi^l > \pi^h\)) from treatment decisions made by people on low income.

Generally speaking, oop payments tend to take two forms that we want to capture: a coinsurance rate, which we denote \(\xi \in [0,1]\), and a maximum expenditure, which we denote \(D\) (for deductible). Some systems have a combination of the two.

Conditional on falling ill, there is a probability \(\zeta_i \in [0,1]\) that the patient is advised to get treatment \(i\) at cost \(x_i\) for \(i\) in the set of “illnesses” \(I\) with \(\sum_{i \in I} \zeta_i =1\).² We define \(I_{\xi}\) as the subset of \(I\) where \(\xi x_i < D\) and \(oop_i = \xi x_i\) and set \(I_D\) where \(\xi x_i \geq D\) and \(oop_i = D\). To keep things simple, we assume that \(\zeta_i\) is exogenous to the patient. We model the treatment decision on the extensive margin only: an agent accepts or rejects the treatment proposed by a physician.³ A pure coinsurance system has \(\xi < 1\) and \(I_{\xi}=I\). A pure deductible system \(\xi=1\) and \(I_D\) non-empty. A combination of the two has \(\xi<1\) and there is a maximum on the oop payment. Health insurance systems in Europe tend to have such maximum oop expenditure.⁴ An increase in either \(\xi\) or \(D\) is interpreted as making health insurance less generous.

Whereas with individual level data one can determine whether an individual faces a positive treatment price at the margin (E.g. using the end-of-year price as in Keeler, Newhouse, and Phelps 1977; Ellis 1986), this is not possible with the aggregate data that we use here. Hence, we rely on an aggregate summary variable, denoted OOP, measured as oop payments over total healthcare expenditure. That is, the fraction of healthcare expenditure paid by patients oop. We interpret this variable as capturing the generosity of the health insurance system. To illustrate, if healthcare is free at point of service, OOP equals zero; if there is no health insurance at all, OOP equals 1. In a pure coinsurance system with rate \(\xi\) applying to all treatments, OOP equals \(\xi\). It is the cap on oop expenditure (like a deductible) that makes the relation between OOP and healthcare use non-linear. The challenge then is to capture changes in \(\xi\) and \(D\) although we do not directly observe these variables in the data. This is what the model sets out to do.

If an ill patient receives treatment, we denote her (expected) health \(\sigma\), while without treatment (expected) health equals \(\sigma_0\) with \(0 < \sigma_0 < \sigma < 1\).⁵ Health is normalized at value one for a patient who does not fall ill. The trade off between health and oop is captured by \(\sigma_0/\sigma <1\) and a simple assumption that utility is multiplicative in health and consumption. That is, consumption yields higher utility if you are healthier. To put it bluntly, if you are healthy and can travel, go skiing etc. consumption yields higher utility than when you are ill, lying in bed all day. We model the patient’s treatment decision as:

where utility \(u(.)\) is determined by how much money can be spent on other goods: income \(y^j\) minus oop in case of treatment and \(y^j\) if no treatment is chosen. The utility function \(u(.)\) is increasing and concave in consumption: \(u(.), u'(.) >0\) and \(u''(.) < 0\). Further, parameter \(\nu\) captures other factors than pure financial ones affecting a patient’s treatment choice.⁶ If the inequality holds, the patient accepts treatment \(i\).

In our data, we have a variable “unmet medical needs” based on a number of motivations: treatment is too far away to travel to, there is a long waiting list, the patient is scared to undergo treatment etc. To make our point, it is enough to assume that such factors affect utility in a multiplicative way. To illustrate, if the patient has to travel far for treatment, utility is reduced by multiplying it with a value of \(\nu < 1\). Agents differ in \(\nu\) and the cumulative distribution function of \(\nu\) is given by \(G(\nu)\), its density function by \(g(\nu)\). Other factors can include waiting time till treatment, belief that the condition will resolve itself without intervention, poor decision making e.g. with a focus on the short term thereby undervaluing the benefit of treatment.

The probability that a patient with income \(y^{j}\) accepts treatment \(i\) offered by a physician equals

that is, \(\nu\) is big enough that inequality \eqref{orgff195e9} holds. With probability \(G\left( \frac{\sigma_0}{\sigma} \frac{u(y^{j})}{u(y^{j}-oop_i)} \right)\) the patient decides to postpone or forgo treatment \(i\).

The probability that a patient postpones or skips a treatment because it is too expensive is given by

These are agents \(\nu\) that would have chosen treatment \(i\) if it were free (\(oop_{i}=0\) and \(u(y^j)/u(y^j-oop_i)=1\)) but who forgo treatment now that it costs \(oop_{i}>0\). The probability \(G(\sigma_{0}/\sigma)\) captures factors like waiting lists or the patient hoping that the health problems resolve themselves without treatment. That is, reasons for postponing treatment not related to oop payments.

In the proof of the lemma below, we show that the probability of accepting treatment, \(\delta_i^j\), is increasing in income \(y^j\) and decreasing in \(oop_{i}\), as one would expect.

Note that this model differs from a standard Rothschild and Stiglitz –R&S– health insurance model (Rothschild and Stiglitz 1976) in the following way. In an R&S model income plays no role and people with low health status have generous insurance coverage. Hence, they would not postpone valuable care. In our model, people with low health tend to have low income and may skip valuable treatment if the oop expense is high. This negatively affects their health.

In Appendix 9 we specify how (unobserved) health and treatment decisions translate into (observed) mortality for each age/gender class.

In our data, the variable Unmet varies with NUTS 2 region and year. In terms of our model, we define this variable with subscript \(2\) for region and \(t\) for calendar year as follows:

with treatment probability \(\delta^j_i\) for illness \(i \in I\) and income class \(j \in \{l,h\}\) varying with country \(c\) and year \(t\) because oop varies with countries over time. In words, for people on low [high] income –fraction $α_2t [1-α_2t]$– there is a probability \(\zeta_i \pi^l [\zeta_i \pi^h]\) of falling ill with disease \(i\) where they forgo treatment with probability \(1-\delta_{ict}^l [1-\delta_{ict}^h]\).

Further, in our data we have the variable OOP defined as oop payments as a percentage of healthcare expenditure. In terms of our model, we write this –ignoring subscripts– as

where \(\zeta_i (\alpha \pi^l \delta^l_i + (1-\alpha) \pi^h \delta^h_i)\) denotes the fraction of people accepting treatment \(i\). If people do not accept treatment, there is no oop and no expenditure. The numerator of OOP contains the oop payments \(oop_{i}\) and the denominator expenditures \(x_i\). If \(I_{\xi} = I\), it is clear that \(\text{OOP} = \xi\). Because \(I_D\) is non-empty (European countries have a maximum oop payment), the expression for OOP is actually non-trivial. We can also write OOP as the ratio of average oop per head and average healthcare expenditure per head:

In our data these variables vary by country \(c\) and year \(t\).

The following lemma summarizes the main results from the model and presents the equations for mortality \(m_{ag2t}\) varying with age, gender, NUTS 2 region, calendar year and the fraction of people forgoing treatment because it is too expensive, \(\text{TooExp}_{2t}\) that we estimate below. Further, our simulation results are presented in terms of the relative increase in deaths due to the increase in oop, \(dm_{ag2t}/m_{ag2t}\). In the lemma we use the following indices: age \(a\), gender \(g \in \{f,m\}\), country \(c\), NUTS 2 region \(2\), calendar year \(t\)

As derived in the appendix, mortality can be written as the multiplication of an age/gender effect with a factor depending on the situation in the NUTS 2 region. We think of the age/gender effect as biology that is the same across regions. This is modeled as a sigmoid of age and gender fixed effects, \(\beta_{ag}\), which makes sure the probability of death is between 0 and 1. We multiply this baseline probability with a multiplier capturing the other effects.

First, NUTS 2 region fixed effects, \(\mu_2\), which capture regional variation in the probability of falling ill. One can think of lifestyle habits that vary by region, external factors affecting health like clean air, road safety, travel distance to closest medical facilities that tend to be longer in rural areas etc.

Second, whether this age \(a\) cohort experienced a health shock in the previous period \(t-1\) when aged \(a-1\). If there was such a negative health shock that increased mortality, we expect that part of this shock spills over in the current period further increasing mortality. We measure the health shock as mortality for this group in the previous period compared to average mortality for this group (across regions and time).

Third, the region’s poverty level and the fraction of people with unmet medical needs in the region in year \(t\) affect mortality. As shown in the proof of the lemma, \(\beta_{poverty} = (\pi^l - \pi^h)(1-\sigma) > 0\): if there are no unmet needs, poverty still raises mortality as poor people fall ill more often (\(\pi^l-\pi^h >0\)) and treatment only partially recovers their health (\(1-\sigma>0\)). This is the health effect of low income. Further, \(\beta_{unmet}=(\sigma-\sigma_0)>0\): with unmet medical needs, patients end up with lower health than they would have gotten with treatment (\(\sigma_0 < \sigma\)). In words, we expect mortality (in a region in year \(t\)) to be higher (for a given age/gender category) if poverty is higher and there are more people with unmet medical needs.

If the sum of these three terms is negative, the multiplier is less than 1 and mortality for this age/gender/region/year combination is reduced compared to the baseline probability given by the sigmoid. If the sum of the terms is positive, mortality for this observation is higher than the baseline probability.

For the second equation in the lemma, we use a linear expansion of TooExp in terms of OOP. The appendix shows how we derive this relation using the policy variables \(\xi\) and \(D\) which affect OOP and TooExp simultaneously. It turns out that there is a direct effect of OOP on TooExp and an interaction effect with the fraction of people below the poverty line in a region. We show that \(b_{oop,c},b_{interaction,c} > 0\): a region that lies in a country with high OOP tends to have high unmet needs (as medical care is expensive) and especially so if the region features a high poverty rate. On the other hand, if OOP equals 0 (healthcare is free at point-of-service) one does not expect poverty to affect TooExp (it will still affect health and mortality through lifestyle choices).

The third equation shows how a 500 euro increase in oop affects mortality for an age/gender category in a NUTS 2 region in year \(t\). In our simulations we present this as the increase in mortality per 1000 deaths. The mortality effect of demand-side cost-sharing goes via unmet medical needs. People fall ill and cannot afford the treatments recommended to them by a physician. This reduces their health status and affects the probability of dying.

As shown below, in our data \(TooExp\) is smaller than 0.5 in all regions. Hence, we see that the mortality effect \(dm/m\) increases with the fraction of people that forgo treatment because it is too expensive. This suggests the following non-linearity in the effect of oop: if demand-side cost-sharing is initially low, \(TooExp\) is low and an increase in oop hardly affects mortality, but in countries where oop is already significant and \(TooExp\) is high a further increase in oop has serious mortality consequences. Finally, an increase in oop increases the fraction of people with unmet medical needs because treatment is too expensive and especially so in regions where poverty is high. The policy implications of this are clear: in a country where oop is already high and where poverty is a serious problem, the government should be careful increasing oop further.

Note that the increase in the number of deaths \(dm_{ag2t}/m_{ag2t}\) is independent of age. This is due to our formulation of mortality where we have a baseline mortality depending on age/gender only and a deviation from this baseline based on poverty and unmet medical needs in the region.

The model sets the stage for the empirical analysis in two ways: (i) it helps us specify functional forms, (ii) it helps us to avoid “causal salad” (McElreath 2020). Because the model is clear on the mechanisms that are covered, we can also identify potential mechanisms that are missing which can confound our estimates. This is illustrated with a Directional Acyclic Graph, DAG (Pearl 2009; Hernán and Robins 2023). The arrow points from the node that has a causal effect to the node that is affected.

The grey nodes capture the model equation where poverty and OOP (and their interaction which is not separately depicted in a DAG) affect the fraction of people that forgo treatment because it is too expensive. The blue nodes capture the mortality equation where TooExp affects (because it is part of) people reporting unmet medical needs, poverty affects health and previous period health (health -1) affects today’s health. Finally, health determines mortality.

The DAG clearly shows two causal paths from poverty to mortality. First, poverty is directly associated with lower health, say due to the financial stress involved of living on low income and due to less healthy lifestyle choices, e.g. because fresh fruit and vegetables tend to be expensive. Second, poverty makes it more likely that people forgo valuable treatments leading to unmet medical needs; especially when OOP is high.

Some variables causing differences between regions that are more or less constant over time –think of lifestyle (smoking habits), pollution in the region etc.– are not explicitly mentioned in the DAG but are captured by region fixed effects in the model.

Effects are potentially confounding if they differ between regions, vary over time (not captured by fixed effects) and simultaneously affect both unmet needs and mortality (in the mortality equation) or both TooExp and the interaction OOP \(\times\) poverty (in the TooExp equation). In the robustness analysis, we focus on two plausible mechanisms that can lead to these effects over time: shocks to healthcare resources and to health itself.

First, consider a shock that reduces government resources available to finance healthcare in the region or country. If this increases waiting lists due to reduced healthcare capacity, Unmet is likely to rise. At the same time, this can also reduce the quality of care, say because equipment maintenance is reduced, equipment is replaced less often or due to wage cuts high quality physicians leave and are replaced by lower quality staff. This reduction in care quality can affect health and mortality introducing a different mechanism from the one we focus on; that is, we have a causal path (indicated with dotted lines) via the red node healthcare quality that goes outside the model.

Note that a shock to government resources which raises poverty (due to reduced social assistance) and forces the government to raise out-of-pocket expenditure (say, by raising the deductible) is not a confounding effect. Indeed, the model captures that due to this shock the fraction of people that forgo treatment because it is too expensive goes up. Further, this increase in TooExp raises the fraction of people with unmet medical needs which will affect health and ultimately mortality. Unlike the healthcare quality example above, here all causal paths are within the model (going through gray and blue nodes).

Second, a health shock (causing a higher fraction of people with bad health) can increase both unmet medical needs (as more people fall ill, more people can have unmet needs) and increase mortality. This is partly captured by previous year mortality in the equation for \(m_{ag2t}\) in the lemma but in the current year this effect is depicted by the dotted arrow from (red) bad health to unmet needs and mortality (via health). This (current year) effect goes outside the model and can potentially confound the effect that we find. In the same vein, a health shock can increase the fraction of people that forgo medical treatment because it is too expensive and can raise poverty through reduced productivity. Again this is a potentially confounding effect outside of the model. Although our data do not include the Covid years, we cannot exclude the possibility of other shocks that tend to either reduce health and raise mortality or increase poverty and the fraction of people forgoing treatment because it is too expensive.

In the robustness section we introduce variables to control for these potentially confounding healthcare quality and health effects. Then we compare our baseline estimate of the OOP effect with the effect that follows from these equations with an extended variable set.

In this section we use the graphviz library to plot the DAG for our model.

g = gr.Digraph()
g.attr('node',color='lightgrey', style='filled')
g.edge("OOP", "TooExp")
g.edge("Poverty", "TooExp")


g.attr('node',color='lightblue2', style='filled')
g.edge("Unmet", "health")
g.edge("health -1", "health")
g.edge("health", "mortality")
g.edge("TooExp", "Unmet")
g.edge("Poverty", "health")



g.attr('node',color='pink', style='filled')
g.edge("bad health", "TooExp", style="dotted")
g.edge("bad health", "Unmet", style="dotted")
g.edge("bad health", "health", style="dotted")
g.edge("bad health", "Poverty", style="dotted")
g.edge("healthcare resources", "Unmet", style="dotted")
g.edge("healthcare resources", "OOP", style="dotted")
g.edge("healthcare resources", "healthcare quality", style="dotted")
g.edge("healthcare quality", "health", style="dotted")
g.render("./figures/DAG",format="png")

figures/DAG.png

The following python code generates Table 1 with summary statistics. The noweb notation <<code-preamble>> runs the code from Listing 2 into the python cell without using copy/paste.

<<code-preamble>>
<<code-data>>
headers = ['count','mean','std','min','median','max']
variables = df[['population','deaths','mortality',\
                'poverty',\
                'deprivation',\
                'TOOEXP','unmet',\
                'HF3_PC_CHE','HF2_PC_CHE',\
           'health expenditure per capita',\
           'infant mortality', 'bad_self_perceived_health']]\
           .describe().T[['count','mean','std','min','50%','max']]
variables = variables.round(1)
variables['count'] = variables['count'].astype(int)
variables.rename({'poverty':'poverty (%)',\
                  'mortality':'mortality (%)',\
                  'deprivation':\
                  'deprivation (%)', 'HF2_PC_CHE':'voluntary (%)',\
                  'HF3_PC_CHE':'out-of-pocket (%)',\
                  'TOOEXP':'too exp. (%)',\
                  'health expenditure per capita':\
                  'expend. per head',\
                  'unmet':'unmet (%)',\
                  'bad_self_perceived_health':'bad health (%)',\
                  'infant mortality':'infant mortality (\\textperthousand)'},inplace=True)
print(tabulate(variables,headers,tablefmt="orgtbl"))

The first equation estimates a binomial model with population size as the number of draws and deaths as the number of events. We do this for every combination of age, gender, NUTS 2 region and calendar year in our data. The probability of \(k \leq n\) deaths out of a population \(n\) is then given by

where \(m\) denotes mortality: the probability of death. The equation that we estimate for mortality \(m_{ag2t}\) is given in the lemma above. The coefficient we are especially interested in is \(\beta_{unmet}\). This is the coefficient through which an increase in unmet medical needs because of financial problems affects mortality.

Figure 4 illustrates that without the multiplicative specification for \(m_{ag2t}\) in the lemma, the coefficients for \(\beta_{unmet}, \beta_{poverty}\) would have to vary with age. Indeed, for the young mortality is low even in regions with high poverty or high unmet needs. Specifying coefficients that vary with age would considerably increase the number of parameters that we need to estimate.

The second equation captures how an increase in OOP affects the fraction of people in a region that postpone or skip treatment because it is too expensive. This fraction TooExp is based on (EU-SILC) survey data where we do not know the number of people interviewed. Hence, we cannot model this as a binomial distribution. In our estimation we want to ensure that TooExp is between 0 and 1. For this we assume that TooExp in the lemma above has a logit-normal distribution. That is, the log-odds of TooExp is normally distributed.

Details of the estimation can be found in Appendix 11.1 and the python code is in the online appendix.

We use Markov Chain Monte Carlo (MCMC), in particular the NUTS sampler to explore the posterior distributions of our parameters. For this sampler, we have the guarantee that the whole posterior distribution is captured as long as we have enough samples. Although this is an asymptotic result, we are confident that drawing four chains of 2000 samples (1000 samples of which are used for tuning) is enough to cover the posterior distribution. In the appendix we discuss a number of checks on this convergence.

It is not straightforward to put priors on the coefficients of the two equations in Lemma 1. To illustrate, how strong is the reaction of mortality to a \(0.1\) increase in the fraction of people reporting unmet medical needs? We are not aware of previous studies looking into this and have no a priori information on the strength of this effect. We use three principles when setting priors. First, we use regularizing priors (“seat belt priors”): priors close to zero with small standard deviations. Hence, a coefficient differs from zero only if there is clear evidence for this in the data. This reduces the risk of over-fitting. Second, we use a hierarchical model to determine the parameters of the prior distributions. Finally, if the theory suggests a parameter is positive, the prior distribution reflects this (e.g. using a half-normal instead of a normal distribution). Details on the priors can be found in the appendix.

As it is hard to judge how sensible a prior for one particular coefficient is, the online appendix to this section shows the prior predictive distributions. That is, the predictions for mortality and TooExp that the model generates without having seen the data. Comparing the prior predictive distributions with the observed distributions, we show that our priors do not exclude relevant possible outcomes nor do they put (much) weight on unlikely outcomes (say, mortality close to 0.9).

Finally, as shown in Table 1, we have fewer observations for poverty than for deprivation. But in the robustness analysis where we use the variable poverty we do not want to change the sample. Appendix 11.2 explains how Bayesian estimation deals with missing values without imputation or dropping observations.

Figure 5 gives an idea of the fit of the model in terms of predicting deaths per age/gender/region/year category and the fraction of people postponing treatment because it is too expensive.

The left panel shows observed number of deaths per category on the horizontal axis and the posterior predictive for this on the vertical axis. For each row in our data, we have observed number of deaths and a distribution of predictions of this number. In the figure, we show the average prediction of deaths across the posterior samples. The predictions are not perfect but do follow the 45-degree line closely.

The right panel shows the (log odds of the) fraction of people per region/year indicating they went without treatment (for a while) because it was too expensive. The difference between this panel compared to the one on the left is that this fraction does not vary by gender and age. Hence, we do not have a prediction for each “row in our data”. The right panel shows the observed and predicted fraction for TooExp per region/year. The dots indicate the average posterior prediction of this log-odds ratio. For small observed values of TooExp (log-odds below \(-5\) in the figure) there is a range of predicted values. Although this range seems wide in log-odds space, both the observed and predicted values are equal to zero. To illustrate, for practical purposes it does not matter if a probability equals 0.0001 (log-odds of \(-9\)) or 0.002 (log-odds of \(-6\)): both values are basically zero. Moreover, given our log-odds specification, the model cannot predict an exact zero probability.

A related observation is that in the data TooExp equals 0 for a number of region/year combinations. To handle this numerically, we use a lower bound for the log-odds. This corresponds to a probability of 0.0001 which is close enough to zero for our purposes. The right panel shows this bunching for a number of observations slightly below \(-9\). The bunching for other values of observed log-odds between \(-5\) and \(-7\) corresponds to regions reporting rounded fractions of \(0.001,0.002\) etc.

Compared to the observed number of deaths, the predictions for TooExp seem less accurate. This is to be expected as there are (a lot) fewer observations for this variable compared to mortality. But all in all the fit does not seem unreasonable as the points cluster around the 45-degree line.

Another way to check how well the model fits, is to see how well it captures the age profile of mortality. This we present in Figure 6. The left panel shows the age profile \(e^{\beta_{ag}}/(1+e^{\beta_{ag}})\). If the other terms in equation \eqref{orgc7f7f41} equal 0, \(e^{\beta_{ag}}/(1+e^{\beta_{ag}})\) gives the probability of death for age/gender category \(ag\). The right panel includes for every region and calendar year the correction on \(e^{\beta_{ag}}/(1+e^{\beta_{ag}})\) to yield mortality for that combination of age/gender/region/year. On average, the model captures the age profile perfectly.

The appendix presents two further checks of the model. Figure 11 shows the trace plots for the parameters of interest. The figures in the left panel show the posterior distribution of the parameters. The coefficients b_oop, b_interaction vary by country and hence we have different colors for the country specific distributions in these graphs. The beta parameters do not vary with country (or another index) and hence there is one color only. In the beta figures it is easy to see that there are four distributions per parameter. These correspond to the four chains that are sampled by the NUTS algorithm.

The right panels show the same samples but now ordered across the horizontal axis as they were drawn. We check these plots for the following three features. First, the plot should be stationary; that is, not trending upward or downward. This implies that the posterior mean of the coefficient is (more or less) constant as we sample. Second, there should be good mixing which translates in condensed zig-zagging. In other words, the algorithm manages to draw values across the whole domain of the posterior quickly one after the other. Finally, the four chains cover the same regions. All three features are satisfied for the coefficients in the right panel of the figure.

Another check on the convergence of the algorithm are the r-hat values in Table 5 in the appendix. This table summarizes the posterior distribution for the slopes that we are interested in. It provides the mean and standard deviation for each of these parameters, the 95% probability/credibility intervals and the number of effective samples for each parameter. As the number of these effective samples (ess_bulk column) is roughly above 500 for all and above 1000 for most parameters, this looks fine. The final column presents the values for r-hat for each parameter. Since these are all equal (close) to one, we can be confident that the NUTS algorithm converged for these parameters.

Table 5 in the appendix presents the posterior values for each of the parameters. Here we focus on the effect we are interested in: what is the increase in mortality due to an increase in oop? Lemma 1 shows the effect \(dm/m\) of a 500 euro increase in oop. Figure 2 reports the expression in this equation multiplied by 1000. That is, we report the increase in deaths due to the oop increase per 1000 deaths. Note that the 500 euro change in OOP enters multiplicatively. In other words, dividing the effect in Figure 2 by ten gives the effect of a 50 euro increase in OOP for each country. In this sense, the choice of 500 euro is a matter of presentation.⁸

As the expression for \(dm/m\) varies with country, year and NUTS 2 region, Figure 2 summarizes our main findings in the following way. For each country we focus on the region where deprivation is highest. This is the region where we expect the mortality effect of an oop increase to be highest as many people could have problems paying medical bills. Table 2 presents this region for each country in our data together with the value of deprivation, the fraction of people with unmet medical needs due to financial constraints and the country’s value for OOP. As the table illustrates, the fraction of people indicating that treatment was too expensive tends to be high when both deprivation and OOP are high.

Substituting these values from the table into the expression for \(dm/m\) we get the numbers in Figure 2. As mentioned, the blue bars give the average effect of the 500 euro increase in oop on mortality. As we have the posterior distributions for each of the parameters, we also have the posterior distribution for the mortality effects per country (taking the uncertainty for all parameters into account). The black horizontal lines present the 95% intervals around the mean effect.

The first observation is that for Bulgaria, Greece, Hungary and Romania the 95% probability interval is bounded away from zero. For these countries we can clearly see that an increase in oop negatively affects health and increases mortality.

Why are the effects smaller for the other countries? The effects are basically zero for the Scandinavian countries, Slovenia and Switzerland. As shown in Table 2, for these countries both deprivation and the fraction of people indicating unmet medical needs because treatment is too expensive are small. For the Scandinavian countries in the region with highest deprivation, TooExp is basically zero. It then follows from the equation for \(dm/m\) in the lemma that the effect on mortality is (close to) zero.

Another reason why the effects are small for some countries is that the underlying parameters b_oop, b_interaction are small for these countries. This can be seen in Table 5 in the appendix. If countries have policies to subsidize healthcare for poor families, the effect of country wide OOP on these families’ unmet medical needs is small as they actually pay a lower fraction (than the national average) of their treatments’ costs oop.

Summarizing, we can identify in our data the effect that an oop increase, raises the number of people with unmet medical needs due to financial constraints and hence increases mortality. This is especially the case in regions with high poverty and high initial OOP. Documenting this effect was the main objective of the paper.

A follow up question is: how big is this effect? In order to interpret the size of the oop effect, Table 3 presents the number of people dying from a particular cause per 1000 dead. If we would consider all causes and add them up, the sum of the second column in Table 3 would equal 1000. The table focuses on causes of death with an order of magnitude comparable to the effects in Figure 2. The table is based on EU wide data in 2017 for ages 35-85.

Note that the comparison of the numbers in the figure with the numbers in the table is just to get an idea of the order of magnitude. But –strictly speaking– the causes are not comparable. Nobody dies of an increase in oop the way people die from pneumonia. Due to an increase in oop, people may have gone without treatment which can then lead to death from, say, lung cancer. Hence, one should be careful in comparing the simulation results with the numbers in Table 3. But the table does provide some context in interpreting the size of the simulated effects.

The average mortality effect due to a 500 euro increase in oop in Romania is approximately 33 (per 1000 dead). This exceeds deaths due to each of the causes in the table. The average effects in Bulgaria and Greece are around 15 and 22 resp. which places them between deaths due to Alzheimer disease and diabetes. In Hungary the order of magnitude is comparable to deaths due to transport accidents.

However, these are effects aggregated at the regional level (of the regions with highest poverty levels). Suppose we are willing to assume that the incidence of the increase in mortality due to the 500 euro increase in oop falls mainly in the group of people who live in material deprivation. Table 2 shows the relative size of this group is around 30% for the relevant regions in Greece, Hungary and Romania. To get these effects at the region level, the effects among this specific group is an order of magnitude bigger (roughly speaking, multiply by 3).

Finally, there is also the following dynamic effect. As oop increases, 35 year olds postpone treatments thereby lowering their health status. Part of this reduced health leads to higher mortality among 35 year olds but some of these people survive this year. Next year, they start with lower than average health which can then raise mortality among 36 year olds. These dynamic feedback effects are captured by the parameter \(\gamma\) in Lemma 1. As shown in Table 5 in the appendix, the estimated value for \(\gamma\) is approximately 0.5 (coefficient beta_lagged_log_mortality). As effects accumulate across age and time, the effect for 85 year olds almost doubles (\(1+\gamma+...+\gamma^{50} \approx 2\)). To illustrate, the long run effect of a 500 euro increase in OOP leads to 66 deaths per 1000 dead for 85 year old Romanians in its poorest region.

One of the advantages of doing a Bayesian analysis is that we can easily show the uncertainty surrounding our estimated effects. This is illustrated in Figure 7 where we show for eight Romanian regions the probability that the mortality effect exceeds a certain value. Region RO22 tends to have the biggest effect (reported in Table 2 and Figure 2), while the effect per 1000 dead is smallest in NUTS 2 region RO42. We are pretty sure (probability close to 1) that in RO22 the effect is at least 15 per 1000 dead. While in RO42 this probability is less than 40%. In RO42 we are 80% sure that the effect exceeds 10 per 1000 dead. This is due to the fact that both deprivation and too expensive are substantially lower in RO42 compared to RO22.

Summarizing the discussion on the size of the effect, we find the following. In countries where poverty and OOP are high, a 500 euro (further) increase in oop leads to an increase in mortality (per 1000 dead) that is comparable to causes varying from Alzheimer disease to diabetes or breast cancer.

Here we provide the code that generates the results on model fit and size of the effects discussed in this section.

As discussed in section we focus on two mechanisms that may confound our estimated effect of oop on health and mortality: a shock to healthcare resources and a health shock.

As explained in section , a plausible path for the shock to resources goes via the quality of care. Controlling for healthcare quality closes this backdoor path (Cinelli, Forney, and Pearl 2022). We use infant mortality to control for the quality of healthcare. If due to the shock there is a reduction in healthcare resources (e.g. high quality physicians leave and are replaced by lower quality staff), we expect infant mortality to reflect this while not being directly related to our outcome variable: adult mortality.

Note that we cannot use quality measures like number of physicians or MRI scanners per 100k population. If due to high oop and poverty many people forgo treatments, these measures will be low and mortality is high. These are mediators of the mechanism that we are analyzing and controlling for a mediator introduces a bias instead of resolving one (Cinelli, Forney, and Pearl 2022).

Recall that there is no confounding if reduced resources lead to longer waiting lists (unmet needs and thus lower health status) with no direct effect of resources on health (say, through healthcare quality): such a pathway is captured by the model.

Health shocks are partly captured by the term \(\gamma \ln(m_{a-1,g,2,t-1}/\bar m_{a-1,g})\). To further close the health shock path, we use the fraction of people who state their (perceived) health is “bad or very bad” to control for health shocks over time. If the estimated effect of unmet medical needs on mortality is caused by health shocks, including the self reported health variable should reduce the effect thereby signaling that our baseline estimates overestimate the true effect.

Hence we re-estimate our two equations by adding to the TooExp equation the variable self-perceived health and to the mortality equation we add both infant mortality and self-perceived health. With the newly estimated parameters we calculate the percentage increase in mortality due to a 500 euro increase in oop using the equation for \(dm/m\). Then we divide this increase in mortality by the baseline estimated mortality increase. If this ratio is clearly smaller than one, this is an indication that the confounding effects are important: including variables to control for healthcare quality and health shocks reduces the estimated effects considerably.

As shown in Figure 8, the distribution of the relative effects (across traces and countries) has a mode (slightly) above one and quite some variation around it. We know from Figure 2 that there is substantial uncertainty about this effect. This is reflected in the dispersion of the histogram in Figure 8.

As the mass of the effect in the histogram is at one and above, there are no strong reasons to believe that the estimated effect in the baseline model is biased due to shocks to healthcare resources or to health itself. Once we control for such shocks we find effects that are on average in line with the baseline effects.

We present three robustness checks in terms of variables. First, we use the at-risk-of-poverty variable, instead of deprivation as our variable to capture the fraction of people on low income. Second, we extend our definition of oop costs with expenditures on voluntary health insurance. Third, instead of focusing on the coefficient \(\beta_{unmet}\) in the mortality equation, we work with TooExp in this equation. The idea here is that other reasons for unmet needs (than too expensive) may have bigger effects on mortality. If this would be the case, the baseline model overestimates the effect of too expensive (and hence of oop) on mortality. Although effect size for some countries can change with some of the checks, the overall point that an increase in oop expenditure raises mortality is robust.

We summarize the robustness checks in Figure 9. The figure has the baseline average effects of Figure 2 on the horizontal axis. For a number of countries these effects are smaller than 5 (per 1000 dead). In all of the robustness checks, these effects remain small.

The four countries with baseline effects exceeding five are explicitly named in the figure. For a given baseline effect, there are three (vertically aligned per country) robustness effects. Ideally, all three points would lie on the 45-degree line. But, obviously, there is some variation in the effect size for different specifications. We know from Figure 2 that there is (considerable) uncertainty about the effect size in the baseline model. Figure 9 illustrates this uncertainty per country by the dashed vertical lines indicating the 95% interval in the baseline outcome.

Including voluntary health insurance payments in our oop measure has the smallest effect on outcomes. The points are very close to the 45-degree line: the effects are basically the same as with the baseline specification.

Using too expensive instead of unmet in the mortality equation yields smaller effects for Hungary and Bulgaria, but the differences with the baseline are relatively small. The ranking of the effect size across the four countries (with significant effects) remains the same.

The differences are bigger if we use the at-risk-of-poverty rate instead of deprivation. Especially in Bulgaria and Hungary, the effects are substantially smaller and outside the 95% interval of the baseline effect. As at-risk-of-poverty is a relative poverty measure, one would expect the effect of this variable on too expensive and hence on unmet needs to be smaller than with deprivation. Being relatively poor in a country does not imply that people lack the funds to pay for treatments in the way deprivation does. This considerably reduces the effect of oop via poverty on mortality.

Overall the picture is that the baseline result that oop expenditure affects mortality is fairly robust to the use of different variables and different specifications.

First we generate the posterior samples for each robustness check. Then we use these samples to create Figure 9.

This section illustrates with two figures the relation between the fraction of people with unmet medical needs due to financial constraints, TooExp, and our variable measuring how generous a health insurance system is, OOP. We illustrate this using data simulated from the model and for one country, Romania, in our data set.

As we assume that TooExp has a logit-normal distribution, we plot the log-odds of TooExp both for Romania and for our simulated data.

We plot these log-odds on the vertical axis and the interaction term OOP*Poverty on the horizontal axis.

df_RO = df[df.country=='Romania']
df_RO['OOP'] = df_RO['HF3_PC_CHE']/100
df_RO['interaction'] = df_RO['HF3_PC_CHE']*\
    df_RO['deprivation']/(100*100)
df_RO['tooexp'] = df_RO['TOOEXP']/100
df_RO['tooexp_lo'] = np.log((df_RO['TOOEXP']/100)/\
                            (1-(df_RO['TOOEXP']/100)))

We use our model to simulate this parametric relation between TooExp and OOP*Poverty by varying the underlying exogenous parameters \(\xi\) and \(D\). As this is just used as illustration, we specify simple/straightforward functions for utility \(u\) and cumulative distribution function \(G\).

α = 0.2
ζ = 0.25
x1 = 5000
x0 = 400
σ_x = 0.8
σ_0 = 0.3
π_l = 0.3
π_h = 0.1
y_h = 20000
y_l = 3000

def u(y):
    return np.sqrt(y)
def G(y,oop):
    return 1-np.exp(-1*(σ_0/σ_x * u(y)/u(y-oop)))
def oop(ξ,D,y):
    return (ζ*ξ*(1-G(y,ξ*x0))*x0+(1-ζ)*(1-G(y,D))*D)
def expend(ξ,D,y):
    return (ζ*(1-G(y,ξ*x0))*x0+(1-ζ)*(1-G(y,D))*x1)
def Expend(ξ,D,α,π_l):
    return α*π_l*expend(ξ,D,y_l)+(1-α)*π_h*expend(ξ,D,y_h)
def OOP(ξ,D,α,π_l):
    return (α*π_l*oop(ξ,D,y_l)+(1-α)*π_h*oop(ξ,D,y_h))/\
        Expend(ξ,D,α,π_l)
def TooExp(ξ,D,α,π_l):
    return α*π_l*(ζ*(G(y_l,ξ*x0)-G(y_l,0))+\
                  (1-ζ)*(G(y_l,D)-G(y_l,0)))+\
                  (1-α)*π_h*(ζ*(G(y_h,ξ*x0)-G(y_h,0))+\
                             (1-ζ)*(G(y_h,D)-G(y_h,0)))

range_D = np.arange(500, 2500,75)[:,None]
range_ξ = (np.arange(0.4,1.00,0.1))[None,:]

Y1 = TooExp(range_ξ,range_D,α,π_l).flatten()+\
    np.random.normal(0,0.0015,size=range_D.shape[0]*range_ξ.shape[1])
Y1_lo = np.clip(np.log(Y1*10/(1-Y1*10)),-5.0,1)
oop_1 = OOP(range_ξ,range_D,α,π_l).flatten()
X1 = oop_1*α
Z1 = Expend(range_ξ,range_D,α,π_l).flatten()*X1

plt.style.use('Solarize_Light2')
fig, (ax1,ax2) = plt.subplots(1, 2, sharex=True,\
                              sharey=True,dpi=140,figsize=(14,6))

fig.suptitle('Relation between OOP $\\times$ Poverty and TooExp')
ax1.scatter(X1,Y1_lo,s=oop_1*90)
ax1.set_xlabel('OOP $\\times$ Poverty')
ax1.set_ylabel('log odds TooExp')
ax1.set_title('Simulated data')
ax1.set_ylim(-4.5,-1.5)
fig2 = sns.scatterplot(ax=ax2,data = df_RO, \
                       x = 'interaction',y='tooexp_lo',\
                       hue='nuts2',size='OOP')
fig2.set(xlabel = 'OOP $\\times$ Poverty', \
         ylabel = 'TooExp', title = 'Romania');

/tmp/ipykernel_136527/144876270.py:36: RuntimeWarning: invalid value encountered in log
  Y1_lo = np.clip(np.log(Y1*10/(1-Y1*10)),-5.0,1)

The mortality equation is specified as follows:

where \(k_{ag2t},n_{ag2t}\) denote the number of deaths, population size resp. in an age/gender/region/year category. Mortality is given by

where we specify the prior as

with \(e^{-3}/(1+e^{-3}) = 0.05\) implying that people above 35 have –on average– still 20 years to live. This is an under-estimation of life expectancy that will be easily adjusted by the data that we have. Both here and for the TooExp equation, we chose relatively narrow priors for the fixed effects. The prior for the regional fixed effects is given by:

We allow \(\gamma\) to be positive or negative with a relatively small standard deviation (determined in a hierarchical way):

Note that \(\gamma\) in the main text is denoted beta_lagged_log_mortality in the code. The standard deviation of this prior is determined by the following prior distribution (this is the hierarchy):

Finally, we know from the theory that the following two parameters are non-negative:

In words, sd_prior_beta captures the extent to which (social) external factors (like poverty, previous period health etc.) can affect mortality at all, instead of it being mainly determined by biology (age/gender).

Next, we turn to the equation for TooExp. We estimate TooExp with a logit-normal distribution.

where the parameter \(\mu\) of the Normal distribution of the log-odds of TooExp is given by the equation in Lemma 1 and the parameter \(\sigma\) has a prior

which allows for a range of values for \(\sigma\) that are positive. For the regional fixed effects, we have

where the \(-5.0\) is approximately equal to the mean log-odds of TooExp (which equals \(-5.17\)) and we chose a relatively narrow standard deviation for the fixed effects. For the time fixed effects we have

Since we know that the following two parameters are non-negative, we specify a half-normal distribution

where the prior of sd_prior_b is given by (hierarchical model):

As can be seen in Table 1, the variable poverty has fewer observations than deprivation. In Bayesian analysis there is a natural way to deal with missing values which is an improvement on two standard ways of dealing with this: (i) dropping observations (rows) with missing values (sometimes called complete case analysis) and (ii) interpolating the missing values. The former would change our sample when comparing the results with poverty and with deprivation. Interpolating data, say by replacing a missing value with the mean value of the variable makes the estimation method “too confident” about this value, thereby negatively affecting the quality of the inference.

We use the following method to deal with missing values (As desribed in McElreath 2020). The uncertainty surrounding the value of a missing observation is taken into account in the posterior distributions of our parameters. When sampling the posterior, if we encounter a missing value in a variable, this value is drawn from its distribution. We work with 4000 samples for the posterior and hence we draw 4000 different values for each missing value. In this way, the uncertainty about the (missing) value translates into posterior uncertainty of the parameters and predictions. The web-appendix provides the details on how this is implemented.

Figure 11 gives the trace plots for the parameters that we are interested in. That is, we leave out the traces for the (age, calendar year and region) fixed effects.

As explained in the main text, we are interested in three features in the plots on the right. First, the plots are stationary; second, condensed zig-zagging and third, the four chains cover the same regions of the parameter space. These features are satisfied for our coefficients of interest.

Table 5 provides summary statistics for the posteriors of the coefficients we are interested in. For some countries the hdi_3% lower bound for the b_oop or b_interaction equals zero. This is compatible with the OOP effect on mortality being bounded away from zero as the coefficients can be correlated: say, low b_oop going together with high b_interaction leading to an overall strictly positive effect.

Hence, to understand the mortality effects of an increase in oop, we use the equation for \(dm/m\) in Lemma 1 with the posterior distributions substituted in for all the parameters. This gives us Figure 2.