Geospatial Analysis and Representation for Data Science Project

Introduction

This notebook aims to perform a spatial statistics analysis concerning ski pass prices in the area (the Dolomiti area in the province of Bolzano), investigating whether there are any relationships between them. Since we are unable to find a complete resource containing the costs of the ski pass for the Dolomiti area, the data used in this notebook was created in the file ski_area_cost.ipynb by consulting the price lists on the websites of the various resorts. The aim is therefore to perform a spatial statistics analysis to understand whether there is a relationship between the location of resort and their prices. In particular, we consider the representative points of each ski resort as coordinates of the ski area.

Set up

if(!require(rgdal)) {
    install.packages("rgdal")
    library(rgdal)
}

if(!require(spdep)) {
    install.packages("spdep")
    library(spdep)
}

if(!require(boot)) {
    install.packages("boot")
    library(boot)
}

Import data

Import South Tirol’s Dolomiti area

The data is build and saved in the ski_area_cost.ipynb file

area_dolomiti_BZ = readOGR(dsn=path.expand("data/area_dolomiti_BZ"), layer="area_dolomiti_BZ")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\elypa\OneDrive\Desktop\Geospatial Project\geospatial_project(DEFINITIVO)\data\area_dolomiti_BZ", layer: "area_dolomiti_BZ"
## with 1 features
## It has 1 fields

head(area_dolomiti_BZ@data)

area_dolomiti_BZ is the shape representing the Dolomiti municipalities in province of Bolzano, which is the location of the ski areas we will investigating in this section.

We can plot this Dolomiti area:

plot(area_dolomiti_BZ)

Import ski areas data

These data represent the coordinates and some information about the representative points of each ski resort and the prices of their ski passes. We computed the representative points in the dataset preprocessing in Python. The various ski areas originally had a Polygon shape, therefore it would have been necessary to calculate a reference point (centroid or representative point) to perform the spatial statistics analyses. We preferred the calculation of the representative point since, due to their shape, the centroids were prone to fall outside the resort area in some cases.

Obviously, these resorts are composed of several slopes and lifts, in fact the price of the ski pass includes access to the entire area to which it refers.

BZ_dolomiti_ski_areas = readOGR(dsn=path.expand("data/ski_slopes_BZ_dissolved"), layer="ski_slopes_BZ_dissolved")

## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\elypa\OneDrive\Desktop\Geospatial Project\geospatial_project(DEFINITIVO)\data\ski_slopes_BZ_dissolved", layer: "ski_slopes_BZ_dissolved"
## with 24 features
## It has 5 fields

BZ_dolomiti_ski_areas@data

We can plot the ski areas representative points:

par(mar=c(0,0,0,0))

plot(area_dolomiti_BZ)
points(BZ_dolomiti_ski_areas, pch=20, col="blue", cex=2)
text(
  coordinates(BZ_dolomiti_ski_areas), 
  labels = BZ_dolomiti_ski_areas$NAME_IT, 
  cex = 0.9, 
  col = "blue",
  pos = 3
  )

The columns of the dataset are:

names(BZ_dolomiti_ski_areas@data)

## [1] "NAME_IT"  "OBJECTID" "NAME_DE"  "AREA"     "prices"

• OBJECTID refers to the ski area’s id;
• NAME_IT refers to the italian name of the ski area;
• NAME_DE refers to the german name of the ski area;
• AREA refers to the area of each complete ski resort original shape;
• prices refers to the price of the ski pass required to access the ski lifts of the area.

Let’s look at the columns types:

str(BZ_dolomiti_ski_areas@data)

## 'data.frame':    24 obs. of  5 variables:
##  $ NAME_IT : chr  "Alpe di Siusi" "Antermoia" "Baranci" "Braies Vecchia" ...
##  $ OBJECTID: num  19831 19750 19749 19759 19781 ...
##  $ NAME_DE : chr  "Seiseralm" "Untermoi" "Haunold" "Altprags" ...
##  $ AREA    : num  2661013 83494 288836 273267 1518392 ...
##  $ prices  : num  62 50 61 61 52 50 62 62 50 62 ...

What is the price range for ski passes?

unique(BZ_dolomiti_ski_areas$prices)

## [1] 62.0 50.0 61.0 52.0 55.0 20.0 37.5

summary(BZ_dolomiti_ski_areas$prices)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   20.00   51.50   61.00   55.77   62.00   62.00

We have 7 unique prices, this means that some ski areas have the same ski pass prices.

How many observations we have for each price?

barplot(table(BZ_dolomiti_ski_areas$prices))

Where are the most expensive ski areas?

par(mar=c(0,0,0,0))

subset <- BZ_dolomiti_ski_areas[BZ_dolomiti_ski_areas$prices== max(BZ_dolomiti_ski_areas$prices),] 
head(subset@data)

plot(area_dolomiti_BZ)
points(subset, pch=20, col="red", cex=2)
text(
  coordinates(subset), 
  labels = subset$NAME_IT, 
  cex = 0.9, 
  col = "red",
  pos = 3
  )

Which ski area is the less expensive and where is it?

par(mar=c(0,0,0,0))

subset <- BZ_dolomiti_ski_areas[BZ_dolomiti_ski_areas$prices== min(BZ_dolomiti_ski_areas$prices),] 
head(subset@data)

plot(area_dolomiti_BZ)
points(subset, pch=20, col="green", cex=2)
text(
  coordinates(subset), 
  labels = subset$NAME_IT, 
  cex = 0.9, 
  col = "green",
  pos = 3
  )

Which are the most extensive resorts (in terms of area)?

par(mar=c(0,0,0,0))

sorted <- BZ_dolomiti_ski_areas[order(-BZ_dolomiti_ski_areas$AREA),]
head(sorted@data)

plot(area_dolomiti_BZ)
points(sorted, pch=20, col="red", cex=sorted$AREA/mean(sorted$AREA))
text(
  coordinates(sorted), 
  labels = sorted$NAME_IT, 
  cex = 0.7, 
  col = "blue",
  pos = 3
  )

There are clearly much larger ski areas than others.

---

Defining neighbours

By referring to the distances between the reference points (for us the representative points) we can then define the neighbourhood relations between the spatial units. Several definitions of proximity are possible: in our case k-nearest neighbours and critical neighbourhood cut-off will be applied; we will not deal with contiguity-based neighbourhood as our units do not share borders.

K-Nearest neighbours

The k-nearest neighbours criterion assumes that two spatial units are considered neighbours if their distance is equal, or less than equal, to the smallest possible distance that can be found between all observations.

This definition of neighbourhood ensures that each spatial unit has precisely the same number k of neighbours.

Let’s start with 1 neighbour (k = 1):

par(mar=c(0,0,0,0))

knn1BZ = knn2nb(knearneigh(BZ_dolomiti_ski_areas, k = 1, longlat = T)) #longlat = T takes into account the fact that the Earth is not flat

plot(area_dolomiti_BZ, border = "grey")
plot(knn1BZ, BZ_dolomiti_ski_areas, add = T)

Let’s now try with 2 neighbour (k = 2):

knn2BZ = knn2nb(knearneigh(BZ_dolomiti_ski_areas, k = 2, longlat = T))

plot(area_dolomiti_BZ, border = "grey")
plot(knn2BZ, BZ_dolomiti_ski_areas, add = T)

With k = 4:

knn4BZ = knn2nb(knearneigh(BZ_dolomiti_ski_areas, k = 4, longlat = T))

plot(area_dolomiti_BZ, border = "grey")
plot(knn4BZ, BZ_dolomiti_ski_areas, add = T)

With k = 6:

knn6BZ = knn2nb(knearneigh(BZ_dolomiti_ski_areas, k = 6, longlat = T))

plot(area_dolomiti_BZ, border = "grey")
plot(knn4BZ, BZ_dolomiti_ski_areas, add = T)

With k = 8:

knn8BZ = knn2nb(knearneigh(BZ_dolomiti_ski_areas, k = 8, longlat = T))

plot(area_dolomiti_BZ, border = "grey")
plot(knn8BZ, BZ_dolomiti_ski_areas, add = T)

Critical cut-off neighbourhood

The neighbourhood criterion of the critical cut-off implies that two spatial units are considered neighbours if their distance is equal, or less than equal, to a certain fixed distance representing a critical cut-off. This cut-off distance should not be less than a minimum value that ensures that each spatial unit has at least one neighbour. For this reason, we must first calculate the minimum threshold distance that allows all regions to have at least one neighbour. As we have seen in class, we can do this by using the list of k nearest neighbours with k = 1 and calculating the maximum distance that separates all pairs of spatial units.

all.linkedT = max(unlist(nbdists(knn1BZ, BZ_dolomiti_ski_areas, longlat = T))) 
all.linkedT

## [1] 0.07901053

We found that the minimum threshold distance is 0.07901, which is actually very small, but it is OK since we are working with a fairly small area. After this calculation we now know that the cut-off distance has to be greater than 0.07901.

We can try different values of the cut-off distance for different neighbourhood definitions:

dnb079 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.0791, longlat=T)
dnb079

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 42 
## Percentage nonzero weights: 7.291667 
## Average number of links: 1.75

plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.079") 
plot(dnb079, BZ_dolomiti_ski_areas, add = TRUE, col = "blue")

dnb089 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.0891, longlat=T)
plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.089") 
plot(dnb089, BZ_dolomiti_ski_areas, add = TRUE, col = "green")

dnb089

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 50 
## Percentage nonzero weights: 8.680556 
## Average number of links: 2.083333

dnb099 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.0991, longlat=T)
plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.099") 
plot(dnb099, BZ_dolomiti_ski_areas, add = TRUE, col = "green")

dnb099

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 50 
## Percentage nonzero weights: 8.680556 
## Average number of links: 2.083333

We have the same results with d = 0.089 and d = 0.99, thus we can try to increase the values further:

dnb109 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.1091, longlat=T)
plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.1091")
plot(dnb109, BZ_dolomiti_ski_areas, add = TRUE, col = "yellow")

dnb109

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 60 
## Percentage nonzero weights: 10.41667 
## Average number of links: 2.5

dnb159 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.1591, longlat=T)
plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.1591") 
plot(dnb159, BZ_dolomiti_ski_areas, add = TRUE, col = "orange")

dnb159

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 114 
## Percentage nonzero weights: 19.79167 
## Average number of links: 4.75

dnb199 = dnearneigh(BZ_dolomiti_ski_areas, 0, 0.1991, longlat=T)
plot(area_dolomiti_BZ, border = "grey", xlab = "", ylab = "", xlim = NULL)
title(main = "d nearest neighbours, d = 0.1991") 
plot(dnb199, BZ_dolomiti_ski_areas, add = TRUE, col = "black")

dnb199

## Neighbour list object:
## Number of regions: 24 
## Number of nonzero links: 148 
## Percentage nonzero weights: 25.69444 
## Average number of links: 6.166667

Obviously with the critical cut-off the links, and therefore the number of neighbours, can only teds to grow rapidly with increasing distance.

Defining spatial weights

Once the neighbourhood relationships between observations have been defined, the spatial weights matrix can be created. In our case we will calculate a standardised spatial weight matrix for each list of critical and k-nearest neighbours created earlier.

# knn
knn1BZ.listw = nb2listw(knn1BZ, style = "W")
knn2BZ.listw = nb2listw(knn2BZ, style = "W")
knn4BZ.listw = nb2listw(knn4BZ, style = "W")
knn6BZ.listw = nb2listw(knn6BZ, style = "W")
knn8BZ.listw = nb2listw(knn8BZ, style = "W")

# critical cut-off
dnb079.listw = nb2listw(dnb079,style="W")
dnb089.listw = nb2listw(dnb089,style="W")
dnb099.listw = nb2listw(dnb099,style="W")
dnb109.listw = nb2listw(dnb109,style="W")
dnb159.listw = nb2listw(dnb159,style="W")
dnb199.listw = nb2listw(dnb199,style="W")


# Example of output
str(dnb079.listw)

## List of 3
##  $ style     : chr "W"
##  $ neighbours:List of 24
##   ..$ : int 6
##   ..$ : int 22
##   ..$ : int [1:2] 21 24
##   ..$ : int 21
##   ..$ : int 14
##   ..$ : int 1
##   ..$ : int [1:3] 8 10 18
##   ..$ : int [1:2] 7 10
##   ..$ : int 23
##   ..$ : int [1:3] 7 8 18
##   ..$ : int 20
##   ..$ : int 23
##   ..$ : int [1:3] 15 17 19
##   ..$ : int 5
##   ..$ : int [1:2] 13 19
##   ..$ : int 24
##   ..$ : int [1:2] 13 19
##   ..$ : int [1:2] 7 10
##   ..$ : int [1:3] 13 15 17
##   ..$ : int [1:2] 11 22
##   ..$ : int [1:2] 3 4
##   ..$ : int [1:2] 2 20
##   ..$ : int [1:2] 9 12
##   ..$ : int [1:2] 3 16
##   ..- attr(*, "class")= chr "nb"
##   ..- attr(*, "region.id")= chr [1:24] "1" "2" "3" "4" ...
##   ..- attr(*, "call")= language dnearneigh(x = BZ_dolomiti_ski_areas, d1 = 0, d2 = 0.0791, longlat = T)
##   ..- attr(*, "dnn")= num [1:2] 0 0.0791
##   ..- attr(*, "bounds")= chr [1:2] "GE" "LE"
##   ..- attr(*, "nbtype")= chr "distance"
##   ..- attr(*, "sym")= logi TRUE
##  $ weights   :List of 24
##   ..$ : num 1
##   ..$ : num 1
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num 1
##   ..$ : num 1
##   ..$ : num 1
##   ..$ : num [1:3] 0.333 0.333 0.333
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num 1
##   ..$ : num [1:3] 0.333 0.333 0.333
##   ..$ : num 1
##   ..$ : num 1
##   ..$ : num [1:3] 0.333 0.333 0.333
##   ..$ : num 1
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num 1
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:3] 0.333 0.333 0.333
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:2] 0.5 0.5
##   ..$ : num [1:2] 0.5 0.5
##   ..- attr(*, "mode")= chr "binary"
##   ..- attr(*, "W")= logi TRUE
##   ..- attr(*, "comp")=List of 1
##   .. ..$ d: num [1:24] 1 1 2 1 1 1 3 2 1 3 ...
##  - attr(*, "class")= chr [1:2] "listw" "nb"
##  - attr(*, "region.id")= chr [1:24] "1" "2" "3" "4" ...
##  - attr(*, "call")= language nb2listw(neighbours = dnb079, style = "W")

The Moran’s I test of spatial autocorrelation

We now perform Moran’s I test of spatial autocorrelation with reference to the variable corresponding to the price of ski passes (BZ_dolomiti_ski_areas$prices). The inspection of these variable may in fact suggest the presence of some form of spatial dependence.

Let’s first plot the quantile distribution of the prices variable:

brks = round(quantile(BZ_dolomiti_ski_areas$prices), digits = 3)
colours = grey((length(brks):2)/length(brks))
plot(area_dolomiti_BZ, col = "#E29B50")

points(BZ_dolomiti_ski_areas, col = "black", pch = 19, cex=3.3) # add a layer to obtain a contour effect
points(BZ_dolomiti_ski_areas, col = colours[findInterval(BZ_dolomiti_ski_areas$prices, brks, all.inside = TRUE)], pch = 19, cex=3)
title(main = "Ski resrorts' skipass prices")

Finally we perform the global Moran’s I test to the prices variable with different specifications of the spatial weights matrix:

moran.test(BZ_dolomiti_ski_areas$prices, knn1BZ.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn1BZ.listw    
## 
## Moran I statistic standard deviate = 0.075958, p-value = 0.4697
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.02412437       -0.04347826        0.06492124

moran.test(BZ_dolomiti_ski_areas$prices, knn2BZ.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn2BZ.listw    
## 
## Moran I statistic standard deviate = 0.78107, p-value = 0.2174
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09565863       -0.04347826        0.03173283

moran.test(BZ_dolomiti_ski_areas$prices, knn4BZ.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn4BZ.listw    
## 
## Moran I statistic standard deviate = 1.135, p-value = 0.1282
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09820404       -0.04347826        0.01558247

moran.test(BZ_dolomiti_ski_areas$prices, knn6BZ.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn6BZ.listw    
## 
## Moran I statistic standard deviate = 1.7592, p-value = 0.03927
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.123520266      -0.043478261       0.009011412

moran.test(BZ_dolomiti_ski_areas$prices, knn8BZ.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn8BZ.listw    
## 
## Moran I statistic standard deviate = 1.6444, p-value = 0.05005
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.079633133      -0.043478261       0.005605112

moran.test(BZ_dolomiti_ski_areas$prices, dnb079.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb079.listw    
## 
## Moran I statistic standard deviate = 0.52792, p-value = 0.2988
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07441268       -0.04347826        0.04986890

moran.test(BZ_dolomiti_ski_areas$prices, dnb089.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb089.listw    
## 
## Moran I statistic standard deviate = 0.53756, p-value = 0.2954
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.06929266       -0.04347826        0.04400940

moran.test(BZ_dolomiti_ski_areas$prices, dnb099.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb099.listw    
## 
## Moran I statistic standard deviate = 0.53756, p-value = 0.2954
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.06929266       -0.04347826        0.04400940

moran.test(BZ_dolomiti_ski_areas$prices, dnb109.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb109.listw    
## 
## Moran I statistic standard deviate = 0.50384, p-value = 0.3072
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.05552170       -0.04347826        0.03860823

moran.test(BZ_dolomiti_ski_areas$prices, dnb159.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb159.listw    
## 
## Moran I statistic standard deviate = 1.3508, p-value = 0.08838
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12347987       -0.04347826        0.01527742

moran.test(BZ_dolomiti_ski_areas$prices, dnb199.listw, randomisation = FALSE)

## 
##  Moran I test under normality
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb199.listw    
## 
## Moran I statistic standard deviate = 1.3193, p-value = 0.09353
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09776138       -0.04347826        0.01146052

The test was performed under the assumption of normality since randomisation parameter is set on False. Let’s try under the assumption of randomization:

moran.test(BZ_dolomiti_ski_areas$prices, knn1BZ.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn1BZ.listw    
## 
## Moran I statistic standard deviate = 0.089743, p-value = 0.4642
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       -0.02412437       -0.04347826        0.04650862

moran.test(BZ_dolomiti_ski_areas$prices, knn2BZ.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn2BZ.listw    
## 
## Moran I statistic standard deviate = 0.92413, p-value = 0.1777
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09565863       -0.04347826        0.02266833

moran.test(BZ_dolomiti_ski_areas$prices, knn4BZ.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn4BZ.listw    
## 
## Moran I statistic standard deviate = 1.343, p-value = 0.08963
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.09820404       -0.04347826        0.01112930

moran.test(BZ_dolomiti_ski_areas$prices, knn6BZ.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn6BZ.listw    
## 
## Moran I statistic standard deviate = 2.0799, p-value = 0.01877
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.123520266      -0.043478261       0.006446833

moran.test(BZ_dolomiti_ski_areas$prices, knn8BZ.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: knn8BZ.listw    
## 
## Moran I statistic standard deviate = 1.9327, p-value = 0.02664
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.079633133      -0.043478261       0.004057547

moran.test(BZ_dolomiti_ski_areas$prices, dnb079.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb079.listw    
## 
## Moran I statistic standard deviate = 0.62499, p-value = 0.266
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.07441268       -0.04347826        0.03558070

moran.test(BZ_dolomiti_ski_areas$prices, dnb089.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb089.listw    
## 
## Moran I statistic standard deviate = 0.63624, p-value = 0.2623
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.06929266       -0.04347826        0.03141618

moran.test(BZ_dolomiti_ski_areas$prices, dnb099.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb099.listw    
## 
## Moran I statistic standard deviate = 0.63624, p-value = 0.2623
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.06929266       -0.04347826        0.03141618

moran.test(BZ_dolomiti_ski_areas$prices, dnb109.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb109.listw    
## 
## Moran I statistic standard deviate = 0.5965, p-value = 0.2754
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.05552170       -0.04347826        0.02754582

moran.test(BZ_dolomiti_ski_areas$prices, dnb159.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb159.listw    
## 
## Moran I statistic standard deviate = 1.5998, p-value = 0.05483
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.12347987       -0.04347826        0.01089191

moran.test(BZ_dolomiti_ski_areas$prices, dnb199.listw, randomisation = TRUE)

## 
##  Moran I test under randomisation
## 
## data:  BZ_dolomiti_ski_areas$prices  
## weights: dnb199.listw    
## 
## Moran I statistic standard deviate = 1.5621, p-value = 0.05913
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.097761383      -0.043478261       0.008175323

Lastly, we can apply the test with different specifications of the spatial weight matrix under permutation bootstrap:

moran.mc(BZ_dolomiti_ski_areas$prices, knn1BZ.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: knn1BZ.listw  
## number of simulations + 1: 1000 
## 
## statistic = -0.024124, observed rank = 590, p-value = 0.41
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, knn2BZ.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: knn2BZ.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.095659, observed rank = 813, p-value = 0.187
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, knn4BZ.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: knn4BZ.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.098204, observed rank = 897, p-value = 0.103
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, knn6BZ.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: knn6BZ.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.12352, observed rank = 969, p-value = 0.031
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, knn8BZ.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: knn8BZ.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.079633, observed rank = 953, p-value = 0.047
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb079.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb079.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.074413, observed rank = 792, p-value = 0.208
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb089.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb089.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.069293, observed rank = 767, p-value = 0.233
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb099.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb099.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.069293, observed rank = 742, p-value = 0.258
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb109.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb109.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.055522, observed rank = 758, p-value = 0.242
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb159.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb159.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.12348, observed rank = 921, p-value = 0.079
## alternative hypothesis: greater

moran.mc(BZ_dolomiti_ski_areas$prices, dnb199.listw, nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  BZ_dolomiti_ski_areas$prices 
## weights: dnb199.listw  
## number of simulations + 1: 1000 
## 
## statistic = 0.097761, observed rank = 921, p-value = 0.079
## alternative hypothesis: greater

The values are only slightly different according to the 3 approaches used (especially between normality and randomization assumptions). When more neighbours are considered (both in the knn and in the critical cut-off) there seems to be a very small positive global spatial autocorrelation. In the most significant case of the knn, with 6 and 8 neighbours and in the case of the highest distances with the critical cut-off, it seems to reach at most a positive correlation of 0.12. This is actually very low and makes us think that there is no effective global correlation.

The Moran’s I test of spatial autocorrelation in OLS residuals

As shown in class, Moran’s I test can also be used as a diagnostic tool to detect the presence of spatial autocorrelation in the residuals of a linear regression model.

Let’s try to apply these to our data considering AREA column:

LinearSolow = lm(prices ~ AREA, BZ_dolomiti_ski_areas)
summary(LinearSolow) 

## 
## Call:
## lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -33.662  -2.951   1.246   7.267   9.072 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 5.281e+01  2.517e+00  20.980 4.89e-16 ***
## AREA        3.257e-06  1.778e-06   1.832   0.0806 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.441 on 22 degrees of freedom
## Multiple R-squared:  0.1323, Adjusted R-squared:  0.0929 
## F-statistic: 3.356 on 1 and 22 DF,  p-value: 0.08055

The results of linear regression are slightly significant, thus we can plot the studentized residuals, which can give a hint about the presence of spatial dependence in the residuals:

studres = rstudent(LinearSolow)
resdistr = quantile(studres) 
colours = grey((length(resdistr):2)/length(resdistr))

plot(area_dolomiti_BZ, col = "#E29B50")
points(BZ_dolomiti_ski_areas, col = colours[findInterval(studres, resdistr, all.inside = T)], cex = 2.5, pch = 18)

We can now apply the test to the studentized residuals of the LinearSolow model, for different specifications of the spatial weights matrix:

lm.morantest(LinearSolow, knn1BZ.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: knn1BZ.listw
## 
## Moran I statistic standard deviate = -0.19491, p-value = 0.5773
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.09610513      -0.04608356       0.06586701

lm.morantest(LinearSolow, knn2BZ.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: knn2BZ.listw
## 
## Moran I statistic standard deviate = 0.20114, p-value = 0.4203
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.01308634      -0.04845623       0.03092287

lm.morantest(LinearSolow, knn4BZ.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: knn4BZ.listw
## 
## Moran I statistic standard deviate = 0.38687, p-value = 0.3494
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##     -0.002202292     -0.049709129      0.015079519

lm.morantest(LinearSolow, knn6BZ.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: knn6BZ.listw
## 
## Moran I statistic standard deviate = 0.95051, p-value = 0.1709
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.040858045     -0.048541067      0.008846062

lm.morantest(LinearSolow, knn8BZ.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: knn8BZ.listw
## 
## Moran I statistic standard deviate = 1.1631, p-value = 0.1224
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      0.040511104     -0.046564518      0.005605025

lm.morantest(LinearSolow, dnb079.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb079.listw
## 
## Moran I statistic standard deviate = -0.073309, p-value = 0.5292
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.05832482      -0.04192036       0.05007348

lm.morantest(LinearSolow, dnb089.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb089.listw
## 
## Moran I statistic standard deviate = -0.068302, p-value = 0.5272
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.05595552      -0.04161646       0.04407365

lm.morantest(LinearSolow, dnb099.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb099.listw
## 
## Moran I statistic standard deviate = -0.068302, p-value = 0.5272
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.05595552      -0.04161646       0.04407365

lm.morantest(LinearSolow, dnb109.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb109.listw
## 
## Moran I statistic standard deviate = -0.082074, p-value = 0.5327
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##      -0.05943414      -0.04337113       0.03830422

lm.morantest(LinearSolow, dnb159.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb159.listw
## 
## Moran I statistic standard deviate = 0.73844, p-value = 0.2301
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.04573131      -0.04434254       0.01487877

lm.morantest(LinearSolow, dnb199.listw, resfun = rstudent)

## 
##  Global Moran I for regression residuals
## 
## data:  
## model: lm(formula = prices ~ AREA, data = BZ_dolomiti_ski_areas)
## weights: dnb199.listw
## 
## Moran I statistic standard deviate = 0.74513, p-value = 0.2281
## alternative hypothesis: greater
## sample estimates:
## Observed Moran I      Expectation         Variance 
##       0.03166855      -0.04734508       0.01124454

Again, as in the previous approaches, we obtain the lowest p-values in the specifications that consider more neighbours. However, the p-values are quite high in this case.

Let’s apply the test under permutation bootstrap as final step:

LinearSolow.lmx = lm(
  prices ~ AREA, 
  data = BZ_dolomiti_ski_areas, 
  x = TRUE
  )

MoraneI.boot = function(var, i, ...) {
  var = var[i]
  lmres = lm(var ~ LinearSolow.lmx$x - 1)
  return(moran(x = residuals(lmres), ...)$I)      
}


boot1 = boot(
  residuals(LinearSolow.lmx),
  statistic = MoraneI.boot,  
  R = 999, 
  sim = "permutation", 
  listw = knn8BZ.listw,  
  n = length(knn8BZ.listw$neighbours), 
  S0 = Szero(knn8BZ.listw)
  )

ti = (boot1$t0 - mean(boot1$t))/sqrt(var(boot1$t))  

plot(boot1)

---

Descriptive spatial statistics for areal data (local spatial spillover)

Moran’s I-statistic is a global measure and therefore does not allow to identify local patterns of spatial autocorrelation. As we saw in class in many circumstances, however, it may also be of interest to assess the presence of local spatial clusters and to check the specific contribution of any particular region to the global pattern of spatial dependence.

We will investigate the local spatial autocorrelation by the means of the Moran scatterplot and the local Moran’s I.

Moran scatterplot

In the Moran scatterplot we can divide the plane into four quadrants that connote the four possible kinds of spatial association between each region and the other neighbouring regions:

• The regions which fall into the High-High and Low-Low quadrants are characterized by positive spatial autocorrelation and are surrounded by regions with similar values;
• High-Low and Low-High identify local patterns of negative spatial autocorrelation since they collect the regions with a high (respectively low) value of x and, in opposition, a low (respectively high) value of wx.

NB: x axis represents the variable of interest and y axis represents the corresponding spatially lagged values (wx).

Let’s plot the Moran scatterplot specifying our variable of interest (prices):

mplot_knn1BZ = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=knn1BZ.listw, 
  main="Moran scatterplot - knn1BZ", 
  return_df=F
  )

mplot_knn2BZ = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=knn2BZ.listw, 
  main="Moran scatterplot - knn2BZ", 
  return_df=F
  )

mplot_knn4BZ = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=knn4BZ.listw, 
  main="Moran scatterplot - knn4BZ", 
  return_df=F
  )

mplot_knn6BZ = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=knn6BZ.listw, 
  main="Moran scatterplot - knn6BZ", 
  return_df=F
  )

mplot_knn8BZ = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=knn8BZ.listw, 
  main="Moran scatterplot - knn8BZ", 
  return_df=F
  )

mplot_dnb079 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb079.listw, 
  main="Moran scatterplot - dnb079", 
  return_df=F
  )

mplot_dnb089 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb089.listw, 
  main="Moran scatterplot - dnb089", 
  return_df=F
  )

mplot_dnb099 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb099.listw, 
  main="Moran scatterplot - dnb099", 
  return_df=F
  )

mplot_dnb109 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb109.listw, 
  main="Moran scatterplot - dnb109", 
  return_df=F
  )

mplot_dnb159 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb159.listw, 
  main="Moran scatterplot - dnb159", 
  return_df=F
  )

mplot_dnb199 = moran.plot(
  BZ_dolomiti_ski_areas$prices, 
  listw=dnb199.listw, 
  main="Moran scatterplot - dnb199", 
  return_df=F
  )

The points (ski areas) that are relatively more influential than others in determining the observed value of global Moran’s I are represented by marked points. These points are identified as having the greatest influence on the regression line based on standard criteria such as Cook’s distance and linkages.

For this reason it may be useful to map regions of significant influence encoded by their quadrant in Moran’s scatterplot. First, we need to identify the influential regions (which was already visible form the plots):

hotspot = as.numeric(row.names(as.data.frame(summary(mplot_dnb159))))

## Potentially influential observations of
##   lm(formula = wx ~ x) :
## 
##    dfb.1_ dfb.x dffit cov.r   cook.d hat    
## 17  0.75  -0.72  0.75  2.70_*  0.29   0.61_*

hotspot = as.numeric(row.names(as.data.frame(summary(mplot_knn6BZ))))

## Potentially influential observations of
##   lm(formula = wx ~ x) :
## 
##    dfb.1_ dfb.x dffit cov.r   cook.d hat    
## 17  0.46  -0.44  0.46  2.76_*  0.11   0.61_*

We can also obtain the spatially lagged values (wx) of the variable of interest for the influential regions:

BZ_dolomiti_ski_areas$wx = lag.listw(knn6BZ.listw, BZ_dolomiti_ski_areas$prices)

Moreover we can assign each influential region to the proper Moran scatterplot quadrant:

BZ_dolomiti_ski_areas$quadrant = rep("None", length(BZ_dolomiti_ski_areas$prices))

for(i in 1:length(hotspot))  {
  if (BZ_dolomiti_ski_areas$prices[hotspot[i]] > mean(BZ_dolomiti_ski_areas$prices) & BZ_dolomiti_ski_areas$wx[hotspot[i]] > mean(BZ_dolomiti_ski_areas$wx)) 
        BZ_dolomiti_ski_areas$quadrant[hotspot[i]] = "HH" 
  if (BZ_dolomiti_ski_areas$prices[hotspot[i]] > mean(BZ_dolomiti_ski_areas$prices) & BZ_dolomiti_ski_areas$wx[hotspot[i]] < mean(BZ_dolomiti_ski_areas$wx)) 
        BZ_dolomiti_ski_areas$quadrant[hotspot[i]] = "HL" 
  if (BZ_dolomiti_ski_areas$prices[hotspot[i]] < mean(BZ_dolomiti_ski_areas$prices) & BZ_dolomiti_ski_areas$wx[hotspot[i]] < mean(BZ_dolomiti_ski_areas$wx)) 
        BZ_dolomiti_ski_areas$quadrant[hotspot[i]] = "LL" 
  if (BZ_dolomiti_ski_areas$prices[hotspot[i]] < mean(BZ_dolomiti_ski_areas$prices) & BZ_dolomiti_ski_areas$wx[hotspot[i]] > mean(BZ_dolomiti_ski_areas$wx)) 
        BZ_dolomiti_ski_areas$quadrant[hotspot[i]] = "LH" 
}

table(BZ_dolomiti_ski_areas$quadrant)

## 
##   LL None 
##    1   23

And we can plot the results, where is our seignificant point on the map:

# which allows us to plot the map of the regions with influence by typing
BZ_dolomiti_ski_areas$colours[BZ_dolomiti_ski_areas$quadrant == "None"] = "#EBEBEB"
BZ_dolomiti_ski_areas$colours[BZ_dolomiti_ski_areas$quadrant == "HH"] = "#745296" 
BZ_dolomiti_ski_areas$colours[BZ_dolomiti_ski_areas$quadrant == "LL"] = "#8B9EB7" # 1 low low 
BZ_dolomiti_ski_areas$colours[BZ_dolomiti_ski_areas$quadrant == "LH"] = "#E85F5C"
BZ_dolomiti_ski_areas$colours[BZ_dolomiti_ski_areas$quadrant == "HL"] = "#FFF07C"

plot(area_dolomiti_BZ, col = "#E29B50")
points(BZ_dolomiti_ski_areas, col = "black", pch=15, cex=2.25)
points(BZ_dolomiti_ski_areas, col = BZ_dolomiti_ski_areas$colours, pch=15, cex=2)

legend(
  x = -10, 
  y = 73, 
  legend = c("None", "Low-Low", "High-Low", "Low-High", "High-High"),
  fill = c("#EBEBEB", "#8B9EB7", "#FFF07C", "#E85F5C", "#745296"), 
  bty = "n", 
  cex = 1
  )

title(main = "Points with influence")

The Local Moran’s I

The Moran scatterplot represents a useful and intuitive visual representation of local patterns of spatial association, but it cannot provide statistical significance of the results. For this reason, to assess the significance of the revealed pattern, we can rely on the so-called local Moran I index:

lmI = localmoran(BZ_dolomiti_ski_areas$prices, dnb159.listw)
head(lmI)

##            Ii         E.Ii     Var.Ii        Z.Ii Pr(z != E(Ii))
## 1 -0.03776105 -0.017915827 0.06909978 -0.07549491      0.9398209
## 2 -0.27668756 -0.015376377 0.03775155 -1.34490190      0.1786569
## 3 -0.03586392 -0.012625307 0.06459607 -0.09143388      0.9271478
## 4 -0.12610217 -0.012625307 0.09066115 -0.37687445      0.7062669
## 5  0.29516084 -0.006565256 0.02561427  1.88526256      0.0593944
## 6 -0.04468586 -0.015376377 0.07845244 -0.10464174      0.9166601

lmI = localmoran(BZ_dolomiti_ski_areas$prices, knn6BZ.listw)
head(lmI)

##            Ii         E.Ii     Var.Ii       Z.Ii Pr(z != E(Ii))
## 1 -0.03996608 -0.017915827 0.05438408 -0.0945535     0.92466949
## 2 -0.25917800 -0.015376377 0.04679619 -1.1270193     0.25973434
## 3  0.09139515 -0.012625307 0.03853099  0.5299241     0.59616456
## 4  0.10065035 -0.012625307 0.03853099  0.5770740     0.56388946
## 5  0.28448232 -0.006565256 0.02015938  2.0498656     0.04037755
## 6 -0.08554151 -0.015376377 0.04679619 -0.3243516     0.74567184

We can also plot the distribution of local Moran’s I-index values:

brks = sort(as.numeric(lmI[,1]))
colours = grey((0:length(lmI[,1]))/length(lmI[,1]))

plot(area_dolomiti_BZ, col = "#E29B50")
points(BZ_dolomiti_ski_areas, col = "black", cex = 2.75, pch = 17)
points(BZ_dolomiti_ski_areas, col = colours[findInterval(lmI[,1], brks, all.inside = TRUE)], cex = 2.5, pch = 17)
title(main = "Local Moran's I values")

With this representation we can identify which ski areas have a positive local spatial correlation, in white (or light grey), and which have a negative one, in dark. However, it should be borne in mind that this index does not provide us with information on the significance of results.

As previously done with the global index, Moran’s local I-statistic can also be tested for deviations using the null hypothesis of no local spatial autocorrelation and thus can provide the statistical significance of the local spatial patterns detected by Moran’s scatterplot. Thus, we can map the corresponding p-values:

pval = as.numeric(lmI[,5])
BZ_dolomiti_ski_areas$colpval[pval>0.05] = "white" 
BZ_dolomiti_ski_areas$colpval[pval<=0.05 & pval>0.01] = "#A5D5F3"
BZ_dolomiti_ski_areas$colpval[pval<=0.01 & pval>0.001] = "#6FBCEB"
BZ_dolomiti_ski_areas$colpval[pval<=0.001 & pval>0.0001] = "#3AA3E4"
BZ_dolomiti_ski_areas$colpval[pval<=0.0001] = "#1878B4"

plot(area_dolomiti_BZ, col = "#E29B50")
points(BZ_dolomiti_ski_areas, col="black", pch=18, cex=3.5)
points(BZ_dolomiti_ski_areas, col=BZ_dolomiti_ski_areas$colpval, pch=18, cex=3)
legend(x=-10, y=73, legend=c("Not significant", 
       "p-value = 0.05", "p-value = 0.01", "p-value = 0.001", 
       "p-value = 0.0001"), fill=c("white", gray(0.9), gray(0.7),    
       gray(0.4), "black"), bty="n", cex=0.8)
title(main="Local Moran's I significance map")

Only a few of the ski areas have a significant spatial correlation. These are at the bottom left of the map and are: Carezza, Nova Ponente, Panorama, Pietralba e Passo Olcini (which also correspond to point 17, the one with the most influence)

Conclusion

We did not obtain significant results with most of the chosen neighbouroods, if we consider the global spatial autocorrelation. The p-value was significant only when we considered many neighbours for each ski area; however, the Moran I statistic value was just above 0, indicating an extremely slight positive correlation.

On the other hand, considering the local index we find (again, using quite high number of neighbours or a quite high distance) that there are some ski areas that have significant positive or negative spatial autocorrelation.