Dynamic effects of health insurance reform

In order to have sufficient data for all our variables we analyze data from 2008-2018. We start with the data from 2008-2020. In the code below we make some further selections like positive population levels and a population that exceeds the number of deaths in a category. At the time of our analysis this narrows down the years to 2008-2018. The categories we are analyzing are indexed by age/gender/region/income as explained below.

We analyze mortality for ages 45-85. At lower ages, mortality is very close to zero; for ages above 85 the number of observations per age group drops rapidly.

low_year = 2008
high_year = 2020
low_age = 45 
high_age = 85
ds = xr.open_dataset('./data/oop_health_data.nc')
ds

<xarray.Dataset> Size: 22MB
Dimensions:                           (nuts2: 313, year: 11, sex: 2, age: 99)
Coordinates:
  * nuts2                             (nuts2) <U4 5kB 'AT11' 'AT12' ... 'UKN0'
  * year                              (year) int64 88B 2008 2009 ... 2017 2018
  * sex                               (sex) <U1 8B 'F' 'M'
  * age                               (age) float64 792B 1.0 2.0 ... 98.0 99.0
    country                           (nuts2) <U14 18kB ...
Data variables:
    population                        (age, sex, year, nuts2) float64 5MB ...
    country_code                      (age, sex, year, nuts2) <U2 5MB ...
    HF3_PC_CHE                        (age, sex, year, nuts2) float64 5MB ...
    deaths                            (age, sex, year, nuts2) float64 5MB ...
    percentage_material_deprivation   (year, nuts2) float64 28kB ...
    TOOEXP                            (year, nuts2) float64 28kB ...
    UNMET                             (year, nuts2) float64 28kB ...
    infant mortality                  (year, nuts2) float64 28kB ...
    number_physicians per inhabitant  (year, nuts2) float64 28kB ...

Following python code generates the list my_regions where we have data on population,deaths and where we have at least one year (per region) with information on UNMET or where UNMET is not equal to zero for all years.

# my_regions = []
# for r in np.unique(ds.nuts2.values):
#     A = np.sum(ds.where(ds.nuts2==r,drop=True).population.isnull().values)
#     B = np.sum(ds.where(ds.nuts2==r,drop=True).deaths.isnull().values)
#     C = np.sum(np.isnan(np.nanmean(ds.where(ds.nuts2==r,drop=True).UNMET.values,axis=0)))
#     D = np.sum(np.nanmean(ds.where(ds.nuts2==r,drop=True).UNMET.values,axis=0)==0)
#     if A + B + C + D == 0:
#         my_regions.append(r)
# my_regions

cohorts = 30
age_matrix = sliding_window_view(np.arange(low_age,high_age+1), cohorts)[:-1]
age_indices = age_matrix-low_age
my_regions = ['BG31', 'BG32', 'BG33', 'BG34', 'BG41', 'BG42', 'CH01', 'CH02', 'CH03', 'CH04', 'CH05', 'CH06', 'CH07', 'DK01', 'DK02', 'DK03', 'DK04', 'DK05', 'EL30', 'EL41', 'EL42', 'EL43', 'EL51', 'EL52', 'EL53', 'EL54', 'EL61', 'EL62', 'EL63', 'EL64', 'EL65', 'FI1B', 'FI1C', 'FI1D', 'HU11', 'HU12', 'HU21', 'HU22', 'HU23', 'HU31', 'HU32', 'HU33', 'LT01', 'LT02', 'NO01', 'NO02', 'NO03', 'NO04', 'NO05', 'NO06', 'NO07', 'SE11', 'SE12', 'SE21', 'SE22', 'SE23', 'SE31', 'SE32', 'SE33', 'SI03', 'SI04', 'SK01', 'SK02', 'SK03', 'SK04']
ds = ds.where((low_age <= ds.age) & (ds.age <= high_age) &\
              (low_year <= ds.year) & (ds.year <= high_year) &\
              (ds.nuts2.isin(my_regions)) &\
              (ds.country_code != "") &\
              (ds.population > 1.0) & (ds.population > ds.deaths), drop=True)
ds["TOOEXP"] = ds["TOOEXP"].isel(age=0,sex=0)
ds["number_physicians per inhabitant"] = ds["number_physicians per inhabitant"].isel(age=0,sex=0)
ds["infant mortality"] = ds["infant mortality"].isel(age=0,sex=0)
ds["UNMET"] = ds["UNMET"].isel(age=0,sex=0)
ds["UNMET_other"] = ds["UNMET"]-ds["TOOEXP"]
ds["percentage_material_deprivation"] = ds["percentage_material_deprivation"].isel(age=0,sex=0)
ds

<xarray.Dataset> Size: 2MB
Dimensions:                           (age: 41, sex: 2, year: 11, nuts2: 65)
Coordinates:
  * nuts2                             (nuts2) <U4 1kB 'BG31' 'BG32' ... 'SK04'
  * year                              (year) int64 88B 2008 2009 ... 2017 2018
  * sex                               (sex) <U1 8B 'F' 'M'
  * age                               (age) float64 328B 45.0 46.0 ... 84.0 85.0
    country                           (nuts2) <U14 4kB 'Bulgaria' ... 'Slovakia'
Data variables:
    population                        (age, sex, year, nuts2) float64 469kB 5...
    country_code                      (age, sex, year, nuts2) object 469kB 'B...
    HF3_PC_CHE                        (age, sex, year, nuts2) float64 469kB n...
    deaths                            (age, sex, year, nuts2) float64 469kB 1...
    percentage_material_deprivation   (year, nuts2) float64 6kB nan nan ... 10.0
    TOOEXP                            (year, nuts2) float64 6kB 12.8 ... 0.4
    UNMET                             (year, nuts2) float64 6kB 27.4 ... 6.5
    infant mortality                  (year, nuts2) float64 6kB 9.1 8.0 ... 8.8
    number_physicians per inhabitant  (year, nuts2) float64 6kB 373.1 ... 344.6
    UNMET_other                       (year, nuts2) float64 6kB 14.6 9.5 ... 6.1

The age_matrix defined above has the following format. In our first year, 2008, we start with people aged 45-74. We follow these age-cohorts over time till 2018 when they are aged 55-84. We summarize this in the paper in Table 1.

print(tabulate(age_matrix, tablefmt="orgtbl"))

| 45 | 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 |
| 46 | 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 |
| 47 | 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 |
| 48 | 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 |
| 49 | 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 |
| 50 | 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 |
| 51 | 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 |
| 52 | 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 |
| 53 | 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 |
| 54 | 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 |
| 55 | 56 | 57 | 58 | 59 | 60 | 61 | 62 | 63 | 64 | 65 | 66 | 67 | 68 | 69 | 70 | 71 | 72 | 73 | 74 | 75 | 76 | 77 | 78 | 79 | 80 | 81 | 82 | 83 | 84 |

country = ds.country.values
n_regions = ds.nuts2.values.shape[0] # number of regions 
n_countries = len(np.unique(ds.country.values)) # number of countries
n_ages = ds.age.values.shape[0] # number of ages

# Mapping of regions to countries
# Example: regions 0, 1 -> country 0; regions 2, 3 -> country 1; regions 4, 5 -> country 2
region_to_country = pd.factorize(country)[0]

# Create the matrix (countries x regions)
country_region_matrix_01 = np.zeros((n_countries, n_regions)) # 0-1 matrix with a 1 at position cr if region r belongs to country c
for r,c in enumerate(region_to_country):
     country_region_matrix_01[c, r] = 1

# Normalize rows to compute averages per country (across regions)
country_region_matrix = country_region_matrix_01 / country_region_matrix_01.sum(axis=1, keepdims=True)

For some country/year combinations we see in the data that OOP (HF3_PC_CHE) equals zero. As in no European country healthcare is completely free, we interpret 0 as a missing observation. With xarray we can turn a 0 into nan by setting da = da.where(da != 0) where da denotes the data array. This is the code we use below for OOP, infant mortality and Unmet. No infant mortality in a region or no unmet medical needs is also highly unlikely; hence a zero value is turned into nan. If Unmet equals zero, the same is true for UNMET_other and TooExp because UNMET = TooExp + UNMET_other.

OOP does not vary with age or gender and hence we can take the average across dimensions (0,1) without losing information. The result is a variable varying over years and across regions. But out-of-pocket varies by country (not by region within a country) and hence we multiply OOP with the country_region_matrix (transformed). This makes sure that OOP_c varies by year and country.

In our empirical analysis we use the standardized version of the number of physicians per 100k inhabitants. Variables like OOP, material deprivation and unmet medical needs are turned from percentages into fractions for the analysis. We also turn infant mortality (per 1000 life-births) into a variable between 0 and 1 by dividing it by 100 (effectively making it infant mortality per 100k life-births).

def standardize(x):
    return (x-np.nanmean(x))/np.nanstd(x)
Population = ds.population.values.astype(int)
physicians = standardize(ds['number_physicians per inhabitant'].values)
OOP = ds.HF3_PC_CHE.where(ds.HF3_PC_CHE!=0)/100.0  # replaces OOP equal to 0 with nan
m_infant = ds['infant mortality'].where(ds['infant mortality']!=0)/100.0  # replaces infant mortality equal to 0 with nan
OOP_c = np.dot(OOP.mean(axis=(0,1)),country_region_matrix.T)
M = ds.deaths.values.astype(int)
Deprivation = ds.percentage_material_deprivation/100.0 
Unmet = ds.UNMET.where(ds.UNMET!=0)/100.0  # replaces UNMET equal to 0 with nan
TooExp = ds.TOOEXP.where(ds.UNMET!=0)/100.0
Unmet_o = ds.UNMET_other.where(ds.UNMET!=0)/100.0 
U_log_odds = np.log((Unmet)/(1-Unmet))
U_o_log_odds = np.log((Unmet_o)/(1-Unmet_o))

We analyze population and mortality in tensor-format. The dimensions of the tensors P_cohort, M_cohort are calendar year, cohort size (first row in Table 1), gender and region.

c_dim, g_dim, t_dim, r_dim = cohorts, 2, len(ds.year), n_regions
n_dim = n_countries
P_cohort = np.zeros((t_dim, c_dim, g_dim, r_dim))
M_cohort = np.zeros((t_dim, c_dim, g_dim, r_dim))

for offset in range(c_dim):
  for g in range(g_dim):
     for r in range(r_dim):
        P_cohort[:,offset,g,r] = np.diagonal(Population, axis1=0, axis2=2, offset=-offset)[g,r,:]
        M_cohort[:,offset,g,r] = np.diagonal(M, axis1=0, axis2=2, offset=-offset)[g,r,:]

We use the following moments to specify our priors for missing values in these variables. To deal with missing values in Deprivation we model it as a distribution from which we draw missing values. Such a distribution cannot have a standard deviation of zero and hence we add a small number to the standard deviation to avoid this. See section for a discussion of missing values.

There are four regions in Norway on which we have no information at all about the number of physicians. For these regions we set mean_physician equal to the average number of physicians in Norway and std_physician equal to the average standard deviation for Norway.

Finally, to define the log-odds of mortality (number of deaths M over population P) we add 1 to M to avoid taking the logarithm of zero.

mean_depr = Deprivation.mean(axis=(0)).values
std_depr = Deprivation.std(axis=(0)).values + 0.0001
mean_oop_c = np.nanmean(OOP_c,axis=0)
std_oop_c = np.nanstd(OOP_c,axis=0)
mean_m_infant = np.nanmean(m_infant,axis=0)
std_m_infant = np.nanstd(m_infant,axis=0)
mean_physician = np.nanmean(physicians,axis=0) # mean over years 
std_physician = np.nanstd(physicians,axis=0)   # std over years
norway = ['NO' in r for r in ds.nuts2.values]
mean_norway = np.nanmean(mean_physician[norway])
mean_physician = np.nan_to_num(mean_physician,nan=mean_norway)
std_norway = np.nanmean(std_physician[norway])
std_physician = np.nan_to_num(std_physician,nan=std_norway)
log_odds_mort = np.log(((1+M)/Population)/(1-((1+M)/Population)))
mort_c = np.log(np.tensordot(((1+M_cohort)/P_cohort), country_region_matrix.T, axes=([3],[0]))/(1- np.tensordot(((1+M_cohort)/P_cohort),\
                                              country_region_matrix.T, axes=([3],[0])))) # varies by age, cohorts, gender, country
gender_effect_π = (mort_c.mean(axis=(0,1,3))-mort_c.mean(axis=(0,1,2,3))) # varies by gender
country_effect_π = (mort_c.mean(axis=(0,1,2))-mort_c.mean(axis=(0,1,2,3))) # varies by country

/var/folders/g3/0hgyvgbn6x3105r_s_rgnrx40000gn/T/ipykernel_11128/260719900.py:7: RuntimeWarning: Mean of empty slice
  mean_physician = np.nanmean(physicians,axis=0) # mean over years
/Users/janboone/anaconda3/envs/pymc_env/lib/python3.12/site-packages/numpy/lib/nanfunctions.py:1879: RuntimeWarning: Degrees of freedom <= 0 for slice.
  var = nanvar(a, axis=axis, dtype=dtype, out=out, ddof=ddof,

The following python code generates Table 2 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>>
<<code-selection>>
stat_list = []
def summary_statistics(var,name,stat_list=stat_list):
   return  stat_list.append([name,np.sum(~np.isnan(var)),np.nanmean(var).round(decimals=1),np.nanstd(var).round(decimals=1),np.nanmin(var).round(decimals=1),np.nanmedian(var).round(decimals=1),np.nanmax(var).round(decimals=1)])

summary_statistics(P_cohort,'population')
summary_statistics(M_cohort,'deaths')
summary_statistics(M_cohort/P_cohort*100,'mortality (%)')
variables = [Deprivation.values*100,TooExp.values*100,Unmet.values*100,m_infant.values*100,ds['number_physicians per inhabitant'].values,OOP_c*100]
names = ['deprivation (%)', 'too expensive (%)','unmet (%)','infant mortality (\\textperthousand)','physicians per 100k inhab.','out-of-pocket (%)']
for i in range(len(variables)):
   summary_statistics(variables[i],names[i])

df = pd.DataFrame(stat_list,columns= ['','count','mean','std','min','median','max'])

headers = ['count','mean','std','min','median','max']

print(tabulate(df,headers,tablefmt="orgtbl",\
               colalign=("left","left", "right", "right", "right", "right", "right", "right")))

Our identification strategy is grounded in a structural modeling approach rather than exogenous quasi-experimental variation. We estimate a dynamic Markov model of health status transitions using cohort-level Eurostat data, where the key policy variable, \(oop\), varies across countries and over time. We link these variations in cost-sharing to observed differences in recovery rates and mortality, conditional on region, age, gender, income, time fixed effects and variables capturing healthcare quality.

The strength of this approach lies in its ability to simulate long-run effects by modeling underlying health dynamics. Unlike DiD estimators, which compare limited time points before and after a reform, our framework accumulates the effects of delayed or forgone treatment over multiple years. It also allows us to examine heterogeneity by income and gender and to propagate parameter uncertainty through Bayesian inference.

However, this comes with trade-offs. Our identification of the effects of cost-sharing relies on structural assumptions, notably the specification of the Markov transition probabilities. Further, \(oop\) variation is not derived from a policy discontinuity or random assignment. As such, our estimates could be biased by omitted variables that simultaneously influence health and insurance \(oop\) (e.g. austerity-induced cuts to both healthcare quality and insurance generosity). We address this concern through robustness checks: (i) the inclusion of unemployment and education to capture lifestyle differences between regions and GDP per capita and imaging device availability as additional controls for healthcare quality and (ii) a placebo test using traffic mortality which should not have the dynamic structure of mortality that we analyse but can be negatively affected by fiscal austerity.

While these checks support the plausibility of our identification, we acknowledge that the approach does not deliver the same level of internal validity as natural experiments. Our findings should therefore be interpreted as credible simulations under a structural model, rather than precise causal point estimates in the experimental sense.

In estimating the model we follow the age cohorts 45-74 over time from 2008 to 2018 when they are aged 55-84. For each age/gender/income/region category we split the population into low and high health status. We use pymc’s scan function to follow the cohort over time according to the difference equation \eqref{org25cccbc}. The function transition takes as inputs \(\pi,\sigma,\delta\) and this period’s fraction of ill. The function returns next period’s fraction of ill.

The dimensions for observed variables (“data”) and parameters are: age, gender, year, region and income. Not all variables or parameters have all these dimensions. E.g. some only vary by country (not region), some parameters are constant or only vary by year etc.

For our first year (2008) we specify the initial fraction of ill, \(\iota_0\). This parameter needs to be estimated and varies by cohort (45-74), gender, region and income. Most probabilities/fractions (like \(\iota_0\) and \(\delta\)) are specified in log-odds terms.

The probability of falling ill is specified as the addition of age, gender, country and income effects (additive in log-odds space).

Then the code specifies the distributions of a number of variables for which we have missing values. This is the case for the deprivation variable, OOP, unmet medical needs in total and for reasons other than the treatment is too expensive. Finally, there are also missing values for infant mortality and the number of physicians per 100k inhabitants. For variables that are between 0 and 1, we either use the clip function or specify a Beta distribution. For the Beta distribution we use the parametrization with mean effect mu and number of observations nu.

The expected number of people with unmet medical needs equals the sum of unmet needs for other reasons plus the fraction of people forgoing treatment because it is too expensive, \(\tau\). We observe \(\tau\) and split it up into \(\tau^l,\tau^h\) such that \(\tau = \alpha \tau^l + (1-\alpha) \tau^h\). Similarly, We observe \(U_o\) and split it up into \(U_o^l, U_o^h\) such that \(U_o = \alpha U_o^l + (1-\alpha) U_o^h\). We do this as follows: \(U_o^j = \xi^j U_o\) with \(\alpha \xi^l + (1-\alpha) \xi^h = 1\) with \(\xi^l = \xi + \varepsilon \geq \xi = \xi^h\). Hence it follows from \(\alpha (\xi + \epsilon) + (1-\alpha) \xi = 1\) that \(\xi = 1-\alpha \varepsilon\) as we model it below in the code.

We know that an increase in oop makes it more likely that people skip treatment because it is too expensive. Hence we specify that \(\zeta^h \geq 0\) (in terms of slope_ζ) and that \(\zeta^l \geq \zeta^h\) using delta_ζ. Both parameters are modelled as HalfNormal to make sure they are non-negative.

We model the log-odds of \(\lambda\) as decreasing in infant mortality and increasing in the number of physicians.

The tensor \(I\) keeps track of the fraction ill per year, cohort, gender, region and income category. Its first “row” is given by \(\iota_0\) which we discussed above.

To determine the probability of death per year, cohort, gender and region, we sum over the income dimension as our data does not specify deaths by income category.

Given this, the number of deaths, M_cohort is a Binomial distribution with parameter \(n\) equal to the population size P_cohort and \(p\) equal to \(\mu\).

def transition(π,σ,I,δ):
    fit = I*σ+(1-π)*(1-I)
    ill = π*(1-I)+pt.clip(1-σ-δ,0.0,1.0)*I
    fraction_ill = ill/(fit+ill)
    return fraction_ill

n_y = 2 # low/high income
n_g = 2 # female/male
countries = np.array(['Bulgaria', 'Switzerland', 'Denmark', 'Greece', 'Finland', 'Hungary', 'Lithuania', 'Norway', 'Sweden', 'Slovenia', 'Slovakia'])
# index: age, gender, year, region, income
# for some variables we use country instead of region

coords = {
    "region" : ds.nuts2.values,\
    "ages" : ds.age.values,\
    "cohort" : age_matrix[0],\
    "gender" : ds.sex.values,\
    "year" : ds.year.values,\
    "income" : ['y_l','y_h'],\
    "country": countries,\
     # we define dummy dimensions for parameters that do not
     # vary by this dimension
     "dummy0": np.arange(1),\
     "dummy1": np.arange(1),\
     "dummy2": np.arange(1)}

with pm.Model(coords=coords) as AR_model:
    logodds_ι = pm.Normal("logodds_ι",-4.0,1,\
              dims=("cohort","gender","region","income"),\
              shape=(cohorts,n_g,n_regions,n_y)) 
    ι_0 = pm.Deterministic("ι_0",pm.math.invlogit(logodds_ι),dims=("cohort","gender","region","income"))
    # ι_0 is used to start the sequence in scan =>
    # no need to add age dimension: we build up the matrix I by following each cohort over the years
    logodds_δ = pm.Normal("logodds_δ",-2.5,0.05,dims=("country")) # conditional on being ill, same prob. of death
    # δ is a non_sequence in the scan loop; no need to add age as a dimension
    δ = pm.Deterministic("δ",pt.dot(pm.math.invlogit(logodds_δ), country_region_matrix_01)[None,None,None,:,None])
    # we build up the prior for π
    logodds_age_π = pm.Normal("logodds_age_π",log_odds_mort.mean(axis=(1,2,3)),log_odds_mort.std(axis=(1,2,3)),dims=("ages"))
    logodds_sex_π = pm.Normal("logodds_sex_π", gender_effect_π, 0.05, dims=("gender"))
    logodds_country_π = pm.Normal("logodds_country_π",country_effect_π,0.05, dims=("country"))
    logodds_π = pm.Deterministic("logodds_π",logodds_age_π[:,None,None] + logodds_sex_π[None,:,None] + pt.dot(logodds_country_π, country_region_matrix_01)[None,None,:]) # prob. of falling ill varies by age, gender and country
    # the priors do not impose that women tend to be healthier than men
    # but we do impose that people on low income tend to have low health status (high π)
    delta_π = pm.HalfNormal("delta_π",sigma=0.05,\
          dims=("ages","dummy0","dummy1"),\
          shape=(n_ages,1,1))
    π = pm.Deterministic("π",pt.stack([pm.math.invlogit(\
                                    logodds_π+ delta_π),\
                   pm.math.invlogit(logodds_π)],axis=-1), dims=("ages","gender","region","income"))
    # when we draw missing values for OOP, we should draw the same value for all regions/age/gender in the same country
    # this we do via oop_c which we later bring back to region dimension with country_region_matrix_01
    pov = pm.Normal("pov", mean_depr[None,:], std_depr[None,:], observed=Deprivation,dims=("year","region")) # fraction of people on low income
    poverty = pt.stack([pt.clip(pov,0,1),1-pt.clip(pov,0,1)],axis=-1) # vector with fraction of low/high income people
    oop_c = pm.Normal("oop_c", mean_oop_c[None,:], std_oop_c[None,:], observed=OOP_c,dims=("year","country"))
    oop = pt.dot(oop_c,country_region_matrix_01)

    slope_ζ = pm.HalfNormal("slope_ζ",sigma=0.05)
    delta_ζ = pm.HalfNormal("delta_ζ",sigma=0.05)
    ζ = pm.Deterministic("ζ",pt.stack([(slope_ζ+ delta_ζ),\
                   (slope_ζ)],axis=-1),dims=("income")) # varies with income
    τ = pm.Deterministic("τ", ζ[None,:]*oop[:,:,None])
    mu_U_o = Unmet_o.mean(axis=0).values[None,:]
    nu_U_o = pm.Exponential("nu_U_o", 0.05)
    U_o = pm.Beta("U_o", alpha = mu_U_o*nu_U_o, beta = (1-mu_U_o)*nu_U_o,dims=("year","region"), observed=Unmet_o)
    mu_U = pt.clip(U_o+pt.sum(poverty*τ,axis=2),0.000001,0.9999)
    nu_U = pm.Exponential("nu_U", 0.05)
    U = pm.Beta("U",alpha=mu_U*nu_U, beta=(1-mu_U)*nu_U,dims=("year","region"), observed=Unmet)
    # λ measures prob. patient heals if treated (i.e. one minus unmet)
    M_infant = pm.Normal("M_infant",mean_m_infant[None,:],std_m_infant[None,:],dims=("year","region"),observed=m_infant)
    Physicians = pm.Normal("Physicians",mean_physician[None,:],std_physician[None,:],dims=("year","region"),observed=physicians)
    lambda_1 = pm.HalfNormal("lambda_1",sigma=0.01)
    lambda_2 = pm.HalfNormal("lambda_2",sigma=0.01)
    lamba_0 = pm.Normal("lambda_0",-4.1,0.05,dims=("year"))
    logodds_λ = pm.Normal("logodds_λ",lamba_0[:,None]-lambda_1*M_infant+lambda_2*Physicians,0.1)
    λ = pm.Deterministic("λ",pm.math.invlogit(logodds_λ))# λ is between 0 and 1
    ε = pm.Beta("ε",1,3) # ε = ξ_l - ξ_h > 0
    ξ = 1-poverty[:,:,0]*ε
    υ = pm.Deterministic("υ", pt.stack([(ξ+ε)*U_o,ξ*U_o],axis=-1) + τ)
    σ_0 = pm.Beta("σ_0",1,3)
    σ = pm.Deterministic("σ",pt.clip(σ_0+λ[:,None,None,:,None]*(1-υ)[:,None,None,:,:],0,1))
    I = pt.zeros((t_dim,cohorts,n_g,n_regions,n_y))
    I = pt.set_subtensor(I[0],ι_0)

    outputs, updates = pytensor.scan(
        transition,
     # age/year is the dynamic part of the model
     # i.e. first dimension of π
        sequences=[dict(input=π[age_indices]),dict(input=σ)], # check sequences with order of arguments in =transition=
     # we loop over the first dimension: age
     # the initial value for I (at low_age)
     # outputs_info has no age dimension
        outputs_info=[dict(initial=ι_0)],
     # variables that have no age dimension:
        non_sequences=[δ[0]],
    )

    I = pt.set_subtensor(I[1:],outputs[:t_dim-1])
    # μ is probability of death for age/gender/year/region
    ill = pm.Deterministic("ill",I)
    μ = pm.Deterministic("μ",(pt.sum(δ*I*poverty[:,None,None,:,:],axis=4)))# aggregate over income
    mortality = pm.Binomial("mortality",n=P_cohort,p=μ,observed=M_cohort,dims=("year","cohort","gender","region"))

A first check of the model is to see whether it fits the data. Figure 2 shows for all our observations the observed number of deaths on the horizontal axis and the predicted number on the vertical axis. Further, for each prediction we add a dashed vertical line representing the 95% prediction interval (of the posterior predictive distribution). The predictions clearly follow the 45-degree line and almost all 95% prediction intervals intersect the 45-degree line. In this sense the model fit seems quite good.

Figure 7 in the appendix presents the trace plots and Table 8 r-hat values for estimated parameters. The r-hat values are close to one and the traceplots satisfy the following three criteria. First, the plots in the right column of the figure are 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 is 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. This is most easily checked in the left column of the trace plot. The three features are satisfied in the figure. All this indicates that the MCMC algorithm manages to sample the posterior distribution correctly.

To gain an intuition on how the model works, Figure 3 shows a number of estimated outcomes. The top-left panel shows the probability of falling ill, \(\pi\), as a function of age for women and men.⁵ This probability is quite low around age 45 and increases to 40% for men in their eighties. As shown in the panel below this one, the probability of death –conditional on being ill– is above 10% per year. In other words, falling ill in our model does not refer to a regular cough or breaking a leg. The observation that the average probability of falling ill is higher for men than for women is in line with the observation that women tend to live longer than men. The 95% posterior probability interval shows that there is quite some variation in the probability of falling ill at older ages but that uncertainty is small till age 75. In order not to crowd the figure we do not report the 95% intervals for males but they are of similar magnitude.

The top-right panel shows the country offsets (in log-odds terms) for four countries. There is a clear dichotomy between Switzerland and Norway on the one hand and Lithuania and Hungary on the other. But even between these countries there is hardly overlap between the distributions. Of these four countries, people in Switzerland are least likely to fall ill (at any age) and in Hungary most likely.

Similarly, the probability of dying conditional on being ill, \(\delta_{c}\), varies by country. Of the four countries presented, it is lowest for Denmark and Finland and highest for Hungary and Bulgaria.

Finally, the bottom-right panel shows the fraction of people ill for four categories: people aged 55 and 60 with low/high income. The fraction of people who are ill is higher for low than for high incomes, as one would expect. Also the fraction increases with age. The posterior distributions for these fractions do not overlap.

The analysis in this section makes the following points. First, the expected qaly’s for a 45 year old vary a lot between countries and within countries between gender and income categories. The concept of qaly’s is much discussed in the literature; here we simply use it as a summary measure of the fractions of people in good health, poor health and deceased. Second, when we increase demand-side cost-sharing by ten percentage points (\(\Delta oop = 0.1\)), there is a clear but modest reduction in qaly’s. Third, the effect is almost zero for people on high income but far higher for people on low income. In this sense, effect heterogeneity is high and considering the average effect (across incomes) is not informative. Finally, instead of comparing qaly’s lost over, say, a thirty year time period (for which one needs a dynamic model), we can compare the loss in qaly’s for the same age groups from one year before the increase in \(oop\) to one year after. The latter is close to a (static) measure used in most event study analyses. We show that the dynamic effect is an order of magnitude bigger than the estimated static effect.

Figure 4 summarizes the two applications that we analyse with the model. First, since our estimation method identifies the fraction of people who have low health status, we can derive qaly’s using equation \eqref{org7f9eb9b}. In the figure we use \(\theta=0.6\): the quality of life for someone who is ill (and has a 10% probability of dying per year) is 60% of the quality of life for someone in full health. This is an arbitrary choice and obviously different choices can be made here. Different values for \(\theta\) change the precise numerical values but not the main (qualitative) conclusions of the analysis; as we show in the robustness analysis.

The figure was made as follows. For each group (female/male, low/high income) of 45 year olds we create a population of 1000 individuals. At 45 we use our estimates of \(\iota_0\) to split the groups into low/high health status. With our estimated parameters \(\pi,\sigma,\delta\) we follow these individuals as they age over time till they are 85. Hence, for each group there are at max. 40,000 life-years (in full health) to be had over this period. For each subgroup gender/income we calculate the realized quality adjusted life-years and divide this by 40,000. The left panel of Figure 4 shows the qaly’s realized for low income men and for high income women. The panel shows that in, say Switzerland, low income males get 80% of the max. number of qaly’s and high income females get 10% of 40,000 qaly’s more. That is, the latter get 12.5% (\(0.10/0.80\)) more qaly’s than the former. This gender/income gap ranges all the way up to 20% in Hungary, Bulgaria and Lithuania.

The panel on the right shows what happens if \(oop\) is increased by 0.1: people pay an additional 10% of their healthcare expenditure out-of-pocket for each country in our data. It turns out that the average effect is around 3.0:⁶ among 1000 people, a region/gender/income group loses 3 life-years for 45-85 year olds due to \(\Delta oop=0.1\). This effect is not huge but clearly bounded away from 0.

What is more surprising is the effect heterogeneity for high and low incomes. The effect is almost ten times bigger for people on low income. People on high income hardly react in terms of healthcare expenditure to the increase in \(oop\) as they can afford it anyway. But for low incomes, the increase reduces their demand for treatment. This raises unmet medical needs because it is too expensive and hence, for this group, reduces the probability \(\sigma\) of recovering to the healthy state. This has two effects on their qaly’s. First, the quality of life in the ill state is lower than in the healthy state. Second, staying longer in the ill state, increases the probability of death (compared to being in the healthy state).

Finally, we compare the dynamic (across years) qaly effect derived above with a year-on-year effect for the same ages. To avoid making additional assumptions we do this for the ages 45-74 in the first row of Table 1. This is the cohort of ages that we use to estimate the model. Hence we have a well defined \(\iota_0\) fraction of people that are ill per age/gender/region/income category in the year before the reform. Then we use the equations for \(\tilde Q_0, \tilde Q_1\) defined in the theory section to find the static effect on qaly’s. We compare this to the dynamic estimate in equation \eqref{org7f9eb9b} for the same ages till \(A=75\).

As shown in Table 3, the dynamic low income estimates are lower than in Figure 4 (right panel): the dynamic estimates across ages till 75 are lower than over the full range till 85.

Comparing static and dynamic estimates shows that for this type of health model the dynamic estimates of qaly’s lost \(Q_0-Q_1\) are a factor twenty higher for low incomes than the static estimates \(\tilde Q_0 - \tilde Q_1\). For high incomes the factor is around seven. As noted in section , we expect the dynamic effect to be bigger than the static effect because the former takes the effect on the stock of ill/healthy into account. However, the size of the difference is larger than one might have expected.

This suggests that for health models where an underlying (latent) health variable is important, short-term estimates (as often with a DiD analysis) severely under-estimate the overall effect of a policy change.

One can see the following trade-off. A DiD approach leads to a cleaner causal identification. It often allows for a convincing way to keep other factors constant. In contrast, identification in our model (which allows for a dynamic analysis) is based on underlying theory but suffers from potential missing variable bias. We show that for reasonable parameter values, the static estimate seriously under-estimates the effect compared to a dynamic one taking account of the cumulative effects over a life-time. In particular, we simulate the data using our dynamic model as true data generating process. Then we use static estimators to see what effects they identify.

Figure 5 illustrates how a static estimate like DiD ignores the dynamics of health. The figure is an abstract representation of the DiD estimate of Finkelstein and McKnight (2008) analysing the introduction of Medicare in 1966 providing universal public health insurance to individuals aged 65 and over. To find the health effect of Medicare, the paper compares the average mortality rate of 65-74 year olds with the rate for 55-64 year olds both before and after the introduction of Medicare; the years 1952-1975.

The logic of comparing age category 55-64 with 65-74 is the underlying assumption that the two categories do not differ except for the introduction of Medicare. This assumption would be harder to motivate when comparing 55-64 with 82 year olds. The idea of the paper is that if Medicare has substantial health effects one expects that the mortality rate for 65-74 year olds falls compared to 55-64 year olds after 1966. The paper does not find such an effect and the figure illustrates one mechanism why this can be the case.

To place this analysis in the context of our paper, we simulate the effect of a 10% point increase in oop (\(\Delta oop = 0.1\)) in year 0 for people aged 65 and over. Hence, in our model this implies that mortality will increase for the 65+ category (not decrease as in the Medicare example). We can simulate both the reform and the counterfactual and hence can make the comparison within age category. This simplifies the analysis as there is no need to compare the effect to 55-64 year olds. For this comparison between the DiD and dynamic estimates we only need the estimated parameters to be reasonable, not perfectly accurate. As shown above, this assumption seems to be satisfied with our estimated coefficients.

We do the simulations for low income men in one region in Bulgaria aged 65 and above who face \(\Delta oop = 0.1\) in year 0. On the horizontal axis we plot age for this group. The vertical axis shows the years after the introduction in year 0. The color intensity of a point indicates the number of years an agent has faced this treatment (higher oop). In the first year everyone aged 65 and above faced the higher out-of-pocket for one year (bottom row in the figure). In the second year, people aged 66 and over faced this for two years and the 65 year olds only for one year etc.

The horizontal black line in the left panel indicates the comparison group used in the paper (65-74) for year 10 (last year in the Finkelstein and McKnight (2008) data). The figure illustrates how –from a dynamic point of view– the effect gets diluted as the 65 year olds in this group have only experienced the reform for one year, the 66 year olds for two years etc. for all years after the reform. In contrast the blue diagonal line gives the true dynamic effect when we follow 65 year olds in year 0 till they become 83 with 19 years of increased oop. The diagonal includes points with a darkness intensity that the horizontal line does not capture for any of the years 0-18. In this sense, the static estimate under-estimates the true dynamic effect.

The right panel shows the consequences of this. The vertical axis shows the increase in mortality by age normalized on the average mortality effect for 65-74 year olds. Hence, the static estimate for 65-74 year olds in this figure is one by construction. But due to the dynamics of health in our model, the effect is first below 1 (till age 69) and then increases to six times the static estimate for 80 year olds. In order to identify this (far bigger) effect, one has to estimate a dynamic model capturing the age effects instead of relying on the DiD estimate. It will depend on the application whether a static precise causal identification is preferable or an estimate taking the dynamics of the problem into account but lacking a convincing causal foundation.

As sampling takes some time, we generate the posterior samples once and save these. Here we read them in.

idata = az.from_netcdf("./trace/baseline_model_lambda.nc")

for the parameters where we do not have (too) many indices, we can summarize the results both in a table and with trace plots.

pd.set_option('display.max_rows', 75)

headers = ['mean', 'sd', 'hdi_3%', 'hdi_97%',\
           'ess_bulk', 'r_hat']
variables = ["logodds_sex_π","ζ","nu_U","nu_U_o","ε","lambda_1","lambda_2"]
df_summary = az.summary(idata,var_names=variables)[headers]
print(tabulate(df_summary,\
               headers,tablefmt='orgtbl',floatfmt=".2f"))

|                  |   mean |    sd |   hdi_3% |   hdi_97% |   ess_bulk |   r_hat |
|------------------+--------+-------+----------+-----------+------------+---------|
| logodds_sex_π[F] |  -0.19 |  0.03 |    -0.25 |     -0.13 |     349.00 |    1.01 |
| logodds_sex_π[M] |   0.50 |  0.03 |     0.44 |      0.56 |     354.00 |    1.01 |
| ζ[y_l]           |   0.27 |  0.01 |     0.25 |      0.30 |    8429.00 |    1.00 |
| ζ[y_h]           |   0.03 |  0.00 |     0.03 |      0.04 |    6048.00 |    1.00 |
| nu_U             | 299.81 | 18.94 |   264.86 |    336.40 |    7950.00 |    1.00 |
| nu_U_o           | 109.35 |  6.59 |    97.78 |    122.34 |    5584.00 |    1.00 |
| ε                |   0.07 |  0.07 |     0.00 |      0.19 |   10331.00 |    1.00 |
| lambda_1         |   0.01 |  0.01 |     0.00 |      0.02 |    6389.00 |    1.00 |
| lambda_2         |   0.06 |  0.01 |     0.04 |      0.07 |    1704.00 |    1.01 |

plt.rcParams['figure.constrained_layout.use'] = True
az.plot_trace(idata,var_names=variables);

The following plot shows the observed number of deaths per age/gender/year/region categories against the model prediction.

predictions = idata.posterior_predictive.mortality.\
              mean(dim=("chain","draw")).values
perc_predictions = np.percentile(idata.posterior_predictive.mortality.values,q=[2.5,97.5],axis=(0,1))
plt.vlines(M_cohort,perc_predictions[0],perc_predictions[1],color='grey',linestyles="dashed",linewidth=1,label="95% interval",alpha=0.4)
plt.scatter((M_cohort),(predictions))
# plt.scatter((M/Population).flatten(),(predictions/Population).flatten())
# plt.plot([0,0.2],[0,0.2],c='k')
plt.plot([0,860],[0,860],c='k')
plt.legend()
plt.xlabel("number of deaths")
plt.ylabel("predicted number of deaths");

We plot the age profile for the probability of falling ill, \(\pi\). As one would expect, this probability increases with age and is higher for men than for women. We plot the country offsets (in log-odds space) for four countries and the probability of death \(\delta_c\) for four (other) countries. Finally, we plot the fraction of people ill by income for two age categories.

plt.style.use('Solarize_Light2')
fig, ((ax1,ax2),(ax3,ax4)) = plt.subplots(2, 2, dpi=140,figsize=(14,14))

pi_avg = np.mean(idata.posterior.π.values,axis=(0,1,4,5))
perc_pi = np.percentile(idata.posterior.π.values,[2.5,97.5],axis=(0,1,4,5))
# fig 1
ax1.plot(np.arange(low_age,high_age-1),pi_avg[:-2,1],label="male")
ax1.plot(np.arange(low_age,high_age-1),pi_avg[:-2,0],label="female")
ax1.vlines(np.arange(low_age,high_age-1),perc_pi[0,:-2,0],perc_pi[1,:-2,0],color='black',linestyles="dashed",linewidth=1,label="95% interval")
ax1.set_title("Age profile of the probability of falling ill, $\\pi$,\n and the 95% interval for females")
ax1.set_xlabel('Age')
ax1.set_ylabel('Probability')
ax1.legend();
# fig 2
post_logodds_pi = idata.posterior.logodds_country_π.values
sns.kdeplot(ax=ax2,data=post_logodds_pi[:,:,1].flatten(),label="Switzerland")
sns.kdeplot(ax=ax2,data=post_logodds_pi[:,:,7].flatten(),label="Norway")
sns.kdeplot(ax=ax2,data=post_logodds_pi[:,:,6].flatten(),label="Lithuania")
sns.kdeplot(ax=ax2,data=post_logodds_pi[:,:,5].flatten(),label="Hungary")
ax2.set_xlabel("Log-odds falling ill")
ax2.set_title("Country offsets falling ill, log-odds $\\pi$,\nfor selected countries")
ax2.legend();
# fig 3
post_logodds_delta = idata.posterior.logodds_δ.values
sns.kdeplot(ax=ax3,data=pm.math.invlogit(post_logodds_delta[:,:,2].flatten()).eval(),label="Denmark")
sns.kdeplot(ax=ax3,data=pm.math.invlogit(post_logodds_delta[:,:,4].flatten()).eval(),label="Finland")
sns.kdeplot(ax=ax3,data=pm.math.invlogit(post_logodds_delta[:,:,5].flatten()).eval(),label="Hungary")
sns.kdeplot(ax=ax3,data=pm.math.invlogit(post_logodds_delta[:,:,0].flatten()).eval(),label="Bulgaria")
ax3.set_xlabel("Probability of dying")
ax3.set_title("Probability of death conditional on being ill, $\\delta_c$")
ax3.legend();
# fig 4
ill_by_income = idata.posterior.ill.values.mean(axis=(4,5))
ax4.hist(ill_by_income[:,:,10,0,1].flatten(),bins=70,density=True,label="55 year old, high income")
ax4.hist(ill_by_income[:,:,10,0,0].flatten(),bins=70,density=True,label="55 year old, low income")
ax4.hist(ill_by_income[:,:,10,5,1].flatten(),bins=70,density=True,label="60 year old, high income")
ax4.hist(ill_by_income[:,:,10,5,0].flatten(),bins=70,density=True,label="60 year old, low income")
ax4.set_xlabel("Fraction ill")
ax4.set_ylabel("Density")
ax4.set_title("Fraction of people ill for\n55 and 60 year olds per income category")
ax4.legend();

An advantage of doing Bayesian analysis in pymc is that it has implemented the “causal do” operator. In the code below we first run a simulation with the baseline value for unmet, \(\upsilon\), which equals the sum of too expensive and other stated reasons why patient decided to forego treatment. The model estimates the fraction of people that forego treatment because it is too expensive as \(\zeta oop\). Hence, if \(oop\) is increased by \(\Delta oop\), the fraction of people skipping treatment because it is too expensive increases by \(\zeta \Delta oop\). This we add to the baseline value of \(\upsilon\) for the counterfactual simulation.

υ_baseline = idata.posterior.υ.mean(axis=(0,1)).values
pov_no_missing = idata.posterior.pov.mean(axis=(0,1)).values
zeta = idata.posterior.ζ.mean(axis=(0,1)).values
delta_oop = 0.1
delta_τ = zeta.reshape(1,2)*delta_oop
iota_0 = idata.posterior.ι_0.mean(axis=(0,1))[0].values
AR_model_0 = do(AR_model, {"υ": υ_baseline})
AR_model_1 = do(AR_model, {"υ": υ_baseline + delta_τ})

Having defined the two models AR_model_0, AR_model_1, we sample from these models with the new values for unmet medical needs \(\upsilon\).

SEED=44
# Sample σ from baseline model
idata_0 = pm.sample_posterior_predictive(
    idata,
    model=AR_model_0,
    predictions=True,
    var_names=["σ","ill"],
    random_seed=SEED,
)
# Sample new σ with higher oop
idata_1 = pm.sample_posterior_predictive(
    idata,
    model=AR_model_1,
    predictions=True,
    var_names=["σ","ill"],
    random_seed=SEED,
)

We can save the posterior distributions of these (counter-factual) simulations.

idata_0.to_netcdf("./trace/baseline_model_do_0.nc")
idata_1.to_netcdf("./trace/baseline_model_do_1.nc")

./trace/baseline_model_do_1.nc

And read in these trace files when we need them:

idata_0 = az.from_netcdf("./trace/baseline_model_do_0.nc")
idata_1 = az.from_netcdf("./trace/baseline_model_do_1.nc")

In order to find the fractions of people who are ill and fit in these simulations, we need the parameters \(\pi,\delta\) and \(\sigma\) from the model. From these parameters, only \(\sigma\) differs between the baseline and the counterfactual simulation.

We want to project the effects for 40 years into the future and hence cannot use the year dimension in the model for \(\sigma\) (which is only 11 years). Therefore we take the latest year for \(\sigma\) as the most recent information for our simulations.

pi = idata.posterior.π.values.transpose((2,3,4,5,0,1))[:,:,:,:,:,:]
delta = idata.posterior.δ.values.transpose((2,3,4,5,6,0,1))[0]
sigma_0 = idata_0.predictions.σ.values.transpose((2,3,4,5,6,0,1))[-1] # use latest year
sigma_1 = idata_1.predictions.σ.values.transpose((2,3,4,5,6,0,1))[-1]

Since we are interested in the dynamic effects of the change in oop, we consider the longest run that we have in the model: people who are 45 at the start of the simulation and we follow this cohort until they are 84.

In order to calculate/normalize the life years, we start with 1000 forty-five year olds per category gender/income and region. We use our estimate of \(\iota_0\) to split the 1000 individuals into fit and ill. Then we run through our Markov model keeping track of who switches health state and who dies.

iota_0 = idata.posterior.ι_0.values.transpose((2,3,4,5,0,1))[0,None,:,:,:,:,:]
# outcomes baseline:
shape = (2,n_ages,n_g,n_regions,n_y,4,2000) # fit/ill, age, gender, region, income, chain, samples
N_0 = np.zeros(shape) 
N_0[0,0] = (1-iota_0)*1000 # fit
N_0[1,0] = iota_0*1000     # ill
# outcomes counterfactual:
N_1 = np.zeros(shape) # fit/ill, age, gender, region, income, samples
N_1[0,0] = (1-iota_0)*1000 # fit
N_1[1,0] = iota_0*1000     # ill
# follow these 45 year olds till they are 85
for j in range(1,n_ages):
   N_0[0,j] = (1-pi[j-1])*N_0[0,j-1] + sigma_0 * N_0[1,j-1]
   N_0[1,j] = pi[j-1]*N_0[0,j-1] + (1-sigma_0 - delta) * N_0[1,j-1] 
   N_1[0,j] = (1-pi[j-1])*N_1[0,j-1] + sigma_1 * N_1[1,j-1]
   N_1[1,j] = pi[j-1]*N_1[0,j-1] + (1-sigma_1 - delta) * N_1[1,j-1] 

To calculate life years over the lifetime (from 45-84), we count a healthy life year as 1.0 and a life year while ill at qaly_ill and add these together. When someone dies, there are no life years left for this person (weight 0.0). The variable delta_lifeyears calculates the difference between life years in the baseline simulation and the counterfactual. We adjust the ’qaly’ for the ill to a value smaller than 1.0 to capture the idea that with a probability of dying above 10% per year (conditional on being ill) suggests that patients do face significant restrictions in daily life.

qaly_ill = 0.6
lifeyears_0 = N_0[0].sum(axis=0) + qaly_ill*N_0[1].sum(axis=0)
lifeyears_1 = N_1[0].sum(axis=0) + qaly_ill*N_1[1].sum(axis=0)
delta_lifeyears = lifeyears_0 - lifeyears_1

In order to summarize the results we present two qaly outcomes per country. The following code presents the unweighted average of regions into countries.

lifeyears_per_country = np.tensordot(country_region_matrix,lifeyears_0,axes=([1,1]))
# index: country/gender/income/chain/samples

Plot qaly’s for low income men and high income women.

category_names = ['qaly low income men', 'qaly high income women']
labels = countries

data = np.stack((lifeyears_per_country.mean(axis=(3,4))[:,1,0],lifeyears_per_country.mean(axis=(3,4))[:,0,1]),axis=1)/40000
category_colors = plt.get_cmap('viridis')(
    np.linspace(0.15, 0.85, data.shape[1]))

fig, (ax1,ax2) = plt.subplots(1,2, figsize=(16, 6))
ax1.invert_yaxis()
# ax.xaxis.set_visible(False)
# ax.set_xlim(0, 1)
ax1.set_title("Fraction of max. qaly's realized for 45 year olds per country\n")

for i, (colname, color) in enumerate(zip(category_names, category_colors)):
    starts = i * data[:,0]
    widths = data[:, i] - starts
    ax1.barh(labels, widths, left = starts, height=0.5,
            label=colname, color=color)
    xcenters = starts + widths / 2

    r, g, b, _ = color
    text_color = 'white' if r * g * b < 0.5 else 'darkgrey'
    for y, (x, c) in enumerate(zip(xcenters, widths)):
        ax1.text(x, y, "{:.2f}".format(c), ha='center', va='center',
                color=text_color)
ax1.legend(ncol=len(category_names), bbox_to_anchor=(0, 1),
          loc='lower left', fontsize='small')


ax2.hist(delta_lifeyears[0,0,1].flatten(),bins=20,label="high income, female",density=True,alpha=0.5)
ax2.hist(delta_lifeyears[1,0,1].flatten(),bins=20,label="high income, male",density=True,alpha=0.5)
ax2.hist(delta_lifeyears[0,0,0].flatten(),bins=20,label="low income, female",density=True,alpha=0.5)
ax2.hist(delta_lifeyears[1,0,0].flatten(),bins=20,label="low income, male",density=True,alpha=0.5)
ax2.legend()
ax2.set_title("Qaly's lost (per 1000 population) due to $oop$ increase by $0.1$");

In order to compare the dynamic and static approach, use the ages 45-74 (first cohort at the start of our data; see the first row in Table 1). Then we compare this with the dynamic estimate delta_lifeyears over the same age-range (that is, not the range 45-84 that we consider above).

The following code calculates the dynamic loss in life years from 45 to 74 and returns the data in the first row of Table 3:

lifeyears_0_d = N_0[0,:30].sum(axis=0) + qaly_ill*N_0[1,:30].sum(axis=0)
lifeyears_1_d = N_1[0,:30].sum(axis=0) + qaly_ill*N_1[1,:30].sum(axis=0)
delta_lifeyears_d = lifeyears_0_d - lifeyears_1_d
delta_lifeyears_d.mean(axis=(1,3,4))

array([[5.22867775, 0.45827571],
       [8.00927005, 0.82308271]])

To find the static estimate, we do the following. For \(\iota_0\) we are only interested in the first year; hence, index 0 in the slicing of ill to get ill_0. We will compare the qaly’s after 1 year of the introduction of \(\Delta oop\) with the counterfactual after 1 year without the increase in \(oop\). The latter is denoted lifeyears_0_s (with s for static) and the former lifeyears_1_s. The code returns the second row of Table 3 (and repeats the first row for ease of reference).

ill_0 = idata_0.predictions.ill.values[:,:,0,:,:,:,:].transpose((2,3,4,5,0,1))*1000
fit_0 = 1000-ill_0
fit_1_0 = (1-pi[:30])*fit_0 + sigma_0 * ill_0
ill_1_0 = pi[:30]*fit_0 + (1-sigma_0-delta) * ill_0
deaths_1_0 = delta * ill_1_0
lifeyears_0_s = fit_1_0 + ill_1_0 * (1-delta) * qaly_ill
fit_1_1 = (1-pi[:30])*fit_0 + sigma_1 * ill_0
ill_1_1 = pi[:30]*fit_0 + (1-sigma_1-delta) * ill_0
deaths_1_1 = delta * ill_1_1
lifeyears_1_s = fit_1_1 + ill_1_1 * (1-delta) * qaly_ill
delta_lifeyears_s = (lifeyears_0_s - lifeyears_1_s).sum(axis=0)
print("The qaly effect based on a one year (static) comparison:")
print(delta_lifeyears_s.mean(axis=(1,3,4)))
print("The qaly effect based on a 30 year (dynamic) comparison:")
print(delta_lifeyears_d.mean(axis=(1,3,4)))

The qaly effect based on a one year (static) comparison:
[[0.24927238 0.05729696]
 [0.31104821 0.11820345]]
The qaly effect based on a 30 year (dynamic) comparison:
[[5.22867775 0.45827571]
 [8.00927005 0.82308271]]

For the extensions we use the variables in Table 4 in addition to the variables in Table 2. We have data on the number of people that die in traffic accidents. The unit is number of deaths per 1 million inhabitants of the region. We have data on GDP per capita (purchasing power standard) and the percentage of people unemployed in the region. We have the percentage of people (aged 25-64) with tertiary education in the region. We only have this from 2013 onward. Finally, we have data on the number of devices for medical imaging (MRI’s, PET scanners etc.) per 100k inhabitants at the country level.⁷

On average 66 people per 1 million inhabitants die in traffics accidents per year. This ranges from six in one region/year combination to a maximum of 254 deaths. GDP per capita is on average 24k but ranges per region from 6600 to almost 60k. Unemployment is on average almost 10% and varies from less than 2% to more than 30%. Almost one third of the working population has tertiary education on average. This varies from 16% to more than half of the population. Finally, the number of imaging devises is on average 40 per 100k inhabitants but varies from no devices in one year for a country to more than 200.

In the extension of the model, we adjust two parameters of the baseline model. First, in the equation for \(\pi\) we also include education and unemployment in the following way:

That is, we multiply the expression for \(\pi\) in the baseline model by a multiplier depending on education and unemployment. With more people having finished tertiary education one may expect more healthy lifestyle choices. Unemployment with the implication of lower income may lead to stress and/or worse lifestyle choices as fresh vegetables, fruit and, say, a gym membership are expensive. With the multiplier we capture that the effect of education and unemployment is small if the probability of falling ill is small to start with, say at age 45. At age 60 when the probability of falling ill is higher, the effect of unemployment and education is higher.

Second, in the equation for the probability of recovery conditional on treatment, \(\lambda\), we add GDP per capita and the number of imaging devices as additional variables capturing healthcare quality. These variables we add in the same way as infant mortality and number of physicians in the baseline specification of \(\lambda\).

With the estimated extended model we can do the same analysis as with the baseline model. Our main result is that the dynamic estimates are an order of magnitude higher than the static estimates. By comparing Tables 5 and 3 we see that this result is robust to the addition of the variables. The mean values of qaly’s lost for the different categories are similar for the base and extended models. To illustrate that the differences in average qaly’s lost are not “significant”, Figure 6 shows the posterior distributions for qaly’s lost for the female, low income category. Indeed, these distributions are almost the same for the base and extended models. The same is true for the other gender/income categories.

As a final check on our proposed causal mechanism, we add a placebo treatment. Instead of analyzing mortality linked to health, we analyze mortality due to traffic accidents. In the latter case there is neither the dynamic element of a health stock developing over time nor a role for out-of-pocket expenditure causing forgone treatments which may cause mortality. The details of the model are in the on-line appendix.

The idea is that if out-of-pocket expenditure has a clear positive effect on traffic mortality, an alternative causal mechanism presents itself: the government’s fiscal space has narrowed (e.g. due to a recession) which simultaneously causes the government to increase demand-side cost-sharing and reduce investments in healthcare. The reduction in these investments, reduces healthcare quality and this is the cause of increased mortality. An increase in traffic accidents is then caused by reduced investment in road safety.

We do not find such a clear relation between our \(oop\) variable and traffic mortality. In particular, we find a posterior mean for the effect of \(oop\) on traffic mortality of 0.35 with a (posterior) standard deviation of 0.60 which is almost twice as high. This suggests that the path from an increase in \(oop\) to forgone treatments because they are too expensive (for low incomes) to reduced health status and higher mortality is plausible.

When evaluating the effect of an increase in oop, we use \(\theta=0.6\) in the qaly calculation –see equation \eqref{org7f9eb9b}. Table 6 illustrates that our main results are not sensitive to this choice. The table reports qaly’s lost (due to \(\Delta oop =0.1\)) for \(\theta=0.3\) and \(0.8\). The middle rows reproduce Table 3 for ease of reference.

For each value of \(\theta\) we see that the effects are far bigger for low than for high incomes and the dynamic estimates are order of magnitude bigger than the static ones. Further, part of the effect on qaly’s lost is –next to higher mortality– that more people are in the low health state compared to the baseline situation with lower \(oop\). When higher weight \(\theta\) is given to the low health state, the effect of \(\Delta oop\) is lower: the qaly penalty for being in the low state is reduced. Indeed, all entries in the table fall with \(\theta\).

Although the size of the effect varies, the main results are robust to the choice of \(\theta\).

We extend the model with 4 variables: per capita GDP per NUTS 2 region, number of devices for medical imaging per 100k inhabitants at the country level, fraction of NUTS 2 inhabitants aged 25-64 who completed tertiary education and unemployment at the NUTS 2 level. The extended model is meant as a robustness check on the causal mechanism where OOP is increased, more people forgo treatment because it is too expensive and this has a negative effect on health.

An alternative interpretation of the negative effect of OOP on health is via the fiscal position of the government: as the government faces deficits, it increases OOP and reduces investment in healthcare quality. This reduction in quality is what is causing the drop in health not so much the forgone treatments. Unemployment and GDP per capita captures the fiscal position. The medical imaging devices (next to infant mortality and number of doctors) captures the quality of healthcare. Education captures knowledge about healthy lifestyles but also the wealth of the region. If more people with tertiary education immigrate to the region, its income will increase and the fiscal position will improve.

We model this as follows. We add education and unemployment as a multiplier to \(\pi\):

where the variables are normalized with expectation zero and standard deviation one. Hence, if education and unemployment are at their expected value, then the multiplier is equal to one. It is only when we variables are not equals to their expectation, the value of \(\pi\) is modified. If for, say, an age \(\pi\) is close to zero, then the multiplier has hardly an effect. Only when \(\pi\) is relatively big, does the multiplier have an impact. The variables GDP per capita and the number of imaging devices is added to the other “health quality variables” (infant mortality, number of physicians) in the expression for \(\lambda\).

Another advantage of standardizing the new variables is that when modelling the missing values they have mean zero and standard deviation one. Further, the variables are based on addition (number of unemployed, people with tertiary education and number of devices; sum of value added in the region) and the central limit theorem suggests a normal distribution for these variables.

def transition(π,σ,I,δ):
    fit = I*σ+(1-π)*(1-I)
    ill = π*(1-I)+pt.clip(1-σ-δ,0.0,1.0)*I
    fraction_ill = ill/(fit+ill)
    return fraction_ill

n_y = 2 # low/high income
n_g = 2 # female/male
countries = np.array(['Bulgaria', 'Switzerland', 'Denmark', 'Greece', 'Finland', 'Hungary', 'Lithuania', 'Norway', 'Sweden', 'Slovenia', 'Slovakia'])
# index: age, gender, year, region, income
# for some variables we use country instead of region

coords = {
    "region" : ds.nuts2.values,\
    "ages" : ds.age.values,\
    "cohort" : age_matrix[0],\
    "gender" : ds.sex.values,\
    "year" : ds.year.values,\
    "income" : ['y_l','y_h'],\
    "country": countries,\
     # we define dummy dimensions for parameters that do not
     # vary by this dimension
     "dummy0": np.arange(1),\
     "dummy1": np.arange(1),\
     "dummy2": np.arange(1)}

with pm.Model(coords=coords) as extension_model:
    logodds_ι = pm.Normal("logodds_ι",-4.0,1,\
              dims=("cohort","gender","region","income"),\
              shape=(cohorts,n_g,n_regions,n_y)) 
    ι_0 = pm.Deterministic("ι_0",pm.math.invlogit(logodds_ι),dims=("cohort","gender","region","income"))
    # ι_0 is used to start the sequence in scan =>
    # no need to add age dimension: we build up the matrix I by following each cohort over the years
    logodds_δ = pm.Normal("logodds_δ",-2.5,0.05,dims=("country")) # conditional on being ill, same prob. of death
    # δ is a non_sequence in the scan loop; no need to add age as a dimension
    δ = pm.Deterministic("δ",pt.dot(pm.math.invlogit(logodds_δ), country_region_matrix_01)[None,None,None,:,None])
    # new variables for the extension
    gdp = pm.Normal('gdp',observed=GDP,dims=("year","region"))
    devices = pm.Normal('devices',observed=Devices,dims=("year","region"))
    education = pm.Normal('education',observed=Education,dims=("year","region"))
    unemployment = pm.Normal('unemployment',observed=Unemployment,dims=("year","region"))
        
    # we build up the prior for π
    ## this forms the multiplier used in the scan-method below
    a1 = pm.Normal("a1",0,0.01)
    a2 = pm.Normal("a2",0,0.01)
    mult_pi = np.exp(a1 * education + a2 * unemployment)
    ## building π from the logodds
    logodds_age_π = pm.Normal("logodds_age_π",log_odds_mort.mean(axis=(1,2,3)),log_odds_mort.std(axis=(1,2,3)),dims=("ages"))
    logodds_sex_π = pm.Normal("logodds_sex_π", gender_effect_π, 0.05, dims=("gender"))
    logodds_country_π = pm.Normal("logodds_country_π",country_effect_π,0.05, dims=("country"))
    logodds_π = pm.Deterministic("logodds_π",logodds_age_π[:,None,None] + logodds_sex_π[None,:,None] +\
                                 pt.dot(logodds_country_π, country_region_matrix_01)[None,None,:])
    # the priors do not impose that women tend to be healthier than men
    # but we do impose that people on low income tend to have low health status (high π)
    delta_π = pm.HalfNormal("delta_π",sigma=0.05,\
          dims=("ages","dummy0","dummy1"),\
          shape=(n_ages,1,1))
    π = pm.Deterministic("π",pt.stack([pm.math.invlogit(\
                                    logodds_π+ delta_π),\
                   pm.math.invlogit(logodds_π)],axis=-1), dims=("ages","gender","region","income"))
    # when we draw missing values for OOP, we should draw the same value for all regions/age/gender in the same country
    # this we do via oop_c which we later bring back to region dimension with country_region_matrix_01
    pov = pm.Normal("pov", mean_depr[None,:], std_depr[None,:], observed=Deprivation,dims=("year","region")) # fraction of people on low income
    poverty = pt.stack([pt.clip(pov,0,1),1-pt.clip(pov,0,1)],axis=-1) # vector with fraction of low/high income people
    oop_c = pm.Normal("oop_c", mean_oop_c[None,:], std_oop_c[None,:], observed=OOP_c,dims=("year","country"))
    oop = pt.dot(oop_c,country_region_matrix_01)

    slope_ζ = pm.HalfNormal("slope_ζ",sigma=0.05)
    delta_ζ = pm.HalfNormal("delta_ζ",sigma=0.05)
    ζ = pm.Deterministic("ζ",pt.stack([(slope_ζ+ delta_ζ),\
                   (slope_ζ)],axis=-1),dims=("income")) # varies with income
    τ = pm.Deterministic("τ", ζ[None,:]*oop[:,:,None])
    mu_U_o = Unmet_o.mean(axis=0).values[None,:]
    nu_U_o = pm.Exponential("nu_U_o", 0.05)
    U_o = pm.Beta("U_o", alpha = mu_U_o*nu_U_o, beta = (1-mu_U_o)*nu_U_o,dims=("year","region"), observed=Unmet_o)
    mu_U = pt.clip(U_o+pt.sum(poverty*τ,axis=2),0.000001,0.9999)
    nu_U = pm.Exponential("nu_U", 0.05)
    U = pm.Beta("U",alpha=mu_U*nu_U, beta=(1-mu_U)*nu_U,dims=("year","region"), observed=Unmet)
    # λ measures prob. patient heals if treated (i.e. one minus unmet)
    M_infant = pm.Normal("M_infant",mean_m_infant[None,:],std_m_infant[None,:],dims=("year","region"),observed=m_infant)
    Physicians = pm.Normal("Physicians",mean_physician[None,:],std_physician[None,:],dims=("year","region"),observed=physicians)
    lambda_1 = pm.HalfNormal("lambda_1",sigma=0.01)
    lambda_2 = pm.HalfNormal("lambda_2",sigma=0.01)
    lambda_3 = pm.HalfNormal("lambda_3",sigma=0.01)
    lambda_4 = pm.HalfNormal("lambda_4",sigma=0.01)
    lamba_0 = pm.Normal("lambda_0",-4.1,0.05,dims=("year"))
    logodds_λ = pm.Normal("logodds_λ",lamba_0[:,None] - lambda_1 * M_infant + lambda_2 * Physicians + lambda_3 * gdp + lambda_4 * devices,0.1)
    λ = pm.Deterministic("λ",pm.math.invlogit(logodds_λ))# λ is between 0 and 1
    ε = pm.Beta("ε",1,3) # ε = ξ_l - ξ_h > 0
    ξ = 1-poverty[:,:,0]*ε
    υ = pm.Deterministic("υ", pt.stack([(ξ+ε)*U_o,ξ*U_o],axis=-1) + τ)
    σ_0 = pm.Beta("σ_0",1,3)
    σ = pm.Deterministic("σ",pt.clip(σ_0+λ[:,None,None,:,None]*(1-υ)[:,None,None,:,:],0,1))
    I = pt.zeros((t_dim,cohorts,n_g,n_regions,n_y))
    I = pt.set_subtensor(I[0],ι_0)

    outputs, updates = pytensor.scan(
        transition,
     # age/year is the dynamic part of the model
     # i.e. first dimension of π; not multiplied by the muliplier mult_pi defined above
        sequences=[dict(input=π[age_indices]*mult_pi[:,None,None,:,None]),dict(input=σ)], # check sequences with order of arguments in =transition=
     # we loop over the first dimension: age
     # the initial value for I (at low_age)
     # outputs_info has no age dimension
        outputs_info=[dict(initial=ι_0)],
     # variables that have no age dimension:
        non_sequences=[δ[0]],
    )

    I = pt.set_subtensor(I[1:],outputs[:t_dim-1])
    # μ is probability of death for age/gender/year/region
    ill = pm.Deterministic("ill",I)
    μ = pm.Deterministic("μ",(pt.sum(δ*I*poverty[:,None,None,:,:],axis=4)))# aggregate over income
    mortality = pm.Binomial("mortality",n=P_cohort,p=μ,observed=M_cohort,dims=("year","cohort","gender","region"))

Death to due traffic accidents does not have this stock nature. Also our eurostat data does not distinguish deaths by age or sex. Hence, we test whether OOP and Unmet needs because treatment is too expensive does correlate with traffic death. This would be evidence of the causal mechanism: government has a fiscal problem and increased OOP but also reduces road quality (e.g. roads are cleaned less frequently –also in winter– and are maintained less frequently)

• poverty is an effect on the roads and perhaps quality of the cars
• the question is whether OOP feeds through unmet needs into health; the other causal interpretation is that OOP is a measure of government means and affect quality of healthcare and quality of roads; hence we add OOP in the regression for the probability of traffic mortality
• nuts 2 fixed effects

Model Traffic_mortality as a Negative Binomial distribution:

• has be integer
• Poisson distribution given its parameter; parameter has Gamma distribution