-
Notifications
You must be signed in to change notification settings - Fork 0
/
SPA.test().R
204 lines (182 loc) · 9.74 KB
/
SPA.test().R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
#This is an example of how the function works.
#Within the example, I use Census data for the population of Camden, London.
#The area of Camden is formatted in the British National grid system.
#Load Packages
setwd('mdecara_pols0010')
library(tmap)
library(leaflet)
library(rgdal)
library(rgeos)
library(sp)
library(RColorBrewer)
library(spdep)
#Load data
Census.Data <-read.csv("worksheet_data/camden/practical_data.csv")
Output.Areas <- readOGR("worksheet_data/camden", "Camden_oa11")
#Merge data
data <- merge(Output.Areas, Census.Data, by.x="OA11CD", by.y="OA")
SPA.test <- function(data, variable, variable_name, neighbour, d1, d2, k, style, method){
#Step 1 - Plot basic template of area and chosen variable to give an idea
#to the user of map and variable presence.
print(tm_shape(data) +
tm_fill(variable_name, palette = "Reds",
style = "quantile", title = "% of chosen variable") +
tm_borders(alpha=.4))
#Step 2 - Function identifies neighbours by chosen measure and adds chosen
#weight style to then plot and summarize the results.
#In this case: ‘queens’ measure.
if(neighbour == 'queens'){
#create neighbour ‘nb’ type object
N <- poly2nb(data, queen = TRUE)
#summary of object
summary.nb(N)
#plot of neighbour object
plot(data, border = 'lightgrey')
plot(N, coordinates(data), add=TRUE, col='blue')
title(main='Neighbour Map by Queens Case')
#create list of neighbour object with chosen weight. This case: ‘W’
listw <- nb2listw(N, style = style)
#summary of list
summary.listw(listw)
} else if(neighbour == 'rooks'){
N <- poly2nb(data, queen = FALSE)
#Plot and Print Neighbour Outcome
summary.nb(N)
plot(data, border = 'lightgrey')
plot(N, coordinates(data), add=TRUE, col='red')
title(main='Neighbour Map by Rooks Case')
listw <- nb2listw(N, style = style)
summary.listw(listw)
} else if(neighbour == 'distance'){
N <- dnearneigh(coordinates(data), d1, d2)
#Plot and Print Neighbour Outcome
summary.nb(N)
plot(data, border = 'lightgrey')
plot(N, coordinates(data), add=TRUE, col='green')
title(main='Neighbour Map by Distance')
listw <- nb2listw(N, style = style)
summary.listw(listw)
} else if(neighbour == 'knn'){
N <- knearneigh(coordinates(data), k = k)
#Plot and Print Neighbour Outcome
plot(data, border = 'lightgrey')
plot(knn2nb(N), coordinates(data), add=TRUE, col='orange')
title(main='Neighbour Map by K Nearest Neighbours')
N <- knn2nb(N)
summary.nb(N)
listw <- nb2listw(N, style = style)
summary.listw(listw)
}
#Step 3 - Function then runs the chosen spatial autocorrelation test,
#prints global and local results and plots certain graphs.
#In this case: everything under ‘Morans I’ method.
if(method == 'Morans I'){
test_scores <- moran.test(variable, listw)
print(test_scores)
if(test_scores$estimate['Moran I statistic'] > 0){
print(paste("It can be said there is POSITIVE spatial autocorrelation with a Global Moran I statistic of: ",test_scores$estimate['Moran I statistic'],sep=""))
} else if(test_scores$estimate['Moran I statistic'] == 0) {
print(paste("It can be said there is NO spatial autocorrelation with a Global Moran I statistic of: ",test_scores$estimate['Moran I statistic'],sep=""))
} else if(test_scores$estimate['Moran I statistic'] < 0) {
print(paste("It can be said there is NEGATIVE spatial autocorrelation with a Global Moran I statistic of: ",test_scores$estimate['Moran I statistic'],sep=""))
}
if(test_scores$statistic['Moran I statistic standard deviate'] >= 1.96){
print(paste('It can be said the data is significantly CLUSTERED with a z-score of:', test_scores$statistic['Moran I statistic standard deviate'], sep=''))
} else if(test_scores$statistic['Moran I statistic standard deviate'] <= -1.96){
print(paste('It can be said the data is significantly DISPERSED with a z-score of:', test_scores$statistic['Moran I statistic standard deviate'], sep=''))
} else if(test_scores$statistic['Moran I statistic standard deviate'] < 1.96 & test_scores$statistic['Moran I statistic standard deviate'] > -1.96){
print(paste('It can be said the data is spatially RANDOM with a z-score of:', test_scores$statistic['Moran I statistic standard deviate'], sep=''))
}
if(test_scores$p.value <= 0.10){
print(paste("It can also be said the model IS statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
} else {
print(paste("It can also be said the model IS NOT statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
}
#creates moran scatter plot
moran <- moran.plot(variable, listw = nb2listw(N, style = style), asp = T)
#creates local moran results plot
local <- localmoran(x = variable, listw = nb2listw(N, style = style))
moran.map <- data
moran.map@data <- cbind(data@data, local)
print(tm_shape(moran.map) + tm_fill(col = "Ii", style = "quantile", title = "Local Moran Statistic", midpoint = NA))
#Step 4 - Creates a LISA cluster map
quadrant <- vector(mode="numeric",length=nrow(local))
#centers the variable of interest around its mean
m.variable <- variable - mean(variable)
#centers the local Moran's around the mean
m.local <- local[,1] - mean(local[,1])
#significance threshold
signif <- 0.1
#builds a data quadrant
quadrant[m.variable < 0 & m.local< 0] <- 1
quadrant[m.variable < 0 & m.local> 0] <- 2
quadrant[m.variable > 0 & m.local< 0] <- 3
quadrant[m.variable > 0 & m.local> 0] <- 4
quadrant[local[,5] > signif] <- 0
#plot
breaks <- c(0,1,2,3,4)
colors <- c("white", "blue", rgb(0,0,1, alpha=0.4), rgb(1,0,0, alpha=0.4), "red")
plot(data, border = "lightgray", col = colors[findInterval(quadrant, breaks, all.inside = FALSE)])
box()
legend("bottomleft", legend = c("insignificant", "low-low", "low-high", "high-low", "high-high"),
fill = colors, bty = "n")
legend('topright', legend = 'LISA cluster map', bty = 'n')
}
else if (method == 'Gearys C'){
test_scores <- geary.test(variable, listw)
print(test_scores)
if(test_scores$estimate['Geary C statistic'] == 1){
print(paste("It can be said there is NO spatial autocorrelation with a Global Geary C statistic of: ",test_scores$estimate['Geary C statistic'],sep=""))
} else if(test_scores$estimate['Geary C statistic'] > 1) {
print(paste("It can be said the variable is NEGATIVELY spatially autocorrelated with a Global Geary C statistic of: ",test_scores$estimate['Geary C statistic'],sep=""))
} else if(test_scores$estimate['Geary C statistic'] < 1) {
print(paste("It can be said the variable is POSITIVELY spatially autocorrelated with a Global Geary C statistic of: ",test_scores$estimate['Geary C statistic'],sep=""))
}
if(test_scores$p.value <= 0.05){
print(paste("It can also be said the model IS statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
} else {
print(paste("It can also be said the model IS NOT statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
}
}
else if(method == 'Getis Ord'){
test_scores <- globalG.test(variable, listw)
print(test_scores)
if(test_scores$estimate['Global G statistic'] == 0){
print(paste("It can be said the variable has NO spatial autocorrelation with a Global Getis Ord G statistic of: ",test_scores$estimate['Global G statistic'],sep=""))
} else if(test_scores$estimate['Global G statistic'] > 0) {
print(paste("It can be said the variable is POSITIVELY spatially autocorrelated with a Global Getis Ord G statistic of: ",test_scores$estimate['Global G statistic'],sep=""))
} else if(test_scores$estimate['Global G statistic'] < 0) {
print(paste("It can be said the variable is NEGATIVELY spatially autocorrelated with a Global Getis Ord G statistic of: ",test_scores$estimate['Global G statistic'],sep=""))
}
if(test_scores$statistic['standard deviate'] > 0){
print(paste('It can be said the HIGH VALUES are clustered within the data with a z-score of:', test_scores$statistic['standard deviate'], sep=''))
} else if(test_scores$statistic['standard deviate'] < 0){
print(paste('It can be said the LOW VALUES are clustered within the data with a z-score of:', test_scores$statistic['standard deviate'], sep=''))
}
if(test_scores$p.value <= 0.05){
print(paste("It can also be said the model IS statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
} else {
print(paste("It can also be said the model IS NOT statistically significant with a Global P-value of: ",test_scores$p.value, sep=""))
}
local_g <- localG(variable, listw)
print('This is the summary of the local G statistic:')
print(summary(local_g))
local_g_sp <- data
local_g_sp@data <- cbind(data@data, as.matrix(local_g))
names(local_g_sp)[6] <- "gstat"
#plot
print(tm_shape(local_g_sp) + tm_fill("gstat", palette = "RdBu", style = "pretty", title = 'gstat') + tm_borders(alpha=.4))
}
}
#This will be parameters imputed into the function for demonstration purposes.
#Therefore, every result displayed by the function in the code below will
#be a result of only these parameters.
SPA.test(data = data,
variable = data$Qualification,
variable_name = 'Qualification',
neighbour = 'queens',
d1 = 0, d2 = 800, k = 4,
style = 'W',
method = 'Morans I')
#Note: parameters ‘d1,d2 & k’ will simply not be run as they are dependent on
#other inputs such as distance or knn measure of neighbours.