Simplified code to produce basic analyses of health spending within OECD

This is a quickly pared down set of a much longer project I started working several months ago (before covid-19). It’s not particularly pretty, but I wanted to make this relatively succinct and easy to follow.

requirements and initialization

library(tidyverse) # dplyr, ggplot, etc
library(janitor)   # column renaming
library(reshape2)  # restructuring data via melt() and dcast()
library(rsdmx)     # API to fetch data from OECD via SDMX (XML)
library(ggpmisc)   # printing formula and fit statistics on plots 
library(ggrepel)   # smarter plot labels
library(hrbrthemes) # theme for ggplot
require(Hmisc)      # pretty printing


fetchSDMXUrl=function(url)
{
  as.data.frame(readSDMX(url))
}


exclude_usFn=function(df) {
  filter(df,location!='USA')
}


simpleData=list()

Make “multiplier” table to convert values expressed in current LCU to various other measures, especially PPPs.

I do this by fetching transactions corresponding to GDP and AIC from OECD’s Main SNA table expressed with measures for both in:

1. Current Dollars (C)
2. constant PPPs, per capita (HVPVOB)
3. current PPPs per capita (HCPC)

The ratios between these measures, specifically Current LCUs to the PPPs define the multiplier to be used for each type.

OECD also publishes PPPs for AIC and GDP, but only in current dollar terms, insofar as I am aware. This method makes it much easier.

getPPPMultipliers=function(url) {
  


  rawData=fetchSDMXUrl(url) %>%
    clean_names() %>%
    filter(
      time_format=='P1Y', # annual
      transact %in% c('B1_GE','P41'), # GDP (expenditure method) and AIC
      measure %in% c('C','HVPVOB','HCPC') # current dollars (nominal LCU) and per capita PPP adjusted in constant and current terms
    ) %>%
    transmute(
      location,
      year=as.integer(obs_time),
      transact,
      measure,
      value=obs_value
    )  


  rawMult=rawData %>%
    group_by(location,year,transact) %>%
    mutate(
      # extract the figure expressed in current LCU for each location, year, and transaction set.
      ref_value=mean(ifelse(measure=='C',value,NA),na.rm=T) 
    ) %>%
    ungroup() %>%
    transmute(
      location,
      year,
      multiplier=value/ref_value,
      multiplier_type=case_when(
        measure=='C' ~ 'CURRENT_LCU_MILLIONS',
        T ~ sprintf(
          "%s_%s_PER_CAPITA",
          ifelse(transact=='P41','AIC','GDP'),
          ifelse(measure=='HVPVOB','CONSTANT','CURRENT')
        )
      )
    ) %>%
    unique() # we dont need multiple results for Current/LCU!
  
    
  return(rawMult)

}


url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SNA_TABLE1/AUS+AUT+BEL+CAN+CHL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+EA19+EU28+EU15+OECDE+OECD+OTF+NMEC+ARG+BRA+BGR+CPV+CHN+COL+CRI+HRV+CYP+HKG+IND+IDN+MDG+MLT+MAR+MKD+PER+ROU+RUS+SAU+SRB+ZAF+ZMB+FRME+DEW.B1_GE+P41.C+HCPC+HVPVOB/all?startTime=1970&endTime=2019'

simpleData$ppp_conversion=getPPPMultipliers(url)

Get population figures from OECD

getPops=function(url) {
  
  fetchSDMXUrl(url) %>%
    clean_names() %>%
    filter( # not necessary to filter this dataset, but just to be clear....
      time_format == 'P1Y', 
      transact == 'POPNC'
    ) %>% 
    transmute(
      location,
      year=as.integer(obs_time),
      pop=obs_value # thousands
    )
  
}

url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SNA_TABLE3/AUS+AUT+BEL+CAN+CHL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+EA19+EU28+EU15+OECDE+OECD+OTF+NMEC+ARG+BRA+BGR+CPV+CHN+COL+CRI+HRV+CYP+IND+IDN+MDG+MLT+MAR+MKD+PER+ROU+RUS+SAU+SRB+ZAF+ZMB+FRME+DEW.POPNC.PER/all?startTime=1970&endTime=2018'
simpleData$pop=getPops(url)

Get transactions from detailed non-financial accounts (SNA)

Use OECD’s SNA table for non-financial accounts by sector to obtain more granular SNA data for different sectors and transaction types. Although I use this for several other purposes (hence the extra transactions and sectors), this provides the transactions in LCU for gross adjusted disposable income and consumption of fixed capital for the household sector.

HH gross adjusted disposable income less their consumption of fixed capital equals net adjusted household disposable income.

getSectoral=function(url,priceCW) {
  sectoral_raw=fetchSDMXUrl(url) %>%
    clean_names()
    
    nfa_sector_cw=tribble(
      ~sector,~sec_short,~sec_long,
      'S1','TOT','Total economy',
      'S11','NFC','Non-Financial corporations',
      'S12','FIN','Financial corporations',
      'S13','GOV','Geneal Government',
      'S14_S15','HH','Households and non-profit institutions serving households'
    )
  
  
  nfa_transact_cw=tribble(
    ~transact,~tx_short,~tx_long,
    'NFB1GR','GVA','Gross Domestic Product / Gross Value Add',
    'NFB5GR','GPI','Gross National Income / Balance of Primary Income',
    'NFB6GP','GDI','Gross Disposable Income',
    'NFB7GP','AGDI','Adjusted Gross Disposable Income',
    'NFB8GP','GS','Gross Saving',
    'NFK1P','CFC','Consumption of Fixed Capital',  
    'NFP5P','GCF','Gross Capital Formation',
    'NFB9P','NL','Net Lending (+), Net Borrowing (-)'
  )
  
  sectoral_long=sectoral_raw %>%
    merge(nfa_sector_cw,by='sector') %>%
    merge(nfa_transact_cw,by='transact') %>%
    transmute(
      location,
      year=as.integer(obs_time),
      sector=sec_short,
      transaction=tx_short,
      value=obs_value
    ) %>%
    merge(priceCW,by=c('location','year')) %>%
    mutate(
      raw_value=value,
      value=raw_value*multiplier
    )
  
  
  # subtract consumption of fixed capital from disposable income, standard and "adjusted" (+STiK)
  # to derive net * disposable income
  sectoral_net_disposable=sectoral_long %>%
    filter(
      transaction %in% c('CFC','GDI','AGDI')
    ) %>%
    dcast(location+year+sector+multiplier_type~transaction) %>%
    mutate(
      ANDI=AGDI-CFC,
      NDI=GDI-CFC
    ) %>%
    melt(measure.vars=c('ANDI','NDI')) %>%
    mutate(
      transaction=as.character(variable)
    ) %>%
    select(-variable)
  
  sectoral_wide=sectoral_long %>%
    bind_rows(sectoral_net_disposable) %>%
    mutate(
      variable=str_to_lower(sprintf("nfa_%s_%s",transaction,sector))
    ) %>%
    dcast(location+year+multiplier_type~variable)
  
  return(sectoral_wide)

}

url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SNA_TABLE14A/AUS+AUT+BEL+CAN+CHL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+EA19+EU28+NMEC+BRA+BGR+CPV+CHN+COL+CRI+HRV+CYP+HKG+IND+IDN+MDG+MLT+MAR+MKD+PER+ROU+RUS+SAU+SRB+ZAF+ZMB+FRME+DEW.NFB5GR+NFB7GP+NFB6GP+NFB8GP+NFP5P+NFK1P+NFB9P+NFB1GR.S1+S11+S12+S13+S14_S15.C/all?startTime=1970&endTime=2018'


simpleData$sectoral=getSectoral(url,simpleData$ppp_conversion)

Get current health expenditures from OECD’s SHA tables

Get SHA health expenditures from OECD. I am fetching this in current LCU and converting them locally to be consistent with the other figures used. OECD does provide measures in per capita AIC PPPs, constant and current. I’ve found the results to be highly comparable, with some rare exceptions, however this ensures the population and PPP values are consistent throughout my analysis. Using their PPP-adjusted figures instead will not substantively change my analysis.

getSHAExpend=function(url,priceCW) {
  
  hce=fetchSDMXUrl(url) %>%
    clean_names() %>%
    filter(
      hf=='HFTOT',
      hc=='HCTOT',
      hp=='HPTOT'
    ) %>%
    transmute(
        location,
        year=as.integer(obs_time),
        variable=sprintf("sha_hce_%s",measure),
        value=obs_value
    ) %>%
    dcast(location+year~variable) %>%
    merge(priceCW,by=c('location','year')) %>%
    mutate(
      sha_hce_adj=sha_hce_MLLNCU*multiplier
    )
  
  return(hce)

}


url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SHA/HFTOT.HCTOT.HPTOT.MLLNCU.AUS+AUT+BEL+CAN+CHL+COL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+BRA+CHN+CRI+IND+IDN+RUS+ZAF/all?startTime=1970&endTime=2019'


simpleData$hce=getSHAExpend(url,simpleData$ppp_conversion)

Get GDP, AIC, and some other NA aggregates for analysis from OECD’s SNA Table 1

Fetch several transactions corresponding to more widely understood SNA aggregates in Current LCU from OECD’s Main SNA table.

getSNA1Basic=function(url,priceCW) {
  
  gdp_disagg_cw=tribble(
    ~transact,~short_lab,~name,
    'B1_GE','gdp','GDP (expenditure approach)',
    'P41','aic','Actual Individual Consumption', 
    'P3','final_consumption','Final Consumption',
    'P5','gross_capital_formation','Gross Capital Formation',
    'B11','net_exports','Net Exports',
    'P31S13','gov_consumption_ind','Government Consumption: Individual',
    'P32S13','gov_consumption_col','Government Consumption: Collective',
    'P3S13','gov_consumption_all','Government Consumption (combined)',
    'P31S14_S15','hh_consumption','Household and NPISH consumption',
    'P6','exports','Exports',
    'P7','imports','Imports',
    'P51','gross_fixed_capital_formation','Gross Fixed Capital Formation',
    'P52_P53','changes_in_inventories','Changes in inventories (stocks)',
    'P3_P5','domestic_demand','Domestic Demand'
  )
  
  gdp_disagg_raw=fetchSDMXUrl(url) %>%
    clean_names() %>%
    transmute(
      location,
      year=as.integer(obs_time),
      transact,
      value=obs_value
    )

  gdp_disagg_long=merge(gdp_disagg_raw,gdp_disagg_cw,by='transact',all.x=T) %>%
    merge(priceCW,by=c('location','year')) %>%
    mutate(
      raw_value=value,
      value=raw_value*multiplier
    )
  
  
  gdp_disagg_wide= gdp_disagg_long %>%
    mutate(
      variable=sprintf("sna1_%s",short_lab)
    ) %>%
    dcast(location+year+multiplier_type~variable)  
  
  return(gdp_disagg_wide)  
}


url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SNA_TABLE1/AUS+AUT+BEL+CAN+CHL+COL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+EA19+EU28+EU15+OECDE+OECD+OTF+NMEC+ALB+ARG+BRA+BGR+CPV+CHN+CRI+HRV+CYP+HKG+IND+IDN+MDG+MLT+MAR+MKD+PER+ROU+RUS+SAU+SRB+SGP+ZAF+ZMB+FRME+DEW.B1_GE+P3_P5+P3+P31S14_S15+P3S13+P41.C/all?startTime=1970&endTime=2019'

simpleData$key_aggregates=getSNA1Basic(url,simpleData$ppp_conversion)

bring the basic data together

simpleData$merged=simpleData$hce %>%
  merge(simpleData$sectoral,by=c('location','year','multiplier_type'),all=T) %>%
  merge(simpleData$key_aggregates,by=c('location','year','multiplier_type'),all=T) %>%
  merge(simpleData$pop,by=c('location','year'),all=T)

simpleData$often_referenced=simpleData$merged %>%
  transmute(
    year,
    location,
    multiplier_type,
    pop,         # population
    sha_hce_adj, # current health expenditures (SHA)
    nfa_andi_hh, # net adjusted household disposable income
    nfa_agdi_hh, # gross adjusted household disposable income
    sna1_gdp,    # GDP
    sna1_aic     # Actual Individual Consumption
  ) %>%
  filter( 
    !is.na(multiplier_type)  
  )




# available PPP-related conversions
cat(unique(simpleData$often_referenced$multiplier_type))

GDP_CONSTANT_PER_CAPITA GDP_CURRENT_PER_CAPITA CURRENT_LCU_MILLIONS AIC_CONSTANT_PER_CAPITA AIC_CURRENT_PER_CAPITA

You will want to be sure to filter these data by the appropriate field in multiplier_type!!

sample usage, number one

exclude_usFn=function(df) {
  filter(df,location!='USA')
}

simpleData$often_referenced %>%
  filter(
    year==2017,
    !is.na(nfa_andi_hh),
    !is.na(sha_hce_adj),
    
    # filter by the desired multiplier or you will get weird results!!!
    multiplier_type=='AIC_CONSTANT_PER_CAPITA', 
    
    # exclude Luxembourg and South Africa
    !location %in% c('ZAF','LUX'), 
    T
  ) %>%
  ggplot(aes(log(nfa_andi_hh),log(sha_hce_adj) )) +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    scale_y_continuous(
      sec.axis = sec_axis(
        ~exp(.),
        labels=dollar_format(),
        breaks=exp(seq(1,20,.5))
      )
    ) +
    scale_x_continuous(
      sec.axis = sec_axis(
        ~exp(.),
        labels=dollar_format(),
        breaks=exp(seq(1,20,.2))
      )
    ) +

    geom_text_repel(aes(label=location)) +
    theme_ipsum() +
    geom_smooth(
      method=lm,
      data=exclude_usFn,
      fullrange=T
    ) +
    labs(
      x='log 2017 net adjusted household disposable income per capita\nat PPPs for AIC (constant dollars)',
      y='log 2017 current heatlh expendiures per capita\nat PPPs for AIC (constant dollars)'
      
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE,
      data=exclude_usFn
    ) 

NA
NA
NA

sample usage, number two

simpleData$often_referenced %>%
  filter(
    year==2017,
    !is.na(nfa_andi_hh),
    !is.na(sna1_aic),
    multiplier_type=='AIC_CONSTANT_PER_CAPITA',
    T
  ) %>%
  ggplot(aes(log(nfa_andi_hh),log(sna1_aic) )) +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    geom_abline(slope=1) +
    scale_y_continuous(
      sec.axis = sec_axis(
        ~exp(.),
        labels=dollar_format(accuracy=1),
        breaks=exp(seq(1,20,.5))
      )
    ) +
    scale_x_continuous(
      sec.axis = sec_axis(
        ~exp(.),
        labels=dollar_format(accuracy=1),
        breaks=exp(seq(1,20,.2))
      )
    ) +

    geom_text_repel(aes(label=location)) +
    theme_ipsum() +
    geom_smooth(
      method=lm,
      data=exclude_usFn,
      fullrange=T
    ) +
    labs(
      x='log 2017 net adjusted household disposable income per capita\nat PPPs for AIC (constant dollars)',
      y='log 2017 AIC per capita\nat PPPs for AIC (constant dollars)'
      
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE,
      data=exclude_usFn
    ) 

NA
NA
NA

compare to BLI Net AHDI figures

Extract Net Adjusted Household Disposable Income per capita from the Better Life Index and compare. This figure is presumably for 2016 and in current PPPs. I’d hope it’d be closer (last time I looked it was!). Still, I think it’s close enough. Statistical discrepancies can arise because of revisions to the underlying data set. It’s more likely the BLI data are a bit outdated than sources I am pulling this.

url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/BLI/AUS+AUT+BEL+CAN+CHL+COL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+OECD+NMEC+BRA+RUS+ZAF.HO_BASE+HO_HISH+HO_NUMR+IW_HADI+IW_HNFW.L.TOT/all?'

bliData=fetchSDMXUrl(url) %>%
    clean_names() %>%
    dcast(location~indicator,value.var='obs_value') %>%
    mutate(
      year=2016 
    )



simpleData$often_referenced %>%
  filter(
    multiplier_type == 'AIC_CURRENT_PER_CAPITA',
    T
  ) %>%
  merge(bliData,by=c('location','year')) %>%
  filter(
    !is.na(IW_HADI),
    !is.na(nfa_andi_hh)
  ) %>%
  ggplot(aes(log(nfa_andi_hh),log(IW_HADI) )) +
    #geom_abline(slope=1) +
    geom_point(color='red',aes(size=pop),alpha=0.5) +
    geom_text_repel(aes(label=location)) +  
    geom_smooth(
      method=lm,
      fullrange=T
    ) +
    labs(
      x='log 2016 net adjusted household disposable income per capita\nat PPPs for AIC (current dollars)\nper my calculations',
      y='log 2016 net adjusted household disposable income per capita\nat PPPs for AIC (current dollars)\nper OECD\'s 2019 Better Life Index'
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE
    ) +
    theme_ipsum()

NA

crosscheck with OECD NAAG

Use estimates published in OECD’s National Accounts At A Glance for gross adjusted household disposable income in current prices, current PPPs (AIC) to cross-check my figures. These correspond exactly.

crosscheck against GDP and AIC

These match exactly (+/- negligible floating point precision error)

url='https://stats.oecd.org/restsdmx/sdmx.ashx/GetData/SNA_TABLE1/AUS+AUT+BEL+CAN+CHL+COL+CZE+DNK+EST+FIN+FRA+DEU+GRC+HUN+ISL+IRL+ISR+ITA+JPN+KOR+LVA+LTU+LUX+MEX+NLD+NZL+NOR+POL+PRT+SVK+SVN+ESP+SWE+CHE+TUR+GBR+USA+EA19+EU28+EU15+OECDE+OECD+OTF+NMEC+ALB+ARG+BRA+BGR+CPV+CHN+CRI+HRV+CYP+HKG+IND+IDN+MDG+MLT+MAR+MKD+PER+ROU+RUS+SAU+SRB+SGP+ZAF+ZMB+FRME+DEW.B1_GE+P41.HCPC+HVPVOB/all?startTime=1970&endTime=2019'

basicAggs=fetchSDMXUrl(url) %>%
  clean_names() %>%
  mutate(
    year=as.integer(obs_time),
    variable=sprintf("%s_%s",transact,measure)
  ) %>%
  dcast(location+year~variable,value.var='obs_value')




basicAggs %>%
  merge(simpleData$often_referenced,by=c('location','year')) %>%
  filter(
    year %in% c(2005,2010,2015,2018),
    multiplier_type=='AIC_CONSTANT_PER_CAPITA'
  ) %>%
  ggplot(aes(log(sna1_aic),log(P41_HVPVOB))) +
    geom_abline(slope=1) +
    facet_wrap(~year,scale='free_y') +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    geom_smooth(
      method=lm,
      fullrange=T
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE
    ) +
    theme_ipsum() +
    guides(size='none') +
    labs(
      x='log AIC per capita at constant PPPs for AIC (my side)',
      y='log AIC per capita at constant PPPs for AIC (OECD side)'
    )

basicAggs %>%
  merge(simpleData$often_referenced,by=c('location','year')) %>%
  filter(
    year %in% c(2005,2010,2015,2018),
    multiplier_type=='AIC_CURRENT_PER_CAPITA'
  ) %>%
  ggplot(aes(log(sna1_aic),log(P41_HCPC))) +
    geom_abline(slope=1) +
    facet_wrap(~year,scale='free_y') +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    geom_smooth(
      method=lm,
      fullrange=T
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE
    ) +
    theme_ipsum() +
    guides(size='none') +
    labs(
      x='log AIC per capita at current PPPs for AIC (my side)',
      y='log AIC per capita at current PPPs for AIC (OECD side)'
    )

basicAggs %>%
  merge(simpleData$often_referenced,by=c('location','year')) %>%
  filter(
    year %in% c(2005,2010,2015,2018),
    multiplier_type=='GDP_CONSTANT_PER_CAPITA'
  ) %>%
  ggplot(aes(log(sna1_gdp),log(B1_GE_HVPVOB))) +
    geom_abline(slope=1) +
    facet_wrap(~year,scale='free_y') +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    geom_smooth(
      method=lm,
      fullrange=T
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE
    ) +
    theme_ipsum() +
    guides(size='none') +
    labs(
      x='log GDP per capita at constant PPPs for GDP (my side)',
      y='log GDP per capita at constant PPPs for GDP (OECD side)'
    )

basicAggs %>%
  merge(simpleData$often_referenced,by=c('location','year')) %>%
  filter(
    year %in% c(2005,2010,2015,2018),
    multiplier_type=='GDP_CURRENT_PER_CAPITA'
  ) %>%
  ggplot(aes(log(sna1_gdp),log(B1_GE_HCPC))) +
    geom_abline(slope=1) +
    facet_wrap(~year,scale='free_y') +
    geom_point(color='red',aes(size=pop),alpha=0.7) +
    geom_smooth(
      method=lm,
      fullrange=T
    ) +
    stat_poly_eq(
      aes(label =  paste(stat(eq.label),stat(rr.label), sep = "*\", \"*")),
      formula = 'y~x',
      parse = TRUE
    ) +
    theme_ipsum() +
    guides(size='none') +
    labs(
      x='log GDP per capita at current PPPs for GDP (my side)',
      y='log GDP per capita at current PPPs for GDP (OECD side)'
    )