Standard approaches
In this section, we focus on classical approaches that consider either the multivariate or the spatial aspect of the data.
Multivariate analysis
Here we consider \(p=6\) variables measured for \(n=85\) individuals (départements of France). As only quantitative variables have been recorded, principal component analysis (PCA, Hotelling 1933) is well adapted. PCA summarizes the data by maximizing simultaneously the variance of the projection of the individuals onto the principal axes and the sum of the squared correlations between the principal component and the variables.
The biplot is simply obtained by
The first two PCA dimensions account for 35.7% and 20% ,respectively, of the total variance.
Correlations between variables and principal components can be represented on a correlation circle. The first axis is negatively correlated to literacy and positively correlated to property crime, suicides and illegitimate births. The second axis is aligned mainly with personal crime and donations to the poor.
Spatial information can be added on the factorial map representing the projections of départements on principal axes by coloring according to colors representing the different regions of France.
s.Spatial(france.map, col = col.region[region.names], plabel.cex = 0)
s.class(xy, region.names, col = col.region, add = TRUE, ellipseSize = 0, starSize = 0)
For the first axis, the North and East are characterized by negative scores, corresponding to high levels of literacy and high rates of suicides, crimes against property and illegitimate births. The second axis mainly contrasts the West (high donations to the the poor and low levels of crime against persons) to the South.
Spatial autocorrelation
Spatial autocorrelation statistics, such as Moran (1948) Coefficient (MC) and Geary (1954) Ratio, aim to measure and analyze the degree of dependency among observations in a geographical context (Cliff and Ord 1973).
The spatial weighting matrix
The first step of spatial autocorrelation analysis is to define a
spatial weighting matrix \(\mathbf{W}=[w_{ij}]\) . In the case of
Guerry’s data, we simply defined a binary neighborhood where two
départements are considered as neighbors if they share a common border.
The spatial weighting matrix is then obtained after row-standardization
(style = "W"):
We can represent this neighborhood on the geographical map:
Moran’s Coefficient
Once the spatial weights have been defined, the spatial autocorrelation statistics can then be computed. Let us consider the \(n\)-by-1 vector \(\mathbf{x}=\left[ {x_1 \cdots x_n } \right]^\textrm{T}\) containing measurements of a quantitative variable for \(n\) spatial units. The usual formulation for Moran’s coefficient of spatial autocorrelation (Cliff and Ord 1973; Moran 1948) is \[\begin{equation} \label{eq1} MC(\mathbf{x})=\frac{n\sum\nolimits_{\left( 2 \right)} {w_{ij} (x_i -\bar {x})(x_j -\bar {x})} }{\sum\nolimits_{\left( 2 \right)} {w_{ij} } \sum\nolimits_{i=1}^n {(x_i -\bar {x})^2} }\mbox{ where }\sum\nolimits_{\left( 2 \right)} =\sum\limits_{i=1}^n {\sum\limits_{j=1}^n } \mbox{ with }i\ne j. \end{equation}\]
MC can be rewritten using matrix notation: \[\begin{equation} \label{eq2} MC(\mathbf{x})=\frac{n}{\mathbf{1}^\textrm{T}\mathbf{W1}}\frac{\mathbf{z}^\textrm{T}{\mathbf{Wz}}}{\mathbf{z}^\textrm{T}\mathbf{z}}, \end{equation}\] where \(\mathbf{z}=\left ( \mathbf{I}_n-\mathbf{1}_n\mathbf{1}_n^\textrm{T} /n \right )\mathbf{x}\) is the vector of centered values (\(z_i=x_i-\bar{x}\)) and \(\mathbf{1}_n\) is a vector of ones (of length \(n\)).
The significance of the observed value of MC can be tested by a Monte-Carlo procedure, in which locations are permuted to obtain a distribution of MC under the null hypothesis of random distribution. An observed value of MC that is greater than that expected at random indicates the clustering of similar values across space (positive spatial autocorrelation), while a significant negative value of MC indicates that neighboring values are more dissimilar than expected by chance (negative spatial autocorrelation).
We computed MC for the Guerry’s dataset. A positive and significant autocorrelation is identified for each of the six variables. Thus, the values of literacy are the most covariant in adjacent departments, while illegitimate births (Infants) covary least.
moran.randtest(df, lw)
#> class: krandtest lightkrandtest
#> Monte-Carlo tests
#> Call: moran.randtest(x = df, listw = lw)
#>
#> Number of tests: 6
#>
#> Adjustment method for multiple comparisons: none
#> Permutation number: 999
#> Test Obs Std.Obs Alter Pvalue
#> 1 Crime_pers 0.4114597 5.797137 greater 0.001
#> 2 Crime_prop 0.2635533 4.160819 greater 0.002
#> 3 Literacy 0.7176053 10.357515 greater 0.001
#> 4 Donations 0.3533613 5.204476 greater 0.001
#> 5 Infants 0.2287241 3.628118 greater 0.002
#> 6 Suicides 0.4016812 6.081180 greater 0.001Moran scatterplot
If the spatial weighting matrix is row-standardized, we can define the lag vector \(\mathbf{\tilde{z}} = \mathbf{Wz}\) (i.e., \(\tilde{z}_i = \sum\limits_{j=1}^n{w_{ij}x_j}\)) composed of the weighted (by the spatial weighting matrix) averages of the neighboring values. Thus, we have: \[\begin{equation} \label{eq3} MC(\mathbf{x})=\frac{\mathbf{z}^\textrm{T}{\mathbf{\tilde{z}}}}{\mathbf{z}^\textrm{T}\mathbf{z}}, \end{equation}\] since in this case \(\mathbf{1}^\textrm{T}\mathbf{W1}=n\). This shows clearly that MC measures the autocorrelation by giving an indication of the intensity of the linear association between the vector of observed values \(\mathbf{z}\) and the vector of weighted averages of neighboring values \(\mathbf{\tilde{z}}\). Anselin (1996) proposed to visualize MC in the form of a bivariate scatterplot of \(\mathbf{\tilde{z}}\) against \(\mathbf{z}\). A linear regression can be added to this Moran scatterplot, with slope equal to MC.
Considering the Literacy variable of Guerry’s data, the Moran scatterplot clearly shows strong autocorrelation. It also shows that the Hautes-Alpes département has a slightly outlying position characterized by a high value of Literacy compared to its neighbors.
x <- df[, 3]
x.lag <- lag.listw(lw, df[, 3])
moran.plot(x, lw)
text(x[5], x.lag[5], dep.names[5], pos = 1, cex = 0.8)
Indirect integration of multivariate and geographical aspects
The simplest approach considered a two-step procedure where the data are first summarized with multivariate analysis such as PCA. In a second step, univariate spatial statistics or mapping techniques are applied to PCA scores for each axis separately. One can also test for the presence of spatial autocorrelation for the first few scores of the analysis, with univariate autocorrelation statistics such as MC. We mapped scores of the départements for the first two axes of the PCA of Guerry’s data. Even if PCA maximizes only the variance of these scores, there is also a clear spatial structure, as the scores are highly autocorrelated. The map for the first axis corresponds closely to the split between la France éclairée (North-East characterized by an higher level of Literacy) and la France obscure.
moran.randtest(pca$li, lw)
#> class: krandtest lightkrandtest
#> Monte-Carlo tests
#> Call: moran.randtest(x = pca$li, listw = lw)
#>
#> Number of tests: 3
#>
#> Adjustment method for multiple comparisons: none
#> Permutation number: 999
#> Test Obs Std.Obs Alter Pvalue
#> 1 Axis1 0.5506343 7.704990 greater 0.001
#> 2 Axis2 0.5614030 8.237265 greater 0.001
#> 3 Axis3 0.1805569 2.701054 greater 0.006
s.value(xy, pca$li[, 1:2], Sp = france.map, pSp.border = "white", symbol = "circle", pgrid.draw = FALSE)
