Take-home Exercise 1: Geographic Analysis of the Supply and Demand of Childcare Services in Singapore

1 Introduction

1.1 Objectives

• To gain understanding on the supply and demand of childcare services in 2017 and 2020 at planning subzone level.
• To gain understanding on the geographic distribution of childcare in 2017 and 2020 in Singapore.

2 Input

2.1 Loading libraries

packages = c('tmap', 'tidyverse', 'sf', 'rgdal','spatstat', 'raster', 'maptools', "rgeos", "dplyr", "spatstat", "OpenStreetMap", "tmaptools")
for (p in packages){
  if(!require(p, character.only = T)){
    install.packages(p)
  }
  library(p,character.only = T)
}

Read data

• The childcare services data in 2017. (You can use the version provided in hands-on exercise)
• The latest childcare data set from data.gov.sg.
• URA Master Plan 2014 Planning Subzone GIS data. This data set is also available at data.gov.sg.
• Singapore Residents by Planning AreaSubzone, Age Group, Sex and Type of Dwelling, June 2011-2019. This data set is available at Singapore Department of Statistics.

## Reading layer `CHILDCARE' from data source `/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1312 features and 18 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
## projected CRS:  SVY21

## Reading layer `CHILDCARE_2020' from data source `/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 1545 features and 25 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 11203.01 ymin: 25667.6 xmax: 45404.24 ymax: 49300.88
## z_range:        zmin: 0 zmax: 0
## projected CRS:  SVY21 / Singapore TM

## Reading layer `MP14_SUBZONE_WEB_PL' from data source `/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial' using driver `ESRI Shapefile'
## Simple feature collection with 323 features and 15 fields
## geometry type:  MULTIPOLYGON
## dimension:      XY
## bbox:           xmin: 2667.538 ymin: 15748.72 xmax: 56396.44 ymax: 50256.33
## projected CRS:  SVY21

2.2 Clean and format data

2.2.1 Geospatial data

2.2.1.1 Checking data

st_crs(sf_childcare_2017)

## Coordinate Reference System:
##   User input: SVY21 
##   wkt:
## PROJCRS["SVY21",
##     BASEGEOGCRS["SVY21[WGS84]",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

st_crs(sf_childcare_2020)

## Coordinate Reference System:
##   User input: SVY21 / Singapore TM 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

st_crs(mpsz)

## Coordinate Reference System:
##   User input: SVY21 
##   wkt:
## PROJCRS["SVY21",
##     BASEGEOGCRS["SVY21[WGS84]",
##         DATUM["World Geodetic System 1984",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]],
##             ID["EPSG",6326]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["Degree",0.0174532925199433]]],
##     CONVERSION["unnamed",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["Degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["(E)",east,
##             ORDER[1],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]],
##         AXIS["(N)",north,
##             ORDER[2],
##             LENGTHUNIT["metre",1,
##                 ID["EPSG",9001]]]]

2.2.1.2 Format data

sf_childcare_2017 <- st_transform(sf_childcare_2017, 3414)
mpsz <- st_transform(mpsz, 3414)

Check again

st_crs(sf_childcare_2017)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

st_crs(sf_childcare_2020)

## Coordinate Reference System:
##   User input: SVY21 / Singapore TM 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

st_crs(mpsz)

## Coordinate Reference System:
##   User input: EPSG:3414 
##   wkt:
## PROJCRS["SVY21 / Singapore TM",
##     BASEGEOGCRS["SVY21",
##         DATUM["SVY21",
##             ELLIPSOID["WGS 84",6378137,298.257223563,
##                 LENGTHUNIT["metre",1]]],
##         PRIMEM["Greenwich",0,
##             ANGLEUNIT["degree",0.0174532925199433]],
##         ID["EPSG",4757]],
##     CONVERSION["Singapore Transverse Mercator",
##         METHOD["Transverse Mercator",
##             ID["EPSG",9807]],
##         PARAMETER["Latitude of natural origin",1.36666666666667,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8801]],
##         PARAMETER["Longitude of natural origin",103.833333333333,
##             ANGLEUNIT["degree",0.0174532925199433],
##             ID["EPSG",8802]],
##         PARAMETER["Scale factor at natural origin",1,
##             SCALEUNIT["unity",1],
##             ID["EPSG",8805]],
##         PARAMETER["False easting",28001.642,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8806]],
##         PARAMETER["False northing",38744.572,
##             LENGTHUNIT["metre",1],
##             ID["EPSG",8807]]],
##     CS[Cartesian,2],
##         AXIS["northing (N)",north,
##             ORDER[1],
##             LENGTHUNIT["metre",1]],
##         AXIS["easting (E)",east,
##             ORDER[2],
##             LENGTHUNIT["metre",1]],
##     USAGE[
##         SCOPE["unknown"],
##         AREA["Singapore"],
##         BBOX[1.13,103.59,1.47,104.07]],
##     ID["EPSG",3414]]

2.2.1.3 Select fields are needed

sf_childcare_2017 <- sf_childcare_2017 %>% 
                    dplyr::select(OBJECTID, ADDRESSPOS, ADDRESSSTR, NAME, geometry)
sf_childcare_2020 <- sf_childcare_2020 %>% 
                    dplyr::select(ADDRESSPOS, ADDRESSSTR, Name, geometry)

2.2.1.4 Check the existence of NA values

sf_childcare_2017[rowSums(is.na(sf_childcare_2017))!=0,]

## Simple feature collection with 0 features and 4 fields
## bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS:  SVY21 / Singapore TM
## [1] OBJECTID   ADDRESSPOS ADDRESSSTR NAME       geometry  
## <0 rows> (or 0-length row.names)

sf_childcare_2020[rowSums(is.na(sf_childcare_2020))!=0,]

## Simple feature collection with 0 features and 3 fields
## bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS:  SVY21 / Singapore TM
## [1] ADDRESSPOS ADDRESSSTR Name       geometry  
## <0 rows> (or 0-length row.names)

mpsz[rowSums(is.na(mpsz))!=0,]

## Simple feature collection with 0 features and 15 fields
## bbox:           xmin: NA ymin: NA xmax: NA ymax: NA
## projected CRS:  SVY21 / Singapore TM
##  [1] OBJECTID   SUBZONE_NO SUBZONE_N  SUBZONE_C  CA_IND     PLN_AREA_N
##  [7] PLN_AREA_C REGION_N   REGION_C   INC_CRC    FMEL_UPD_D X_ADDR    
## [13] Y_ADDR     SHAPE_Leng SHAPE_Area geometry  
## <0 rows> (or 0-length row.names)

2.2.1.5 Check validity of geometry

length(which(st_is_valid(sf_childcare_2017)==FALSE))

## [1] 0

length(which(st_is_valid(sf_childcare_2020)==FALSE))

## [1] 0

length(which(st_is_valid(mpsz)==FALSE))

## [1] 9

Make mpsz geometry valid

mpsz <- st_make_valid(mpsz)
st_is_valid(mpsz)

##   [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [16] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [31] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [46] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [61] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [76] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [91] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [106] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [121] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [136] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [151] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [166] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [181] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [196] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [226] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [241] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [256] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [271] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [286] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [301] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [316] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Check dulplicates

str(sf_childcare_2017[duplicated(sf_childcare_2017[,c("ADDRESSSTR", "NAME")]), ]) 

## Classes 'sf' and 'data.frame':   0 obs. of  5 variables:
##  $ OBJECTID  : int 
##  $ ADDRESSPOS: chr 
##  $ ADDRESSSTR: chr 
##  $ NAME      : chr 
##  $ geometry  :sfc_GEOMETRY of length 0 - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA NA
##   ..- attr(*, "names")= chr [1:4] "OBJECTID" "ADDRESSPOS" "ADDRESSSTR" "NAME"

str(sf_childcare_2020[duplicated(sf_childcare_2020[,c("ADDRESSSTR", "Name")]), ]) 

## Classes 'sf' and 'data.frame':   0 obs. of  4 variables:
##  $ ADDRESSPOS: chr 
##  $ ADDRESSSTR: chr 
##  $ Name      : chr 
##  $ geometry  :sfc_GEOMETRY of length 0 - attr(*, "sf_column")= chr "geometry"
##  - attr(*, "agr")= Factor w/ 3 levels "constant","aggregate",..: NA NA NA
##   ..- attr(*, "names")= chr [1:3] "ADDRESSPOS" "ADDRESSSTR" "Name"

2.2.2 Aspatial data

2.2.2.1 Check the existence of NA values

popagsex[rowSums(is.na(popagsex))!=0,]

## # A tibble: 0 x 7
## # … with 7 variables: planning_area <chr>, subzone <chr>, age_group <chr>,
## #   sex <chr>, type_of_dwelling <chr>, resident_count <dbl>, year <dbl>

2.2.2.2 Check for duplicates

str(popagsex[duplicated(popagsex[,c("planning_area", "subzone", "age_group", "sex", "type_of_dwelling", "year")]), ]) 

## tibble [0 × 7] (S3: tbl_df/tbl/data.frame)
##  $ planning_area   : chr(0) 
##  $ subzone         : chr(0) 
##  $ age_group       : chr(0) 
##  $ sex             : chr(0) 
##  $ type_of_dwelling: chr(0) 
##  $ resident_count  : num(0) 
##  $ year            : num(0)

3 Data Preparation

As the objective is to understand the supply and demand of childcare services in 2017 and 2020 at planning subzone level. We need to extract both 2017 and 2020 records from popagesex data set. However, the latest data set of Singapore residents by planning area, subzone, age group, sex and type of dwelling only provided until year = 2019. Therefore, we use year = 2019 as 2020 data.

Referring to the description of child care centres and kindergartens in Singapore Early Childhood Development Agency (ECDA). Child care centres provide child care services and pre-school developmental programmes for children aged between 18 months and below 7 years old¹. As the data set only provides age with range. We need to use age group 0_to_4 and 5_to_9 to calculate it.

3.1 Data wrangling

3.1.1 Supply of child care centre

Calculate the number of child care centres in each subzone and save different year in different columns.

childcare_supply_mpsz <- mpsz
childcare_supply_mpsz$`CHILDCARE 2017` <- lengths(st_intersects(childcare_supply_mpsz, sf_childcare_2017))
childcare_supply_mpsz$`CHILDCARE 2020` <- lengths(st_intersects(childcare_supply_mpsz, sf_childcare_2020))
childcare_supply_mpsz <- childcare_supply_mpsz %>%
  mutate(`CHILDCARE DIFF 2017-20` = `CHILDCARE 2020` - `CHILDCARE 2017`)

3.1.2 Demand of child care centre

summary(popagsex)

##  planning_area        subzone           age_group             sex           
##  Length:883728      Length:883728      Length:883728      Length:883728     
##  Class :character   Class :character   Class :character   Class :character  
##  Mode  :character   Mode  :character   Mode  :character   Mode  :character  
##                                                                             
##                                                                             
##                                                                             
##  type_of_dwelling   resident_count         year     
##  Length:883728      Min.   :   0.00   Min.   :2011  
##  Class :character   1st Qu.:   0.00   1st Qu.:2013  
##  Mode  :character   Median :   0.00   Median :2015  
##                     Mean   :  39.83   Mean   :2015  
##                     3rd Qu.:  10.00   3rd Qu.:2017  
##                     Max.   :2860.00   Max.   :2019

popagsex_2017 <- popagsex %>% 
  filter(year == 2017) %>% 
  group_by(`planning_area`, `subzone`,`age_group`) %>%
  summarize(`POP` = sum(`resident_count`)) %>%
  ungroup() %>%
  spread(`age_group`, POP) %>% 
  dplyr::select(planning_area, subzone, `0_to_4`, `5_to_9`) %>% 
  mutate(`1_to_4` = `0_to_4`/5*4) %>%
  mutate(`5_to_6` = `5_to_9`/5*2) %>%
  mutate(`1_to_6` = `1_to_4` + `5_to_6`) %>%
  dplyr::select(`planning_area`, `subzone`, `1_to_6`) %>%
  mutate_at(.vars = vars(planning_area, subzone), .funs = funs(toupper))

## `summarise()` regrouping output by 'planning_area', 'subzone' (override with `.groups` argument)

## Warning: `funs()` is deprecated as of dplyr 0.8.0.
## Please use a list of either functions or lambdas: 
## 
##   # Simple named list: 
##   list(mean = mean, median = median)
## 
##   # Auto named with `tibble::lst()`: 
##   tibble::lst(mean, median)
## 
##   # Using lambdas
##   list(~ mean(., trim = .2), ~ median(., na.rm = TRUE))
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.

popagsex_2020 <- popagsex %>%
  filter(year == 2019) %>% 
  group_by(`planning_area`, `subzone`,`age_group`) %>%
  summarize(`POP` = sum(`resident_count`)) %>%
  ungroup() %>%
  spread(`age_group`, POP) %>% 
  dplyr::select(planning_area, subzone, `0_to_4`, `5_to_9`) %>% 
  mutate(`1_to_4` = `0_to_4`/5*4) %>%
  mutate(`5_to_6` = `5_to_9`/5*2) %>%
  mutate(`1_to_6` = `1_to_4` + `5_to_6`) %>%
  dplyr::select(`planning_area`, `subzone`, `1_to_6`) %>%
  mutate_at(.vars = vars(planning_area, subzone), .funs = funs(toupper))

## `summarise()` regrouping output by 'planning_area', 'subzone' (override with `.groups` argument)

3.2 Joining the attribute data and geospatial data

mpsz_child_2017 <- left_join(mpsz, popagsex_2017, 
                              by = c("SUBZONE_N" = "subzone"))
mpsz_child_2020 <- left_join(mpsz, popagsex_2020, 
                              by = c("SUBZONE_N" = "subzone"))

4 Analysis

4.1 Section A

The supply and demand of childcare and kindergarten services by planning subzone.

4.1.1 Exploratory Spatial Data Analysis

4.1.1.1 Supply of child care centre

tm_shape(childcare_supply_mpsz) + 
  tm_fill("CHILDCARE 2017", 
          n = 4,
          palette = "Blues", 
          title = "Number of Childcare Centres") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(main.title = "Number of Childcare Centres by Planning Subzone in 2017", 
            main.title.position = "center", 
            main.title.size = 1, 
            frame = FALSE) +
  tmap_style("white")

## tmap style set to "white"

## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

tm_shape(childcare_supply_mpsz) + 
  tm_fill("CHILDCARE 2020", 
          n = 4,
          palette = "Blues", 
          title = "Number of Childcare Centres") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(main.title = "Number of Childcare Centres by Planning Subzone in 2020", 
            main.title.position = "center", 
            main.title.size = 1, 
            frame = FALSE) +
  tmap_style("white")

## tmap style set to "white"
## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

4.1.1.2 Demand of child care centre

tm_shape(mpsz_child_2017) + 
  tm_fill("1_to_6", 
          palette = "Oranges", 
          title = "Number of children aged between 1 to 6") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(main.title = "Number of children aged between 1 to 6 by Planning Subzone in 2017", 
            main.title.position = "center", 
            main.title.size = 1, 
            frame = FALSE) +
  tmap_style("white")

## tmap style set to "white"

## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

tm_shape(mpsz_child_2020) + 
  tm_fill("1_to_6", 
          palette = "Oranges",  
          title = "Number of children aged between 1 to 6") + 
  tm_borders(alpha = 0.5) + 
  tm_layout(main.title = "Number of children aged between 1 to 6 by Planning Subzone in 2020", 
            main.title.position = "center", 
            main.title.size = 1, 
            frame = FALSE) +
  tmap_style("white")

## tmap style set to "white"
## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

The figure above shows that there are five clusters inside this map. West region, north east region, and east region contain one cluster in each. And north region contains two clusters.

The demand from west region has decreased as the color of the cluster become lighter. The subzone with the largest number of children between 1 to 6 has changed from range 4,000 to 5,000 in 2017 to 3,000 to 4,000 in 2020.

The demand from east region has decreased as the color of some subzones become lighter in 2020. And the demand from north and north east regions has a significant increase with showing the colors of the subzones become darker than 2017 in 2020.

4.1.2 Analytics mapping

Using appropriate EDA and choropleth mapping techniques to reveal the supply and demand of childcare services in 2017 and 2020 at the planning subzone level. Describe the spatial patterns observed.

4.1.2.1 Difference between 2017 and 2020

4.1.2.1.1 Childcare centre

tm_shape(childcare_supply_mpsz) + 
  tm_fill("CHILDCARE DIFF 2017-20", 
          n = 5,
          palette = "RdBu",
          title = "Number of Childcare Centres Increased") + 
  tm_borders(alpha = 0.5) +   
  tm_layout(main.title = "Number of Childcare Centres Increased between 2017-2020 by Planning Subzone", 
            main.title.position = "center", 
            main.title.size = 1, 
            frame = FALSE) +
  tmap_style("white")

## tmap style set to "white"

## other available styles are: "gray", "natural", "cobalt", "col_blind", "albatross", "beaver", "bw", "classic", "watercolor"

## Variable(s) "CHILDCARE DIFF 2017-20" contains positive and negative values, so midpoint is set to 0. Set midpoint = NA to show the full spectrum of the color palette.

4.1.3 Geocommunication

The figure above have shown the increase of childcare centre from 2017 to 2020. The demand and supply between childcare centre is not obvious. For example, the child care center still increase in the west region, but the supply is has decreased.

4.2 Section B

4.2.1 Exploratory Spatial Data Analysis

Using point mapping techniques, display the location of childcare services in 2019 and 2020 at the national level. Describe the spatial patterns reveal by their respective distribution. #### Data preparation

Importing the spatial data

ogr_childcare_2017 <- readOGR(dsn = "../data/geospatial", layer="CHILDCARE")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial", layer: "CHILDCARE"
## with 1312 features
## It has 18 fields

ogr_childcare_2020 <- readOGR(dsn = "../data/geospatial", layer="CHILDCARE_2020")

## Warning in OGRSpatialRef(dsn, layer, morphFromESRI = morphFromESRI,
## dumpSRS = dumpSRS, : Discarded datum SVY21 in CRS definition: +proj=tmerc
## +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642
## +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial", layer: "CHILDCARE_2020"
## with 1545 features
## It has 25 fields
## Integer64 fields read as strings:  tessellate extrude visibility drawOrder

ogr_mpsz <- readOGR(dsn = "../data/geospatial", layer = "MP14_SUBZONE_WEB_PL")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial", layer: "MP14_SUBZONE_WEB_PL"
## with 323 features
## It has 15 fields

ogr_sg <- readOGR(dsn = "../data/geospatial" , layer = "CostalOutline")

## OGR data source with driver: ESRI Shapefile 
## Source: "/Users/admin/Documents/GitHub/IS415-Exercises/take home exercises/take_home_ex_01/data/geospatial", layer: "CostalOutline"
## with 60 features
## It has 4 fields

Ensure that they are projected in same projection system

crs(ogr_childcare_2017)

## CRS arguments:
##  +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
## +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

crs(ogr_childcare_2020)

## CRS arguments:
##  +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
## +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0
## +units=m +no_defs

crs(ogr_mpsz)

## CRS arguments:
##  +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
## +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

crs(ogr_sg)

## CRS arguments:
##  +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1
## +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

Examine the imported geospatial data by using plot().

par(mfrow=c(1,3))
plot(ogr_childcare_2017)
plot(ogr_childcare_2020)
plot(ogr_mpsz)

plot(ogr_sg)

Converting the spatial point data frame into generic sp format.

sp_childcare_2017 <- as(ogr_childcare_2017, "SpatialPoints")
sp_childcare_2020 <- as(ogr_childcare_2020, "SpatialPoints")
sp_mpsz <- as(ogr_mpsz, "SpatialPolygons")
sp_sg <- as(ogr_sg, "SpatialPolygons")

sp_childcare_2017

## class       : SpatialPoints 
## features    : 1312 
## extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

sp_childcare_2020

## class       : SpatialPoints 
## features    : 1545 
## extent      : 11203.01, 45404.24, 25667.6, 49300.88  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs

sp_mpsz

## class       : SpatialPolygons 
## features    : 323 
## extent      : 2667.538, 56396.44, 15748.72, 50256.33  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

sp_sg

## class       : SpatialPolygons 
## features    : 60 
## extent      : 2663.926, 56047.79, 16357.98, 50244.03  (xmin, xmax, ymin, ymax)
## crs         : +proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +datum=WGS84 +units=m +no_defs

Converting the generic sp format into spatstat’s ppp format

childcare_2017_ppp <- as(sp_childcare_2017, "ppp")
childcare_2020_ppp <- as(sp_childcare_2020, "ppp")
plot(childcare_2017_ppp)

plot(childcare_2020_ppp)

summary(childcare_2017_ppp)

## Planar point pattern:  1312 points
## Average intensity 1.623186e-06 points per square unit
## 
## *Pattern contains duplicated points*
## 
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
## 
## Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
##                     (34200 x 23630 units)
## Window area = 808287000 square units

summary(childcare_2020_ppp)

## Planar point pattern:  1545 points
## Average intensity 1.91145e-06 points per square unit
## 
## *Pattern contains duplicated points*
## 
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
## 
## Window: rectangle = [11203.01, 45404.24] x [25667.6, 49300.88] units
##                     (34200 x 23630 units)
## Window area = 808287000 square units

Handling duplicated points

We can check the duplication in a ppp object by using the code chunk below.

any(duplicated(childcare_2017_ppp))

## [1] TRUE

any(duplicated(childcare_2020_ppp))

## [1] TRUE

Use jittering to add a small perturbation to the duplicate points so that they do not occupy the exact same space.

childcare_2017_ppp_jit <- rjitter(childcare_2017_ppp, retry = TRUE, nsim=1, drop = TRUE)
childcare_2020_ppp_jit <- rjitter(childcare_2020_ppp, retry = TRUE, nsim=1, drop = TRUE)
any(duplicated(childcare_2017_ppp_jit))

## [1] FALSE

any(duplicated(childcare_2020_ppp_jit))

## [1] FALSE

Creating owin

sg_owin <- as(sp_sg, "owin")
plot(sg_owin)

summary(sg_owin)

## Window: polygonal boundary
## 60 separate polygons (no holes)
##             vertices        area relative.area
## polygon 1         38 1.56140e+04      2.09e-05
## polygon 2        735 4.69093e+06      6.27e-03
## polygon 3         49 1.66986e+04      2.23e-05
## polygon 4         76 3.12332e+05      4.17e-04
## polygon 5       5141 6.36179e+08      8.50e-01
## polygon 6         42 5.58317e+04      7.46e-05
## polygon 7         67 1.31354e+06      1.75e-03
## polygon 8         15 4.46420e+03      5.96e-06
## polygon 9         14 5.46674e+03      7.30e-06
## polygon 10        37 5.26194e+03      7.03e-06
## polygon 11        53 3.44003e+04      4.59e-05
## polygon 12        74 5.82234e+04      7.78e-05
## polygon 13        69 5.63134e+04      7.52e-05
## polygon 14       143 1.45139e+05      1.94e-04
## polygon 15       165 3.38736e+05      4.52e-04
## polygon 16       130 9.40465e+04      1.26e-04
## polygon 17        19 1.80977e+03      2.42e-06
## polygon 18        16 2.01046e+03      2.69e-06
## polygon 19        93 4.30642e+05      5.75e-04
## polygon 20        90 4.15092e+05      5.54e-04
## polygon 21       721 1.92795e+06      2.57e-03
## polygon 22       330 1.11896e+06      1.49e-03
## polygon 23       115 9.28394e+05      1.24e-03
## polygon 24        37 1.01705e+04      1.36e-05
## polygon 25        25 1.66227e+04      2.22e-05
## polygon 26        10 2.14507e+03      2.86e-06
## polygon 27       190 2.02489e+05      2.70e-04
## polygon 28       175 9.25904e+05      1.24e-03
## polygon 29      1993 9.99217e+06      1.33e-02
## polygon 30        38 2.42492e+04      3.24e-05
## polygon 31        24 6.35239e+03      8.48e-06
## polygon 32        53 6.35791e+05      8.49e-04
## polygon 33        41 1.60161e+04      2.14e-05
## polygon 34        22 2.54368e+03      3.40e-06
## polygon 35        30 1.08382e+04      1.45e-05
## polygon 36       327 2.16921e+06      2.90e-03
## polygon 37       111 6.62927e+05      8.85e-04
## polygon 38        90 1.15991e+05      1.55e-04
## polygon 39        98 6.26829e+04      8.37e-05
## polygon 40       415 3.25384e+06      4.35e-03
## polygon 41       222 1.51142e+06      2.02e-03
## polygon 42       107 6.33039e+05      8.45e-04
## polygon 43         7 2.48299e+03      3.32e-06
## polygon 44        17 3.28303e+04      4.38e-05
## polygon 45        26 8.34758e+03      1.11e-05
## polygon 46       177 4.67446e+05      6.24e-04
## polygon 47        16 3.19460e+03      4.27e-06
## polygon 48        15 4.87296e+03      6.51e-06
## polygon 49        66 1.61841e+04      2.16e-05
## polygon 50       149 5.63430e+06      7.53e-03
## polygon 51       609 2.62570e+07      3.51e-02
## polygon 52         8 7.82256e+03      1.04e-05
## polygon 53       976 2.33447e+07      3.12e-02
## polygon 54        55 8.25379e+04      1.10e-04
## polygon 55       976 2.33447e+07      3.12e-02
## polygon 56        61 3.33449e+05      4.45e-04
## polygon 57         6 1.68410e+04      2.25e-05
## polygon 58         4 9.45963e+03      1.26e-05
## polygon 59        46 6.99702e+05      9.35e-04
## polygon 60        13 7.00873e+04      9.36e-05
## enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
##                      (53380 x 33890 units)
## Window area = 748741000 square units
## Fraction of frame area: 0.414

Combining childcare points and the study area

childcare_SG_ppp_2017 = childcare_2017_ppp_jit[sg_owin]
childcare_SG_ppp_2020 = childcare_2020_ppp_jit[sg_owin]
plot(childcare_SG_ppp_2017)

plot(childcare_SG_ppp_2020)

summary(childcare_SG_ppp_2017)

## Planar point pattern:  1312 points
## Average intensity 1.752274e-06 points per square unit
## 
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
## 
## Window: polygonal boundary
## 60 separate polygons (no holes)
##             vertices        area relative.area
## polygon 1         38 1.56140e+04      2.09e-05
## polygon 2        735 4.69093e+06      6.27e-03
## polygon 3         49 1.66986e+04      2.23e-05
## polygon 4         76 3.12332e+05      4.17e-04
## polygon 5       5141 6.36179e+08      8.50e-01
## polygon 6         42 5.58317e+04      7.46e-05
## polygon 7         67 1.31354e+06      1.75e-03
## polygon 8         15 4.46420e+03      5.96e-06
## polygon 9         14 5.46674e+03      7.30e-06
## polygon 10        37 5.26194e+03      7.03e-06
## polygon 11        53 3.44003e+04      4.59e-05
## polygon 12        74 5.82234e+04      7.78e-05
## polygon 13        69 5.63134e+04      7.52e-05
## polygon 14       143 1.45139e+05      1.94e-04
## polygon 15       165 3.38736e+05      4.52e-04
## polygon 16       130 9.40465e+04      1.26e-04
## polygon 17        19 1.80977e+03      2.42e-06
## polygon 18        16 2.01046e+03      2.69e-06
## polygon 19        93 4.30642e+05      5.75e-04
## polygon 20        90 4.15092e+05      5.54e-04
## polygon 21       721 1.92795e+06      2.57e-03
## polygon 22       330 1.11896e+06      1.49e-03
## polygon 23       115 9.28394e+05      1.24e-03
## polygon 24        37 1.01705e+04      1.36e-05
## polygon 25        25 1.66227e+04      2.22e-05
## polygon 26        10 2.14507e+03      2.86e-06
## polygon 27       190 2.02489e+05      2.70e-04
## polygon 28       175 9.25904e+05      1.24e-03
## polygon 29      1993 9.99217e+06      1.33e-02
## polygon 30        38 2.42492e+04      3.24e-05
## polygon 31        24 6.35239e+03      8.48e-06
## polygon 32        53 6.35791e+05      8.49e-04
## polygon 33        41 1.60161e+04      2.14e-05
## polygon 34        22 2.54368e+03      3.40e-06
## polygon 35        30 1.08382e+04      1.45e-05
## polygon 36       327 2.16921e+06      2.90e-03
## polygon 37       111 6.62927e+05      8.85e-04
## polygon 38        90 1.15991e+05      1.55e-04
## polygon 39        98 6.26829e+04      8.37e-05
## polygon 40       415 3.25384e+06      4.35e-03
## polygon 41       222 1.51142e+06      2.02e-03
## polygon 42       107 6.33039e+05      8.45e-04
## polygon 43         7 2.48299e+03      3.32e-06
## polygon 44        17 3.28303e+04      4.38e-05
## polygon 45        26 8.34758e+03      1.11e-05
## polygon 46       177 4.67446e+05      6.24e-04
## polygon 47        16 3.19460e+03      4.27e-06
## polygon 48        15 4.87296e+03      6.51e-06
## polygon 49        66 1.61841e+04      2.16e-05
## polygon 50       149 5.63430e+06      7.53e-03
## polygon 51       609 2.62570e+07      3.51e-02
## polygon 52         8 7.82256e+03      1.04e-05
## polygon 53       976 2.33447e+07      3.12e-02
## polygon 54        55 8.25379e+04      1.10e-04
## polygon 55       976 2.33447e+07      3.12e-02
## polygon 56        61 3.33449e+05      4.45e-04
## polygon 57         6 1.68410e+04      2.25e-05
## polygon 58         4 9.45963e+03      1.26e-05
## polygon 59        46 6.99702e+05      9.35e-04
## polygon 60        13 7.00873e+04      9.36e-05
## enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
##                      (53380 x 33890 units)
## Window area = 748741000 square units
## Fraction of frame area: 0.414

summary(childcare_SG_ppp_2020)

## Planar point pattern:  1545 points
## Average intensity 2.063463e-06 points per square unit
## 
## Coordinates are given to 3 decimal places
## i.e. rounded to the nearest multiple of 0.001 units
## 
## Window: polygonal boundary
## 60 separate polygons (no holes)
##             vertices        area relative.area
## polygon 1         38 1.56140e+04      2.09e-05
## polygon 2        735 4.69093e+06      6.27e-03
## polygon 3         49 1.66986e+04      2.23e-05
## polygon 4         76 3.12332e+05      4.17e-04
## polygon 5       5141 6.36179e+08      8.50e-01
## polygon 6         42 5.58317e+04      7.46e-05
## polygon 7         67 1.31354e+06      1.75e-03
## polygon 8         15 4.46420e+03      5.96e-06
## polygon 9         14 5.46674e+03      7.30e-06
## polygon 10        37 5.26194e+03      7.03e-06
## polygon 11        53 3.44003e+04      4.59e-05
## polygon 12        74 5.82234e+04      7.78e-05
## polygon 13        69 5.63134e+04      7.52e-05
## polygon 14       143 1.45139e+05      1.94e-04
## polygon 15       165 3.38736e+05      4.52e-04
## polygon 16       130 9.40465e+04      1.26e-04
## polygon 17        19 1.80977e+03      2.42e-06
## polygon 18        16 2.01046e+03      2.69e-06
## polygon 19        93 4.30642e+05      5.75e-04
## polygon 20        90 4.15092e+05      5.54e-04
## polygon 21       721 1.92795e+06      2.57e-03
## polygon 22       330 1.11896e+06      1.49e-03
## polygon 23       115 9.28394e+05      1.24e-03
## polygon 24        37 1.01705e+04      1.36e-05
## polygon 25        25 1.66227e+04      2.22e-05
## polygon 26        10 2.14507e+03      2.86e-06
## polygon 27       190 2.02489e+05      2.70e-04
## polygon 28       175 9.25904e+05      1.24e-03
## polygon 29      1993 9.99217e+06      1.33e-02
## polygon 30        38 2.42492e+04      3.24e-05
## polygon 31        24 6.35239e+03      8.48e-06
## polygon 32        53 6.35791e+05      8.49e-04
## polygon 33        41 1.60161e+04      2.14e-05
## polygon 34        22 2.54368e+03      3.40e-06
## polygon 35        30 1.08382e+04      1.45e-05
## polygon 36       327 2.16921e+06      2.90e-03
## polygon 37       111 6.62927e+05      8.85e-04
## polygon 38        90 1.15991e+05      1.55e-04
## polygon 39        98 6.26829e+04      8.37e-05
## polygon 40       415 3.25384e+06      4.35e-03
## polygon 41       222 1.51142e+06      2.02e-03
## polygon 42       107 6.33039e+05      8.45e-04
## polygon 43         7 2.48299e+03      3.32e-06
## polygon 44        17 3.28303e+04      4.38e-05
## polygon 45        26 8.34758e+03      1.11e-05
## polygon 46       177 4.67446e+05      6.24e-04
## polygon 47        16 3.19460e+03      4.27e-06
## polygon 48        15 4.87296e+03      6.51e-06
## polygon 49        66 1.61841e+04      2.16e-05
## polygon 50       149 5.63430e+06      7.53e-03
## polygon 51       609 2.62570e+07      3.51e-02
## polygon 52         8 7.82256e+03      1.04e-05
## polygon 53       976 2.33447e+07      3.12e-02
## polygon 54        55 8.25379e+04      1.10e-04
## polygon 55       976 2.33447e+07      3.12e-02
## polygon 56        61 3.33449e+05      4.45e-04
## polygon 57         6 1.68410e+04      2.25e-05
## polygon 58         4 9.45963e+03      1.26e-05
## polygon 59        46 6.99702e+05      9.35e-04
## polygon 60        13 7.00873e+04      9.36e-05
## enclosing rectangle: [2663.93, 56047.79] x [16357.98, 50244.03] units
##                      (53380 x 33890 units)
## Window area = 748741000 square units
## Fraction of frame area: 0.414

4.2.2 With reference to the spatial point patterns observed

4.2.2.1 Formulate the null hypothesis and alternative hypothesis and select the confidence level.

Extracting selected study area due to limited computation power

sp_mpsz <- as_Spatial(mpsz)

## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO"): Discarded datum Unknown based on WGS84 ellipsoid in CRS definition,
##  but +towgs84= values preserved

sk = sp_mpsz[sp_mpsz@data$PLN_AREA_N == "SENGKANG",]
bd = sp_mpsz[sp_mpsz@data$PLN_AREA_N == "BEDOK",]
bb = sp_mpsz[sp_mpsz@data$PLN_AREA_N == "BUKIT BATOK",]
hg = sp_mpsz[sp_mpsz@data$PLN_AREA_N == "HOUGANG",]
plot(sk)

plot(bb)

plot(bd)

plot(hg)

Convert SpatialPolygonDataFrame into generic spatialpolygons objects

sp_sk = as(sk, "SpatialPolygons")
sp_bb = as(bb, "SpatialPolygons")
sp_bd = as(bd, "SpatialPolygons")
sp_hg = as(hg, "SpatialPolygons")

Convert these SpatialPolygons objects into owin objects that is required by spatstat.

owin_sk = as(sp_sk, "owin")
owin_bb = as(sp_bb, "owin")
owin_bd = as(sp_bd, "owin")
owin_hg = as(sp_hg, "owin")

Combine points and study area

childcare_sk_ppp_2017 = childcare_2017_ppp_jit[owin_sk]
childcare_bb_ppp_2017 = childcare_2017_ppp_jit[owin_bb]
childcare_bd_ppp_2017 = childcare_2017_ppp_jit[owin_bd]
childcare_hg_ppp_2017 = childcare_2017_ppp_jit[owin_hg]

childcare_sk_ppp_2020 = childcare_2020_ppp_jit[owin_sk]
childcare_bb_ppp_2020 = childcare_2020_ppp_jit[owin_bb]
childcare_bd_ppp_2020 = childcare_2020_ppp_jit[owin_bd]
childcare_hg_ppp_2020 = childcare_2020_ppp_jit[owin_hg]

Plot to check

plot(childcare_sk_ppp_2017)

plot(childcare_bb_ppp_2017)

plot(childcare_bd_ppp_2017)

plot(childcare_hg_ppp_2017)

plot(childcare_sk_ppp_2020)

plot(childcare_bb_ppp_2020)

plot(childcare_bd_ppp_2020)

plot(childcare_hg_ppp_2020)

4.2.2.2 Perform the test by using appropriate 2nd order spatial point patterns analysis technique.

Analysing Spatial Point Process Using L-Function

4.2.2.2.1 Year 2017

4.2.2.2.1.1 Computing L Fucntion estimation

L_ck = Lest(childcare_sk_ppp_2017, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_bb_ppp_2017, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_bd_ppp_2017, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_hg_ppp_2017, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

###### Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at 4 planning areas are randomly distributed in year 2017.

H1= The distribution of childcare services at 4 planning areas are not randomly distributed in year 2017.

The null hypothesis will be rejected if p-value if smaller than alpha value of 0.001.

As the alpha value of this hypothesis test is 0.01, we run 2/0.01 = 200 simulations for our Monte Carlo simulations.

The code chunk below is used to perform the hypothesis testing.

L_sk_2017.csr <- envelope(childcare_sk_ppp_2017, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_bb_2017.csr <- envelope(childcare_bb_ppp_2017, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_bd_2017.csr <- envelope(childcare_bd_ppp_2017, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_hg_2017.csr <- envelope(childcare_hg_ppp_2017, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

plot(L_sk_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Sengkang 2017")

plot(L_bb_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Bukit Batok 2017")

plot(L_bd_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Bedok 2017")

plot(L_hg_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Hougang 2017")

Figure L-funtion in Sengkang in 2017 largely lies above the envelope and has the signs of clusters forming after 100m. Therefore spatial clustering is statistically significant and we reject the Ho.

Figure L-funtion in Bukit Batok 2017 largely lies within the envelope. Therefore spatial clustering is not statistically significant and we fail to reject the Ho.

Figure L-funtion in Bedock 2017 largely lies within the envelope. Therefore spatial clustering is not statistically significant and we fail to reject the Ho.

Figure L-funtion in Hougang 2017 shows some signs of clustering at certain distances but it largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.

4.2.2.2.2 Year 2020

4.2.2.2.2.1 Computing L Fucntion estimation

L_ck = Lest(childcare_sk_ppp_2020, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_bb_ppp_2020, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_bd_ppp_2020, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

L_ck = Lest(childcare_hg_ppp_2020, correction = "Ripley")
plot(L_ck, . -r ~ r, 
     ylab= "L(d)-r", xlab = "d(m)")

###### Performing Complete Spatial Randomness Test

To confirm the observed spatial patterns above, a hypothesis test will be conducted. The hypothesis and test are as follows:

Ho = The distribution of childcare services at 4 planning areas are randomly distributed in year 2020.

H1= The distribution of childcare services at 4 planning areas are not randomly distributed in year 2020.

The null hypothesis will be rejected if p-value if smaller than alpha value of 0.001.

As the alpha value of this hypothesis test is 0.01, we run 2/0.01 = 200 simulations for our Monte Carlo simulations.

The code chunk below is used to perform the hypothesis testing.

L_sk_2017.csr <- envelope(childcare_sk_ppp_2020, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_bb_2017.csr <- envelope(childcare_bb_ppp_2020, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_bd_2017.csr <- envelope(childcare_bd_ppp_2020, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

L_hg_2017.csr <- envelope(childcare_hg_ppp_2020, Lest, nsim = 199, rank = 1, glocal=TRUE)

## Generating 199 simulations of CSR  ...
## 1, 2, 3, 4.6.8.10.12.14.16.18.20.22.24.26.28.30.32.34.36.38.40
## .42.44.46.48.50.52.54.56.58.60.62.64.66.68.70.72.74.76.78.80
## .82.84.86.88.90.92.94.96.98.100.102.104.106.108.110.112.114.116.118.120
## .122.124.126.128.130.132.134.136.138.140.142.144.146.148.150.152.154.156.158.160
## .162.164.166.168.170.172.174.176.178.180.182.184.186.188.190.192.194.196.198 199.
## 
## Done.

plot(L_sk_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Sengkang 2020")

plot(L_bb_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Bukit Batok 2020")

plot(L_bd_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Bedok 2020")

plot(L_hg_2017.csr, . - r ~ r, 
     xlab="d", ylab="L(d)-r", xlim=c(0,500), main = "L-function in Hougang 2020")

Figure L-funtion in Sengkang 2020 largely lies above the envelope and has the signs of clusters forming after around 120m. Therefore spatial clustering is statistically significant and we reject the Ho.

Figure L-funtion in Bukit Batok 2020 shows some signs of clustering at certain distances but it largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.

Figure L-funtion in Bedok 2020 shows some signs of clustering at certain distances but it largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.

Figure L-funtion in Hougang 2020 shows some signs of clustering at certain distances but it largely lies within the envelope, which means it is not statistically significant and we fail to reject the Ho.

4.2.3 Kernel density

4.2.3.1 Derive kernel density maps of childcare services in 2017 and 2020.

As the default unit of svy21 is in meters and density values computed is using “number of points per square meter”.

Therefore, we need to use rescale is used to covert the unit of measurement from meters to kilometers.

childcare_sk_ppp_2017.km <- rescale(childcare_sk_ppp_2017, 1000, "km")
childcare_bb_ppp_2017.km <- rescale(childcare_bb_ppp_2017, 1000, "km")
childcare_bd_ppp_2017.km <- rescale(childcare_bd_ppp_2017, 1000, "km")
childcare_hg_ppp_2017.km <- rescale(childcare_hg_ppp_2017, 1000, "km")

childcare_sk_ppp_2020.km <- rescale(childcare_sk_ppp_2020, 1000, "km")
childcare_bb_ppp_2020.km <- rescale(childcare_bb_ppp_2020, 1000, "km")
childcare_bd_ppp_2020.km <- rescale(childcare_bd_ppp_2020, 1000, "km")
childcare_hg_ppp_2020.km <- rescale(childcare_hg_ppp_2020, 1000, "km")

Generation

kde_childcare_sk_2017 <- density(childcare_sk_ppp_2017.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcare_bb_2017 <- density(childcare_bb_ppp_2017.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcare_bd_2017 <- density(childcare_bd_ppp_2017.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcare_hg_2017 <- density(childcare_hg_ppp_2017.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")

kde_childcare_sk_2020 <- density(childcare_sk_ppp_2020.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcare_bb_2020 <- density(childcare_bb_ppp_2020.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")

## Warning: Berman-Diggle Cross-Validation criterion was minimised at right-hand
## end of interval [0, 0.235]; use argument 'hmax' to specify a wider interval for
## bandwidth 'sigma'

kde_childcare_bd_2020 <- density(childcare_bd_ppp_2020.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")
kde_childcare_hg_2020 <- density(childcare_hg_ppp_2020.km, sigma= bw.diggle, edge=TRUE, kernel="gaussian")

Plot

plot(kde_childcare_sk_2017)

plot(kde_childcare_bb_2017)

plot(kde_childcare_bd_2017)

plot(kde_childcare_hg_2017)

plot(kde_childcare_sk_2020)

plot(kde_childcare_bb_2020)

plot(kde_childcare_bd_2020)

plot(kde_childcare_hg_2020)

Create kernel density estimation map Function

plot_kd_map <- function(ppp, osm, sz, title) {
  # convert from m to km
  ppp_km <- rescale(ppp, 1, "km")
  
  # kernel density
  kde <- density(ppp_km, sigma = bw.diggle, edge=TRUE, kernel="gaussian")
  
  # convert kde to grid object
  kde_grid <- as.SpatialGridDataFrame.im(kde)
  
  # convert grid object into raster
  kde_raster <- raster(kde_grid)
  
  projection(kde_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
  
  # plot on openstreetmap
  tm_shape(osm) + 
    tm_rgb() +
  tm_shape(sz) +
    tm_borders(col = "black", lwd = 2, lty="longdash") +
  tm_shape(kde_raster) + 
    tm_raster("v", alpha=0.5, palette = "BuPu") +
    tm_layout(legend.outside = TRUE, title=title)
}

4.2.3.2 Display the kernel density maps on openstreetmap of Singapore

sk_bbox <- st_bbox(mpsz_child_2017 %>% filter(PLN_AREA_N == "SENGKANG"))
bd_bbox <- st_bbox(mpsz_child_2017 %>% filter(PLN_AREA_N == "BEDOK"))
bb_bbox <- st_bbox(mpsz_child_2017 %>% filter(PLN_AREA_N == "BUKIT BATOK"))
hg_bbox <- st_bbox(mpsz_child_2017 %>% filter(PLN_AREA_N == "HOUGANG"))

sk_osm <- read_osm(sk_bbox, ext=1.1)
bd_osm <- read_osm(bd_bbox, ext=1.1)
bb_osm <- read_osm(bb_bbox, ext=1.1)
hg_osm <- read_osm(hg_bbox, ext=1.1)

Due to the limitation of size on the document that can be publish on RPubs, we only show the kernel density on OpenStreetMap at Senkang. But the code of other parts still being provided.

Sengkang

plot_kd_map(childcare_sk_ppp_2017, sk_osm, sk, "Sengkeng in year 2017")

## stars object downsampled to 1596 by 626 cells. See tm_shape manual (argument raster.downsample)

plot_kd_map(childcare_sk_ppp_2020, sk_osm, sk, "Sengkeng in year 2020")

## stars object downsampled to 1596 by 626 cells. See tm_shape manual (argument raster.downsample)

There is not much difference between 2017 and 2020.

Bedok

#plot_kd_map(childcare_bd_ppp_2017, bd_osm, bd, "Bedok in year 2017")
#plot_kd_map(childcare_bd_ppp_2020, bd_osm, bd, "Bedok in year 2020")

The center cluster become more obvious in 2020 comparing to 2017.

Bukit Batok

#plot_kd_map(childcare_bb_ppp_2017, bb_osm, bb, "Bukit Batok in year 2017")
#plot_kd_map(childcare_bb_ppp_2020, bb_osm, bb, "Bukit Batok in year 2020")

The size of the south west cluster has increased from 2017 to 2020.

Hougang

#plot_kd_map(childcare_hg_ppp_2017, hg_osm, hg, "Hougang in year 2017")
#plot_kd_map(childcare_hg_ppp_2020, hg_osm, hg, "Hougang in year 2020")

There is a huge increase of size for all clusters from 2017 to 2020.

4.2.3.3 Compare the advantage of kernel density maps over point maps

Kernel density map gives us a clearer idea of the area having higher density as comparing to the point maps.However, it does not plot the exact location on the map since it have group the points together with the subzone. Hence, we cannot really point out the actual location of the clusters.

5 References

---

1. Early Childhood Development Agency. (2019). Retrieved from: https://www.ecda.gov.sg/pages/aboutus.aspx#:~:text=Child%20care%20centres%20provide%20child,2%20and%2018%20months%20old.↩︎