Spatial Analysis: 10

Hotspot Analysis: FDR

0. Loading R packages and data

rm(list = ls())
library(spdep)
library(sf)
library(tmap)

# setwd("C:/Data")
TWPOP_sf <- st_read("./Data/Popn_TWN2.shp",options = "encoding=Big5" )

options:        encoding=Big5 
Reading layer `Popn_TWN2' from data source 
  `C:\Users\wenth\Desktop\WK10\Part2\Data\Popn_TWN2.shp' using driver `ESRI Shapefile'
Simple feature collection with 368 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -197572.5 ymin: 2295201 xmax: 606875.6 ymax: 2919551
Projected CRS: TWD97 / TM2 zone 121

ID<- TWPOP_sf$COUNTY_ID
Sel<- ID == "65000" | ID == "63000" 
NorthTW_sf<- TWPOP_sf[Sel,]

Popn<-NorthTW_sf$A15A64_CNT
NorthTW_sf$Density <- Popn * 10^6 / st_area(NorthTW_sf)

NorthTW_lyr <- tm_shape(NorthTW_sf) + 
  tm_polygons("Density", palette = "OrRd", style = "jenks", title = "Popn Density")

1. Spatial Neighbors

TWN_nb<-poly2nb(NorthTW_sf) #QUEEN = TRUE
TWN_nb_w<- nb2listw(TWN_nb, zero.policy=T) # default: style = "W" (row standardised)

2. Global Moran I

Density <- NorthTW_sf$Density
Density <- as.vector(Density)
moran.test(Density, listw=TWN_nb_w, zero.policy=T)

    Moran I test under randomisation

data:  Density  
weights: TWN_nb_w    

Moran I statistic standard deviate = 6.3433, p-value = 1.125e-10
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
      0.596361561      -0.025000000       0.009595443 

3. Local Moran

LISA.Popn <- localmoran(Density, TWN_nb_w, zero.policy=T, conditional = TRUE)

NorthTW_sf$z.li <- LISA.Popn[,4]
NorthTW_sf$pvalue <- LISA.Popn[,5]

tm_shape(NorthTW_sf) + 
  tm_polygons("z.li", palette = "OrRd", style = "jenks", title = "Local Moran")

[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual
variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named "brewer.or_rd"Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

tm_shape(NorthTW_sf) + tm_polygons("pvalue", palette = "-OrRd", style = "jenks", title = "p-value")

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual
variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-OrRd" is named "or_rd" (in
long format "brewer.or_rd")

4. Mapping Local Moran

quadr <- attr(LISA.Popn, "quadr")$mean
quadr <- factor(quadr, levels = c(levels(quadr), "NoSig"))
NorthTW_sf$Type <- quadr
head(NorthTW_sf)

Simple feature collection with 6 features and 18 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 300874.7 ymin: 2766756 xmax: 309745.8 ymax: 2776127
Projected CRS: TWD97 / TM2 zone 121
     TOWN_ID   TOWN COUNTY_ID COUNTY A0A14_CNT A0A14_M A0A14_F A15A64_CNT A15A64_M A15A64_F
221 63000010 松山區     63000 臺北市     30591   15780   14811     138248    64322    73926
222 63000020 信義區     63000 臺北市     27314   14197   13117     156271    74712    81559
223 63000030 大安區     63000 臺北市     46775   24357   22418     202471    93198   109273
224 63000040 中山區     63000 臺北市     28732   14741   13991     160976    74475    86501
225 63000050 中正區     63000 臺北市     25997   13552   12445     105921    49926    55995
226 63000060 大同區     63000 臺北市     17082    8964    8118      89519    43692    45827
    A65UP_CNT A65UP_M A65UP_F INFO_TIME                       geometry          Density
221     37724   16800   20924   107Y06M MULTIPOLYGON (((307703.1 27... 15921.74 [1/m^2]
222     40712   18085   22627   107Y06M MULTIPOLYGON (((307788.7 27... 13916.57 [1/m^2]
223     60297   26917   33380   107Y06M MULTIPOLYGON (((304591.5 27... 17868.67 [1/m^2]
224     40148   17643   22505   107Y06M MULTIPOLYGON (((305699 2776... 11514.05 [1/m^2]
225     27834   12696   15138   107Y06M MULTIPOLYGON (((302203.6 27... 14198.26 [1/m^2]
226     21977   10018   11959   107Y06M MULTIPOLYGON (((302217.9 27... 18861.47 [1/m^2]
        z.li       pvalue      Type
221 1.365617 0.1720592912 High-High
222 1.244017 0.2134934671 High-High
223 1.898417 0.0576411094 High-High
224 2.153648 0.0312678035 High-High
225 3.434886 0.0005928042 High-High
226 1.995359 0.0460037835 High-High

table(NorthTW_sf$Type)

  Low-Low  High-Low  Low-High High-High     NoSig 
       21         0         7        13         0 

signif <- 0.05
quadr[LISA.Popn[, 5]> signif] <- "NoSig"
NorthTW_sf$Type <- quadr

colors<- c( 'High-High' ='red', 'Low-Low'='blue', 
            'High-Low'='lightpink', 'Low-High'='skyblue2', 'NoSig'='grey')
tm_shape(NorthTW_sf) + tm_polygons("Type", palette = colors)

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of the visual
variable `fill` namely 'palette' (rename to 'values') to fill.scale = tm_scale(<HERE>).

5. FDR correction

LISA.Popn <- localmoran(Density, TWN_nb_w, zero.policy=T)

NorthTW_sf$pvalue.adj <- p.adjust(LISA.Popn[,5], method="fdr") 

tm_shape(NorthTW_sf) + 
  tm_polygons("pvalue.adj", palette = "-OrRd", style = "jenks", title = "p-value(FDR correction)")

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual
variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-OrRd" is named "or_rd" (in
long format "brewer.or_rd")

min(NorthTW_sf$pvalue.adj)

[1] 0.001073422

min(NorthTW_sf$pvalue) * 41

[1] 0.001073422

max(NorthTW_sf$pvalue.adj)

[1] 0.9288335

max(NorthTW_sf$pvalue)

[1] 0.9288335

p.adjustSP

NorthTW_sf$pvalue.adj2.1 <- p.adjust(LISA.Popn[,5], method="fdr")
NorthTW_sf$pvalue.adj2.2 <- p.adjustSP(LISA.Popn[,5], TWN_nb, method="fdr") 

head(NorthTW_sf)

Simple feature collection with 6 features and 21 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 300874.7 ymin: 2766756 xmax: 309745.8 ymax: 2776127
Projected CRS: TWD97 / TM2 zone 121
     TOWN_ID   TOWN COUNTY_ID COUNTY A0A14_CNT A0A14_M A0A14_F A15A64_CNT A15A64_M A15A64_F
221 63000010 松山區     63000 臺北市     30591   15780   14811     138248    64322    73926
222 63000020 信義區     63000 臺北市     27314   14197   13117     156271    74712    81559
223 63000030 大安區     63000 臺北市     46775   24357   22418     202471    93198   109273
224 63000040 中山區     63000 臺北市     28732   14741   13991     160976    74475    86501
225 63000050 中正區     63000 臺北市     25997   13552   12445     105921    49926    55995
226 63000060 大同區     63000 臺北市     17082    8964    8118      89519    43692    45827
    A65UP_CNT A65UP_M A65UP_F INFO_TIME                       geometry          Density
221     37724   16800   20924   107Y06M MULTIPOLYGON (((307703.1 27... 15921.74 [1/m^2]
222     40712   18085   22627   107Y06M MULTIPOLYGON (((307788.7 27... 13916.57 [1/m^2]
223     60297   26917   33380   107Y06M MULTIPOLYGON (((304591.5 27... 17868.67 [1/m^2]
224     40148   17643   22505   107Y06M MULTIPOLYGON (((305699 2776... 11514.05 [1/m^2]
225     27834   12696   15138   107Y06M MULTIPOLYGON (((302203.6 27... 14198.26 [1/m^2]
226     21977   10018   11959   107Y06M MULTIPOLYGON (((302217.9 27... 18861.47 [1/m^2]
        z.li       pvalue      Type pvalue.adj pvalue.adj2.1 pvalue.adj2.2
221 1.365617 0.1720592912     NoSig 0.32065595    0.32065595   1.000000000
222 1.244017 0.2134934671     NoSig 0.33666277    0.33666277   1.000000000
223 1.898417 0.0576411094     NoSig 0.19798039    0.19798039   0.345846656
224 2.153648 0.0312678035 High-High 0.18313999    0.18313999   0.218874624
225 3.434886 0.0005928042 High-High 0.01215249    0.01215249   0.004149629
226 1.995359 0.0460037835 High-High 0.19212781    0.19212781   0.276022701

tm_shape(NorthTW_sf) + 
  tm_polygons("pvalue.adj2.1", palette = "-OrRd", style = "jenks", title = "p-value(FDR correction)")

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual
variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.
[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "-OrRd" is named "or_rd" (in
long format "brewer.or_rd")

6. Mapping LISA Map (FDR)

signif <- 0.05
quadr[ NorthTW_sf$pvalue.adj > signif ] <- "NoSig"
NorthTW_sf$Type <- quadr

colors <- c('High-High' ='red', 'Low-Low'='blue', 
            'High-Low'='lightpink', 'Low-High'='skyblue2', NoSig='grey')
tm_shape(NorthTW_sf) + tm_polygons("Type", palette = colors)

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of the visual
variable `fill` namely 'palette' (rename to 'values') to fill.scale = tm_scale(<HERE>).

7. Local Gi*

TWN_nb_in<-include.self(TWN_nb); summary(TWN_nb_in)

Neighbour list object:
Number of regions: 41 
Number of nonzero links: 243 
Percentage nonzero weights: 14.45568 
Average number of links: 5.926829 
Link number distribution:

 3  4  5  6  7  8  9 10 12 
 3  7  6 12  7  3  1  1  1 
3 least connected regions:
278 292 296 with 3 links
1 most connected region:
231 with 12 links

TWN_nb_in_w<- nb2listw(TWN_nb_in, zero.policy=T)

LG <-localG(Density, TWN_nb_in_w); LG

 [1]  1.7233446  1.5358246  2.2955770  2.2575385  3.5680055  2.4243832  4.4002452  1.5935843
 [9] -0.2011367  0.1873412  0.2039434 -1.7879950  2.6015773  3.0486377  2.4917352  2.9034051
[17]  1.1389338 -0.8962848  0.4607859 -0.9741356 -1.4210311 -1.4498460 -1.6708723 -1.8162000
[25]  0.4273101  1.0063072  0.4332790 -0.3465971 -1.3485541 -1.0850525 -2.2394723 -2.0938802
[33] -1.7500732 -1.5259770 -1.6278155 -2.1806525 -2.0660263 -1.5440033 -1.9083605 -1.4421784
[41] -1.9553767
attr(,"internals")
               G*i     E(G*i)       V(G*i)     Z(G*i) Pr(z != E(G*i))
 [1,] 0.0423710270 0.02439024 1.088612e-04  1.7233446    8.482623e-02
 [2,] 0.0421929468 0.02439024 1.343658e-04  1.5358246    1.245814e-01
 [3,] 0.0483415035 0.02439024 1.088612e-04  2.2955770    2.170008e-02
 [4,] 0.0458835729 0.02439024 9.064358e-05  2.2575385    2.397445e-02
 [5,] 0.0583601259 0.02439024 9.064358e-05  3.5680055    3.597090e-04
 [6,] 0.0496854224 0.02439024 1.088612e-04  2.4243832    1.533442e-02
 [7,] 0.0662836240 0.02439024 9.064358e-05  4.4002452    1.081286e-05
 [8,] 0.0365111192 0.02439024 5.785193e-05  1.5935843    1.110292e-01
 [9,] 0.0226255011 0.02439024 7.698039e-05 -0.2011367    8.405917e-01
[10,] 0.0263448971 0.02439024 1.088612e-04  0.1873412    8.513931e-01
[11,] 0.0257598513 0.02439024 4.509962e-05  0.2039434    8.383977e-01
[12,] 0.0057349205 0.02439024 1.088612e-04 -1.7879950    7.377682e-02
[13,] 0.0491590607 0.02439024 9.064358e-05  2.6015773    9.279615e-03
[14,] 0.0511385254 0.02439024 7.698039e-05  3.0486377    2.298815e-03
[15,] 0.0481132874 0.02439024 9.064358e-05  2.4917352    1.271208e-02
[16,] 0.0580454283 0.02439024 1.343658e-04  2.9034051    3.691288e-03
[17,] 0.0362734865 0.02439024 1.088612e-04  1.1389338    2.547308e-01
[18,] 0.0158569921 0.02439024 9.064358e-05 -0.8962848    3.701007e-01
[19,] 0.0291979250 0.02439024 1.088612e-04  0.4607859    6.449522e-01
[20,] 0.0094131116 0.02439024 2.363842e-04 -0.9741356    3.299892e-01
[21,] 0.0095636952 0.02439024 1.088612e-04 -1.4210311    1.553077e-01
[22,] 0.0053413032 0.02439024 1.726227e-04 -1.4498460    1.471015e-01
[23,] 0.0084823832 0.02439024 9.064358e-05 -1.6708723    9.474691e-02
[24,] 0.0005279264 0.02439024 1.726227e-04 -1.8162000    6.933970e-02
[25,] 0.0288486498 0.02439024 1.088612e-04  0.4273101    6.691535e-01
[26,] 0.0376117070 0.02439024 1.726227e-04  1.0063072    3.142678e-01
[27,] 0.0281917673 0.02439024 7.698039e-05  0.4332790    6.648121e-01
[28,] 0.0198364444 0.02439024 1.726227e-04 -0.3465971    7.288940e-01
[29,] 0.0066721362 0.02439024 1.726227e-04 -1.3485541    1.774802e-01
[30,] 0.0101341783 0.02439024 1.726227e-04 -1.0850525    2.778984e-01
[31,] 0.0061480308 0.02439024 6.635347e-05 -2.2394723    2.512520e-02
[32,] 0.0001187687 0.02439024 1.343658e-04 -2.0938802    3.627065e-02
[33,] 0.0041040510 0.02439024 1.343658e-04 -1.7500732    8.010568e-02
[34,] 0.0009286651 0.02439024 2.363842e-04 -1.5259770    1.270156e-01
[35,] 0.0074061780 0.02439024 1.088612e-04 -1.6278155    1.035640e-01
[36,] 0.0016380678 0.02439024 1.088612e-04 -2.1806525    2.920913e-02
[37,] 0.0004416416 0.02439024 1.343658e-04 -2.0660263    3.882600e-02
[38,] 0.0006515132 0.02439024 2.363842e-04 -1.5440033    1.225875e-01
[39,] 0.0044790683 0.02439024 1.088612e-04 -1.9083605    5.634464e-02
[40,] 0.0054420449 0.02439024 1.726227e-04 -1.4421784    1.492521e-01
[41,] 0.0017242502 0.02439024 1.343658e-04 -1.9553767    5.053863e-02
attr(,"cluster")
 [1] High High High High High High High Low  Low  Low  Low  Low  High High High High High
[18] Low  Low  Low  Low  Low  Low  Low  Low  High Low  Low  Low  Low  Low  Low  Low  Low 
[35] Low  Low  Low  Low  Low  Low  Low 
Levels: Low High
attr(,"gstari")
[1] TRUE
attr(,"call")
localG(x = Density, listw = TWN_nb_in_w)
attr(,"class")
[1] "localG"

NorthTW_sf$LG <- LG

tm_shape(NorthTW_sf) + 
  tm_polygons("LG", palette = "OrRd", style = "jenks", title = "Local Gi*")

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_polygons()`: instead of `style = "jenks"`, use fill.scale =
`tm_scale_intervals()`.
ℹ Migrate the argument(s) 'style', 'palette' (rename to 'values') to
  'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the visual
variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'[cols4all] color palettes: use palettes from the R package cols4all. Run
`cols4all::c4a_gui()` to explore them. The old palette name "OrRd" is named "brewer.or_rd"Multiple palettes called "or_rd" found: "brewer.or_rd", "matplotlib.or_rd". The first one, "brewer.or_rd", is returned.

Hotspots (using Local Gi*, p < 0.05)

cluster_type <- attr(LG, "cluster")
cluster_type <- factor(cluster_type, 
                       levels = c(levels(cluster_type), "NoSig"))
NorthTW_sf$CLUSTER <- cluster_type
head(NorthTW_sf)

Simple feature collection with 6 features and 23 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: 300874.7 ymin: 2766756 xmax: 309745.8 ymax: 2776127
Projected CRS: TWD97 / TM2 zone 121
     TOWN_ID   TOWN COUNTY_ID COUNTY A0A14_CNT A0A14_M A0A14_F A15A64_CNT A15A64_M A15A64_F
221 63000010 松山區     63000 臺北市     30591   15780   14811     138248    64322    73926
222 63000020 信義區     63000 臺北市     27314   14197   13117     156271    74712    81559
223 63000030 大安區     63000 臺北市     46775   24357   22418     202471    93198   109273
224 63000040 中山區     63000 臺北市     28732   14741   13991     160976    74475    86501
225 63000050 中正區     63000 臺北市     25997   13552   12445     105921    49926    55995
226 63000060 大同區     63000 臺北市     17082    8964    8118      89519    43692    45827
    A65UP_CNT A65UP_M A65UP_F INFO_TIME                       geometry          Density
221     37724   16800   20924   107Y06M MULTIPOLYGON (((307703.1 27... 15921.74 [1/m^2]
222     40712   18085   22627   107Y06M MULTIPOLYGON (((307788.7 27... 13916.57 [1/m^2]
223     60297   26917   33380   107Y06M MULTIPOLYGON (((304591.5 27... 17868.67 [1/m^2]
224     40148   17643   22505   107Y06M MULTIPOLYGON (((305699 2776... 11514.05 [1/m^2]
225     27834   12696   15138   107Y06M MULTIPOLYGON (((302203.6 27... 14198.26 [1/m^2]
226     21977   10018   11959   107Y06M MULTIPOLYGON (((302217.9 27... 18861.47 [1/m^2]
        z.li       pvalue      Type pvalue.adj pvalue.adj2.1 pvalue.adj2.2       LG CLUSTER
221 1.365617 0.1720592912     NoSig 0.32065595    0.32065595   1.000000000 1.723345    High
222 1.244017 0.2134934671     NoSig 0.33666277    0.33666277   1.000000000 1.535825    High
223 1.898417 0.0576411094     NoSig 0.19798039    0.19798039   0.345846656 2.295577    High
224 2.153648 0.0312678035     NoSig 0.18313999    0.18313999   0.218874624 2.257538    High
225 3.434886 0.0005928042 High-High 0.01215249    0.01215249   0.004149629 3.568005    High
226 1.995359 0.0460037835     NoSig 0.19212781    0.19212781   0.276022701 2.424383    High

table(NorthTW_sf$CLUSTER)

  Low  High NoSig 
   28    13     0 

signif <- 0.05

pvalue <- attr(LG,"internals")[,5]

cluster_type[ pvalue > signif] <- "NoSig"
NorthTW_sf$CLUSTER <- cluster_type


colors <- c('High' ='red','Low' ='blue','NoSig' = 'lightgray')
tm_shape(NorthTW_sf) + tm_polygons("CLUSTER", palette = colors)

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of the visual
variable `fill` namely 'palette' (rename to 'values') to fill.scale = tm_scale(<HERE>).

Mapping Gi* (Bonferroni)

LGV <- NorthTW_sf$LG
qnorm(1-0.05/41, 0,1)

[1] 3.030805

chk <- LGV-3.03

quadrant <- vector(length=nrow(LISA.Popn))
levels(quadrant) <- c("Y","N")
quadrant[chk<0] <- "N" # Non-cluster
quadrant[chk>0] <- "Y" # Cluster
NorthTW_sf$cluster <- quadrant

colors <- c("lightgray","red")
tm_shape(NorthTW_sf) + tm_polygons("cluster", palette = colors)

── tmap v3 code detected ────────────────────────────────────────────────────────────────────
[v3->v4] `tm_tm_polygons()`: migrate the argument(s) related to the scale of the visual
variable `fill` namely 'palette' (rename to 'values') to fill.scale = tm_scale(<HERE>).