# environment prepare
library(tidyverse)

## ── Attaching packages ─────────────────────────────────────── tidyverse 1.3.1 ──

## ✓ ggplot2 3.3.5     ✓ purrr   0.3.4
## ✓ tibble  3.1.3     ✓ dplyr   1.0.7
## ✓ tidyr   1.1.3     ✓ stringr 1.4.0
## ✓ readr   2.0.1     ✓ forcats 0.5.1

## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()

library(data.table)

## 
## Attaching package: 'data.table'

## The following objects are masked from 'package:dplyr':
## 
##     between, first, last

## The following object is masked from 'package:purrr':
## 
##     transpose

library(sp)
library(sf)

## Linking to GEOS 3.7.2, GDAL 2.4.2, PROJ 5.2.0

library(table1)

## 
## Attaching package: 'table1'

## The following objects are masked from 'package:base':
## 
##     units, units<-

library(tm)

## Loading required package: NLP

## 
## Attaching package: 'NLP'

## The following object is masked from 'package:ggplot2':
## 
##     annotate

library(spatstat)

## Loading required package: spatstat.data

## Loading required package: spatstat.geom

## Registered S3 method overwritten by 'spatstat.geom':
##   method     from
##   print.boxx cli

## spatstat.geom 2.2-2

## 
## Attaching package: 'spatstat.geom'

## The following object is masked from 'package:data.table':
## 
##     shift

## Loading required package: spatstat.core

## Loading required package: nlme

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

## Loading required package: rpart

## spatstat.core 2.3-0

## Loading required package: spatstat.linnet

## spatstat.linnet 2.3-0

## 
## spatstat 2.2-0       (nickname: 'That's not important right now') 
## For an introduction to spatstat, type 'beginner'

library(here)

## here() starts at /Users/fangzeqiang/Github/Master-Dissertation

library(sp)
library(rgeos)

## rgeos version: 0.5-5, (SVN revision 640)
##  GEOS runtime version: 3.7.2-CAPI-1.11.2 
##  Linking to sp version: 1.4-2 
##  Polygon checking: TRUE

# library(maptools)
library(tmap)
library(sf)
library(geojson)

## 
## Attaching package: 'geojson'

## The following object is masked from 'package:graphics':
## 
##     polygon

library(geojsonio)

## Registered S3 method overwritten by 'geojsonsf':
##   method        from   
##   print.geojson geojson

## 
## Attaching package: 'geojsonio'

## The following object is masked from 'package:base':
## 
##     pretty

library(tmaptools)
library(RColorBrewer)
library(spdep)

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`

library(lubridate)

## 
## Attaching package: 'lubridate'

## The following objects are masked from 'package:rgeos':
## 
##     intersect, setdiff, union

## The following objects are masked from 'package:data.table':
## 
##     hour, isoweek, mday, minute, month, quarter, second, wday, week,
##     yday, year

## The following objects are masked from 'package:base':
## 
##     date, intersect, setdiff, union

Data Read

1.the index dataset

# read data in R
raw_data = read.csv(here::here("Dataset","df_hh_EnR_ttwaName_asset_not_drop.csv"))

# you can have a overview of this dataset
print("The number of rows is: ")

## [1] "The number of rows is: "

nrow(raw_data)

## [1] 3089

print("The number of columns is: ")

## [1] "The number of columns is: "

ncol(raw_data)

## [1] 9

print("70 of all varriables are:")

## [1] "70 of all varriables are:"

head(names(raw_data),n = 70)

## [1] "X"              "ttwa"           "ttwa_code"      "year"          
## [5] "firms"          "total_firms"    "entry_rate"     "hh"            
## [9] "average_assets"

spatial dataset

# import spatial dataset

eng_ttwa = st_read(here::here("Dataset","Spatial","eng_ttwa_boundary.geojson"))

## Reading layer `eng_ttwa_boundary' from data source 
##   `/Users/fangzeqiang/Github/Master-Dissertation/Dataset/Spatial/eng_ttwa_boundary.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 149 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -6.418524 ymin: 49.86474 xmax: 1.762942 ymax: 55.54107
## CRS:           4326

# plot the map
plot(st_geometry(eng_ttwa))

try to merge, but not a spatial type

dfm = merge(raw_data,eng_ttwa,by.x="ttwa_code",by.y="TTWA11CD",all = TRUE)

Data Prepare

# select year
df_1998 = raw_data %>% filter(year==1998)
df_2008 = raw_data %>% filter(year==2008)
df_2018 = raw_data %>% filter(year==2018)
# df_2017_2018 = raw_data %>% filter(year==2017 | year==2018)

# merge
dfm_1998 = merge(eng_ttwa,df_1998,by.x="TTWA11CD",by.y="ttwa_code",all.x = TRUE)
dfm_2008 = merge(eng_ttwa,df_2008,by.x="TTWA11CD",by.y="ttwa_code",all.x = TRUE)
dfm_2018 = merge(eng_ttwa,df_2018,by.x="TTWA11CD",by.y="ttwa_code",all.x = TRUE)
# dfm_2017_2018 = merge(eng_ttwa,df_2017_2018,by.x="TTWA11CD",by.y="ttwa_code",all.x = TRUE)
dfm_all = merge(eng_ttwa,raw_data,by.x="TTWA11CD",by.y="ttwa_code",all.x = TRUE)

# n/a -> 0
dfm_1998 = dfm_1998 %>% replace_na(list(entry_rate=0,hh=0,average_assets=0))
dfm_2008 = dfm_2008 %>% replace_na(list(entry_rate=0,hh=0,average_assets=0))
dfm_2018 = dfm_2018 %>% replace_na(list(entry_rate=0,hh=0,average_assets=0))
# dfm_2017_2018 = dfm_2017_2018 %>% replace_na(list(entry_rate=0,hh=0,average_assets=0))
dfm_all_dna = dfm_all %>% replace_na(list(entry_rate=0,hh=0,average_assets=0))

# projection
dfm_1998 = dfm_1998 %>% st_transform(.,27700)
dfm_2008 = dfm_2008 %>% st_transform(.,27700)
dfm_2018 = dfm_2018 %>% st_transform(.,27700)

cus_break = c(0.00,0.007,0.01,0.013,0.03,0.035,0.077,0.12153,0.14,0.16,0.32)
# library(viridis)
tmap_mode("plot")

## tmap mode set to plotting

Distributiom

df = dfm_1998
map_var = "entry_rate"
spatial_var = "ttwa"
legend_name = "Entry Rate (1998)"
title_name = "The Distribution of the Entry Rate of the England Tech Clusters(ttwa) in 1998"

enR_1998_map = tm_shape( df )+ 
        tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("left","center"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
        tm_scale_bar(position=c("left", "center")
                    )+ 
        tm_fill(col = map_var,
              palette = "BuPu", n = 12,
              breaks = cus_break, 
              title=legend_name,
              legend.hist = TRUE,
              legend.hist.position = c("left", "center")
              )+ 
        tm_layout(main.title = "",main.title.size = .85,
                 frame=FALSE,
                 legend.outside =TRUE,
                 legend.hist.width = 1.2,
                 legend.hist.height = .25 )+
        tm_borders(col = "#D2D2D2", lwd = .5, lty = "solid", alpha = NA)
        
enR_1998_map

df = dfm_2008
map_var = "entry_rate"
spatial_var = "ttwa"
legend_name = "Entry Rate (2008)"
title_name = "The Distribution of the Entry Rate of the England Tech Clusters(ttwa) in 2008"


enR_2008_map = tm_shape( df )+ 
        tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("left","center"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
        tm_scale_bar(position=c("left", "center")
                    )+ 
        tm_fill(col = map_var,
              palette = "BuPu", n = 12,
              breaks = cus_break, 
              title=legend_name,
              legend.hist = TRUE,
              legend.hist.position = c("left", "center")
              )+ 
        tm_layout(main.title = "",main.title.size = .85,
                 frame=FALSE,
                 legend.outside =TRUE,
                 legend.hist.width = 1.2,
                 legend.hist.height = .25 )+
        tm_borders(col = "#D2D2D2", lwd = .5, lty = "solid", alpha = NA)
        
enR_2008_map

df = dfm_2018
map_var = "entry_rate"
spatial_var = "ttwa"
legend_name = "Entry Rate (2018)"
title_name = "The Distribution of the Entry Rate of the England Tech Clusters(ttwa) in 2018"


enR_2018_map = tm_shape( df )+ 
        tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("left","center"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
        tm_scale_bar(position=c("left", "center")
                    )+ 
        tm_fill(col = map_var,
              palette = "BuPu", n = 12,
              breaks = cus_break, 
              title=legend_name,
              legend.hist = TRUE,
              legend.hist.position = c("left", "center")
              )+ 
        tm_layout(main.title = "",main.title.size = .85,
                 frame=FALSE,
                 legend.outside =TRUE,
                 legend.hist.width = 1.2,
                 legend.hist.height = .25 )+
        tm_borders(col = "#D2D2D2", lwd = .5, lty = "solid", alpha = NA)
        
enR_2018_map

enR_tmap = tmap_arrange(enR_1998_map,enR_2008_map,enR_2018_map,ncol = 1)

enR_tmap

# enR_tmap %>% tmap_save(.,here("Img","map_entry_rate_1998_to_2018.png"),dpi = 300)

Add the Density Variables (1998, 2008 & 2018)

dfm_1998 = dfm_1998 %>% 
    mutate(area = st_area(geometry))%>%
    mutate(area_km = units::set_units(area,km^2))%>%
    mutate(density = firms / area_km)%>%
    select(TTWA11CD,ttwa,year,entry_rate,hh,average_assets,density,firms,area_km)

dfm_2008 = dfm_2008 %>% 
    mutate(area = st_area(geometry))%>%
    mutate(area_km = units::set_units(area,km^2))%>%
    mutate(density = firms / area_km)%>%
    select(TTWA11CD,ttwa,year,entry_rate,hh,average_assets,density,firms,area_km)

dfm_2018 = dfm_2018 %>% 
    mutate(area = st_area(geometry))%>%
    mutate(area_km = units::set_units(area,km^2))%>%
    mutate(density = firms / area_km)%>%
    select(TTWA11CD,ttwa,year,entry_rate,hh,average_assets,density,firms,area_km)

dfm_1998 = dfm_1998 %>% replace_na(list(density=0))
dfm_2008 = dfm_2008 %>% replace_na(list(density=0))
dfm_2018 = dfm_2018 %>% replace_na(list(density=0))

Centroids and neighbour list

# projection

dfm_1998 = dfm_1998 %>% st_transform(.,27700)
dfm_2008 = dfm_2008 %>% st_transform(.,27700)
dfm_2018 = dfm_2018 %>% st_transform(.,27700)

coordsW_1998 = dfm_1998%>%
   st_centroid()%>%
   st_geometry()

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

coordsW_2008 = dfm_2008%>%
   st_centroid()%>%
   st_geometry()

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

coordsW_2018 = dfm_2018%>%
   st_centroid()%>%
   st_geometry()

## Warning in st_centroid.sf(.): st_centroid assumes attributes are constant over
## geometries of x

Ettwa_nb_1998 = dfm_1998 %>% 
  poly2nb(., queen=T)

Ettwa_nb_2008 = dfm_2008 %>% 
  poly2nb(., queen=T)

Ettwa_nb_2018 = dfm_2018 %>% 
  poly2nb(., queen=T)

df = dfm_2018
coordsW = coordsW_2018
nb = Ettwa_nb_2018

plot(df$geometry)
plot(nb,coordsW,col="blue",add=T)

# plot(coordsW,axes=TRUE)

#Create a spatial weights object from these weights
Ettwa_1998.wl = nb2listw(Ettwa_nb_1998, style="W",zero.policy = TRUE)

Ettwa_2008.wl = nb2listw(Ettwa_nb_2008, style="W",zero.policy = TRUE)

Ettwa_2018.wl = nb2listw(Ettwa_nb_2018, style="W",zero.policy = TRUE)

Spatial autocorrelation: Global Moran’I Index

Global Moran’s Index

In order to calculate the three time periods global Moran’s I

Morans’I in entry rate in 1998

I_Ettwa_enR_1998 = dfm_1998 %>%
  pull(entry_rate) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_1998.wl,zero.policy = TRUE)

I_Ettwa_enR_2008 = dfm_2008 %>%
  pull(entry_rate) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_2008.wl,zero.policy = TRUE)

I_Ettwa_enR_2018 = dfm_2018 %>%
  pull(entry_rate) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_2018.wl,zero.policy = TRUE)

Morans’I in entry rate in 2008

I_Ettwa_enR_1998 %>%  head()

## $statistic
## Moran I statistic standard deviate 
##                           2.497747 
## 
## $p.value
## [1] 0.006249272
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.129084083      -0.006849315       0.002961799 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_1998.wl  n reduced by no-neighbour observations\n  \n"

print("---------------------------------------------------------------")

## [1] "---------------------------------------------------------------"

head(I_Ettwa_enR_2008)

## $statistic
## Moran I statistic standard deviate 
##                           1.797589 
## 
## $p.value
## [1] 0.03612104
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.090494213      -0.006849315       0.002932467 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_2008.wl  n reduced by no-neighbour observations\n  \n"

print("---------------------------------------------------------------")

## [1] "---------------------------------------------------------------"

head(I_Ettwa_enR_2018)

## $statistic
## Moran I statistic standard deviate 
##                           1.600701 
## 
## $p.value
## [1] 0.05472156
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.079332556      -0.006849315       0.002898754 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_2018.wl  n reduced by no-neighbour observations\n  \n"

I_Ettwa_den_1998 = dfm_1998 %>%
  pull(density) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_1998.wl,zero.policy = TRUE)

I_Ettwa_den_2008 = dfm_2008 %>%
  pull(density) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_2008.wl,zero.policy = TRUE)

I_Ettwa_den_2018 = dfm_2018 %>%
  pull(density) %>% 
  as.vector() %>% 
  moran.test(.,Ettwa_2018.wl,zero.policy = TRUE)

head(I_Ettwa_den_1998)

## $statistic
## Moran I statistic standard deviate 
##                           8.755935 
## 
## $p.value
## [1] 1.012089e-18
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.401144651      -0.006849315       0.002171213 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_1998.wl  n reduced by no-neighbour observations\n  \n"

print("---------------------------------------------------------------")

## [1] "---------------------------------------------------------------"

head(I_Ettwa_den_2008)

## $statistic
## Moran I statistic standard deviate 
##                           6.933101 
## 
## $p.value
## [1] 2.05857e-12
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.316352833      -0.006849315       0.002173169 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_2008.wl  n reduced by no-neighbour observations\n  \n"

print("---------------------------------------------------------------")

## [1] "---------------------------------------------------------------"

head(I_Ettwa_den_2018)

## $statistic
## Moran I statistic standard deviate 
##                           4.475233 
## 
## $p.value
## [1] 3.816395e-06
## 
## $estimate
## Moran I statistic       Expectation          Variance 
##       0.179614715      -0.006849315       0.001736036 
## 
## $alternative
## [1] "greater"
## 
## $method
## [1] "Moran I test under randomisation"
## 
## $data.name
## [1] ".  \nweights: Ettwa_2018.wl  n reduced by no-neighbour observations\n  \n"

Table 2 Spatio-temporal dynamics cluster analysis

Year	Moran’s Index	p value	z-score
1998	0.118	0.00	2.4385
2008	0.062	0.08	1.3483
2018	0.056	0.10	1.2382

Table 3 Spatio-temporal density cluster analysis

Year	Moran’s Index	p value	z-score
1998 0.4359	0.00	10.07783
2008	0.3248	0.00	7.54488
2018	0.1891	0.00	4.982051

Global G test

Getis Ord General G…? This tells us whether high or low values are clustering. If G > Expected = High values clustering; if G < expected = low values clustering

# take 1998 for example
G_Ettwa_Global_enR = dfm_1998 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  globalG.test(., Ettwa_1998.wl,zero.policy = T)

## Warning in globalG.test(., Ettwa_1998.wl, zero.policy = T): Binary weights
## recommended (sepecially for distance bands)

G_Ettwa_Global_enR

## 
##  Getis-Ord global G statistic
## 
## data:  . 
## weights: Ettwa_1998.wl 
## 
## standard deviate = 2.5217, p-value = 0.00584
## alternative hypothesis: greater
## sample estimates:
## Global G statistic        Expectation           Variance 
##       7.239014e-03       6.849315e-03       2.388295e-08

This condition: G > E = high values means that high values tending to clustering in 1998

Local Moran’s Index

Data Process

# entry rate
I_Ettwa_Local_enR_1998 =  dfm_1998 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localmoran(., Ettwa_1998.wl,zero.policy = T)%>%
  as_tibble()

I_Ettwa_Local_enR_2008 =  dfm_2008 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localmoran(., Ettwa_2008.wl,zero.policy = T)%>%
  as_tibble()

I_Ettwa_Local_enR_2018 =  dfm_2018 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localmoran(., Ettwa_2018.wl,zero.policy = T)%>%
  as_tibble()


# density
I_Ettwa_Local_den_1998 <- dfm_1998 %>%
  pull(density) %>%
  as.vector()%>%
  localmoran(., Ettwa_1998.wl,zero.policy = T)%>%
  as_tibble()

I_Ettwa_Local_den_2008 <- dfm_2008 %>%
  pull(density) %>%
  as.vector()%>%
  localmoran(., Ettwa_2008.wl,zero.policy = T)%>%
  as_tibble()

I_Ettwa_Local_den_2018 <- dfm_2018 %>%
  pull(density) %>%
  as.vector()%>%
  localmoran(., Ettwa_2018.wl,zero.policy = T)%>%
  as_tibble()


#what does the output (the localMoran object) look like?
slice_head(I_Ettwa_Local_den_1998, n=5)

## # A tibble: 5 × 5
##   Ii         E.Ii         Var.Ii     Z.Ii       `Pr(z > 0)`
##   <localmrn> <localmrn>   <localmrn> <localmrn> <localmrn> 
## 1 0.01945781 -0.006756757 0.11814043 0.07626818 0.4696029  
## 2 0.01822073 -0.006756757 0.14237834 0.06619524 0.4736112  
## 3 0.08392046 -0.006756757 0.14237834 0.24031242 0.4050440  
## 4 0.09354003 -0.006756757 0.08784304 0.33840234 0.3675300  
## 5 0.12815488 -0.006756757 0.17873521 0.31911290 0.3748205

add Local Moran’s I Indicators

dfm_1998 <- dfm_1998 %>%
  mutate(entry_rate_I = as.numeric(I_Ettwa_Local_enR_1998$Ii))%>%
  mutate(entry_rate_Iz =as.numeric(I_Ettwa_Local_enR_1998$Z.Ii))%>%
  mutate(density_I =as.numeric(I_Ettwa_Local_den_1998$Ii))%>%
  mutate(density_Iz =as.numeric(I_Ettwa_Local_den_1998$Z.Ii))

dfm_1998 = dfm_1998 %>%   
  replace_na(list(entry_rate_I=0,entry_rate_Iz=0,density_I=0,density_Iz=0))

dfm_2008 <- dfm_2008 %>%
  mutate(entry_rate_I = as.numeric(I_Ettwa_Local_enR_2008$Ii))%>%
  mutate(entry_rate_Iz =as.numeric(I_Ettwa_Local_enR_2008$Z.Ii))%>%
  mutate(density_I =as.numeric(I_Ettwa_Local_den_2008$Ii))%>%
  mutate(density_Iz =as.numeric(I_Ettwa_Local_den_2008$Z.Ii))

dfm_2008 = dfm_2008 %>%   
  replace_na(list(entry_rate_I=0,entry_rate_Iz=0,density_I=0,density_Iz=0))

dfm_2018 <- dfm_2018 %>%
  mutate(entry_rate_I = as.numeric(I_Ettwa_Local_enR_2018$Ii))%>%
  mutate(entry_rate_Iz =as.numeric(I_Ettwa_Local_enR_2018$Z.Ii))%>%
  mutate(density_I =as.numeric(I_Ettwa_Local_den_2018$Ii))%>%
  mutate(density_Iz =as.numeric(I_Ettwa_Local_den_2018$Z.Ii))

dfm_2018 = dfm_2018 %>%   
  replace_na(list(entry_rate_I=0,entry_rate_Iz=0,density_I=0,density_Iz=0))

Map making for LISA

data points >2.58 or <-2.58 standard deviations away from the mean are significant at the 99% level (<1% chance that autocorrelation not present);
1.96 - <2.58 or <-1.96 to >-2.58 standard deviations are significant at the 95% level (<5% change that autocorrelation not present).
1.65 = 90% etc.

breaks1 = c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
MoranColours = rev(brewer.pal(8, "RdGy"))

# breaks1 = c(-5,-1,-0.7,-0.5,-0.1,0.1,0.5,0.7,1,5)
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)
tmap_mode("plot")

## tmap mode set to plotting

df = dfm_1998 %>% st_transform(.,27700)
legend_name = "Local Moran's I (1998)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
Iz_map_1998 = tm_shape(df) +
    tm_polygons("density_Iz",
        style="fixed",
        breaks=breaks1,
        palette=MoranColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              size = 2.5,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
    tm_scale_bar(position=c("right", "bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

Iz_map_1998

df = dfm_2008
legend_name = "Local Moran's I (2008)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
Iz_map_2008 = tm_shape(df) +
    tm_polygons("density_Iz",
        style="fixed",
        breaks=breaks1,
        palette=MoranColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              size = 2.5,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
    tm_scale_bar(position=c("right", "bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

Iz_map_2008

df = dfm_2018
legend_name = "Local Moran's I (2018)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
Iz_map_2018 = tm_shape(df) +
    tm_polygons("density_Iz",
        style="fixed",
        breaks=breaks1,
        palette=MoranColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              size = 2.5,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
    tm_scale_bar(position=c("right", "bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

Iz_map_2018

library(grid)
grid.newpage()
pushViewport(viewport(layout = grid.layout(
        nrow = 2,
        ncol = 3,
        heights = c(0.1, 0.9))))

grid.text("Entry Rate in Travel to Work Area in England", vp = viewport(layout.pos.row = 1, layout.pos.col = 1:3))
print(Iz_map_1998, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(Iz_map_2008, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
print(Iz_map_2018, vp = viewport(layout.pos.row = 2, layout.pos.col = 3))

Gi-Score

Data Process

calculate the Gi-Score

# Gi-Score

# entry rate
G_Ettwa_Local_enR_1998 =  dfm_1998 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localG(., Ettwa_1998.wl,zero.policy = T)

G_Ettwa_Local_enR_2008 =  dfm_2008 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localG(., Ettwa_2008.wl,zero.policy = T)

G_Ettwa_Local_enR_2018 =  dfm_2018 %>%
  pull(entry_rate) %>%
  as.vector()%>%
  localG(., Ettwa_2018.wl,zero.policy = T)

# density
G_Ettwa_Local_den_1998 =  dfm_1998 %>%
  pull(density) %>%
  as.vector()%>%
  localG(., Ettwa_1998.wl,zero.policy = T)

G_Ettwa_Local_den_2008 =  dfm_2008 %>%
  pull(density) %>%
  as.vector()%>%
  localG(., Ettwa_2008.wl,zero.policy = T)

G_Ettwa_Local_den_2018 =  dfm_2018 %>%
  pull(density) %>%
  as.vector()%>%
  localG(., Ettwa_2018.wl,zero.policy = T)

add columns in dataset

# entry rate
dfm_1998 = dfm_1998%>%
  mutate(entry_rate_G = as.numeric(G_Ettwa_Local_enR_1998))%>% 
  replace_na(list(entry_rate_G=0))

dfm_2008 = dfm_2008%>%
  mutate(entry_rate_G = as.numeric(G_Ettwa_Local_enR_2008))%>% 
  replace_na(list(entry_rate_G=0))

dfm_2018 = dfm_2018%>%
  mutate(entry_rate_G = as.numeric(G_Ettwa_Local_enR_2018))%>% 
  replace_na(list(entry_rate_G=0))

# density
dfm_1998 = dfm_1998%>%
  mutate(density_G = as.numeric(G_Ettwa_Local_den_1998))%>% 
  replace_na(list(density_G=0))

dfm_2008 = dfm_2008%>%
  mutate(density_G = as.numeric(G_Ettwa_Local_den_2008))%>% 
  replace_na(list(density_G=0))

dfm_2018 = dfm_2018%>%
  mutate(density_G = as.numeric(G_Ettwa_Local_den_2018))%>% 
  replace_na(list(density_G=0))

Map making

This research set the breaks manually based on the rule that entry rate >2.58 or <-2.58 standard deviations away from the mean are significant at the 99% level (<1% chance that autocorrelation not present); >1.96 - <2.58 or <-1.96 to >-2.58 standard deviations are significant at the 95% level (<5% change that autocorrelation not present). >1.65 = 90% etc.

# breaks1 = c(-5,-1,-0.7,-0.5,-0.1,0.1,0.5,0.7,1,5)
breaks1<-c(-1000,-2.58,-1.96,-1.65,1.65,1.96,2.58,1000)

tmap_mode("plot")

## tmap mode set to plotting

GIColours = rev(brewer.pal(8, "RdBu"))
df = dfm_1998  
legend_name = "Gi* (1998)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
G_map_1998 = tm_shape(df) +
    tm_polygons("density_G",
        style="fixed",
        breaks=breaks1,
        palette=GIColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              size = 2.5,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              just = NA)+ 
    tm_scale_bar(position=c("right", "bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

G_map_1998

df = dfm_2008
legend_name = "Gi* (2008)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
G_map_2008 = tm_shape(df) +
    tm_polygons("density_G",
        style="fixed",
        breaks=breaks1,
        palette=GIColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              size = 2.5,
              just = NA)+ 
    tm_scale_bar(position=c("right","bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

G_map_2008

df = dfm_2018
legend_name = "Gi* (2018)"
title_name = ""
# title_name = "Entry Rate in Travel to Work Area in England"
#now plot on an interactive map
G_map_2018 = tm_shape(df) +
    tm_polygons("density_G",
        style="fixed",
        breaks=breaks1,
        palette=GIColours,
        midpoint=NA,
        title=legend_name)+
      tm_compass( north = 0,
              type = "4star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              position = c("right","top"),
              bg.color = NA,
              bg.alpha = NA,
              size = 2.5,
              just = NA)+ 
    tm_scale_bar(position=c("right","bottom"))+
    tm_layout(main.title = title_name,
              main.title.size = .85,
              legend.text.size = 0.5)

G_map_2018

arrange the map

grid.newpage()
pushViewport(viewport(layout = grid.layout(
        nrow = 2,
        ncol = 3,
        heights = c(0.1, 0.9))))
grid.text("Entry Rate in Travel to Work Area in England", vp = viewport(layout.pos.row = 1, layout.pos.col = 1:3))
print(G_map_1998, vp = viewport(layout.pos.row = 2, layout.pos.col = 1))
print(G_map_2008, vp = viewport(layout.pos.row = 2, layout.pos.col = 2))
print(G_map_2018, vp = viewport(layout.pos.row = 2, layout.pos.col = 3))

# tmap_save(here("Img","G_score_enR.png"),dpi=300)
# tmap_arrange(G_map_1998,G_map_2008,G_map_2018,nrow = 1)

`{r} # st_write(obj = dfm_1998, dsn = here("Dataset","Spatial","dfm_1998_output.geojson")) # # st_write(obj = dfm_2008, dsn = here("Dataset","Spatial","dfm_2008_output.geojson")) # # st_write(obj = dfm_2018, dsn = here("Dataset","Spatial","dfm_2018_output.geojson")) #`

dfm_1998 %>% head(
)

## Simple feature collection with 6 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 155558.7 ymin: 11430.46 xmax: 448105.9 ymax: 484656.3
## CRS:           EPSG:27700
##    TTWA11CD                    ttwa year  entry_rate        hh average_assets
## 1 E30000004                Barnsley 1998 0.013351135 0.1200000       19405.25
## 2 E30000018                Bradford 1998 0.006211180 0.3388430       29856.56
## 3 E30000029                 Halifax 1998 0.012698413 0.1250000      167935.50
## 4 E30000039                 Skipton 1998 0.004444444 1.0000000      209782.00
## 5 E30000046 Dorchester and Weymouth 1998 0.009630819 0.3333333       85320.17
## 6 E30000051    Falmouth and Helston 1998 0.009404389 0.3333333      904365.00
##                 density firms          area_km                       geometry
## 1 0.0291974984 [1/km^2]    10  342.4951 [km^2] MULTIPOLYGON (((436695.4 41...
## 2 0.0319872061 [1/km^2]    11  343.8875 [km^2] MULTIPOLYGON (((403773.1 45...
## 3 0.0329725143 [1/km^2]    12  363.9395 [km^2] MULTIPOLYGON (((397952.9 43...
## 4 0.0008482608 [1/km^2]     1 1178.8828 [km^2] MULTIPOLYGON (((380072.8 48...
## 5 0.0083329237 [1/km^2]     6  720.0354 [km^2] MULTIPOLYGON (((384177.4 79...
## 6 0.0084678515 [1/km^2]     3  354.2811 [km^2] MULTIPOLYGON (((166091 1758...
##   entry_rate_I entry_rate_Iz  density_I density_Iz entry_rate_G  density_G
## 1 -0.449882946   -1.11571534 0.01945781 0.07626818  -2.13006137  0.3547993
## 2  0.375249417    0.87501490 0.01822073 0.06619524  -1.02516678  0.1989258
## 3 -0.167621598   -0.36847345 0.08392046 0.24031242  -0.94513033  0.8063152
## 4 -0.037254138   -0.08928316 0.09354003 0.33840234   0.06880698 -0.4572062
## 5 -0.089814155   -0.16958532 0.12815488 0.31911290   0.96483411 -0.6340875
## 6  0.002315084    0.01300948 0.17439996 0.30171023  -0.01642535 -0.6097366

Study Area

# study area

eng_ttwa = st_read(here::here("Dataset","Spatial","eng_ttwa_boundary.geojson"))

## Reading layer `eng_ttwa_boundary' from data source 
##   `/Users/fangzeqiang/Github/Master-Dissertation/Dataset/Spatial/eng_ttwa_boundary.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 149 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -6.418524 ymin: 49.86474 xmax: 1.762942 ymax: 55.54107
## CRS:           4326

eng_ttwa = eng_ttwa

uk_ttwa = st_read(here::here("Dataset","Spatial","Travel_to_Work_Areas__December_2011__Generalised_Clipped_Boundaries_in_United_Kingdom.geojson"))

## Reading layer `Travel_to_Work_Areas__December_2011__Generalised_Clipped_Boundaries_in_United_Kingdom' from data source `/Users/fangzeqiang/Github/Master-Dissertation/Dataset/Spatial/Travel_to_Work_Areas__December_2011__Generalised_Clipped_Boundaries_in_United_Kingdom.geojson' 
##   using driver `GeoJSON'
## Simple feature collection with 228 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -8.649996 ymin: 49.86474 xmax: 1.762942 ymax: 60.86078
## CRS:           4326

title_name =  ""

tmap_options(check.and.fix = TRUE)

study_area = tm_shape(uk_ttwa) +
    tm_borders(alpha = 0.5, lwd = .15, lty = 1
        )+
    tm_fill(
        )+
    tm_shape(eng_ttwa)+
      tm_fill(col = "blue", 
            alpha = 0.3)+
      tm_borders(col = "blue", lwd = .25, lty = 1, alpha = 0.5)+
    tm_layout(main.title = title_name)+
    tm_compass( north = 0,
              type = "8star",
              show.labels = TRUE,
              cardinal.directions = c("N", "E", "S", "W"),
              lwd = 1,
              size = 2.5,
              position = c("left","top"),
              just = NA)+ 
    tm_scale_bar(position=c("right", "bottom"))

study_area

## Warning: The shape uk_ttwa is invalid. See sf::st_is_valid

## Warning: The shape eng_ttwa is invalid. See sf::st_is_valid

Morans-I-and-GWR

Franklin

18/08/2021

Data Read

Data Prepare

Distributiom

Add the Density Variables (1998, 2008 & 2018)

Centroids and neighbour list

Spatial autocorrelation: Global Moran’I Index

Global Moran’s Index

Global G test

Local Moran’s Index

Data Process

Map making for LISA

Gi-Score

Data Process

Map making

`{r} # st_write(obj = dfm_1998, dsn = here("Dataset","Spatial","dfm_1998_output.geojson")) # # st_write(obj = dfm_2008, dsn = here("Dataset","Spatial","dfm_2008_output.geojson")) # # st_write(obj = dfm_2018, dsn = here("Dataset","Spatial","dfm_2018_output.geojson")) #`

Study Area