The package `adespatial`

contains functions for the multiscale analysis of spatial multivariate data. It implements some new functions and reimplements existing functions that were available in packages of the sedaR project hosted on R-Forge (`spacemakeR`

, `packfor`

, `AEM`

, etc.). It can be seen as a bridge between packages dealing with mutltivariate data (e.g., `ade4`

, Dray and Dufour (2007)) and packages that deals with spatial data (`spdep`

). In `adespatial`

, the spatial information is considered as a spatial weighting matrix, object of class `listw`

provided by the `spdep`

package (Figure 1). It allows to build Moran’s Eigenvector Maps (MEM, Dray, Legendre, and Peres-Neto (2006)) that are orthogonal vectors maximizing the spatial autocorrelation (measured by Moran’s index of autocorrelation). These spatial predictors can be used in multivariate statistical methods to provide spatially-explicit multiscale tools (Dray et al. 2012). This document provides a description of the main functionalities of the package.

Figure 1: Schematic representation of the functioning of the `adespatial`

package. Classes are represented in pink frames and functions in blue frames. Classes and functions provided by `adespatial`

are in bold.

To run the different analysis described, several packages are required and are loaded:

```
library(adespatial)
library(ade4)
```

```
##
## Attaching package: 'ade4'
```

```
## The following object is masked from 'package:adespatial':
##
## multispati
```

`library(adegraphics)`

```
##
## Attaching package: 'adegraphics'
```

```
## The following objects are masked from 'package:ade4':
##
## kplotsepan.coa, s.arrow, s.class, s.corcircle, s.distri,
## s.image, s.label, s.logo, s.match, s.traject, s.value,
## table.value, triangle.class
```

`library(spdep)`

`## Loading required package: sp`

`## Loading required package: Matrix`

```
##
## Attaching package: 'spdep'
```

```
## The following object is masked from 'package:ade4':
##
## mstree
```

`library(maptools)`

`## Checking rgeos availability: TRUE`

Spatial neighborhoods are managed in `spdep`

as objects of class `nb`

. It corresponds to the notion of connectivity matrices discussed in Dray, Legendre, and Peres-Neto (2006) and can be represented by an unweighted graph. Various functions are devoted to create `nb`

objects from geographic coordinates of sites. We present different alternatives according to the design of the sampling scheme.

The function `poly2nb`

allows to define neighborhood when the sampling sites are polygons and not points (two regions are neighbors if they share a common boundary).

`columbus <- readShapePoly(system.file("etc/shapes/columbus.shp", package="spdep")[1])`

`## Warning: use rgdal::readOGR or sf::st_read`

```
xx <- poly2nb(columbus)
plot(columbus, border = "grey")
plot(xx, coordinates(columbus), add = TRUE, pch = 20, col = "red")
title(main="Neighborhood for polygons")
```

If the sampling scheme is based on grid of 10 rows and 8 columns, spatial coordinates can be easily generated:

```
xygrid <- expand.grid(x = 1:10, y = 1:8)
plot(xygrid, pch = 20, asp = 1)
```

For a regular grid, spatial neighborhood can be created with the function `cell2nb`

. Two types of neighborhood can be defined. The `queen`

specification considered horizontal, vertical and diagonal edges:

```
nb1 <- cell2nb(10, 8, type = "queen")
plot(nb1, xygrid, col = "red", pch = 20)
title(main = "Queen neighborhood")
```

`nb1`

```
## Neighbour list object:
## Number of regions: 80
## Number of nonzero links: 536
## Percentage nonzero weights: 8.375
## Average number of links: 6.7
```

The `rook`

specification considered only horizontal and vertical edges:

```
nb2 <- cell2nb(10, 8, type = "rook")
plot(nb2, xygrid, col = "red", pch = 20)
title(main = "Rook neighborhood")
```

`nb2`

```
## Neighbour list object:
## Number of regions: 80
## Number of nonzero links: 284
## Percentage nonzero weights: 4.4375
## Average number of links: 3.55
```

The easiest way to deal with transects is to consider them as grids with only one row:

```
xytransect <- expand.grid(1:20, 1)
nb3 <- cell2nb(20, 1)
plot(nb3, xytransect, col = "red", pch = 20)
title(main = "Transect of 20 sites")
```

`summary(nb3)`

```
## Neighbour list object:
## Number of regions: 20
## Number of nonzero links: 38
## Percentage nonzero weights: 9.5
## Average number of links: 1.9
## Link number distribution:
##
## 1 2
## 2 18
## 2 least connected regions:
## 1:1 20:1 with 1 link
## 18 most connected regions:
## 2:1 3:1 4:1 5:1 6:1 7:1 8:1 9:1 10:1 11:1 12:1 13:1 14:1 15:1 16:1 17:1 18:1 19:1 with 2 links
```

All sites have two neighbors except the first and the last one.

There are many ways to define neighborhood in the case of irregular samplings. We consider a random sampling with 10 sites:

```
set.seed(3)
xyir <- matrix(runif(20), 10, 2)
plot(xyir, pch = 20, main = "Irregular sampling with 10 sites")
```

The most intuitive way is to consider that sites are neighbors (or not) according to the distances between them. This definition is provided by the `dnearneigh`

function:

```
nbnear1 <- dnearneigh(xyir, 0, 0.2)
nbnear2 <- dnearneigh(xyir, 0, 0.3)
nbnear3 <- dnearneigh(xyir, 0, 0.5)
nbnear4 <- dnearneigh(xyir, 0, 1.5)
plot(nbnear1, xyir, col = "red", pch = 20)
title(main = "neighbors if 0<d<0.2")
plot(nbnear2, xyir, col = "red", pch = 20)
title(main = "neighbors if 0<d<0.3")
plot(nbnear3, xyir, col = "red", pch = 20)
title(main = "neighbors if 0<d<0.5")
plot(nbnear4, xyir, col = "red", pch = 20)
title(main = "neighbors if 0<d<1.5")
```

Using a distance-based criteria could lead to unbalanced graphs. For instance, if the maximum distance is too low, some points have no neighbors:

`nbnear1`

```
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 14
## Percentage nonzero weights: 14
## Average number of links: 1.4
## 3 regions with no links:
## 2 7 10
```

On the other hand, if the maximum distance is to high, all sites could connected to the 9 others:

`nbnear4`

```
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 90
## Percentage nonzero weights: 90
## Average number of links: 9
```

It is also possible to possible to define neighborhood by a criteria based on nearest neighbors. However, this option can lead to non-symmetric neighborhood: if site A is the nearest neighbor of site B, it does not mean that site B is the nearest neighbor of site A.

The function `knearneigh`

creates an object of class `knn`

. It can be transformed into a `nb`

object with the function `knn2nb`

. This function has an argument `sym`

which can be set to `TRUE`

to force the output neighborhood to symmetry.

```
knn1 <- knearneigh(xyir, k = 1)
nbknn1 <- knn2nb(knn1, sym = TRUE)
knn2 <- knearneigh(xyir, k = 2)
nbknn2 <- knn2nb(knn2, sym = TRUE)
plot(nbknn1, xyir, col = "red", pch = 20)
title(main = "Nearest neighbors (k=1)")
plot(nbknn2, xyir, col = "red", pch = 20)
title(main="Nearest neighbors (k=2)")
```

This definition of neighborhood can lead to unconnected subgraphs. The function `n.comp.nb`

finds the number of disjoint connected subgraphs:

`n.comp.nb(nbknn1)`

```
## $nc
## [1] 3
##
## $comp.id
## [1] 1 2 1 1 3 3 1 1 3 2
```

`n.comp.nb(nbknn2)`

```
## $nc
## [1] 1
##
## $comp.id
## [1] 1 1 1 1 1 1 1 1 1 1
```

More elaborate procedures are available to define neighborhood. For instance, Delaunay triangulation is obtained with the function `tri2nb`

. It requires the package `deldir`

. Other graph-based procedures are also available:

`nbtri <- tri2nb(xyir)`

```
##
## PLEASE NOTE: The components "delsgs" and "summary" of the
## object returned by deldir() are now DATA FRAMES rather than
## matrices (as they were prior to release 0.0-18).
## See help("deldir").
##
## PLEASE NOTE: The process that deldir() uses for determining
## duplicated points has changed from that used in version
## 0.0-9 of this package (and previously). See help("deldir").
```

```
nbgab <- graph2nb(gabrielneigh(xyir), sym = TRUE)
nbrel <- graph2nb(relativeneigh(xyir), sym = TRUE)
nbsoi <- graph2nb(soi.graph(nbtri, xyir), sym = TRUE)
plot(nbtri, xyir, col = "red", pch = 20)
title(main="Delaunay triangulation")
plot(nbgab, xyir, col = "red", pch = 20)
title(main = "Gabriel Graph")
plot(nbrel, xyir, col = "red", pch = 20)
title(main = "Relative Neighbor Graph")
plot(nbsoi, xyir, col = "red", pch = 20)
title(main = "Sphere of Influence Graph")
```

The function `chooseCN`

provides a simple way to build spatial neighborhoods. It is a wrap up to many of the `spdep`

functions presented above. The function `createlistw`

discussed in section XX is an interactive graphical interface that allows to generate R code to build neighborhood objects.

`nb`

objectsA `nb`

object is a list of neighbors. The neighbors of the first site are in the first element of the list:

`nbgab[[1]]`

`## [1] 4 7`

Various tools are provided by `spdep`

to deal with these objects. For instance, it is possible to identify differences between two neighborhoods:

`diffnb(nbsoi,nbrel)`

```
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 16
## Percentage nonzero weights: 16
## Average number of links: 1.6
## 2 regions with no links:
## 4 5
```

Usually, it can be useful to remove some connections due to edge effects. In this case, the function `edit.nb`

provides an interactive tool to add or delete connections.

The function `include.self`

allows to include a site itself in its own list of neighbors:

`str(nbsoi)`

```
## List of 10
## $ : int [1:4] 3 4 7 8
## $ : int 10
## $ : int [1:3] 1 4 8
## $ : int [1:3] 1 3 8
## $ : int [1:2] 6 9
## $ : int [1:2] 5 9
## $ : int [1:2] 1 10
## $ : int [1:3] 1 3 4
## $ : int [1:2] 5 6
## $ : int [1:2] 2 7
## - attr(*, "region.id")= chr [1:10] "1" "2" "3" "4" ...
## - attr(*, "call")= language soi.graph(tri.nb = nbtri, coords = xyir)
## - attr(*, "class")= chr "nb"
## - attr(*, "sym")= logi TRUE
```

`str(include.self(nbsoi))`

```
## List of 10
## $ : int [1:5] 1 3 4 7 8
## $ : int [1:2] 2 10
## $ : int [1:4] 1 3 4 8
## $ : int [1:4] 1 3 4 8
## $ : int [1:3] 5 6 9
## $ : int [1:3] 5 6 9
## $ : int [1:3] 1 7 10
## $ : int [1:4] 1 3 4 8
## $ : int [1:3] 5 6 9
## $ : int [1:3] 2 7 10
## - attr(*, "region.id")= chr [1:10] "1" "2" "3" "4" ...
## - attr(*, "call")= language soi.graph(tri.nb = nbtri, coords = xyir)
## - attr(*, "class")= chr "nb"
## - attr(*, "sym")= logi TRUE
## - attr(*, "self.included")= logi TRUE
```

The `spdep`

package provides many other tools to manipulate `nb`

objects:

```
intersect.nb(nb.obj1, nb.obj2)
union.nb(nb.obj1, nb.obj2)
setdiff.nb(nb.obj1, nb.obj2)
complement.nb(nb.obj)
droplinks(nb, drop, sym = TRUE)
nblag(neighbours, maxlag)
```

Spatial weighting matrices are computed by a transformation of the spatial neighborhood objects. In R, they are not stored as matrices but as objects of the class `listw`

. This format is more efficient than a matrix representation to manage large data sets. An object of class `listw`

can be easily created from an object of class `nb`

with the function `nb2listw`

.

Different objects `listw`

can be obtained from a `nb`

object. The argument `style`

allows to define a transformation of the matrix such as standardization by row sum, by total sum or binary coding, etc. General spatial weights can be introduced by the argument `glist`

. This allows to introduce, for instance, a weighting relative to the distances between the points. For this task, the function `nbdists`

is very useful as it computes Euclidean distance between neighbor sites defined by an `nb`

object.

To obtain a simple row-standardization, the function is simply called by:

`nb2listw(nbgab)`

```
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 26
## Percentage nonzero weights: 26
## Average number of links: 2.6
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 10 100 10 8.513889 41.04167
```

More sophisticated forms of spatial weighting matrices can be defined. For instance, it is possible to weight edges between neighbors as functions of geographic distances. In a fist step, distances between neighbors are obtained by the function :

```
distgab <- nbdists(nbgab, xyir)
str(distgab)
```

```
## List of 10
## $ : num [1:2] 0.166 0.403
## $ : num [1:3] 0.424 0.383 0.286
## $ : num [1:4] 0.4236 0.0617 0.3682 0.3538
## $ : num [1:3] 0.166 0.0617 0.1501
## $ : num [1:2] 0.0383 0.0384
## $ : num [1:4] 0.383 0.3682 0.0383 0.3344
## $ : num [1:2] 0.403 0.534
## $ : num [1:2] 0.15 0.334
## $ : num 0.0384
## $ : num [1:3] 0.286 0.354 0.534
## - attr(*, "class")= chr "nbdist"
## - attr(*, "call")= language nbdists(nb = nbgab, coords = xyir)
```

Then, spatial weights are defined as a function of distance (e.g. \(1-d_{ij}/max(d_{ij})\)):

`fdist <- lapply(distgab, function(x) 1-x/max(dist(xyir)))`

And the spatial weighting matrix is then created:

```
listwgab <- nb2listw(nbgab, glist = fdist, style = "B")
listwgab
```

```
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 10
## Number of nonzero links: 26
## Percentage nonzero weights: 26
## Average number of links: 2.6
##
## Weights style: B
## Weights constants summary:
## n nn S0 S1 S2
## B 10 100 18.19528 27.02395 148.4577
```

`names(listwgab)`

`## [1] "style" "neighbours" "weights"`

`listwgab$neighbours[[1]]`

`## [1] 4 7`

`listwgab$weights[[1]]`

`## [1] 0.8170501 0.5558821`

The matrix representation of a `listw`

object can also be obtained:

`print(listw2mat(listwgab),digits=3)`

```
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## 1 0.000 0.000 0.000 0.817 0.000 0.000 0.556 0.000 0.000 0.000
## 2 0.000 0.000 0.533 0.000 0.000 0.578 0.000 0.000 0.000 0.685
## 3 0.000 0.533 0.000 0.932 0.000 0.594 0.000 0.000 0.000 0.610
## 4 0.817 0.000 0.932 0.000 0.000 0.000 0.000 0.835 0.000 0.000
## 5 0.000 0.000 0.000 0.000 0.000 0.958 0.000 0.000 0.958 0.000
## 6 0.000 0.578 0.594 0.000 0.958 0.000 0.000 0.631 0.000 0.000
## 7 0.556 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.412
## 8 0.000 0.000 0.000 0.835 0.000 0.631 0.000 0.000 0.000 0.000
## 9 0.000 0.000 0.000 0.000 0.958 0.000 0.000 0.000 0.000 0.000
## 10 0.000 0.685 0.610 0.000 0.000 0.000 0.412 0.000 0.000 0.000
```

To facilitate the building of spatial neighborhoods (`nb`

object) and associated spatial weighting matrices (`listw`

object), the package `adespatial`

provides an interactive graphical interface. The interface is launched by the call `createlistw()`

assuming that spatial coordinates are still stored in an object of the R session (Figure 2).

The package `adespatial`

provide different tools to build spatial predictors that can be incorporated in multivariate analysis. They are orthogonal vectors stored in a object of class `orthobasisSp`

. Orthogonal polynomials of geographic coordinates can be computed by the function `orthobasis.poly`

whereas traditional principal coordinates of neighbour matrices (PCNM, Borcard and Legendre (2002)) are obtained by the function `dbmem`

. The more flexible Moran’s eigenvectors maps (MEMs) of a spatial weighting matrix are computed by the functions `scores.listw`

or `mem`

of the `adespatial`

package. These two functions are exactly identical and return an object of class `orthobasisSp`

.

```
mem.gab <- mem(listwgab)
mem.gab
```

```
## Orthobasis with 10 rows and 9 columns
## Only 6 rows and 4 columns are shown
## MEM1 MEM2 MEM3 MEM4
## 1 -0.99150068 1.1963752 -0.93642712 -0.04953977
## 2 -0.03655431 -1.6549057 0.09973653 -0.23657908
## 3 -0.66077128 -0.9284951 0.82861853 1.35172328
## 4 -1.26547947 1.0414066 0.91626372 1.02967682
## 5 1.84724812 0.4858047 -0.09173118 0.18858500
## 6 0.96155231 -0.3553900 1.15183204 -1.28527087
```

This object contains MEMs, stored as a `data.frame`

and other attributes:

`str(mem.gab)`

```
## Classes 'orthobasisSp', 'orthobasis' and 'data.frame': 10 obs. of 9 variables:
## $ MEM1: num -0.9915 -0.0366 -0.6608 -1.2655 1.8472 ...
## $ MEM2: num 1.196 -1.655 -0.928 1.041 0.486 ...
## $ MEM3: num -0.9364 0.0997 0.8286 0.9163 -0.0917 ...
## $ MEM4: num -0.0495 -0.2366 1.3517 1.0297 0.1886 ...
## $ MEM5: num 1.441 0.341 0.636 -0.027 0.558 ...
## $ MEM6: num -1.13 -2.0264 1.4244 0.0542 0.4736 ...
## $ MEM7: num -0.8913 1.0131 0.7105 -0.0329 -1.0514 ...
## $ MEM8: num -1.057 0.947 -0.705 1.404 1.619 ...
## $ MEM9: num -0.664 -0.222 -1.323 1.562 -1.022 ...
## - attr(*, "values")= num 0.12766 0.09508 0.0749 0.00496 -0.02487 ...
## - attr(*, "weights")= num 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1 0.1
## - attr(*, "call")= language mem(listw = listwgab)
```

The eigenvalues associated to MEMs are stored in the attribute called `values`

:

```
barplot(attr(mem.gab, "values"),
main = "Eigenvalues of the spatial weighting matrix", cex.main = 0.7)
```

A `plot`

method is provided to represent MEMs. By default, eigenvectors are represented as a table (sites as rows, MEMs as columns):

`plot(mem.gab)`

The previous representation is not really informative and MEMs can be represented in the geographical space as maps if the argument `SpORcoords`

is documented:

`plot(mem.gab, SpORcoords = xyir, nb = nbgab)`

Moran’s I can be computed and tested for each eigenvector with the `moran.randtest`

function:

```
moranI <- moran.randtest(mem.gab, listwgab, 99)
moranI
```

```
## class: krandtest lightkrandtest
## Monte-Carlo tests
## Call: moran.randtest(x = mem.gab, listw = listwgab, nrepet = 99)
##
## Number of tests: 9
##
## Adjustment method for multiple comparisons: none
## Permutation number: 99
## Test Obs Std.Obs Alter Pvalue
## 1 MEM1 0.7015926 3.4453729 greater 0.01
## 2 MEM2 0.5225529 2.5618940 greater 0.01
## 3 MEM3 0.4116201 2.3076371 greater 0.01
## 4 MEM4 0.0272521 0.5418607 greater 0.29
## 5 MEM5 -0.1366833 0.2527610 greater 0.45
## 6 MEM6 -0.2873389 -0.6194100 greater 0.73
## 7 MEM7 -0.4986655 -1.5993989 greater 0.94
## 8 MEM8 -0.7296825 -2.2195544 greater 1.00
## 9 MEM9 -1.0106474 -3.4273999 greater 1.00
```

By default, the function `moran.randtest`

tests against the alternative hypothesis of positive autocorrelation (`alter = "greater"`

) but this can be modified by setting the argument `alter`

to `"less"`

or `"two-sided"`

. The function is not only devoted to MEMs and can be used to compute spatial autocorrelations for all kind of variables.

As demonstrated in Dray, Legendre, and Peres-Neto (2006), eigenvalues and Moran’s I are equal (post-multiply by a constant):

`attr(mem.gab, "values") / moranI$obs`

```
## MEM1.statistic MEM2.statistic MEM3.statistic MEM4.statistic MEM5.statistic
## 0.1819528 0.1819528 0.1819528 0.1819528 0.1819528
## MEM6.statistic MEM7.statistic MEM8.statistic MEM9.statistic
## 0.1819528 0.1819528 0.1819528 0.1819528
```

Then, it is possible to map only positive significant eigenvectors (i.e., MEMs with significant positive spatial autocorrelation):

```
signi <- which(moranI$p < 0.05)
signi
```

`## integer(0)`

`plot(mem.gab[,signi], SpORcoords = xyir, nb = nbgab)`

The choice of a spatial weighting matrix is an important step and Dray, Legendre, and Peres-Neto (2006) proposed a data-driven procedure of selection based on AICc. The function `ortho.AIC`

orders variables and returns AICc for all models of one, two, …, \(p\) variables. We illustrate its use with the oribatid data-set which is available in the `ade4`

package. Data are Hellinger-transformed and then the linear trend is removed:

```
data(oribatid)
fau <- sqrt(oribatid$fau / outer(apply(oribatid$fau, 1, sum),
rep(1, ncol(oribatid$fau)), "*"))
faudt <- resid(lm(as.matrix(fau) ~ as.matrix(oribatid$xy)))
```

For instance, we consider the binary spatial weighting matrix based on the Delaunay triangulation.

```
nbtri <- tri2nb(as.matrix(oribatid$xy))
sc.tri <- scores.listw(nb2listw(nbtri, style = "B"))
AIC.tri <- ortho.AIC(faudt, sc.tri)
head(AIC.tri)
```

```
## MEM4 MEM3 MEM2 MEM1 MEM52 MEM8
## -89.78627 -91.31327 -92.54627 -93.53458 -93.80113 -93.98751
```

The minimum value and the rank of the corresponding are obtained easily:

`min(AIC.tri, na.rm = TRUE)`

`## [1] -94.18844`

`which.min(AIC.tri)`

```
## MEM5
## 7
```

Note that the order of the variables can also be obtained from the function `ortho.AIC`

by setting the `ord.var`

argument to `TRUE`

. In this case, the returned object is a list of two vectors:

```
AIC.tri <- ortho.AIC(faudt, sc.tri, ord.var = TRUE)
head(AIC.tri$AICc)
```

```
## MEM4 MEM3 MEM2 MEM1 MEM52 MEM8
## -89.78627 -91.31327 -92.54627 -93.53458 -93.80113 -93.98751
```

`head(AIC.tri$ord)`

`## [1] 4 3 2 1 52 8`

The user-friendly function `test.W`

simplifies the procedure of selection of a spatial weighting matrix. It takes at least two arguments: a response matrix and an object of the class `nb`

.

If only two arguments are considered, the function prints the results for the best model. All the results are stored in the element `best`

of the list. It contains eigenvectors and eigenvalues of the spatial weighting matrix considered and the results of the AIC-based procedure.

`tri.res <- test.W(faudt,nbtri)`

```
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM5 -94.18844 7
```

`names(tri.res)`

`## [1] "all" "best"`

`names(tri.res$best)`

`## [1] "MEM" "AIC"`

The function can also be used to estimate the best values of parameters if we consider a function of the distance. This can be illustrated with the function \(f_2=1-(x^\alpha)/dmax^\alpha\) with the connectivity defined by Delaunay triangulation. We considered the sequence of integers between 2 and 10 for \(\alpha\).

```
f2 <- function(x, dmax, y) {
1 - (x ^ y) / (dmax) ^ y
}
maxi <- max(unlist(nbdists(nbtri, as.matrix(oribatid$xy))))
tri.f2 <- test.W(faudt, nbtri, f = f2, y = 2:10, dmax = maxi,
xy = as.matrix(oribatid$xy))
```

```
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 3 -95.44519 6
```

In this case, the element `best`

contains the results for the best values of the parameter \(\alpha\).

`names(tri.f2$best)`

`## [1] "MEM" "AIC"`

Lastly, the function `test.W`

can be used to evaluate different definitions of neighborhood. We illustrate this possibility by the definition of a sequence of neighborhood by distance criteria. Firstly, we choose the range of values to be tested with an empirical multivariate variogram using the function `variogmultiv`

.

The function has been applied to oribatid mites data:

```
mvspec <- variogmultiv(faudt, oribatid$xy, nclass = 20)
plot(mvspec$d, mvspec$var, type = 'b', pch = 20, xlab = "Distance", ylab = "C(distance)")
```

We will construct ten neighborhood matrices with a distance criterion varying along the sequence of 10 evenly distributed values between 1.012 and 4 m:

```
dxy <- seq(give.thresh(dist(oribatid$xy)), 4, length = 10)
nbdnnlist <- lapply(dxy, dnearneigh, x = as.matrix(oribatid$xy), d1 = 0)
```

Then, the function `test.W`

can be applied to this list of neighborhood matrices:

`dnn.bin <- lapply(nbdnnlist, test.W, Y = faudt)`

```
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM10 -96.52487 7
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM55 -97.1354 5
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM10 -99.41955 7
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM53 -100.56 5
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM1 -99.28275 4
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM5 -97.79955 5
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM23 -95.23283 4
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM56 -95.51426 7
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM37 -93.40506 5
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## AICc NbVar
## MEM27 -97.21409 8
```

The object `dnn.bin`

is a list with the results of `test.W`

for each neighborhood matrix:

`length(dnn.bin)`

`## [1] 10`

For each neighborhood matrix, we can find the lowest :

`minAIC <- sapply(dnn.bin, function(x) min(x$best$AIC$AICc, na.rm = TRUE))`

And select the best spatial weighting matrix corresponding to a distance of 2.007 m:

`which.min(minAIC)`

`## [1] 4`

`dxy[which.min(minAIC)]`

`## [1] 2.007458`

A similar approach can be used with a spatial weighting function:

`f2 <- function(x,dmax,y) {1-(x^y)/(dmax)^y}`

It is a little bit more complicate if some parameters (here dmax) vary with the neighborhood matrix:

```
dnn.f2 <- lapply(nbdnnlist, function(x)
test.W(x, Y = faudt, f = f2, y = 2:10,
dmax = max(unlist(nbdists(x, as.matrix(oribatid$xy)))), xy = as.matrix(oribatid$xy)))
```

```
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
## Warning in nb2listw(nb, style = "B", glist = lapply(nbdist, f), zero.policy
## = TRUE): zero sum general weights
```

```
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 1.011187 -101.09 8
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 3 4 1.341641 -98.58201 5
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 4 5 1.662077 -99.31352 9
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 9 10 2.005617 -100.3819 7
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 2.332381 -102.6408 9
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 2 3 2.668333 -102.6977 7
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 3 -102.5816 6
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 3.328663 -100.8116 6
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## MEM1 2 3.65821 -99.30806 3
##
##
## AICc for the null model: -87.47112034786
##
## Best spatial model:
## y dmax AICc NbVar
## 1 2 4 -98.60206 4
```

```
minAIC <- sapply(dnn.f2, function(x) min(x$best$AIC$AICc, na.rm = TRUE))
min(minAIC)
```

`## [1] -102.6977`

`which.min(minAIC)`

`## [1] 6`

`dnn.f2[[which.min(minAIC)]]$all`

```
## y dmax AICc NbVar
## 1 2 2.668333 -102.3561 7
## 2 3 2.668333 -102.6977 7
## 3 4 2.668333 -100.9536 5
## 4 5 2.668333 -100.6003 5
## 5 6 2.668333 -101.4106 6
## 6 7 2.668333 -100.7664 5
## 7 8 2.668333 -100.6898 6
## 8 9 2.668333 -100.4506 5
## 9 10 2.668333 -100.3007 5
```

Lastly, Eigenvectors of the best spatial weighting matrix can be mapped. They are represented by the order given by the selection procedure. The third MEM explains the largest part of the oribatid community, then it is the second and the eighth:

`plot(dnn.f2[[7]]$best$MEM[, dnn.f2[[7]]$best$AIC$ord[1:3]], oribatid$xy)`

Borcard, D, and P Legendre. 2002. “All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices.” *Ecological Modelling* 153: 51–68.

Dray, S, and A B Dufour. 2007. “The ade4 package: implementing the duality diagram for ecologists.” *Journal of Statistical Software* 22 (4): 1–20.

Dray, S, P Legendre, and P R Peres-Neto. 2006. “Spatial modeling: a comprehensive framework for principal coordinate analysis of neighbor matrices (PCNM).” *Ecological Modelling* 196: 483–93.

Dray, S, R Pélissier, P Couteron, M J Fortin, P Legendre, P R Peres-Neto, E Bellier, et al. 2012. “Community ecology in the age of multivariate multiscale spatial analysis.” *Ecological Monographs* 82 (3): 257–75.