1 Overview

The availability and proximity of childcare services are increasingly important to parents in Singapore. Often times today both parents work and as such if the governments wishes to encourage more couple to have children they need to provide childcare facilites. On the most part this is the aim of the government.

1.1 The Problem Statement

As such this study aims to understand if there has been positive changes in the supply of childcare services and if they are able to meet the demand especially. As such it seeks to understand the temporal supply and demand of childcare services from the year 2017 to 2020 through this study.

1.2 The data

To provide answers to the questions above, two data sets will be used. They are:

CHILDCARE, a point feature data providing both location and attribute information of childcare centers as of 20202. It is downloaded from www.data.gov.sg and is in KML format.
CHILDCARE 2017, a point feature data providing both location and attribute information of childcare centers as of 2017 and is in ESRI shapefile format. (data was procured from Prof Kam Tin Seong)
MP14_SUBZONE_WEB_PL, a polygon feature data providing information of URA 2014 Master Plan Planning Subzone boundary data. It is downloaded from www.data.gov.sg and is in ESRI shapefile format.
Singapore Residents by Planning Area/Subzone, Age Group and Sex, June 2011 - 2019 , It is downloaded from www.singstat.gov.sg and is in csv format (i.e. respopagsex2011to2019.csv). This is an aspatial data fie. Although it does not contain any coordinates values, but it’s PA and SZ fields can be used as unique identifiers to georeference to MP14_SUBZONE_WEB_PL shapefile.

2 Setting up the environment

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

## Warning in fun(libname, pkgname): Java home setting is INVALID, it will be ignored.
## Please do NOT set it unless you want to override system settings.

3 Geospatial Data Wrangling

3.1 Importing aspatial data

The latest population statistics are available in csv format, as such the below chunk will parse in the data into R as a tibble data.frame.

popdata <- read_csv("data/aspatial/respopagesextod2011to2019.csv")

Lets take a look at the popdata and the various fields that it contains.

head(popdata)

Lets also ensure that there are no NA values in the data.frame.

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

There are no NA columns, we can proceed on to preparing the table to be utilised in data analysis.

3.1.1 Data Prepeparation

In Singapore, the Early Childhood Development Agency (ECDA) manages child care centers and recommends childcare center services to be offered to children between the ages of 18 months and 7 years. However the imported population statistics aggregates the age ranges of childcare children into “0_to_4” and “5_to_9” for the years 2017-2019. We can deal with this by:

Estimating the proportion splits for the age group “5_to_9 and mutate the data to retrieve only the relevant age groups (e.g splitting the total sum of”5_to_9" into 5 equal proportions and utilising only 2/5 of the total as a 5-to-6 age group for the 2017 data and 1/5 of the total for the 2019 data to estimate for 2020). This is because to utilise the data as is will lead to inaccuracies in demand as children aged 7 and above are not eligible officially for childcare and will need to utilse student care facilities, which is beyond the scope of this analysis.
Utilise the available data as is for the age group “0_to_4”, even though it may contain kids not yet eligible for childcare (e.g. below 18months) but it may still serve as a indicator of demand for childcare as these kids will utilise childcare and additionaly certain childcare centers also offer infant care as young as 2 months as old.

To this we will be utilising the following code chunk

#retrieve popdata of childcare age kids (0-7) for 2017 
popdata2017 <- popdata %>% 
  filter(Time == 2017) %>% 
  group_by(PA,SZ,AG) %>%
  summarize(`POP` = sum(`Pop`)) %>%
  ungroup() %>%
  spread(AG, POP) %>% 
  mutate(`5_to_6` = 2/5*`5_to_9`)%>% 
  dplyr::select(PA, SZ, `0_to_4`, `5_to_6`)%>% 
  mutate(`Childcare_0_to_6` = rowSums(.[3:4])) %>% 
  dplyr::select(PA, SZ, `Childcare_0_to_6`) %>%
  mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper))

#retrieve popdata of childcare age kids (0-7) for 2020 by using 2019 data
popdata2020 <- popdata %>% 
  filter(Time == 2019) %>% 
  group_by(PA,SZ,AG) %>%
  summarize(`POP` = sum(`Pop`)) %>%
  ungroup() %>%
  spread(AG, POP) %>% 
#only select 1/5 of the population in estimate of 2020 data where 6 year olds will become 7
  mutate(`5_to_6` = 1/5*`5_to_9`)%>% 
  dplyr::select(PA, SZ, `0_to_4`, `5_to_6`)%>% 
  mutate(`Childcare_0_to_6` = rowSums(.[3:4])) %>% 
  dplyr::select(PA, SZ, `Childcare_0_to_6`) %>%
  mutate_at(.vars = vars(PA, SZ), .funs = funs(toupper))

Lets do a quick check to ensure there are no NA values introduced through the manipulations

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

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

There are no NA values we can now conduct a georelational join of the datas with the planning subzone layer.

3.2 Importing the spatial data

3.2.1 URA Master Plan 2017 Planning Subzone

The following code chunk imports the Master Plan data into R as a simple feature dataframe.

sf_mpsz <- st_read(dsn = "data/geospatial", layer="MP14_SUBZONE_WEB_PL")

## Reading layer `MP14_SUBZONE_WEB_PL' from data source `E:\Y4S1\gis\in_class_ex\Take-Home_Ex01\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

We can take a look at this data utilising:

head(sf_mpsz)

3.2.1.1 CRS check and Reassignment

After importing, the following chunk allows us to check that the imported data is projected in an appropriate projection system.

st_crs(sf_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]]]]

We will assign EPSG 3414 to sf_mpsz

sf_mpsz3414 <- st_set_crs(sf_mpsz, 3414)

Check the crs again to ensure that assignment is correct

st_crs(sf_mpsz3414)

## 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]]

3.2.1.2 NA and Geometric Validity Checks and Correction

We will also check to see if there are any NA values.

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

## Warning in `[<-.data.frame`(`*tmp*`, is_list, value = list(`16` = "<S3:
## sfc_GEOMETRY>")): replacement element 1 has 1 row to replace 0 rows

There are no NA values detected.Additionally we will also conduct a geometry check to ensure that the data is not corrupted or invalid when imported into R.The following code chunk does this check

test <- st_is_valid(sf_mpsz3414,reason=TRUE)
length(which(test!= "Valid Geometry"))

## [1] 9

It looks like there are 9 invalid geometries, we can use the following code chunk to query the reason for invalidity

test[which(test!= "Valid Geometry")]

## [1] "Ring Self-intersection[27932.3925999999 21982.7971999999]"
## [2] "Ring Self-intersection[26885.4439000003 26668.3121000007]"
## [3] "Ring Self-intersection[26920.1689999998 26978.5440999996]"
## [4] "Ring Self-intersection[15432.4749999996 31319.716]"       
## [5] "Ring Self-intersection[12861.3828999996 32207.4923]"      
## [6] "Ring Self-intersection[19681.2353999997 31294.4521999992]"
## [7] "Ring Self-intersection[41375.108 40432.8588999994]"       
## [8] "Ring Self-intersection[38542.2260999996 44605.4089000002]"
## [9] "Ring Self-intersection[21702.5623000003 48125.1154999994]"

It seems that certain of the polygons are considered as invalid, we will need to make them valid.

sf_mpsz3414<- st_make_valid(sf_mpsz3414)
test <- st_is_valid(sf_mpsz3414)
length(which(test == FALSE))

## [1] 0

We have utilised the st_make_valid function, to make our invalid geometries valid. Now we can quickly plot the layer to view it.

plot(sf_mpsz3414$geometry)

3.2.1.3 Georelational join with population data

By utilising planning subzone as a unique identifier we will join the sf_mpsz3414 to popdata2017 and 2020.

sf_mpszpop2017 <- left_join(sf_mpsz3414, popdata2017, by = c("SUBZONE_N" = "SZ"))

sf_mpszpop2020 <- left_join(sf_mpsz3414, popdata2020, by = c("SUBZONE_N" = "SZ"))

3.2.2 Childcare Centers 2017

The following code chunk imports the 2017 childcare data into R as a as simple feature dataframes.

sf_childcare2017 <- st_read(dsn = "data/geospatial", layer="CHILDCARE")

## Reading layer `CHILDCARE' from data source `E:\Y4S1\gis\in_class_ex\Take-Home_Ex01\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

3.2.2.1 Data Cleaning

We can take a look at this data utilising:

head(sf_childcare2017, n=4)

There are quite a few additional features and columns that do not value add to this analysis, as such the following code chunk selects the columns that are needed.

sf_childcare2017 <- sf_childcare2017 %>%
  dplyr::select(OBJECTID,ADDRESSPOS,ADDRESSSTR,NAME,X_ADDR,Y_ADDR,geometry)

3.2.2.2 CRS check and Reassignment

The following chunk allows us to check that the imported data is projected in an appropriate projection system.

st_crs(sf_childcare2017)

## 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]]]]

We will assign EPSG 3414 to sf_childcare2017

sf_childcare2017_3414 <- st_set_crs(sf_childcare2017, 3414)

Lets check the CRS again.

st_crs(sf_childcare2017_3414)

## 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]]

3.2.2.3 NA and Geometric Validity Checks and Correction

We will also check to see if there are any NA values.

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

## Warning in `[<-.data.frame`(`*tmp*`, is_list, value = list(`7` = "<S3:
## sfc_GEOMETRY>")): replacement element 1 has 1 row to replace 0 rows

There are no NA values detected.Additionally we will also conduct a geometry check to ensure that the data is not corrupted or invalid.

test <- st_is_valid(sf_childcare2017_3414,reason=TRUE)
length(which(test!= "Valid Geometry"))

## [1] 0

All geometries are valid. We can now plot sf_childcare2017_3414.

plot(sf_childcare2017_3414$geometry)

3.2.3 Childcare Centers 2020

The following code chunk imports the 2017 childcare data which is in KML form into R as a as simple feature dataframes.

sf_childcare2020 = st_read("data/geospatial/child-care-services-kml.kml")

## Reading layer `CHILDCARE' from data source `E:\Y4S1\gis\in_class_ex\Take-Home_Ex01\data\geospatial\child-care-services-kml.kml' using driver `KML'
## Simple feature collection with 1545 features and 2 fields
## geometry type:  POINT
## dimension:      XYZ
## bbox:           xmin: 103.6824 ymin: 1.248403 xmax: 103.9897 ymax: 1.462134
## z_range:        zmin: 0 zmax: 0
## geographic CRS: WGS 84

3.2.3.1 CRS check and Transforming

After importing, the following chunk allows us to check that the imported data is projected in an appropriate projection system.

st_crs(sf_childcare2020)

## Coordinate Reference System:
##   User input: WGS 84 
##   wkt:
## GEOGCRS["WGS 84",
##     DATUM["World Geodetic System 1984",
##         ELLIPSOID["WGS 84",6378137,298.257223563,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["geodetic latitude (Lat)",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["geodetic longitude (Lon)",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4326]]

We will need to transform the sf_childcare2020 into the projected coordinate system that we are utilising, which is svy21(i.e. EPSG 3414)

sf_childcare2020_3414 <- st_transform(sf_childcare2020, 3414)
st_crs(sf_childcare2020_3414)

## 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]]

3.2.3.2 Data Cleaning

We can take a look at this data utilising:

head(sf_childcare2020_3414, n=1)

It looks like most of the data on the point is in the description column in a table format,it might be useful to be able to access at least the name, as such we will need to retrieve the data from the table.

getName <- function(input){
  noHtml <- str_replace_all(input,"<.*?>","")
  return(str_match(noHtml, "\\sNAME\\s*(.*?)\\s*PHOTOURL")[,2])
}

sf_childcare2020_3414 <- sf_childcare2020_3414 %>%
  rename(OBJECTID=Name)%>%
  mutate(`NAME` = getName(sf_childcare2020_3414$Description)) 
head(sf_childcare2020_3414, 1)

3.2.3.3 NA and Geometric Validity Checks and Correction

We will also check to see if there are any NA values.

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

## Warning in `[<-.data.frame`(`*tmp*`, is_list, value = list(`3` = "<S3:
## sfc_GEOMETRY>")): replacement element 1 has 1 row to replace 0 rows

There are no NA values detected.Additionally we will also conduct a geometry check to ensure that the data is not corrupted or invalid.

test <- st_is_valid(sf_childcare2020_3414,reason=TRUE)
length(which(test!= "Valid Geometry"))

## [1] 0

All geometries are valid. We can now plot sf_childcare2020_3414.

plot(sf_childcare2020_3414$geometry)

4 Section A: Exploring the change in supply and demand of childcare between year 2017 and 2020

4.1 Exploratory Data Analysis

4.1.1 Chloropleth Map of Childcare Demand by Subzone

To understand the demand for childcare services, we can utilise a chloropleth map to understand the population distribution of children of childcare age, (0 to 6 years old). As only children of childcare age will demand childcare services, and a direct indicator of demand will be the availabilit of the popularion

Lets create the chloropleth map for 2017 and 2020.

pop_map2017<-tm_shape(sf_mpszpop2017) + 
   tm_polygons() +
  #select only values that have atleast one child
   tm_shape(sf_mpszpop2017[sf_mpszpop2017$Childcare_0_to_6>0, ]) +
  tm_polygons(col = "Childcare_0_to_6", style="jenks", palette="Blues",title = "Population of Childcare Age Children") +
   tm_layout("Childcare Demand by Subzone, 2017",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Population data from Department of 
             Statistics (DOS)",position = c("left", "bottom"))

pop_map2020<-tm_shape(sf_mpszpop2020) + 
   tm_polygons() +
  #select only values that have atleast one child
   tm_shape(sf_mpszpop2020[sf_mpszpop2020$Childcare_0_to_6>0, ]) +
  tm_polygons(col = "Childcare_0_to_6", style="jenks", palette="Blues", title = "Population of Childcare Age Children") +
  tm_layout("Childcare Demand by Subzone, 2020",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Population data from Department of 
             Statistics (DOS)",position = c("left", "bottom"))

tmap_arrange(pop_map2017, pop_map2020, asp=1, ncol=2)

Based on a visual analysis of the chloropleth map we can observe that the demand has dropped between 2017 and 2020, across most subzone. This based on the assumption that children of childcare age will attend child care.

4.1.2 Summary Statistics of Childcare Population

The main reason for the drop can be seen when we look at the summary statistics for each year.

summary(sf_mpszpop2017$Childcare_0_to_6)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0   278.0   831.1  1105.0  8070.0

summary(sf_mpszpop2020$Childcare_0_to_6)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##     0.0     0.0   232.0   697.4   896.0  6796.0

You can see that on average the population of children at child care age have dropped from the year 2017 (831 Children) to 2020 (687 children). This will indicate a similar overall drop in demand for childcare services.

4.1.3 Chloropleth Map of Number of Childcare by Subzone

Understanding the number of childcare available in each subzone will assit us to understand the supply of childcare services. As explorative analysis step, lets create Chloropleth to explore this. Firstly we need to identify then number of childcare centres in each subzone, we can do so whith the following code chunk.

sf_mpszpop2017$`COUNT_CHILDCARE` <- lengths(st_intersects(sf_mpszpop2017, sf_childcare2017_3414))

sf_mpszpop2020$`COUNT_CHILDCARE` <- lengths(st_intersects(sf_mpszpop2020, sf_childcare2020_3414))

We can now create the chloropleth map.

care_map2017<-tm_shape(sf_mpszpop2017) +
   tm_polygons() + 
  #select only values that have atleast one child
   tm_shape(sf_mpszpop2017[sf_mpszpop2017$COUNT_CHILDCARE>0, ]) +
  tm_polygons(col = "COUNT_CHILDCARE", style="jenks", palette="Greens",title = "Number of Childcare Services") +
   tm_layout("Childcare Supply by Subzone, 2017",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Childcare data from Early Childhood 
             Development Agency (ECDA)",position = c("left", "bottom"))

care_map2020<-tm_shape(sf_mpszpop2020) + 
   tm_polygons() +
  #select only values that have atleast one child
   tm_shape(sf_mpszpop2020[sf_mpszpop2020$COUNT_CHILDCARE>0, ]) +
  tm_polygons(col = "COUNT_CHILDCARE", style="jenks", palette="Greens", title = "Number of Childcare Services") +
  tm_layout("Childcare Supply by Subzone, 2020",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Childcare data from Early Childhood 
             Development Agency (ECDA)",position = c("left", "bottom"))

tmap_arrange(care_map2017, care_map2020, asp=1, ncol=2)

Through visual analysis, we can identify that there has been an increase in the supply of childcare services from the year 2017 to the year 2020. There are more subzones with darker shades in the year 2020 than 2017.

4.1.4 Summary Statistics of Childcare Services

We can check on this by utilizing summary statistics.

summary(sf_mpszpop2017$COUNT_CHILDCARE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   2.000   4.062   6.000  36.000

summary(sf_mpszpop2020$COUNT_CHILDCARE)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.000   0.000   3.000   4.783   7.500  42.000

There has been and increase in the number of childcare services from the year 2017 to the year 2020, however the increase in childcare per subzone has not been substantial on average when you observe that the mean only increased by 0.721.

4.1.5 Top 3 Subzones with Most Number of Childcare Age Children and Most Number of Childcare services

It will additionally be interesting to identify which subzones have the most number of children as correspondingly these subzones will be where more childcare services will be most needed.

The top 3 subzones for 2017 based on number of childcare age children are:

top_n(sf_mpszpop2017,n=3,Childcare_0_to_6)%>%
  select(SUBZONE_N,Childcare_0_to_6,COUNT_CHILDCARE)

The top 3 subzones for 2020 based on number of childcare age children are:

top_n(sf_mpszpop2020,n=3,Childcare_0_to_6)%>%
  select(SUBZONE_N,Childcare_0_to_6,COUNT_CHILDCARE)

It is interesting to note that the population of childcare age kids in Mathilda and Waterway East has risen, while Fernvale’s has dropped from the year 2017 to 2020. This indicates a shift in demand, in these subzones.

To understand if the supply is meeting demand, we will also need to identify the top subzones with the most number of childcare services. We will be assuming that it is preferable to send your children to a chilcare in your own subzone. As such more childcare in the subzone will indicate more supply of childcare services.

The top 3 subzones for 2017 based on number of childcare services are:

top_n(sf_mpszpop2017,n=3,COUNT_CHILDCARE)%>%
  select(SUBZONE_N,COUNT_CHILDCARE,Childcare_0_to_6)

The top 3 subzones for 2020 based on number of childcare services are:

top_n(sf_mpszpop2020,n=3,COUNT_CHILDCARE)%>%
  select(SUBZONE_N,COUNT_CHILDCARE,Childcare_0_to_6)

Tampines East has consistently had the most number of child care services and children of childcare age for both 2017 and 2020. This is an inital positive indication that at least for this subzone the demand is being met with a appropriate supply of childcare service. While maintaining their top ranks the other 2 subzones with the most childcare services (SENGKANG TOWN CENTRE,WOODLANDS EAST), are not the subzones with the most number of childcare age children.

4.2 Analytics Mapping of Temporal Supply and Demand Changes

Understanding changes in the population and the corresponding changes in the childcare centers, will enable us to understand if across the country supply is increasing with demand.

sf_mpszchange <- st_join(sf_mpszpop2017, sf_mpszpop2020, by = c("SUBZONE_N" = "SUBZONE_N"),suffix=c("_2017","_2020"))

sf_mpszchange<- sf_mpszchange%>%
  mutate(POP_CHANGE = sf_mpszchange$Childcare_0_to_6_2020-sf_mpszchange$Childcare_0_to_6_2017)%>%
  mutate(CARE_CHANGE = sf_mpszchange$COUNT_CHILDCARE_2020-sf_mpszchange$COUNT_CHILDCARE_2017)

popchange_map<-tm_shape(sf_mpszchange) +
   tm_fill(col="POP_CHANGE", style="jenks", palette = "Blues",title = "Change in Childcare Age Population") +
   tm_layout("Childcare Age Population Change by Subzone, 2017 to 2020",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Childcare data from Early Childhood 
             Development Agency (ECDA)",position = c("left", "bottom"))

carechange_map<-tm_shape(sf_mpszchange) + 
  tm_fill(col = "CARE_CHANGE", style="jenks", palette="Blues", title = "Change in Number of Childcare") +
  tm_layout("Childcare Supply Change by Subzone, 2017 to 2020",legend.position = c("right", "bottom")) +
  tm_borders(alpha = 0.5) +
  tmap_style("white")+
  tm_credits("Source: Planning Area Sub-zone boundary MP2014 from Urban Redevelopment 
             Authorithy (URA)\n and Childcare data from Early Childhood 
             Development Agency (ECDA)",position = c("left", "bottom"))

tmap_arrange(popchange_map, carechange_map, asp=1, ncol=2)

From the map you can observe that there has been huge decreases in population in many subzones, however we can also observe that there are increases in populations in a few subzones. When compared with the change of child care services,it looks like there is close to proportional increase in these subzones as well. Lets take a closer look at the map.

Similarly to the earlier maps, the maps shaded blue are, the child care age population changes while the map in green are the childcare services changes. We can see that on the most part the changes do correspond however in some regions (e.g. East Region) an increase in number of children does not always correspond with a similar increase in childcare services.

4.2.1 Supply vs Demand over the Years

To understand if supply is indeed meeting demand we have to see if the available childcare services have the ability to hold all the children in each subzone. Based on the ECDA fact sheet as of May 2020, there was a total of 139,259 childcare slots available for enrollment in 2017 and a total of 169,785 childcare slots available for enrollment in 2020.

As such by taking a simple average of the number of slots by the number of available childcare services per year we can get the average childcare slots per year. As such for the year 2017 there where 106 slots (rounded down) per childcare center and there where 129 slots (rounded down) per childcare center in 2020.

Based on these figures we can categorize the subzones to see whether or not demand is greater than supply or vice versa.

getDemandSupply <- function(population,childcare,no_slots){
  if(population>0 & childcare>0){
    childcarePerSubzone = childcare*no_slots
    if(population>childcarePerSubzone){
      return("Demand > Supply")
    }
    else if(population==childcarePerSubzone){
      return("Demand = Supply")
    }
    else if(population<childcarePerSubzone){
      return("Demand < Supply")
    }
  }
  else{
    if(population>0 & childcare==0){
      return("No Supply")
    }
    else if(population==0 & childcare>0){
      return("No Demand")
    }
    else{
      return("NA (No Supply & Demand)")
    }
  }
}

sf_mpszpop2017 <- sf_mpszpop2017 %>%
  rowwise()%>%
  mutate(DEMAND_SUPPLY = getDemandSupply(Childcare_0_to_6,COUNT_CHILDCARE,round(139259/sum(sf_mpszpop2017$COUNT_CHILDCARE))))

sf_mpszpop2020 <- sf_mpszpop2020 %>%
  rowwise()%>%
  mutate(DEMAND_SUPPLY = getDemandSupply(Childcare_0_to_6,COUNT_CHILDCARE,round(169785/sum(sf_mpszpop2020$COUNT_CHILDCARE))))

We can then utilise these categories to plot them into a map to better visualise the differences over the years.

sf_mpszpop2017 <- st_as_sf(sf_mpszpop2017)
demandsupply_map2017 <- tm_shape(sf_mpszpop2017) +
  tm_polygons(col = "DEMAND_SUPPLY", palette = "Accent",
              title="Supply vs Demand",legend.hist=TRUE) +
   tm_layout(main.title = "Supply vs Demand of Childcare Services, 2017", title.size = 1.5,main.title.position="center",legend.outside = TRUE, legend.outside.position = "right") +
  tm_borders(alpha = 0.5) +
  tmap_style("white")

demandsupply_map2017

Based on the map, for the year 2017 there was generally more demand than supply, with r length(which(sf_mpszpop2017$DEMAND_SUPPLY == "Demand > Supply")) subzones having more demand than supply. This is indicative that more attention will have to be paid in these regions to see if the number of child cares can be increased. Lets take alook if there was increase in 2020.

sf_mpszpop2020 <- st_as_sf(sf_mpszpop2020)
demandsupply_map2020 <- tm_shape(sf_mpszpop2020) +
  tm_polygons(col = "DEMAND_SUPPLY", palette = "Accent",
              title="Supply vs Demand",legend.hist=TRUE) +
   tm_layout(main.title = "Supply vs Demand of Childcare Services, 2020", title.size = 1.5,main.title.position="center",legend.outside = TRUE, legend.outside.position = "right") +
  tm_borders(alpha = 0.5) +
  tmap_style("white")

demandsupply_map2020

The slight increase in supply coupled with the drop in demand helped to lower the number of subzones that had more demand than supply. In 2020 there where r length(which(sf_mpszpop2017$DEMAND_SUPPLY == "Demand > Supply")) subzones with “Demand > Supply”. There are still more subzones that are facing a demand problem, this is indicative that more centers may need to be made available. However it is also interesting to note that based on the ECDA fact sheet not all available slots where taken up in 2020 as well. So it may be that not all parents desire to send their kids to childcare or it maybe that there is no available childcare near them to send their kids. As such more focused on the ground interviews and studies will need to be embarked on

5 Spatial Point Pattern Analyis

Further analysis will be carried out with spatial pattern analysis, this will allow us to understand how various events may interact and affect each other.

5.1 Data Preperation

5.1.1 Converting to PPP format

For the purpose of the analysis as we will be utilising spatstat for our pattern analysis we will have to convert our current data into ppp format.

#sf data.frame to SpatialPointsDataFrame
spdf_childcare2017<-as_Spatial(sf_childcare2017_3414)
spdf_childcare2020<-as_Spatial(sf_childcare2020_3414)

#SpatialPointsDataFrame to generic sp format
sp_childcare2017 <- as(spdf_childcare2017, "SpatialPoints")
sp_childcare2020 <- as(spdf_childcare2020, "SpatialPoints")

#SpatialPoints to spatstat’s ppp format
ppp_childcare2017 <- as(sp_childcare2017, "ppp")
ppp_childcare2020 <- as(sp_childcare2020, "ppp")

5.1.2 Handling duplicated points

Lets check if there are any duplicated point events and if so how many for each each event.

sum(multiplicity(ppp_childcare2017) > 1)

## [1] 85

sum(multiplicity(ppp_childcare2020) > 1)

## [1] 128

It looks like both data set has duplicated points, we will make use of the jittering approach to deal with this. This will add small perturbation to the points to ensure they dont occupy the exact same space.

ppp_childcare2017 <- rjitter(ppp_childcare2017, retry=TRUE, nsim=1, drop=TRUE)
ppp_childcare2020 <- rjitter(ppp_childcare2020, retry=TRUE, nsim=1, drop=TRUE)

Now lets check if the duplicates have been handled.

any(duplicated(ppp_childcare2017))

## [1] FALSE

any(duplicated(ppp_childcare2020))

## [1] FALSE

All clear! We can move on to creating an owin now.

5.1.3 Creating an owin

We will be utilising an owin of the singapore subzone boundaries, so that we can confine our analyis.

# convert to sp data.frame
spdf_mpszpop2017 <- as_Spatial(sf_mpszpop2017) 
spdf_mpszpop2020 <- as_Spatial(sf_mpszpop2020) 

#Convert to generic sp format
sp_mpszpop2017 <- as(spdf_mpszpop2017, "SpatialPolygons")
sp_mpszpop2020 <- as(spdf_mpszpop2020, "SpatialPolygons")

#Convert to owin
owin_mpszpop2017 <- as(sp_mpszpop2017, "owin")
owin_mpszpop2020 <- as(sp_mpszpop2020, "owin")

Lets now combine our childcare point with our owin. We will use the following code chunk to do so.

childcareSG_ppp2017 = ppp_childcare2017[owin_mpszpop2017]
childcareSG_ppp2020 = ppp_childcare2020[owin_mpszpop2020]

5.2 Exploratory Data Analysis

Lets take a look at the spatial patterns in the year 2017.

plot(childcareSG_ppp2017,main="Childcare Services in 2017", border="lightgrey")

Based on a visual analysis of the plot we can observe some clustering in the north, north-east and western regions. Lets take a lookg at child care services in 2020.

plot(childcareSG_ppp2020,main="Childcare Services in 2020", border="lightgrey")

Visually comparing with 2017,there seems to be greater signs of clustering in 2020. With clusters deepening in existing regions. Let see if statistically there is clustering as well.

5.2.1 Quadrat Analaysis of Childcare Services in 2017 and 2020

To better understand if the childcare locations are clustered or random we utilise the following test hypothesis:

Ho = The distribution of childcare services are randomly distributed
H1 = The distribution of childcare services are not randomly distributed
The 95% confident interval will be used.
We will utilise the Conditional monte carlo test of Complete Spatial Randomness

qt_childcare2017 <- quadrat.test(childcareSG_ppp2017, 
                   nx = 20, ny = 15, method="M",
             nsim=999)

qt_childcare2020 <- quadrat.test(childcareSG_ppp2020, 
                   nx = 20, ny = 15, method="M",
             nsim=999)

qt_childcare2017

## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  childcareSG_ppp2017
## X2 = 2199, p-value = 0.002
## alternative hypothesis: two.sided
## 
## Quadrats: 193 tiles (irregular windows)

qt_childcare2020

## 
##  Conditional Monte Carlo test of CSR using quadrat counts
##  Test statistic: Pearson X2 statistic
## 
## data:  childcareSG_ppp2020
## X2 = 2628.6, p-value = 0.002
## alternative hypothesis: two.sided
## 
## Quadrats: 193 tiles (irregular windows)

For both years we can see that the p-value of 0.002 is below our established a-values of 0.05. Therefore we are able to reject the hypothesis that point patterns are randomly distributed for both years. The monte carlo test gives a large x2 value for both years. Therefore it shows that for both years the childcare services spatial pattern exhibits clustering.

5.3 2nd Order Spatial Point Pattern Analysis

We can further try to identify if there is significant clustering in a distance range. For this we will be performing the following test:

L-function

We will be running 200 simulations.

These are the test parameters:

*Ho = The distribution of childcare services are randomly distributed.

*H1= The distribution of childcare services are not randomly distributed.

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

For the purpose of this spatial point analysis, we will be focousing on the following subzones: Seng Kang, Bedok, Bukit Batok and Hougang.

5.3.1 Extrating Subzones

The below code chunk will extract the relevant study area.

sk2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "SENGKANG",]
bk2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "BEDOK",]
bb2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "BUKIT BATOK",]
hg2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "HOUGANG",]

sk2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "SENGKANG",]
bk2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "BEDOK",]
bb2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "BUKIT BATOK",]
hg2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "HOUGANG",]

# convert to generic sp format
sk2017 = as(sk2017, "SpatialPolygons")
bk2017 = as(bk2017, "SpatialPolygons")
bb2017 = as(bb2017, "SpatialPolygons")
hg2017 = as(hg2017, "SpatialPolygons")

sk2020 = as(sk2020, "SpatialPolygons")
bk2020 = as(bk2020, "SpatialPolygons")
bb2020 = as(bb2020, "SpatialPolygons")
hg2020 = as(hg2020, "SpatialPolygons")
# create owin
sk_owin2017 = as(sk2017, "owin")
bk_owin2017 = as(bk2017, "owin")
bb_owin2017 = as(bb2017 , "owin")
hg_owin2017 = as(hg2017 , "owin")

sk_owin2020 = as(sk2020, "owin")
bk_owin2020 = as(bk2020, "owin")
bb_owin2020 = as(bb2020 , "owin")
hg_owin2020 = as(hg2020 , "owin")

Now lets comnine the child care points with the study area.

childcareSG_sk_ppp2017 = ppp_childcare2017[sk_owin2017]
childcareSG_bk_ppp2017 = ppp_childcare2017[bk_owin2017]
childcareSG_bb_ppp2017 = ppp_childcare2017[bb_owin2017]
childcareSG_hg_ppp2017 = ppp_childcare2017[hg_owin2017]

childcareSG_sk_ppp2020 = ppp_childcare2020[sk_owin2020]
childcareSG_bk_ppp2020 = ppp_childcare2020[bb_owin2020]
childcareSG_bb_ppp2020 = ppp_childcare2020[bb_owin2020]
childcareSG_hg_ppp2020 = ppp_childcare2020[hg_owin2020]

Lets quickly visulise the graph and ensure all is okay.

par(mfrow=c(2,4))
plot(childcareSG_sk_ppp2017,main="Seng Kang 2017")
plot(childcareSG_bk_ppp2017,main="Bedok 2017")
plot(childcareSG_bb_ppp2017,main="Bukit Batok 2017")
plot(childcareSG_hg_ppp2017,main="Hougang 2017")
plot(childcareSG_sk_ppp2020,main="Seng Kang 2020")
plot(childcareSG_bk_ppp2020,main="Bedok 2020")
plot(childcareSG_bb_ppp2020,main="Bukit Batok 2020")
plot(childcareSG_hg_ppp2020,main="Hougang 2020")

Lets create the function to conduct the test.

test <- function(input){
  # RUNNING THE L FUNCTION TEST WITH RIPLEY CORRECTION
  L = Lest(input, correction = "Ripley")
  L_test.csr <- envelope(input, Lest, nsim = 199, rank = 1, glocal=TRUE)
  
  # create a plot
  plot(L_test.csr, . - r ~ r, xlab="d", ylab="L(d)-r")
}

5.3.2 Comparing Seng kang 2017 & 2020

test(childcareSG_sk_ppp2017)

## 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.

test(childcareSG_sk_ppp2020)

## 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.

For both years, we can only reject the null hypothesis only beyond 100 meters and since its above the L(theo), it is a sign that child care services in seng kang have signs of clustering at distance above 100 meters

5.3.3 Comparing Bedok 2017 & 2020

test(childcareSG_bk_ppp2017)

## 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.

test(childcareSG_bk_ppp2020)

## 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.

For bedok, in the year 2017 we are not able to reject the null hypothesis as the point pattern of childcare services lie within the envelope, which means it is not statistically significant to reject the null hypothesis.

However in the year 2020 we can see that beyond 400 meters we can reject the null hypothesis, as the pattern exits the envelope, and as it it above the L(theo) line we can conclude that there are clusters. This can indicate that there are more childcare services that are being set up in close proximity to each other in bedok from the year 2017-2020.

5.3.4 Comparing Bukit Batok 2017 & 2020

test(childcareSG_bb_ppp2017)

## 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.

test(childcareSG_bb_ppp2020)

## 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.

For Bukit batok, in the year 2017 and 2020 we can reject the null hypothesis as the point pattern largely lies above the enevelope. However in the year 2017 the pattern exited the envolpe at under 400 meters but in 2020 it exits the envelop at above 400 meters.

This could be indicative that more chilcares where built at a larger distance from other childcare in bukit batok in the year 2020. As the observed L value is greater than the corresponding L(theo) value, spatial clustering is statistically significant.

5.3.5 Comparing Hougang 2017 & 2020

test(childcareSG_hg_ppp2017)

## 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.

test(childcareSG_hg_ppp2020)

## 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.

For the year 2017 we can reject the null hypothesis at around d=400+ while in the year 2020 we can rejuect the null hypothesis at around d=500+. As the observed L value is greater than the corresponding L(theo) value, spatial clustering is statistically significant for these distances.

5.4 Deriving Kernel Density Maps

We will now derive kernel density maps, we will be utilising Open Street Map as our base layer.

sk2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "SENGKANG",]
bk2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "BEDOK",]
bb2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "BUKIT BATOK",]
hg2017 = spdf_mpszpop2017[spdf_mpszpop2017@data$PLN_AREA_N == "HOUGANG",]

sk2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "SENGKANG",]
bk2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "BEDOK",]
bb2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "BUKIT BATOK",]
hg2020 = spdf_mpszpop2020[spdf_mpszpop2020@data$PLN_AREA_N == "HOUGANG",]

sk_bb2017 <- st_bbox(sf_mpszpop2017%>% filter(PLN_AREA_N == "SENGKANG"))
bk_bb2017 <- st_bbox(sf_mpszpop2017 %>% filter(PLN_AREA_N == "BEDOK"))
bb_bb2017 <- st_bbox(sf_mpszpop2017 %>% filter(PLN_AREA_N == "BUKIT BATOK"))
hg_bb2017 <- st_bbox(sf_mpszpop2017 %>% filter(PLN_AREA_N == "HOUGANG"))


sk_osm2017 <- read_osm(sk_bb2017, ext=1.1)
bk_osm2017 <- read_osm(bk_bb2017, ext=1.1)
bb_osm2017 <- read_osm(bb_bb2017, ext=1.1)
hg_osm2017 <- read_osm(hg_bb2017, ext=1.1)


sk_bb2020 <- st_bbox(sf_mpszpop2020%>% filter(PLN_AREA_N == "SENGKANG"))
bk_bb2020 <- st_bbox(sf_mpszpop2020 %>% filter(PLN_AREA_N == "BEDOK"))
bb_bb2020 <- st_bbox(sf_mpszpop2020 %>% filter(PLN_AREA_N == "BUKIT BATOK"))
hg_bb2020 <- st_bbox(sf_mpszpop2020 %>% filter(PLN_AREA_N == "HOUGANG"))


sk_osm2020 <- read_osm(sk_bb2020, ext=1.1)
bk_osm2020 <- read_osm(bk_bb2020, ext=1.1)
bb_osm2020 <- read_osm(bb_bb2020, ext=1.1)
hg_osm2020 <- read_osm(hg_bb2020, ext=1.1)

We will utilse a fixed bandwidth of 250meters.

deriveKDE <- function(title,osm,planningarea,ppp){
  # Rescale to km
  ppp_km <- rescale(ppp, 1000, "km")
  
  # compute kde
  kde <- density(ppp_km, sigma=0.25, edge=TRUE, kernel="gaussian")
  
  # convert to grid object
  kde_grid <- as.SpatialGridDataFrame.im(kde)
  
  # Convert into raster
  kde_raster <- raster(kde_grid)
  
  projection(kde_raster) <- crs("+init=EPSG:3414 +datum=WGS84 +units=km")
  
  # Plot on osm
  tm_shape(osm) + 
    tm_rgb() +
  tm_shape(planningarea) +
    tm_borders(col = "black", lwd = 2, lty="longdash") +
  tm_shape(kde_raster) + 
    tm_raster("v", alpha=0.5, palette = "YlOrRd") +
    tm_layout(legend.outside = TRUE, title=title)
}

5.4.1 Seng Kang

Lets take a look at the KDE for Seng kang.

sk2017KDE <- deriveKDE("Seng Kang, 2017",sk_osm2017, sk2017 , childcareSG_sk_ppp2017)
sk2020KDE <- deriveKDE("Seng Kang, 2020",sk_osm2020, sk2020 , childcareSG_sk_ppp2020)
tmap_arrange(sk2017KDE, sk2020KDE, asp=1, ncol=2)

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

Based on the kernal density map, from 2017 to 2020 there has not been that large of a difference and the greatest density of childcare services are in the north-east portion of sengkang.

5.4.2 Bedok

Lets take a look at the KDE for Bedok.

bk2017KDE <- deriveKDE("Bedok, 2017",bk_osm2017, bk2017 , childcareSG_bk_ppp2017)
bk2020KDE <- deriveKDE("Bedok, 2020",bk_osm2020, bk2020 , childcareSG_bk_ppp2020)
tmap_arrange(bk2017KDE, bk2020KDE, asp=1, ncol=2)

From 2017 to 2020 there is a visible reduction in the intensity of clusters. Instead in 2020 there is one visbily intense location right in the centre of bedok.

5.4.3 Bukit Batok

Lets take a look at the KDE for Bukit batok.

bb2017KDE <- deriveKDE("Bukit batok, 2017",bb_osm2017, bb2017 , childcareSG_bb_ppp2017)
bb2020KDE <- deriveKDE("Bukit batok, 2020",bb_osm2020, bb2020 , childcareSG_bb_ppp2020)
tmap_arrange(bb2017KDE, bb2020KDE, asp=1, ncol=2)

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

There is a visible increase in intensity of childcare services from the year 2017 to 2020. With intensity increasing in the west of bedok in 2020, indicating more childcare services being available in these locations.

5.4.4 Hougang

Lets take a look at the KDE for Hougang.

hg2017KDE <- deriveKDE("Hougang, 2017",hg_osm2017, hg2017 , childcareSG_hg_ppp2017)
hg2020KDE <- deriveKDE("Hougang, 2020",hg_osm2020, hg2020 , childcareSG_hg_ppp2020)
tmap_arrange(hg2017KDE, hg2020KDE, asp=1, ncol=2)

There is only a slight change in intensity from 2017 to 2020 in hougang, with the intense areas spreading out more. Indicating that there are more childcare services being set up in close proximity to other child care services from the year 2017 to 2020 in hougang.

5.5 Advantage of Kernal Density Map over Point Maps

There are few advantages of a Kernal density map over point map, they are:

Kernel Density map allows quick and easy identification of areas with high intensity/density of childcare.
Unlike point map which relies on the interpretation of the viewer which can subjective, kernal density utilises statistical analysis as such the density of an area is expressed interms of statistical confidence

6 Conclusion

In conclusion, from this study of the temporal changes in supply from 2017 to 2020 we have found that:

There is an increase in the number of child care services available over the years, while number childcare age children has dropped
We identified that even though there is an increase overall in the number of childcare services there are still subzones that are under-supplied.
There are signs of clustering across the subzones in-terms of number of childcare services and it is deepening from 2020-2017
When we look at select subzones we can observe that the distribution of subzones are not statistically random but instead they often cluster
By utilising Kernel density map we where able to identify that often the childcare services are most dense in certain regions (near malls, mrts). And the density often deepens from 2017 to 2020.

Take Home Exercise 1

Rajiv Abraham Xavier

9/16/2020 Updated: 24 September, 2020

1 Overview

1.1 The Problem Statement

1.2 The data

2 Setting up the environment

3 Geospatial Data Wrangling

3.1 Importing aspatial data

3.1.1 Data Prepeparation

3.2 Importing the spatial data

3.2.1 URA Master Plan 2017 Planning Subzone

3.2.1.1 CRS check and Reassignment

3.2.1.2 NA and Geometric Validity Checks and Correction

3.2.1.3 Georelational join with population data

3.2.2 Childcare Centers 2017

3.2.2.1 Data Cleaning

3.2.2.2 CRS check and Reassignment

3.2.2.3 NA and Geometric Validity Checks and Correction

3.2.3 Childcare Centers 2020

3.2.3.1 CRS check and Transforming

3.2.3.2 Data Cleaning

3.2.3.3 NA and Geometric Validity Checks and Correction

4 Section A: Exploring the change in supply and demand of childcare between year 2017 and 2020

4.1 Exploratory Data Analysis

4.1.1 Chloropleth Map of Childcare Demand by Subzone

4.1.2 Summary Statistics of Childcare Population

4.1.3 Chloropleth Map of Number of Childcare by Subzone

4.1.4 Summary Statistics of Childcare Services

4.1.5 Top 3 Subzones with Most Number of Childcare Age Children and Most Number of Childcare services

4.2 Analytics Mapping of Temporal Supply and Demand Changes

4.2.1 Supply vs Demand over the Years

5 Spatial Point Pattern Analyis

5.1 Data Preperation

5.1.1 Converting to PPP format

5.1.2 Handling duplicated points

5.1.3 Creating an owin

5.2 Exploratory Data Analysis

5.2.1 Quadrat Analaysis of Childcare Services in 2017 and 2020

5.3 2nd Order Spatial Point Pattern Analysis

5.3.1 Extrating Subzones

5.3.2 Comparing Seng kang 2017 & 2020

5.3.3 Comparing Bedok 2017 & 2020

5.3.4 Comparing Bukit Batok 2017 & 2020

5.3.5 Comparing Hougang 2017 & 2020

5.4 Deriving Kernel Density Maps

5.4.1 Seng Kang

5.4.2 Bedok

5.4.3 Bukit Batok

5.4.4 Hougang

5.5 Advantage of Kernal Density Map over Point Maps

6 Conclusion