Introduction

This short document shows two things:

  1. How to use Rmarkdown to create an interactive html report.
  2. Calculates catchment areas for general practices based on LSOA - GP mapping.

It is inspired by this blog post, and based on this dataset from NHS Digital which provides a table of LSOA populations by sex for each general practice.

The basic workflow is:

The source code for the markdown file is available at…<>

Method

To complete this task we need a few R libraries.

Then import the data:

data <- read_csv("http://www.content.digital.nhs.uk/catalogue/PUB23139/gp-reg-patients-LSOA-alt-tall.csv") %>% janitor::clean_names()     ## The data is long format and ~ 37Mb
## Parsed with column specification:
## cols(
##   PRACTICE_CODE = col_character(),
##   NAME = col_character(),
##   LSOA_CODE = col_character(),
##   MALE_PATIENTS = col_integer(),
##   FEMALE_PATIENTS = col_integer(),
##   ALL_PATIENTS = col_integer()
## )
head(data)
## # A tibble: 6 × 6
##   practice_code                name lsoa_code male_patients
##           <chr>               <chr>     <chr>         <int>
## 1        A81001 THE DENSHAM SURGERY E01012187             0
## 2        A81001 THE DENSHAM SURGERY E01012188            42
## 3        A81001 THE DENSHAM SURGERY E01012189            22
## 4        A81001 THE DENSHAM SURGERY E01012190            63
## 5        A81001 THE DENSHAM SURGERY E01012191            71
## 6        A81001 THE DENSHAM SURGERY E01012192             0
## # ... with 2 more variables: female_patients <int>, all_patients <int>

This shows the first few rows of data in the format we need for further analysis.

Calculate catchments

## Calculate LSOA totals
lsoa_data <- data %>%
  group_by(lsoa_code) %>%
  summarise(m = sum(male_patients), 
            f = sum(female_patients), 
            total = sum(all_patients))

## Attach totals to dataset and calculate practice proportions
data <- data %>%
  left_join(lsoa_data) %>%
  mutate(percent_m = 100 * male_patients/m, 
         percent_f = 100 * female_patients/f, 
         percent_tot = 100 * all_patients/total)
## Joining, by = "lsoa_code"
## Calculate first past the post

fpp <- data %>%
  arrange(lsoa_code, -percent_tot) %>%
  select(lsoa_code,percent_tot, practice_code) %>%
  group_by(lsoa_code) %>%
  slice(1)

fpp 
## Source: local data frame [32,939 x 3]
## Groups: lsoa_code [32,939]
## 
##    lsoa_code percent_tot practice_code
##        <chr>       <dbl>         <chr>
## 1  E01000001    95.93220        F84640
## 2  E01000002    94.62617        F84640
## 3  E01000003    91.57695        F84640
## 4  E01000005    53.42363        F84081
## 5  E01000006    25.90543        F82647
## 6  E01000007    30.76058        Y02583
## 7  E01000008    38.28899        F82034
## 8  E01000009    27.82573        F82647
## 9  E01000010    56.87023        Y02583
## 10 E01000011    36.02484        F82625
## # ... with 32,929 more rows

Obtain LSOA boundary file

The easiest place to find boundary files is the ONS Geography Portal. We can obtain LSOA boundary files in GeoJSON format. There are options - we’ll use super-generalised boundaries for speed.

shape <- geojson_read("http://geoportal.statistics.gov.uk/datasets/da831f80764346889837c72508f046fa_3.geojson", what = "sp")
## using 'what = "sp" ensures the data is a SpatialPolygonsDataFrame
head(shape)
## class       : SpatialPolygonsDataFrame 
## features    : 6 
## extent      : -1.799004, 1.402579, 51.22011, 53.56475  (xmin, xmax, ymin, ymax)
## coord. ref. : +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0 
## variables   : 6
## names       : objectid,  lsoa11cd,     lsoa11nm,    lsoa11nmw, st_areashape, st_lengthshape 
## min values  :        1, E01000722, Bromley 023A, Bromley 023A,     280143.7,       2228.371 
## max values  :        6, E01030056, Swindon 007D, Swindon 007D,   18633214.6,      22361.255

This gives a spatial data frame to which we can link our data.

shape <- subset(shape, substr(lsoa11cd, 1, 1) == "E") ## restrict to English LSOAs only
shape@data <- left_join(shape@data, fpp, by=c("lsoa11cd" = "lsoa_code"))
## Warning in left_join_impl(x, y, by$x, by$y, suffix$x, suffix$y): joining
## character vector and factor, coercing into character vector

We can map the data

# shape
library(viridis)
palette <- rev(plasma(200))
credits <- "Contains ordnance survey data (c) \nCrown copyright and database right 2016"


t <- tm_shape(shape) + 
  tm_fill("practice_code", title = "Test", 
          palette = palette) +
  tm_credits( credits, size = 0.5, align = "right") +
  tm_borders(col = "grey", lwd = 0.2) +
  tm_layout(legend.outside = TRUE, 
            frame = FALSE) +
  tm_compass(position = c("left", "center")) + 
  tm_scale_bar(position = c("left", "center"))

tmap_mode("plot")
## tmap mode set to plotting
t

Making it interative

tmap_mode("view")
## tmap mode set to interactive viewing
t
## Warning: Credits not supported in view mode.
## Warning: Compass not supported in view mode.
## Warning: Scale bar not yet supported in view mode, it will be in the next
## version.