age_cty <- age_cty %>% mutate(Age_gt_37 = Mean_Age > 37)

library(spdep) # A spatial statistics package for R
col_cty <- ifelse(age_cty$Age_gt_37,'Red','Blue') %>% factor()
joincount.test(col_cty, 
               nb2listw(poly2nb(age_cty),style='B'))[[1]]

    Join count test under nonfree sampling

data:  col_cty 
weights: nb2listw(poly2nb(age_cty), style = "B") 

Std. deviate for Blue = 2.6138, p-value = 0.004477
alternative hypothesis: greater
sample estimates:
Same colour statistic           Expectation              Variance 
             8.000000              3.747692              2.646642 

• Based on Blue-Blue counts there is reasonable evidence to reject a hypothesis of no clustering of Blue counties (p < 0.005)

• NOTES:
  • Uses the spdep library
  • col_cty is a factor variable for the colour of each county
  • poly2nb extracts a neighbour list from a polygon sf
  • nb2listw turns this into a listw object - a weighted neighbour list, although here all the weights are 1 (specified by style="B")
  • joincount.test does the test - it returns a list with 2 results (one for each colour), we select the first item (for \(\textrm{BB}\))

set.seed(19760604) # Reproducibility!!
joincount.mc(col_cty, 
               nb2listw(poly2nb(age_cty),style='B'),nsim=10000)[[1]]

    Monte-Carlo simulation of join-count statistic

data:  col_cty 
weights: nb2listw(poly2nb(age_cty), style = "B") 
number of simulations + 1: 10001 

Join-count statistic for Blue = 8, rank of observed statistic =
9868, p-value = 0.0133
alternative hypothesis: greater
sample estimates:
    mean of simulation variance of simulation 
              3.733800               2.664204 

• Based in Blue-Blue counts there is reasonable evidence (but weaker) to reject a hypothesis of no clustering of Blue counties (p < 0.05)

• NOTES:
  • This randomly permutes the Red and Blue colours among counties several times
  • Doesn’t rely on asymptotic normal approximation
  • joincount.mc then compares the actual join count with the simulated distribution
  • Then the experimental p-value is based on the rank of the actual count
  • joincount.mc works pretty much like joincount.test with an extra nsim (number of simulations) parameter.

• The Moran plot
• Works well if \(W_{ij} = J_{ij}/J_{i.}\) (ie rows sum to 1)
• The plot here is of \(x_i\) vs \(\sum_j W_{ij} x_i\)
• ie \(x_i\) vs. mean of neighbours \(x_j\)’s
• R code is below (style='W' for \(W_{ij}\) as above)

moran.plot(age_cty$Mean_Age, 
           nb2listw(poly2nb(age_cty),style='W'))

moran.test(age_cty$Mean_Age, 
               nb2listw(poly2nb(age_cty),style='W'))

    Moran I test under randomisation

data:  age_cty$Mean_Age  
weights: nb2listw(poly2nb(age_cty), style = "W")    

Moran I statistic standard deviate = 4.2282, p-value = 1.178e-05
alternative hypothesis: greater
sample estimates:
Moran I statistic       Expectation          Variance 
       0.49267372       -0.04000000        0.01587125 

• Based on similar analysis to join count
• \(E(I) = -1/(n-1)\)
• \(\textsf{Var}(I)\) very complicated - see Cliff & Ord!
• Again, evidence of spatial pattern

set.seed(19981115)
moran.mc(age_cty$Mean_Age, 
               nb2listw(poly2nb(age_cty),style='W'),nsim=1000)

    Monte-Carlo simulation of Moran I

data:  age_cty$Mean_Age 
weights: nb2listw(poly2nb(age_cty), style = "W")  
number of simulations + 1: 1001 

statistic = 0.49267, observed rank = 1001, p-value = 0.000999
alternative hypothesis: greater

• Monte-Carlo variation
• Avoids aymptotic normality assumption
• Again, evidence of spatial pattern

• Now consider point data - are points randomly distributed?
• New data, thanks to City of Toronto Open Data Portal
  • https://open.toronto.ca
• Here reported broken water main locations are downloaded in R

library(opendatatoronto)
watermain <- search_packages("Watermain Breaks") %>% 
  list_package_resources() %>% .[2,] %>% get_resource()
neighbourhoods <- list_package_resources(
    "https://open.toronto.ca/dataset/neighbourhoods/") %>%
  get_resource() %>% st_transform(3857)
cline <- search_packages("TCL") %>% 
  list_package_resources() %>% .[1,] %>% 
  get_resource() %>%
  st_transform(3857) %>% 
  filter(FCODE_DESC %in% c("Major Arterial",
                           "Minor Arterial"))

ggplot(neighbourhoods) + geom_sf() + 
  geom_sf(data=cline,col='darkred') + 
  geom_sf(data=watermain %>% filter(BREAK_YEAR==2014))

• This measures spatial clustering of the points, the lower this is, the more items are ‘clumped’
• Can compute this via the RANN package in R

nnd <- watermain %>% 
  st_coordinates() %>% 
  RANN::nn2(k=2) %>% .$nn.dists 
nnd <- nnd[,2]
glue::glue(
 "Mean NND is {round(mean(nnd),2)} for {length(nnd)} points.")

Mean NND is 46.79 for 35465 points.

• However we need to know how the observed value (in m) compares to a random process
• For a random process on a 2D plane, the probability density of \(2\pi\lambda d_{i}^2\) is \(\chi^2_2\) Where \(d_{i}\) is the distance from \(p_i\) to its nearest neighbour and \(\lambda\) is the density of points per unit area - (Skellam, 1951)

  ---

  Skellam, JG, 1951, “Random Dispersal in Theoretical Populations”, Biometrika, 38, 196-218

• This needs a lower tail test
• We can use some sf tools to estimate \(\lambda\)

lambda <- length(nnd)/sum(st_area(neighbourhoods))
test_stat <- 2 * pi * lambda * sum(nnd^2)
glue::glue(
 " Test statistic = {round(test_stat,1)}
  D.F. = {2*length(nnd)}
  p-value = {round(pchisq(test_stat,2*length(nnd)),2)}")

Test statistic = 37216.7
D.F. = 70930
p-value = 0

• Although, note ideally we would allow for edge effects
• … and think carefully about how \(\lambda\) is computed