Mapping PHNs in Australia using R

This program uses the esri shapefile for Primary Health Networks (PHNs) in Australia to generate a map and then the MORT data from the Australian Institute of Health and Welfare (AIHW) to show all-cause deaths and age-standardised rates by PHNs for persons, 2009–2013.

#Import shapefile from data directory using rgdal package
phn <- readOGR(dsn = "data", layer = "PHN_boundaries_AUS_Sep2015_V5")
## OGR data source with driver: ESRI Shapefile 
## Source: "data", layer: "PHN_boundaries_AUS_Sep2015_V5"
## with 31 features
## It has 6 fields
head(phn@data, n = 2)
##                       PHN_Name PHN_Code AREA_sqkm
## 0                     Adelaide   PHN401  1553.011
## 1 Australian Capital Territory   PHN801  2351.360
##                          STATE Shape_Leng Shape_Area
## 0              South Australia   3.433834  0.1540940
## 1 Australian Capital Territory   3.026809  0.2342085
#Generate a map of PHNs in Australia
plot(phn)

summary(phn)
## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min        max
## x  96.81694 159.109219
## y -43.74051  -9.142176
## Is projected: FALSE 
## proj4string : [+proj=longlat +ellps=GRS80 +no_defs]
## Data attributes:
##                                          PHN_Name     PHN_Code 
##  Adelaide                                    : 1   PHN101 : 1  
##  Australian Capital Territory                : 1   PHN102 : 1  
##  Brisbane North                              : 1   PHN103 : 1  
##  Brisbane South                              : 1   PHN104 : 1  
##  Central and Eastern Sydney                  : 1   PHN105 : 1  
##  Central Queensland, Wide Bay, Sunshine Coast: 1   PHN106 : 1  
##  (Other)                                     :25   (Other):25  
##    AREA_sqkm                                STATE      Shape_Leng     
##  Min.   :    626   New South Wales             :10   Min.   :  2.184  
##  1st Qu.:   3094   Queensland                  : 7   1st Qu.:  5.131  
##  Median :  32047   Victoria                    : 5   Median : 15.789  
##  Mean   : 244975   Western Australia           : 3   Mean   : 34.175  
##  3rd Qu.: 127530   South Australia             : 2   3rd Qu.: 43.465  
##  Max.   :2477561   Australian Capital Territory: 1   Max.   :220.409  
##                    (Other)                     : 3                    
##    Shape_Area       
##  Min.   :  0.06168  
##  1st Qu.:  0.31103  
##  Median :  3.05774  
##  Mean   : 22.43806  
##  3rd Qu.: 12.43801  
##  Max.   :226.84342  
## 
#Import MORT data for PHNs to merge with shapefile data
#Import data
mort_data <- read.csv("MORTTABLE2datagovau.csv", stringsAsFactors = FALSE)
#Put data into local data frame
mort_df <- tbl_df(mort_data)
# Subset data into PHNs
phn_mort <- filter(mort_df, category=="Primary Health Network (PHN)") %>% 
  select(PHN_Code=mort, category, geography, SEX, rank, cause_of_death, deaths, deaths_percent, age_standardised_rate)
#Subset data into all cause deaths for people in PHNs 2009-2013, excluding unknowns and total Aust
phn_mort_all <- filter(phn_mort, cause_of_death=="All causes", SEX=="Persons", 
                       age_standardised_rate != ".", PHN_Code != "PHN999")
#Make new variable asr and convert age_standardised_rate to numeric
phn_mort_all$asr <- as.numeric(phn_mort_all$age_standardised_rate)

#Subset data into top cause of death
phn_mort_top <- filter(phn_mort, rank==1)

#Use dplyr and left join to merge all cause mort data with phn map data using PHN_Code variable
phn@data <- left_join(phn@data, phn_mort_all)
## Joining, by = "PHN_Code"

Chart of the number of deaths by PHN, persons, 2009–2013

#Bar chart of deaths by PHN 2009-2013
ggplot (phn_mort_all , aes( x = reorder ( geography , -deaths ), y = deaths, fill = deaths )) + 
  geom_bar ( stat = "identity" ) + xlab("PHN") + 
  theme ( axis.text.x = element_text ( angle = 30 , hjust = 1 , vjust = 1 ))

Map of age-standardised all-cause mortality rates per 100,000 population by PHN, persons, 2009–2013

qtm(phn, "asr", title = "Age standardised death rates, Persons, by PHN, Australia, 2009-2013", fill.palette = "Purples")

Age-standardised all-cause mortality rates per 100,000 population facetted by PHN, persons, 2009–2013

#facet ASRs by PHN
tm_shape(phn) + tm_fill("asr", thres.poly = 0, palette = "Purples") +
  tm_facets("PHN_Name", free.coords=TRUE, drop.units=TRUE)

Chart of age-standardised all-cause mortality rates per 100,000 population by PHN, persons, 2009–2013

ggplot (phn_mort_all , aes( x = reorder ( geography , -asr ), y = asr, fill = asr )) + 
  geom_bar ( stat = "identity" ) + xlab("PHN") + 
  theme ( axis.text.x = element_text ( angle = 30 , hjust = 1 , vjust = 1 )) + 
  ggtitle ( "Age standardised death rates per 100,000 population by PHN, Australia 2009-2013" )