General Practitioners and population by Medicare Primary Health Network area

### start necessary libraries

# data import libraries
library(xlsx)
library(foreign)
library(rgdal)

library(tidyverse)
library(dplyr)

# reproducible research library
library(knitr)

# plotting and maps libraries
library(rmapshaper) # for simplifying maps
library(ggplot2)
library(leaflet)
library(ggmap)

### set figure save path to figure/
knitr::opts_chunk$set(fig.path = 'figure/')
knitr::opts_chunk$set(fig.width=12, fig.height=8)

Download and import data files

Population data in Australia 2009 to 2014, grouped by Statistical Area Level 3 (SA3) from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Demographic_Data

Primary Health Network (‘PHN’) map data from http://www.health.gov.au/internet/main/publishing.nsf/content/phn-boundaries

Concordance table between PHN geographical area and Statistical Area Level 3 from http://www.health.gov.au/internet/main/publishing.nsf/Content/PHN-Concordances

General Practice (‘GP’, also known as family medicine) workforce by PHN found from Australian Government Health Workforce Data website (http://hwd.health.gov.au/), direct links for General Practitioner data by PHN is not available without login and setting search terms. Pre-downloaded 2016 data.

# download data files if necessary, and extract if necessary
# minimal processing to rename variables and remove variables not required

# Primary Health Network map boundaires
url.phnboundaries <- 'http://www.health.gov.au/internet/main/publishing.nsf/Content/20FD74ABB14A1297CA257F150001FD3B/$File/PHN_boundaries_AUS_May2017_V7_Shapefile.zip'

if (!file.exists('PHN_boundaries_AUS_May2017_V7.shp')) {
  if (!file.exists('PHNboundaries.zip')) {
    if (download.file(url.phnboundaries,'PHNboundaries.zip')) {
      stop('Unable to download data file')
    }
    unzip('PHNboundaries.zip')
  }
}
oz.phn.map <- readOGR('PHN_boundaries_AUS_May2017_V7.shp')
## OGR data source with driver: ESRI Shapefile 
## Source: "PHN_boundaries_AUS_May2017_V7.shp", layer: "PHN_boundaries_AUS_May2017_V7"
## with 31 features
## It has 8 fields
## Integer64 fields read as strings:  OBJECTID
# map data includes the list of PHNs in a DBF
oz.phn.description <- read.dbf('PHN_boundaries_AUS_May2017_V7.dbf')

# Populations, in Statistical Area 3 (SA3) areas
url.populations <- 'http://www.health.gov.au/internet/main/publishing.nsf/Content/AEADA0341B987748CA257F15000413F6/$File/PopulationReport.xlsx'
if (!file.exists('PopulationReport.xlsx')) {
  if (download.file(url.populations,'PopulationReport.xlsx', method = 'curl')) {
    # method 'curl' required for Windows
    stop('Unable to download data file')
  }
}
population.2014 <- as.tibble(read.xlsx('PopulationReport.xlsx', sheetName = 'Pop2014', startRow = 9, endRow = 359, header = FALSE))
# variable names needs to be added back into population table
names(population.2014) <- c('Region','Region2','Nothing','AllAges','Age0_4','Age5_9','Age10_14','Age15_19','Age20_24','Age25_29','Age30_34','Age35_39','Age40_44','Age45_49','Age50_54','Age55_59','Age60_64','Age65_69','Age70_74','Age75_79','Age80_84','Age85plus')

population.2014 <- population.2014 %>%
  select(-Nothing,-Region) %>%
  rename(SA3_NAME_2011 = Region2)
# the original 'Region' variable had leading spaces/tabs

url.concordance <- 'http://www.health.gov.au/internet/main/publishing.nsf/Content/4FAAB33C3ED3E520CA257F150001FD3A/$File/CG_SA3_2011_PHN_2015.xls'
# The correlation between Primary Health Network (PHN) zones and Statistical Area 3 (SA3) areas
# There may be multiple SA3s in a PHN, and a SA3 might span over the 'boundary' between PHNs
# this is the 2011 SA3 names, which is required for 2009 to 2014 population data
# 
if (!file.exists('CGSA3_2011PHN2015.xlsx')) {
  if (download.file(url.concordance,'CGSA3_2011PHN2015.xlsx', method = 'curl')) {
    stop('Unable to download data file')
  }
}
sa3.phn.concordance <- as.tibble(read.xlsx('CGSA3_2011PHN2015.xlsx',
                                           sheetName = 'Table 3', startRow = 8, header = FALSE))
names(sa3.phn.concordance) <-  c('SA3_CODE_2011','SA3_NAME_2011','PHN_CODE','PHN_NAME','RATIO','PERCENTAGE')

gp.numbers <- as.tibble(read.csv('table_2018-03-12_14-41-09.csv', skip = 11, row.names = NULL))
# number of general practitioners 'GP' in each primary health network (PHN) 2016
names(gp.numbers) <- c('PHN_NAME','gp','X') # rename columns
gp.numbers <- gp.numbers %>%
  select(-X) # remove column not needed (and empty)
gp.numbers[gp.numbers$PHN_NAME == 'Central Queensland','PHN_NAME'] <-
  'Central Queensland, Wide Bay, Sunshine Coast'

Place information into Map file

Places Population, General Practitioner numbers and Population per General practitioner information into map file

Creates a data frame with

# choose all ages statistic
select.age.population <- select(population.2014, SA3_NAME_2011, AllAges)

# calculate population in each PHN
phn.population <- sa3.phn.concordance %>%
  merge(select.age.population, by = 'SA3_NAME_2011') %>%
  # merge population numbers by SA3 areas
  mutate(PHN_Number = as.numeric(AllAges)*as.numeric(RATIO)) %>%
  # calculate population in each SA3 as proportion in each PHN
  # some SA3 are spread across more than one PHN
  group_by(PHN_NAME) %>%
  # PHN can have multiple SA3
  summarize(population = sum(PHN_Number)) 

# add GP information to each PHN
# calculate population per GP ratio in each PHN
combined.oz.phn.description <- oz.phn.description %>%
  merge(gp.numbers, by = 'PHN_NAME') %>%
  merge(phn.population, by = 'PHN_NAME') %>%
  mutate(pop.gp.ratio = population/gp) %>%
  # just select the variables that we need
  select(PHN_NAME, gp, population, pop.gp.ratio)

add extra map data to the map

oz.phn.map@data <-  merge(oz.phn.map@data, combined.oz.phn.description, by = 'PHN_NAME')

Population and General Practitioner numbers in Australia, by Primary Health Network areas

‘Click’ on region to reveal Population and GP numbers

oz_loc <- geocode('Australia', source='dsk')
## Information from URL : http://www.datasciencetoolkit.org/maps/api/geocode/json?address=Australia&sensor=false
phn.map <- ms_simplify(oz.phn.map) %>% # reduce points to about 5% of original
  leaflet() %>% # the map generator
  addTiles() %>% # base map, usually OpenStreetMap
  setView(lng = oz_loc$lon, lat = oz_loc$lat, zoom = 4) %>% # set view
  addPolygons(color="#444444", weight=1, smoothFactor=0.5,
              opacity=0.8, fillOpacity = 0.15,
              popup = ~paste(PHN_NAME, '<br>Population :',
                             sprintf('%.0f',population),
                             '<br>GPs :', gp,
                             '<br>Population to GP :',
                             sprintf('%.0f',pop.gp.ratio)),
              fillColor = ~colorQuantile("YlOrRd", pop.gp.ratio)(pop.gp.ratio),
              highlightOptions =
                highlightOptions(color='green',weight=2,bringToFront = TRUE))

phn.map

‘Click’ on region to reveal Population and GP numbers