Megan's IS415 Journey: Hands-On Exercise 5B

1.0 Overview

In Hands-On Exercise 5A, we learned about Network-constrained Spatial Point Pattern Analysis (NetSPPA). In this exercise, we’ll follow-up on what we’ve learned with Marked Point Pattern Analysis.

Firstly, what are marked point patterns? You might recall what point patterns are: they’re events occurring within a defined study region (usually a country or city) - from traffic incidents to locations of business establishments. Now, these events might have associated continuous or categorical measurements, such as volume of sales (continuous) or types of schools (categorical). These measurements are called marks - and naturally, events with said marks are called marked point patterns.

One more thing to note! Marked point patterns have first-order and second-order properties:

first-order properties: related to the intensity (i.e., density) of events and associated marks
second-order properties: related to the spatial dependence (i.e., spatial arrangement) of the events and associated marks

1.1 Problem Statement

For this exercise, our aim is to discover & investigate the spatial point processes of childecare centres by operators in Singapore. In other words:

Are the locations of childcare centre by different business groups (i.e. NT, PT, RC, ST) spatially independent?
(In relation to the question above) if our answer is NO, is there any phenomena of attraction/ repulsion?

2.0 Setup

2.1 Packages Used

The R packages we’ll be using today are: - rgdal: provides bindings to the Geospatial Data Analysis Library (GDAL) and used for projectoin/transforamtion operations - maptools: a set of tools for manipulating geographic data - raster: reads, writes, manipulates, analyses and models gridded spatial data (i.e. raster-based geographical data) - spatstat: used for point pattern analysis - tmap: used for creating thematic maps, such as choropleth and bubble maps

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

2.2 Data Used

The datasets used for this exercise are:

CHILDCARE, originally a .kml that has been converted to an ESRI shapefile, providing both location and attribute information of childcare centres
MP14_SUBZONE_WEB_PL, a polygon feature data in ESRI shapefile format from data.gov.sg, providing information of URA 2014 Master Plan Planning Subzone boundary data

2.3 Importing Data

Like what we’ve learned in Hands-On Ex 5A, we’ll be using readOGR() to import geospatial data and output a suitable Spatial vector object. In this case, we’ll be importing CHILDCARE and MP14_SUBZONE_WEB_PL as SpatialPointsDataFrame and SpatialPolygonsDataFrame respectively.

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

OGR data source with driver: ESRI Shapefile 
Source: "C:\IS415\IS415_msty\_posts\2021-09-23-hands-on-exercise-5b\data\geospatial", layer: "CHILDCARE"
with 1885 features
It has 1 fields

mpsz = readOGR(dsn = "data/geospatial", layer="MP14_SUBZONE_WEB_PL")

OGR data source with driver: ESRI Shapefile 
Source: "C:\IS415\IS415_msty\_posts\2021-09-23-hands-on-exercise-5b\data\geospatial", layer: "MP14_SUBZONE_WEB_PL"
with 323 features
It has 15 fields

2.4 Data Preparation

We can take a look at our data with str():

str(childcare)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1885 obs. of  1 variable:
  .. ..$ Type: chr [1:1885] "PT" "PT" "ST" "ST" ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:1885, 1:2] 11227 11783 11894 11961 12128 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  ..@ bbox       : num [1:2, 1:2] 11227 25524 44936 49308
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +to"| __truncated__
  .. .. ..$ comment: chr "PROJCRS[\"SVY21 / Singapore TM\",\n    BASEGEOGCRS[\"SVY21\",\n        DATUM[\"SVY21\",\n            ELLIPSOID["| __truncated__

As we are working with marked data, and we know that the values are categorical (different business groups), we need to ensure that the marked field is of factor data type. However, as seen from the output, our Type field is of chr data type, not factor! Let’s rectify that:

childcare@data$Type <- as.factor(childcare@data$Type)

DIY section: Using the skill you learned from previous step, check to ensure that Type field is in factor data type now. Now, let’s check to ensure that Type field is now of the factor data type:

str(childcare)

Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
  ..@ data       :'data.frame': 1885 obs. of  1 variable:
  .. ..$ Type: Factor w/ 4 levels "NT","PT","RC",..: 2 2 4 4 2 2 3 2 2 2 ...
  ..@ coords.nrs : num(0) 
  ..@ coords     : num [1:1885, 1:2] 11227 11783 11894 11961 12128 ...
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : NULL
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  ..@ bbox       : num [1:2, 1:2] 11227 25524 44936 49308
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "coords.x1" "coords.x2"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+proj=tmerc +lat_0=1.36666666666667 +lon_0=103.833333333333 +k=1 +x_0=28001.642 +y_0=38744.572 +ellps=WGS84 +to"| __truncated__
  .. .. ..$ comment: chr "PROJCRS[\"SVY21 / Singapore TM\",\n    BASEGEOGCRS[\"SVY21\",\n        DATUM[\"SVY21\",\n            ELLIPSOID["| __truncated__

2.5 Data Visualisation

Let’s visualise our geospatial data in both static and interactive forms to better help with our analysis:

tmap_mode("view")
tm_shape(mpsz) +
  tm_borders(alpha = 0.5) +
  tmap_options(check.and.fix = TRUE) +
tm_shape(childcare) +
  tm_dots(col = 'Type', size = 0.02)

tmap_mode("plot")

Alternatively, we can also create four small point maps by using tm_facets():

tm_shape(mpsz) +
  tm_borders(alpha = 0.5) +
tm_shape(childcare) +
  tm_dots(col = 'Type', 
          size = 0.5) +
tm_facets(by="Type")

3.0 Spatial Data Wrangling

We’ll be using spatstat functions extensively in this section for wrangling geospatial data, so here’s a handy list of functions:

3.1 Converting the SpatialPointsDataFrame into ppp format

Firstly, let’s convert our SpatialPointsDataFrame into ppp format with as.(x, “ppp”) or as.ppp(x). In doing so, the additional field in x data frame will become the marks of the point pattern z.

childcare_ppp <- as(childcare, "ppp")
plot(childcare_ppp)

From our graph, we can see that there are four sub-types in the marks list: NT, PT, RC and ST. Let’s also take a look at the summary statistics:

summary(childcare_ppp)

Marked planar point pattern:  1885 points
Average intensity 2.351049e-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

Multitype:
   frequency proportion    intensity
NT       138 0.07320955 1.721193e-07
PT      1183 0.62758620 1.475486e-06
RC       200 0.10610080 2.494482e-07
ST       364 0.19310340 4.539957e-07

Window: rectangle = [11226.55, 44936.07] x [25523.51, 49308.17] units
                    (33710 x 23780 units)
Window area = 801770000 square units

Now, we know that PT is the largest childcare operator in Singapore with a market share of 63%. This is followed by ST, RC and NT.

Notice that the summary also tells us, “Pattern contains duplicated points”. We’ll learn how to resolve this issue in the next section.

3.2 Avoiding duplicated spatial point event by using jittering method

Let’s use the jittering method to avoid duplicated spatial point event issues:

childcare_ppp_jit <- rjitter(childcare_ppp, retry=TRUE, nsim=1, drop=TRUE)

Are there any duplicated spatial point events? Let’s check:

any(duplicated(childcare_ppp_jit))

[1] FALSE

Yay! 😄 Our duplicated points issue has been resolved!

3.3 Creating owin: Data Context

When analysing spatial point patterns, it is a good practice to confine the analysis within a geographical area, such as the Singapore boundary. To do this, we use an object called owin (from the spatstat package) which is designed to represent this polygonal region.

“So, what’s wrong with our current geographical area?” Take a look at the map below:

Have you noticed any pecularities? Here’s my annotated version:

The distribution of settlements in Singapore are constrained: either through natural geography, such as the water catchment areas, or strategic locations, such as the areas in the Western region reserved for airports.

As such, it would be wiser for us to narrow down the study area - such as looking at specific planning areas, which would be more appropriate.

3.4 Extracting study area

For this analysis, we’ll focus on the Jurong West planning area:

jw = mpsz[mpsz@data$PLN_AREA_N == "JURONG WEST",]
plot(jw, main = "Jurong West")

3.5 Converting the spatial point data frame into generic sp format

Next, let’s convert these SpatialPolygonsDataFrame layers into generic spatialpolygons layers by using as.SpatialPolygons.tess(x):

jw_sp = as(jw, "SpatialPolygons")

#always a best practice to review the structure with str()
str(jw_sp)

Formal class 'SpatialPolygons' [package "sp"] with 4 slots
  ..@ polygons   :List of 9
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 12930 34858
  .. .. .. .. .. .. ..@ area   : num 2098176
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:105, 1:2] 14212 14156 14122 14085 14067 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 12930 34858
  .. .. .. ..@ ID       : chr "142"
  .. .. .. ..@ area     : num 2098176
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 11290 35200
  .. .. .. .. .. .. ..@ area   : num 1524551
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:93, 1:2] 11769 11774 11901 11925 11934 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 11290 35200
  .. .. .. ..@ ID       : chr "143"
  .. .. .. ..@ area     : num 1524551
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 15522 35189
  .. .. .. .. .. .. ..@ area   : num 1484296
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:95, 1:2] 15941 15935 15924 15922 15921 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 15522 35189
  .. .. .. ..@ ID       : chr "145"
  .. .. .. ..@ area     : num 1484296
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 13397 35812
  .. .. .. .. .. .. ..@ area   : num 1404537
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:65, 1:2] 13673 13657 12785 12776 12776 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 13397 35812
  .. .. .. ..@ ID       : chr "151"
  .. .. .. ..@ area     : num 1404537
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 14632 35241
  .. .. .. .. .. .. ..@ area   : num 1287950
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:57, 1:2] 14106 14097 14085 14070 14065 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 14632 35241
  .. .. .. ..@ ID       : chr "157"
  .. .. .. ..@ area     : num 1287950
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 14404 36445
  .. .. .. .. .. .. ..@ area   : num 906317
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:41, 1:2] 14220 14209 14157 14090 13766 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 14404 36445
  .. .. .. ..@ ID       : chr "171"
  .. .. .. ..@ area     : num 906317
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 15508 36890
  .. .. .. .. .. .. ..@ area   : num 1793464
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:173, 1:2] 16296 16297 16297 16297 16297 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 15508 36890
  .. .. .. ..@ ID       : chr "176"
  .. .. .. ..@ area     : num 1793464
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 13954 37533
  .. .. .. .. .. .. ..@ area   : num 1974943
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:44, 1:2] 13849 13735 13710 13616 13596 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 13954 37533
  .. .. .. ..@ ID       : chr "186"
  .. .. .. ..@ area     : num 1974943
  .. .. .. ..$ comment: chr "0"
  .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots
  .. .. .. ..@ Polygons :List of 1
  .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots
  .. .. .. .. .. .. ..@ labpt  : num [1:2] 12593 36300
  .. .. .. .. .. .. ..@ area   : num 2206305
  .. .. .. .. .. .. ..@ hole   : logi FALSE
  .. .. .. .. .. .. ..@ ringDir: int 1
  .. .. .. .. .. .. ..@ coords : num [1:106, 1:2] 12671 12711 12723 12726 12732 ...
  .. .. .. ..@ plotOrder: int 1
  .. .. .. ..@ labpt    : num [1:2] 12593 36300
  .. .. .. ..@ ID       : chr "192"
  .. .. .. ..@ area     : num 2206305
  .. .. .. ..$ comment: chr "0"
  ..@ plotOrder  : int [1:9] 9 1 8 7 2 3 4 5 6
  ..@ bbox       : num [1:2, 1:2] 10373 33982 16297 38489
  .. ..- attr(*, "dimnames")=List of 2
  .. .. ..$ : chr [1:2] "x" "y"
  .. .. ..$ : chr [1:2] "min" "max"
  ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
  .. .. ..@ projargs: chr "+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"
  .. .. ..$ comment: chr "PROJCRS[\"SVY21\",\n    BASEGEOGCRS[\"SVY21[WGS84]\",\n        DATUM[\"World Geodetic System 1984\",\n         "| __truncated__

3.6 Creating owin object

Now to actually convert these SpatialPolygons objects into owin objects:

jw_owin = as(jw_sp, "owin")
str(jw_owin)

List of 5
 $ type  : chr "polygonal"
 $ xrange: num [1:2] 10373 16297
 $ yrange: num [1:2] 33982 38489
 $ bdry  :List of 9
  ..$ :List of 2
  .. ..$ x: num [1:104] 14212 14240 14250 14250 14243 ...
  .. ..$ y: num [1:104] 35397 35520 35596 35689 35740 ...
  ..$ :List of 2
  .. ..$ x: num [1:92] 11769 11764 11760 11756 11688 ...
  .. ..$ y: num [1:92] 35484 35486 35489 35493 35582 ...
  ..$ :List of 2
  .. ..$ x: num [1:94] 15941 15944 15944 15945 15950 ...
  .. ..$ y: num [1:94] 34916 34979 34995 35012 35161 ...
  ..$ :List of 2
  .. ..$ x: num [1:64] 13673 13689 13704 13733 13747 ...
  .. ..$ y: num [1:64] 35226 35229 35233 35245 35252 ...
  ..$ :List of 2
  .. ..$ x: num [1:56] 14106 14195 14305 14329 14353 ...
  .. ..$ y: num [1:56] 34492 34512 34537 34542 34548 ...
  ..$ :List of 2
  .. ..$ x: num [1:40] 14220 14444 14878 14938 14945 ...
  .. ..$ y: num [1:40] 35850 35933 36095 36113 36145 ...
  ..$ :List of 2
  .. ..$ x: num [1:172] 16296 16296 16296 16296 16296 ...
  .. ..$ y: num [1:172] 36694 36695 36696 36697 36698 ...
  ..$ :List of 2
  .. ..$ x: num [1:43] 13849 14496 14639 14684 14723 ...
  .. ..$ y: num [1:43] 36652 37092 37190 37216 37242 ...
  ..$ :List of 2
  .. ..$ x: num [1:105] 12671 12648 12605 12637 12678 ...
  .. ..$ y: num [1:105] 35650 35702 35777 35827 35882 ...
 $ units :List of 3
  ..$ singular  : chr "unit"
  ..$ plural    : chr "units"
  ..$ multiplier: num 1
  ..- attr(*, "class")= chr "unitname"
 - attr(*, "class")= chr "owin"

3.7 Combining childcare points and the study area

Now, let’s extract the childcare data that is within our specified region (Jurong West):

childcare_jw_ppp = childcare_ppp_jit[jw_owin]

And once again checking for summary statistics:

summary(childcare_jw_ppp)

Marked planar point pattern:  114 points
Average intensity 7.765382e-06 points per square unit

Coordinates are given to 3 decimal places
i.e. rounded to the nearest multiple of 0.001 units

Multitype:
   frequency proportion    intensity
NT        16  0.1403509 1.089878e-06
PT        58  0.5087719 3.950808e-06
RC        12  0.1052632 8.174086e-07
ST        28  0.2456140 1.907287e-06

Window: polygonal boundary
9 separate polygons (no holes)
           vertices    area relative.area
polygon 1       104 2098180        0.1430
polygon 2        92 1524550        0.1040
polygon 3        94 1484300        0.1010
polygon 4        64 1404540        0.0957
polygon 5        56 1287950        0.0877
polygon 6        40  906317        0.0617
polygon 7       172 1793460        0.1220
polygon 8        43 1974940        0.1350
polygon 9       105 2206310        0.1500
enclosing rectangle: [10373.179, 16297.184] x [33981.5, 38488.61] 
units
                     (5924 x 4507 units)
Window area = 14680500 square units
Fraction of frame area: 0.55

DIY section: By referring to previous discussion, describe the content of the output. As we can see from our summary statistics output, PT is once again reigning in the Jurong West area with a 50% market share. However, since that is smaller than their national average of 63% market share, it can be said that their presence in Jurong West is smaller than that of other areas.

3.8 Plotting childcare points and the study area

Lastly, let’s ensure that the spatial point events are indeed contained within the study area by plotting the combined childcare point and the study area:

plot(childcare_jw_ppp)

4.0 Analysing Marked Point Patterns

4.1 First-order Spatial Point Patterns Analysis

There are two things we need to do to get our intiial visualisation: use density() of the spatstat package to compute the kernel density objects, and then plot() it out. We can keep them as two separate lines of code, or we can combine them, like so:

plot(density(split(rescale(childcare_jw_ppp, 1000))))

Question: Can you recall what is the purpose of rescale() and why it is used in our case? rescale() of the spatstat package is used for converting a point pattern dataset to another unit of length. In this context, we’re converting metres to kilometres: we’re telling the function that the new unit length is 1000 metres (or 1km), so the function will divide the old coordinate values by 1000 to obtain coordinates expressed in kilometres ✨

DIY: What observations can you draw from the figure above? At first glance, it might seem that NT has a high intensity (aka density of childcare centres) in Jurong West due to its gradient, which has a high amount of yellow (high-intensity), but notice the legend at the side: its range is from 0 to 2.5. Now compare this to PT’s range of up to 14: that’s a lot more congruent with our earlier analysis in 3.7, where we summarised that PT had an impressive 50% market share in Jurong West.

Surprisingly, a good portion of PT’s childcare centres are nestled in the Southeast tail of Jurong West - though they also have a middling portion in the central region, which seems to be the case for other business groups as well. In addition, NT seems to be the only operator with a ‘unified’ area - most of its centres are in the central region, with a few more expanding to the Eastern region and a consistent gradient throughout, as compared to the other operators where there tends be a a space of blues/purples (low-intensity) between their clusters.

Something interesting to note is that there across all the childcare operators, none of them have a childcare centre in the westmost part. A quick look at Google Maps tells us why: that’s hwere the Singapore Discovery Centre is!

Well, those are some good observations we’ve made! Let’s check if we’re correct - by using the intensity() function, we’ll reveal the density of childcare centres by operators:

intensity(rescale(childcare_jw_ppp, 1000))

       NT        PT        RC        ST 
1.0898782 3.9508083 0.8174086 1.9072868

Hooray 🎉 We were right on the money: the output reveals that childcare centres operated by PT has the highest density of 3.95 units per km square, followed by ST’s 1.91 units per km square, NT’s 1.09 unit per km square and RC’s 0.82 unit per km square.

4.2 Second-order Multi-type Point Patterns Analysis: Cross K-Function

Now, let’s analyse the relationship between PT (i.e. privately-operated childcare centres) and ST (i.e. PAP-run childcare centres) by using Kcross():

childcare_Kcross <- Kcross(childcare_jw_ppp, 
                           i="PT", j="ST",
                           correction='border')
plot(childcare_Kcross)

Hmm… it seems that the marked spatial point events are not spatially independent. We’ll need to confirm this observation statistically through a hypothesis test.

4.3 Performing CSR testing on the Cross K-Function

The hypothesis and test are as follows:

Ho = The distribution of PT childcare centres and ST childcare centres are spatially independent.
H1 = The distribution of PT childcare centres and ST childcare centres are NOT at spatially independent.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001 (i.e. at 99.9% confidence interval).

In order to perform the CSR test, we’ll use the envelope() function of our spatstat package.

childcare_Kcross.csr <- envelope(childcare_jw_ppp, Kcross, i="PT", j="ST", correction='border', nsim=999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.

Done.

plot(childcare_Kcross.csr, xlab="distance(m)", xlim=c(0,500))

Question: Why nsim=999 is used?

We’ll be running through this with 999 simulations of CSR!

Hmm… as we can see, there are signs that the distribution of childcare centres operated by PT and ST are not spatially independent. Unfortunately, we failed to reject the null hypothesis because the empirical k-cross line is within the envelope of the 99.9% confidence interval.

4.4 Second-order Multi-type Point Patterns Analysis: Cross L-Function

Let’s compute Cross L-function:

childcare_Lcross <- Lcross(childcare_jw_ppp, i="PT", j="ST", correction='border')
plot(childcare_Lcross, . -r ~ r, 
     xlab = "distance(m)", 
     xlim=c(0, 500))

DIY: With reference to discussion in the earlier section, what observation(s) can you draw from the plot above?

Once again, there are signs that the marked spatial point events are not spatially independent. We’ll need to perform a hypothesis test to confirm this observation statistically.

4.5 Performing CSR testing on the Cross L-Function

DIY: With reference to the example given in previous section, define the hypothesis null, hypothesis alternative and rejection criterion.

The hypothesis and test are as follows:

Ho = The distribution of PT childcare centres and ST childcare centres are spatially independent.
H1 = The distribution of PT childcare centres and ST childcare centres are NOT at spatially independent.

The null hypothesis will be rejected if p-value is smaller than alpha value of 0.001 (i.e. at 99.9% confidence interval).

Similar to Cross-K-Function, we can perform the CSR test by using the envelope() fucntion of our spatstat package:

childcare_Lcross.csr <- envelope(childcare_jw_ppp, Lcross, i="PT", j="ST", correction='border', nsim=999)

Generating 999 simulations of CSR  ...
1, 2, 3, ......10.........20.........30.........40.........50.........60
.........70.........80.........90.........100.........110.........120
.........130.........140.........150.........160.........170.........180
.........190.........200.........210.........220.........230.........240
.........250.........260.........270.........280.........290.........300
.........310.........320.........330.........340.........350.........360
.........370.........380.........390.........400.........410.........420
.........430.........440.........450.........460.........470.........480
.........490.........500.........510.........520.........530.........540
.........550.........560.........570.........580.........590.........600
.........610.........620.........630.........640.........650.........660
.........670.........680.........690.........700.........710.........720
.........730.........740.........750.........760.........770.........780
.........790.........800.........810.........820.........830.........840
.........850.........860.........870.........880.........890.........900
.........910.........920.........930.........940.........950.........960
.........970.........980.........990........ 999.

Done.

plot(childcare_Lcross.csr, . -r ~ r, xlab="distance(m)", xlim=c(0,500))

DIY: Interpret the analysis result and draw conclusion with reference to the statistical testing result.

Once again, the empirical k-cross line is within the envelope of the 99.9% confidence interval, so we fail to reject the null hypothesis.