-
Notifications
You must be signed in to change notification settings - Fork 8
/
chapter_dependance.qmd
1690 lines (1128 loc) · 80.8 KB
/
chapter_dependance.qmd
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
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
---
title: "Chapter : statistical tests with dependent data"
bibliography: references.bib
execute:
freeze: auto
output:
html_document:
toc: true
toc_float: true
editor:
markdown:
wrap: sentence
---
```{r}
library(tidyverse)
```
# Introduction
### 1. What is autocorrelation ?
Statistical analyses, including many regression methods, often assume that observations are independent of each other.
Most of the time, it seems natural to analyse two different datas, X and Y, assuming their independence.
But there may still be a dependency intrinsic to each of them: this is known as autocorrelation.
Autocorrelation indicates that the independence assumption is not satisfied, which can lead to biased parameter estimates, incorrect confidence intervals and a loss of statistical power.
In statistics, autocorrelation is a measure of the linear dependence between the values of one variable itself, at different times or locations.
Checking for autocorrelation in a variable means analyzing the relationship between an observation and other previous observations in the series.
In other words, this means analyzing the possible link between what is measured at time 1 and what is measured at time 2, to see if the data at time 2 can be predicted with what we know about the data at time 1, if they are close in space or time.
Correlation between two observations is often expressed as a correlation coefficient, such as Pearson's correlation coefficient.
Values are correlated with an offset function in time or space, giving an autocorrelation coefficient which may be:
- positive: indicating that the observations tend to be similar,
- negative: indicating that the observations tend to be opposite to each other.
### 2. What types of autocorrelation ?
Autocorrelation is commonly used to analyze data for which the order or location of observations has an important significance, such as:
- Time series: data taken at the same location at several points in time,
- Spatial data: data taken at a given point in time at several locations in the same geographical area, and therefore close in space.
#### 2.1. Time series
A typical example of autocorrelation in a time series is a seasonal autocorrelation, where values are observed at a same time each year (seasons), and are therefore correlated with each other.
Temporal autocorrelation shows spatial similarities that depend on the scale used.
![**(from P. LeGouar)**](Intro_dependance.jpg)
Temporal autocorrelation in ecology is important for understanding the processes of reproduction, migration, population dispersal and species responses to seasonal environmental changes.
##### Example of temporal autocorrelation
An example of temporal autocorrelation in ecology concerns populations of animal or plant species, particularly when studying seasonal variations.
Let's imagine a study of annual fluctuations in the population of migratory birds in a given region.
Temporal autocorrelation means that the number of birds observed at a given time of the year is correlated with the number of birds observed at the same time of the previous year.
Suppose the data show positive temporal autocorrelation for a migratory bird species.
This would mean that if, for example, in April of the current year, a large number of these birds are observed, it is highly likely that in April of the previous year, a large number of these same birds were also observed.
This temporal dependence may be due to seasonal bird migrations, the availability of food resources or other cyclical environmental factors.
#### 2.2. Spatial dependence
Spatial autocorrelation refers to the correlation between observations in close spatial locations.
It shows temporal similarities that depend on the scale used.
Spatial autocorrelation is important for understanding patterns of species distribution, population dispersal, biotic interactions and ecological processes at different spatial scales.
Spatial autocorrelation can be due to various factors, such as :
- Microclimate effects: local characteristics of terrain, vegetation or topography can influence temperature on a small scale,
- Dispersion phenomena: the propagation of heat, wind, humidity or other environmental factors can cause spatial autocorrelation,
- Biotic interactions: interactions between plant, animal and micro-organism species in an ecosystem can also influence the spatial distribution of variables such as temperature.
##### Example of spatial dependence
An example of spatial autocorrelation is the observation of temperature distribution in a large forest.
Let's imagine a large, dense forest with many weather stations measuring temperatures at different locations.
Spatial autocorrelation would manifest itself in the fact that temperatures recorded at nearby locations in the forest are correlated, i.e. they tend to be similar.
Suppose the data show positive spatial autocorrelation.
This would mean that, if you have two weather stations located close to each other (say, a few meters apart), the temperatures measured at these two stations at any given time are highly positively correlated.
In other words, when one of the stations records a high temperature, the other station will also tend to record a high temperature, and vice versa.
### 3. Why is this important?
Autocorrelation is an essential concept in statistics, as it can have a significant impact on data analysis.
In the presence of pseudoreplications, meaning that several different values are thought to be taken when in fact they are dependent on each other, it may appear that one variable has a significant effect on another, when in reality this effect is biased due to the temporal or spatial dependence of the data.
It can therefore have a significant impact on statistical analysis, and often requires adjustments to ensure that results are reliable and unbiased.
Autocorrelation checking is an essential process in statistics to guarantee the accuracy, reliability and validity of analyses.
It enables appropriate statistical methods to be applied, biases to be corrected and informed modeling and forecasting decisions to be made.
It is also essential for understanding patterns of variation in biological data.
<br />
# I. Time dependence
Temporal dependency means that values sampled at a given time can be related to the value sampled before.
Therefore, the variable can be correlated to itself and create **pseudoreplication**.
Diverse packages on R are useful to manage time series, such as TSA (Time Series Analysis).
```{r, echo = F, error=FALSE}
## Install and load the packages required
library(tidyverse) # Data cleaning, ordering, etc.
library(lubridate) # Helps manage the date format data
library(ade4) # Study ecological data
library(nlme) # glm and gls
library(car) # Type II anova.
#To install older version of a package to use the arimax function.
# require(devtools)
# install_version("TSA", version = "1.2.1", repos = "http://cran.us.r-project.org")
# install_version("TSA", version = "1.2.1", repos = "https://pbil.univ-lyon1.fr/CRAN/")
library(TSA) # Time series analysis (autocorrelation tests)
```
### 1. The dataset
To understand better how to handle time series, data are used to illustrate the process from an experiment on vultures in time.
```{r}
## Load the data
dataT=read.table("vautour.txt",header=T,dec=",", stringsAsFactors = T)
head(dataT) # What are the data
summary(dataT) # Summary of the spread of each variable
```
Two colonies 'C' and 'O' are studied from 1972 to 2003.
In particular, the reproductive success (success), the growth rate (growth), the density of population (denstot) and the resource availability (ratio) of the colonies.
```{r}
## Arrange the dataset
dataT_ordered=dataT %>% arrange(col,year)
## Create 2 sub-datasets for the 2 colonies
C=dataT[dataT$col=="C",]
O=dataT[dataT$col=="O",]
```
### 2. Data exploration
```{r}
## Data representation
# Window size
par(mfrow=c(2,3),mar=c(2,4,2.5,0),oma=c(2,2,2,4))
## Colony C
plot(C[,2],C[,3], xlab="year", ylab="mean reproductive success")
lines(C[,2],C[,3])
plot(C[,2],C[,4], pch=17, xlab="year", ylab="Colony growth rate")
lines(C[,2],C[,4])
plot(C[,2],C[,5],pch=16, xlab="year", ylab="Number of pairs")
lines(C[,2],C[,5])
mtext(text="Parameter for Colony C",side=3,line=-1,outer=TRUE)
## Colony O
plot(O[,2],O[,3], xlab="year", ylab="mean reproductive success")
lines(O[,2],O[,3])
plot(O[,2],O[,4],pch=17, xlab="year", ylab="Colony growth rate")
lines(O[,2],O[,4])
plot(O[,2],O[,5], pch=16, xlab="year", ylab="Number of pairs")
lines(O[,2],O[,5])
mtext(text="Parameter for Colony O",side=1,line=-14.5,outer=TRUE)
```
Through time, all the variables do not have the same behaviour for both of the colonies.
Some variables fluctuate in time like the O colony growth rate, or can evolve gradually, like the number of pairs for the O colony.
The irregularity in time of the variables can mean that the data are correlated to the precedent ones.
A key parameter to know if there really is a correlation of the variable to itself is **autocorrelation**.
### 3. Autocorrelation
Autocorrelation is given by the expression: $$\rho(k) = \frac{E(X_tX_{t+1}) - E(X_t)E(X{t+1})}{\sqrt{E(X_t^2)E(X_{t+1}^2)}}$$
With the discrete sequence $\{X_1, ... X_T\}$ with $X = \{X_t\}_{t\in\mathbb{N}}$ as $E[X_t^2] < \infty$ for $t\in\mathbb{Z}$.
In other words, the values are compared to their neighbors to conclude if they are correlated or not.
The autocorrelation can be detected in a **non-stationary process**.
**Stationnarity** is the constistency of the statistical parameters of a series such as the mean and the variance (Figure).
- In non-stationary processes, these parameters vary in time.
- In stationary processes, the values gravitate around a mean value within the same range.
![Stationary and non-stationary processes](C:/Users%5CUser%5COneDrive%5CBureau%5CMODE%5CM2%5COCR%5CStationarycomparison.png)
Given the precedent plots displayed, one would think there are both stationary and non-stationary processes within the different variables.
Tests will allow us to know with more details if there is a dependency.
On R, autocorrelation can be known thanks to the function acf (Auto- and Cross- Covariance and -Correlation Function Estimation) of the TSA package.
**Important note:** In order for the acf() function to work, the data have to be sorted according to time with one value per time step.
To do so, *NA* or mean values can be added if necessary to fill the dataset accordingly.
The package lubridate can be useful to manage dates data.
```{r}
par(mfrow=c(2,3),mar=c(2,4,4,2))
acf(C[,3], main="Reproductive success for C")
acf(C[,4], main="Growth rate for C")
acf(C[,5], main="Density for C")
acf(O[,3], main="Reproductive success for O")
acf(O[,4], main="Growth rate for O")
acf(O[,5], main="Density for O")
```
The plots display the correlation coefficient at the associated time step with the neighbors (lag).
The blue dotted lines represent the significance threshold of correlation.
These plots allow us to know if the time dependency exists and what profile it follows.
- Growth rate of O and reproductive success of C : **White noise**. There is no significant value and particular pattern. Therefore, no time dependency.
- Density of C and reproductive success of O : **Short positive dependency**. Some values are significantly positively correlated and a pattern emerges (cyclic here).
- Growth rate of C : **Short negative dependency**. Some values are significantly negatively correlated and a pattern emerges (cyclic here).
- Density of O : **Long term dependency**. Several values are significantly correlated and a pattern appears (linear here).
For the following analyses, we will focus on the Density of the O colony.
The time is considered an **autocovariable**.
### 4. Models
Now that the time dependency have been proved to be true thanks to autocorrelation test, it is possible to proceed with adapted statistical tests.
Different models can be used according to the type of dependency found in the data.
The goal will be to apply the model and verify that there is no time dependency left on the **residuals**.
#### 5. Least square ordinary regression
Least square ordinary regression methods are used during the process of selection of the model.
$$Y = \beta_0 + \beta_1X_1 + \beta_2X_2+...+\beta_nX_n + \epsilon$$ - $Y$ is the dependent variable - $X_n$ are the explanatory variables - $\beta_n$ are the coefficients of $X$ - $\epsilon$ is the random error term or **residuals**
To select a relevant model is the same as selecting a model that maximizes **likelihood** and minimize the **number of parameters**.
#### 6. Weak time dependency
If the data display a weak time dependency, it is possible to simply add the time as a random variable to the model.
$$Y\sim\alpha.X + time$$
**Important note:** If the variance of $Y$ increases beforehand, you will need to perform a **logarithmic** transformation.
If the variance of $Y$ decreases beforehand, one would use an **exponential** transformation.
Here, the variance increases.
A useful package on R to manage glm and gls is be nlme.
```{r}
par(mfrow=c(1,2))
hist(O$denstot, main = paste("Histogram of total density", "\n" ,"of O colony "), xlab = "Density")
hist(log(O$denstot), main = paste("Histogram of total density", "\n" ,"of colony O (log transformation)"), xlab = "Log(Density)")
```
The plots indicates that a log transformation could help better manage the data.
A simple lme model can be used (Linear Mixed-Effects Model) with X = log(total density) and Y the ressource availability.
We can check if the model is satisfactory if the residuals are not significant.
```{r}
## Model construction
modlme=lme(log(denstot) ~ ratio, data=O, random=~1|year)
## Residuals
par(mfrow=c(1,1))
acf(modlme$residuals[,1], main = "Residuals of the lme model")
```
```{r}
## Residuals
summary(modlme)
```
The plot shows there is still a time dependency in this kind of model.
The correlation value is too strong, the model is not adapted to this case.
### 7. Medium to strong time dependency
If a medium or a strong time dependency are proven, adding time as a random effect will not be effective enough.
It will be necessary to use **Time Models**.
#### 7.1 Stationary process
To deal with stationary processes, it is possible to use GLS (Generalized Leat Square) models.
Three major GLS models can be used.
- **AR** or **Autoregressive Model** states that a value $Y_{(t)}$ will depend on the values $Y$ at $p$ the precedent times.
The model adds a an autoregressive coefficient $\phi_p$ to each $Y_{(t-1)}$ value.
In summary, this model explains the current value as a linear combination of the precedent values with $\epsilon_t$ the random component or white noise.
$$Y_{(t)}=\mu+\phi_1 \times Y_{(t-1)}+\phi_2 \times Y_{(t-2)} + ...+\phi_q \times Y_{(t-p)}+\epsilon_{(t)}$$
- **MA** or **Moving Average model** functions in the same way as the AR model.
However, it will consider that the current value of a series is the weighted linear combination of the previous $q$ errors.
$$Y_{(t)}= \mu+ \theta_1 \times \epsilon_{(t-1)}+\theta_2 \times \epsilon_{(t-2)}) + ...+\theta_q \times \epsilon_{(t-q)} $$
- **ARMA** is a combination of **AR** and **MA** models and is a linear combination of the $p$ past values and $q$ past residual errors.
$$Y_{(t)}=\mu+\phi_1 \times Y_{(t-1)}+\phi_2 \times Y_{(t-2)} + ...+\phi_q \times Y_{(t-p)}+\epsilon_{(t)} + \theta_1 \times \epsilon_{(t-1)}+\theta_2 \times \epsilon_{(t-2)}) + ...+\theta_q \times \epsilon_{(t-q)} $$ The gls() function covers the three types of model depending on the correlation structure of the data.
```{r}
## Model construction
modgls=gls(log(denstot) ~ year + ratio, data=O, correlation = corAR1(value=0.9,form = ~ year))
acf(modgls$residuals, main = "Residuals of the gls model")
```
```{r}
summary(modgls)
```
The autoregressive coefficient Psi is equal to 0.72.
The autocorrelation is not completely removed.
Still, the variable 'ratio' does not have a significant effect on $Y$ anymore.
We can compare the two previous models thanks to an analysis of their variance (anova).
```{r}
anova(modgls,modlme)
Anova(modgls)
Anova(modlme)
```
```{r}
modgls2=gls(log(denstot) ~ ratio, O, correlation = corAR1(value=0.9,form = ~ year))
anova(modgls2,modlme)
anova(modgls2,modgls)
anova(modgls2)
```
The second models fits the data better than the first.
We can conclude that the autocorrelation made the ratio variable effect wrongly significant.
In order to chose the range of application of the parameters $p$ and $q$, it is possible to use a **partial correlation function** that will display the direct correlation link between the time steps.
```{r}
par(mfrow=c(1,2))
acf(log(O$denstot), main="Density for O")
pacf(log(O$denstot), main="Density for O")
```
So, the order $p$ (AR) is generally chosen based on significant lags in the PACF, and the order $q$ (MA) is generally chosen based on significant lags in the ACF.
#### 7.2 Non-stationary process
To manage non-stationary processes, another model called **ARIMA** is used.
A new parameter $d$ is introduced with $p$ and $q$.
**ARIMA(p,d,q)** models are a combination of **AR** models with $p$, **MA** models with $q$, and finally includes difference of order $d$.
This difference helps make the series stationary by subtracting from the current value the values from the previous $d$ time points.
In general, we will choose an order $d=1$ when the series shows a **constant trend**, and an order $d=2$ when the series shows a **fluctuating trend**.
On R, ARIMA models can be used with the function arima().
The argument 'order'contains the values of $p$, $d$ and $q$, in that order.
```{r}
arima1=stats::arima(log(O$denstot),order=c(1,1,1),xreg=O$ratio)
arima1
```
Different values of these parameters can be tested to find the best model with the AIC criteria displayed in the results.
$d$ must always be equal to 1.
The ARIMA function can help make **predictions**.
- prediction: now we will use this statistical model to make prediction under a specific scenario
- scenario: constant and maximum availability in resource (ratio=1) during 20 years
```{r}
newxreg=data.frame(O.ratio=rep(1,20))
predict(arima1, n.ahead = 20,newxreg)
pred=predict(arima1, n.ahead = 20,newxreg)
par(mfrow=c(1,1))
obspred=rbind (data.frame(denstot=O$denstot), data.frame(denstot=exp(predict(arima1, n.ahead = 20,newxreg)$pred)))
plot(obspred$denstot, main = "Prediction on a 20 years scenario", ylab ="Predicted observations for the total density of O")
lines(obspred$denstot)
```
### 8. To go further
The **ARIMAX** model, which stands for **AutoRegressive Integrated Moving Average with eXogenous inputs**, is an extension of the ARIMA model that incorporates additional exogenous (external) variables.
The ARIMAX model is particularly useful when there are external factors that can influence the time series to model.
The eXogenous (X) Inputs are additional external variables that are not part of the time series being modeled but may influence it.
These variables are included in the model to capture their impact on the dependent variable.
The ARIMAX model allows to account for the influence of external factors when modeling time series data, making it more flexible and potentially improving the accuracy of predictions.
# II. Spatial dependance
Spatial relationships are multidirectional and multilateral.
They are distinct, in this sense, from temporal relationships, which allow only sequential relationships along the past-present-future axis.
These are the steps that will be followed here :
**Define connection criteria**: We establish criteria for connections, determining that point A is connected to point B, which in turn is connected to point C, and so forth.
This linkage needs to be translated into a numerical value, providing a quantitative measure of connectivity.
**Transform information into numerical value:** The established connection criteria are converted into numerical values, creating a quantifiable representation of the relationships between points.
This numerical transformation is crucial for subsequent analyses.
**Calculate autocorrelation coefficient:** Utilizing the numerical values derived from the connection criteria, we compute an autocorrelation coefficient.
This statistical measure helps quantify the degree of spatial autocorrelation, providing insights into the patterns of similarity or dissimilarity among neighbouring points.
**Integrating spatial autocorrelation into analyses:** Upon detecting spatial autocorrelation, the next step involves seamlessly integrating this information into our analytical processes.
This will be discussed later in this chapter.
------------------------------------------------------------------------
## **1 - Codifying the neighbour structure**
------------------------------------------------------------------------
### **1.1 Defining neighbours**
#### **1.1.1 Characteristics of the relationships between spatial objects**
Consider a surface ℜ.
This surface can be divided into n mutually exclusive zones.
Two adjacent zones are separated by a common boundary.
::: callout-note
## Mathematical definition of spatial relationships
Spatial relationships $B$ are a subset of the Cartesian product $ℝ²× ℝ² = {(i, j) : i ∈ ℝ ², j ∈ ℝ²}$ of couples $(i, j)$ of spatial objects, i.e. all couples $(i, j)$ such that $i$ and $j$ are both spatial objects identified by their geographical coordinates, and such that $(i, j)$ is different from $(j,i)$.
A spatial object cannot be linked to itself: $(i,i) ⊈ B$.
Moreover, if $(i, j) ⊆ B$ and $(j,i) ⊆ B$ for all couples of spatial objects, the spatial relationships are said to be symmetrical (Tiefelsdorf, 1998).
:::
*Figure 2.1* illustrates the codifying process of spatial relationships.
This approach makes it possible to systematically transcribe the complexity of geographic space into a final set of data analyzable by a computer.
![**Figure 2.1 :** Codifying spatial correlation (Tiefelsdorf, 1998)](image-1363470329.png){width="517"}
First, the study zone is divided into mutually exclusive areas.
Each area contains a reference point (often its centroid).
Then, the spatial relationships can be specified by a neighbourhood graph connecting the areas considered to be neighbouring, or by a matrix containing the geographical coordinates of the reference points.
The third step consists in coding the graph in a neighbourhood matrix, or transforming the geographic coordinates into a distance matrix.
The neighborhood matrix measures the similarity between observations.
A value greater than zero indicates that the observations are considered neighbors.
For example, in the case of the binary matrix shown in figure 1 :
![](image-232123009.png){width="470" height="59"}
Conversely, the distance matrix measures dissimilarity between zones.
The higher $d_{ij}$ , the more different the zones.
With, if an Euclidian distance is used : $d_{ij}=\sqrt{(x_i-x_j )^2+(y_i-y_j )²}$, $α$ and $β$ being the geographical coordinates of the observations.
The spatial dependence structure may not be geographical.
Any relevant dual relationship may be used to define a neighbourhood graph.
For instance:
- Genetic proximity between families
- Cultural exchanges
- Migratory flows
#### **1.1.2 The "list of neighbours" object in R**
| Package *spdep* makes it possible to define the relationships between spatial objects.
In R, the class of an object defines all its properties and how the statistician can use it.
Neighbourhood relationships are recorded in an object of class nb.
Assume n spatial observations and `neighbours_nb` the spatial object containing the associated neighbourhood relationships.
`neighbours_nb` is a list of length n.
Each element `[i]` of the list contains a vector with the index of the neighbours of the item indexed i.
If `[i]` does not have neighbours, the list contains only 0.
The list also contains a vector of characters corresponding to the attributes of each neighbourhood zone, as well as a logical value indicating whether the relationship is symmetrical.
The main information about the object `neighbours_nb` can be derived using the function: `summary(neighbours_nb)`
Let's load a few packages useful for analysis :
```{r}
# Spatial series analysis : packages
#Libraries for spatial series analyses
library(spdep) #to extracted the neighbors is a spatial data
library(ade4) #used for plotting spatial data
library(spatialreg) #used for spatial data modelling
library(gwrr) #to run geographically weighted regression
#these two libraries (spdep and ade4 also provive function that make object compatibles between library)
```
Let's then begin the analysis of a first dataset dealing with richness of birds from agricultural habitats and of plant species from permanent pastures at the scale of Ille-et-Vilaine, "*div35*".
```{r}
#Spatial series analysis : dataset div35
#For a REGULAR configuration of points
# 0 - Exploring the dataset
#Investigate the data set of plant and bird richness at the scale of Ille-et-Vilaine
#Import dataset
div35=read.table("div35.txt", header=T, dec=",")
#Summary
summary(div35)
#Plot the spatial distribution of the points
par(mfrow=c(1,1))
plot(div35[,3:4])
#This is a map of Ille-et-Vilaine department.
#We can point St-Malo in the North and Rennes in the center for of the map, for example.
#The spatial configuration of the points is regular.
```
Let's also start the analysis of second dataset dealing with chemical measures in a set of ponds, "*mares*".
```{r}
#Spatial series analysis : dataset mares
#For an IRREGULAR configuration of points
# 0 - Exploring the data set
#The data set is issued from a study that analysed several physico-chemical properties of several ponds. each pond is geolocalised
#Import the dataset
mare=read.table("mares.txt",header=TRUE)
#Summary
summary(mare)
row.names(mare)<-seq(1,31,1) #we name the ponds from 1 to 31.
#We can plot the ponds on a map according to their coordinates
#We can see that the configuration of the ponds is not regular over space
par(mfrow=c(1,1))
plot(mare$x,mare$y)
```
### **1.2 Connection among neighbours**
In order to identify points in close proximity, it is imperative to establish specific criteria beyond merely considering distance along the x-axis.
This approach will enable us to clearly define the parameters that determine proximity between points in our analysis.
#### **1.2.1 Neighbourhood criterion**
**a) Defining neighbours based on distance**
Once we have a set of points spread across the territory, we can calculate the distance between them.
These points may be specific locations where the information has been observed, or points representative of each zone, for example their centroid.
In this case, the underlying assumption is that the distribution of the variable within each zone is sufficiently homogeneous to approximate it to a single point.
Neighbourhood graphs materialise the links between the various entities.
They are defined in such a way that they represent the underlying spatial structure as closely as possible.
There exists many different types of neighbourhood graphs.
**a.1) Neighbourhood graphs based on geometric concepts :**
• **Delaunay's Triangulation** Delaunay's Triangulation is a geometric method that connects points into triangles such that the minimum angle of all triangles is maximised.
![**Figure 2.2 :** Delaunay triangle (source : Pascaline Le Gouar)](image-1235549498.png){width="190"}
• **Sphere-of-influence based graph**
The sphere-of-influence based graph links two points if their "circles from the nearest neighbour" overlap.
The "circle of the nearest neighbour" of point $P$ is the largest circle centred in $P$ and that contains no other points than $P$.
All points in the study set are not necessarily interconnected.
![**Figure 2.3 :** The sphere of influence graph (edges) of a set of points, showing their nearest neighbour circles (source : Toussaint & Emirates, 2014)](image-1767308128.png){width="361"}
• **Gabriel's graph**
Gabriel's graph links two points $p_i$ and $p_j$ if and only if all other points are outside the circle with diameter $[p_i, p_j]$.
![**Figure 2.4 :** Gabriel graph (source : Pascaline Le Gouar)](image-1231590752.png){width="227"}
• **Graph of relative neighbours**
The graph of relative neighbours considers that two points $pi$ and $pj$ are neighbours if $$d(p_i,p_j) ≤ max [d(p_i,p_k),d(p_j,p_k)] \,\,\,\,\,\,∀k = 1,..., n k\ne i,j$$ with $d(p_i, p_j)$ the distance between $p_i$ and $p_j$.
![**Figure 2.5 :** Graph of relative neighbours (source : J. O'Rourke, 1982)](image-1901385633.png){width="304"}
For example, here's the application of these geometric concepts to Parisian districts :
![**Figure 2.6 :** Four neighbourhood graphs of Parisian districts based on geometric concepts](image-691203873.png){width="610"}
**Application in R :**
Coming back to our "*mares*" dataset, we can explore these geometric concepts :
```{r}
#Spatial series analysis : dataset "mares"
#Connexions among neighbours
#For an IRREGULAR configuration of points
#Several criteria exist to connect points that are irregularly distributed.
# It is necessary to tests several criteria to test different hypothesis and different scale of the phenomena.
#Criteria are determined when we have an a priori idea of the processes involved.
#For example, may be due to eutrophication (lots of nutrients added, lots of organisms making excrements), heterogeneity in the pond - pond animal and plant communities may have an influence
#Two ponds should look alike because they're on the same soil
#Animals and plants: depending on the dispersal abilities of these species, will influence a pattern in water quality.
#Ponds supposed to be in a specific landscape will look alike.
#Larger-scale phenomena (flooding): if ponds are distributed in a valley, all the ponds will be homogenized.
# 1 - Examine the connection between points ====================
# for first step extract the coordinate
marexy=mare[,11:12]
# A - The gabriel graph method ---------------------------
#Note that for this function, the coordinates should be in a matrix:
mare.gab<-gabrielneigh(as.matrix(marexy))
# We can represent it spatially with this function:
s.label(marexy,clabel=.5, cpoint=.1,neig=nb2neig(graph2nb(mare.gab)))
# the part nb2neig(graph2nb(mare.gab)) is necessary to translate the object from one package to another
# s.label is a function from ade4 and we need to translate our object from tripack->spdep->ade4
# if we want a synthetic view of the neighbor relationship we can use this function:
graph2nb(mare.gab)
#31 pools, 47 links, each point will have an average of 1.5 neighbors, 3 regions have no links (31, 13, 28) = considered as non-neighbors
gab=graph2nb(mare.gab)
#note: the ponds 13, 28 and 31 seems to have no link. This is just because they do not appear
#in the $from object but only in the $to (this is a bug in the function - we would assume)
# B - The Delaunay triangulation graph method -----------
#to compute a neighbor relationship using the delaunay triangulation method
#We should use this function using directly the coordinates:
mare.tri<-tri2nb(marexy)
# graphic representation:
s.label(marexy,clabel=0.6,cpoint=1,neig=nb2neig(mare.tri))
# synthetic view:
mare.tri
#With this connection criterion, ponds have an average of 5.2 neighbors: more connections are created
```
**a.2) Neighbourhood graphs based on nearest neighbours**
A second method consists in selecting the k closest points as neighbours.
This method has the advantage that it leaves no point without a neighbour, which is not required when conducting a spatial analysis, but generally offers a better reflection of reality.
The choice can also be made to keep only the points located at a certain distance.
For example :
- Nearest neighbour
- Two nearest neighbour
- Three nearest neighbour
- Neighbours at a minimum distance
etc.
Still using the example of Parisian districts, here's the nearest neighbour application:
![**Figure 2.7 :** Four graphs based on the nearest neighbours of Parisian districts](image-1634475889.png)
The `nbdists` function of R can be used to calculate the vector of distances between neighbours.
It makes it possible to determine the minimum distance $d_{min}$ above which all points have at least one neighbour, then the `dnearneighb` function allows to keep as neighbours only the points between distances $0$ and $d_{min}$.
**Application in R**
Let's see what it looks like with the "*mares*" dataset with the criterion of neighbours at a minimum distance (30 km) :
```{r}
#Spatial series analysis : dataset "mares"
#Connexions among neighbours
#For an IRREGULAR configuration of points
# 1 - Examine the connection between points ====================
# for first step extract the coordinate
marexy=mare[,11:12]
# C - The neighbours at a minimum distance criterion --------------------------
# to compute a neighbor relationship using a distance criteria we should use this function using directly the coordinates as a matrix and indicating the min and max of distance.
# the max of distance is defined depending on our knowledge of the system or for different hypothesis we want to test:
mare.dnear<-dnearneigh(as.matrix(marexy),0,30) #here we specify that ponds are connected if they are appart from less than 30m
# Graphic representation:
s.label(marexy,clabel=0.6, cpoint=1,neig=nb2neig(mare.dnear))
# Synthetic view:
mare.dnear #in that specific case: the ponds 4 and 31 are not connected to any other pond
```
Let's then compare the current criterions in the *mares* dataset:
```{r}
#Spatial series analysis : dataset "mares"
#Connexions among neighbours
#For an IRREGULAR configuration of points
par(mfrow=c(1,3))
s.label(marexy,clabel=.5, cpoint=.1,neig=nb2neig(graph2nb(mare.gab))) #Gabriel
s.label(marexy,clabel=0.6,cpoint=1,neig=nb2neig(mare.tri)) #Triangulation
s.label(marexy,clabel=0.6, cpoint=1,neig=nb2neig(mare.dnear)) #Neighbour at a min distance
#When comparing gabriel, triangulation and distance results, we observe different connexions.
#We can think of several biological processes according to the different critera :
#Gabriel: we can imagine gradients in relation to the ground - we can also imagine a small hydro network that connects the two.
#Triangular: Flooding in a valley would make a connection on a larger scale
#Neighbour at a min distance: composition of plant and animal communities, dispersion
```
**b) Defining neighbours based on contiguity**
When the areal data consist in a partition of the entire territory, the concept of "distance between observations" can become quite ambiguous.
The example below illustrates the limits of using the distance between centroids to define the notion of neighbourhood.
::: callout-note
##### Example - Ambiguity of the notion of distance between centroids
Let $R_1$, $R_2$, $R_3$ be three distinct zones.
It can be considered that since $R_2$ and $R_3$ are separated in space, but both are adjacent to $R_1$, they are both closer to $R_1$ than to one another.
However, the centroids in these zones are equidistant from each other (see *Figure 2.8*).
Summarising the proximity between zones by the distance between the centroids results in a partial loss of the richness of the spatial relationships.
![**Figure** **2.8** - Left: three zones - Right: distance between centroids (Source: Smith, 2016)](image-2131022957.png){width="579"}
:::
**Types of contiguity :**
- **Rook contiguity**
In the sense of Rook contiguity, neighbours have at least two common boundary points (a segment).
This matches the movement of the Rook in chess.
![**Figure 2.9 :** Rook configuration (source : Pascaline Le Gouar)](image-1114536119.png){width="211"}
- **Queen contiguity**
For two zones to be adjacent in the sense of Queen contiguity, they only need to share one common boundary point.
This matches the movement of the Queen in chess.
![**Figure 2.10 :** Queen configuration (source : Pascaline Le Gouar)](image-594369851.png)
- **Bishop contiguity**
In the sense of bishop contiguity, two regions are considered as spatial neighbors if they meet at one point ;it is based on the existence of common vertices between two spatial units.
This matches the movement of the Bishop in chess.
![**Figure 2.11 :** Bishop configuration](image-2117709373.png){width="177"}
Still using the example of Paris districts, here's an example of how continuity can be applied:
![**Figure 2.12 :** Queen and Rook contiguity in Paris districts](image-351519621.png){width="474"}
**Application in R**
We can return to the first dataset to investigate the configuration of the points :
```{r}
#Spatial series analysis : dataset div35
#Connexion between points
#For a REGULAR configuration of points
# 1 - Examine the connection between points ====================
# To investigate the various criteria of connection between points we need to
# extract the coordinates in a separate object
xy=div35[,3:4]
# A - If we choose that points are connected to their 4th nearest neighbors (the "rook connection") ----
div.knear4<-knearneigh(as.matrix(xy),4) #to extract the neighbour points
div.knear4
knn2nb(div.knear4) #knn2nb convert the object returned by the function knearneigh to a class nb (necessary for the subsequent function)
# The information printed corresponds to all the links that have been created compared to all those that could have been created (272)
knn2nb(div.knear4) # transformation to get as many lists than points/observation
plot.new()
plot(knn2nb(div.knear4), xy, add=TRUE)
# B - If we specify a queen type of connection with points connected to their 8th nearest neighbors (the "queen connection") ----
div.knear8<-knearneigh(as.matrix(xy),8)
div.knear8
knn2nb(div.knear8)
# Plotting the connection between points for each criterion
## MPE this in an intercative command, can not be used with a markodon document.
par(mfrow=c(1,2))
ade4::s.label(xy, ylim=c(min(div35[,4])-10000,max(div35[,4])+10000), xlim=c(min(div35[,3])-10000,max(div35[,3])+10000), clabel=0.6, cpoint=1,neig=nb2neig(knn2nb(div.knear4)))
ade4::s.label(xy, ylim=c(min(div35[,4])-10000,max(div35[,4])+10000), xlim=c(min(div35[,3])-10000,max(div35[,3])+10000), clabel=0.4, cpoint=1,neig=nb2neig(knn2nb(div.knear8)))
#Be careful, the distance between points is not always 10km !
par(mfrow=c(1,1))
```
**c) Defining neighbours based on the optimisation of a trajectory**
Some methods such as spatial sampling require prior data sorting.
When the latter are characterised by two variables (i.e. their geographical coordinates in the plan), how to choose a sorting method becomes a complex theoretical problem.
One solution consists in running a path along all the points, and sorting them by their order of appearance when the path is taken.
The neighbours of a given point are then the points located just before or just after along the path.
Out of the set of possible paths, some have characteristics that are better suited to the desired objectives, such as, for instance, reducing sampling variance.
- **Shortest path (Hamilton path)**
This is the case of the **shortest path**.
It minimises the sum of the distances between two consecutive points.
This path, which does not set any particular constraints on the starting or arrival point, is known in the literature as the **Hamilton path** (*Figure 2.13b*).
- **Hamiltonian cycle**
A particular and well-known case of shortest path is that of the travelling salesman.
It represents the path which a travelling salesman must take to visit all his customers, minimising the distance travelled and managing to return home in the evenings.
Such a path corresponds to a **Hamiltonian cycle** (*Figure 2.13c*).
The `TSP` package in R allows to compute this thanks to the function `insert-dummy` and the `concorde` program.
- **General randomized tessellation stratified method (GRTS)**
The ***general randomized tessellation stratified method*** (GRTS, Stevens Jr et al. 2004) is popular in spatial sampling, as it makes it possible to get a spatially-balanced sample for a finite population of individuals (distinct and identifiable units of dimension 0 of a discrete population, e.g. trees in a forest), a linear population (continuous units of dimension 1, e.g. rivers) or a population of surfaces (continuous units of dimension 2, e.g. forests).
The idea of the method is to project the coordinates on a unit square, then cut this square into four cells, each of which is cut again into four sub-cells, etc.
To each cell, a value is assigned, resulting from the order in which the division was carried out, ultimately making it possible for the units to be placed on the path going through the two-dimensional space.
This method can be implemented with package `spsurvey` in R (Kincaid et al. 2016).
However, with the GRTS method, *large jumps* (*Figures 2.13d*) are created along the paths, which can affect the accuracy of the estimates.
![Figure 2.13 : Looking for paths that cross through all the neighbourhoods of Paris](image-1294038395.png)
### **1.3 Neighborhood weights**
#### **1.3.1 From a list of neighbours to a weight matrix**
Once the neighbourhood graph has been defined and codified into a list of neighbours, the link between points $i$ and $j$ is transformed into the element $w_{ij}$ of the weight matrix $W$.
The weight matrix $W$ is the "formal expression of spatial dependency between observations" (Anselin et al. 1988)\
Formally, the weights express the neighbor structure between the observations as a $n\times n$ matrix $W$ in which the elements $w_{ij}$ of the matrix are the spatial weights:
![](image-672436714.png){width="285"}
**Defining the weight matrix**
Most commonly, the weight matrix is a binary contiguity matrix (see *Figure 2.13*):
> ![](image-562063019.png){width="431"}
![**Figure 2.14 :** Binary weight matrix](image-497770688.png){width="592"}
The weight matrices can also take into account the **distance** between the geographical zones, as relationships becoming smaller with distance: $1$ if $d < d_0 - 0$ otherwise, $\frac{1}{d^α}$ , or $e^{−αd}$ with $α$ an estimated or predetermined parameter.
Using a maximum distance beyond which $w_{ij} = 0$ makes it possible to limit the number of components with a value different from zero.
When the size of the zones is heterogeneous, this method increases the risk of a considerable variability in the number of neighbours.
Lastly, certain matrices take the **strength of relations** between the zones into account.
For example, weight can be defined by $\frac{b^α_{ij}}{d^β_{ij}}$ with $b_{ij}$ a measure of the strength of relationships between zones $i$ and $j$ (which is not necessarily symmetrical), such as the percentage of common boundaries, the total population, the wealth and $d_{ij}$ the distance between the zones.
**Application in R**
Let's go back to our previous data set, *div35*, to see the transformation into a binary matrix :
```{r}
#Spatial series analysis : dataset div35
#Weight matrix
#For a REGULAR configuration of points
# about the code :
# nb2neig(knn2nb()) transform the matrix to subsequently use it is ade4
# neig = define the connection matrix between points
# Transform the list of the names of the neighbors to a 2 by 2 matrix of neighboors with "0" and "1" -----
w2<-knn2nb(div.knear4)
str(w2) # we have as much list as points
nb2neig(knn2nb(div.knear4)) # these lists are transformed in a 2D matrix with 0 or 1.
```
#### **1.3.2 Weight matrix standardisation**
The sum of the weights of the neighbours of a zone is called its *degree of connection*.
If the weight matrix is not standardised ("$B$" coding scheme), the degree of connection will depend on the number of its neighbours, which creates heterogeneity between the zones.
According to Tiefelsdorf (1998), four types of standardisation can be distinguished:
--- **Line standardisation** (*"W" coding scheme*): for a given zone, the weight ascribed to each neighbour is divided by the sum of the weights of its neighbours : $∑^n_{j=1}w_{ij} = 1$.
This standardisation makes the interpretation of the weight matrix easier, because $∑^n_{j=1}w_{ij}x_j$ represents the average of variable $x$ on all neighbours of observation $i$.
Each weight $w_{ij}$ can be interpreted as the fraction of spatial influence on observation $i$ ascribable to $j$.
In contrast, such standardisation implies a certain degree of competition between neighbours: the fewer neighbours a zone has, the greater their weight.
Moreover, when weights are inversely proportional to the distance between the zones, row standardisation makes them difficult to interpret.
--- **Global standardisation** *("C" coding scheme)*: weights are standardised so that the sum of all weights is equal to the total number of entities.
All weights are multiplied by $\frac{n}{∑^n_{j=1}∑^n_{i=1}w_{ij}}$.
--- **Uniform standardisation** (*"U" coding scheme*): weights are standardised so that the sum of all weights equals $1$ : $∑^n_{j=1}∑^n_{i=1}w_{ij} = 1$.
--- **Standardisation by variance stabilisation** (*"S" coding scheme*): let **q** be the vector defined by : **q** = $(\sqrt{∑^n_{j=1}w^2_{1j}}, \sqrt{∑^n_{j=1}w^2_{2j}},...., \sqrt{∑^n_{j=1}w^2_{nj}})^T$.
Let matrix $S^∗$ = $[diag(q)]^{−1}W.^4$ From $S^*$, we calculate $Q = ∑^n_{j=1}∑^n_{i=1} s^∗_{ij}$ from which we deduce the standardised weight matrix: $S = \frac{n}{Q} S^∗$.
#### **1.3.3 Importance of the choice of weight matrix**
When trying to test the importance of economic or social relationships between certain variables, the geographical location of the observations is a key parameter.
First of all, observations in the same geographical zone are subject to the same external parameters (climate, pollution, etc.) Secondly, neighbouring observations mutually influence one another.
Spatial models take these various interactions into account ; these models use neighbourhood specification via weight matrix $W$.
**The "weight list" object in R**
The function `nb2listw` of package *spdep* makes it possible to convert a "list of neighbours" object into a "weight list" object.
It is important to note that the "weight list" object, which corresponds to the weight matrix described above, is not a matrix n×n as represented in theory.
It is a list containing the standardisation style and then for each observation: its attribute, the list of observation numbers of its neighbours, the list of the attributes of its neighbours and the list of the weights of its neighbours.
Reference is often made to sparse matrices.
When a zone has no neighbours, the option `zero.policy=TRUE` makes it possible to generate a list of weights which takes value 'zero' for observations without neighbours (if the option is `FALSE`, an error message is generated).
**Application in R**
In the dataset "*div35*" :
```{r}
#Spatial series analysis : dataset div35
#Weight matrix
#For a REGULAR configuration of points
# 2 - Define the spatial weight matrix
#Note that spatial weight matrix is always associated to a neighborhood matrix (if two points are not considered as connected,
#their spatial weight is of "0")
# A - In the case in which points are connected to their 4th nearest neighbours (the "rook connection") ----
#compare spatial weighted matrix with rook connection (that is with the 4th nearest neighbors)
#but different weight formulation (binary B vs standardized W)
pond4.bin<- nb2listw(knn2nb(div.knear4), style="B",zero.policy=TRUE)
# "B" for binary = all connections have the same weight
pond4.stan<- nb2listw(knn2nb(div.knear4), style="W", zero.policy=TRUE)
# "W" for standardized (line standardisation) = 1/nb of connections
pond4.bin
pond4.stan
pond4.bin$weights
#List of created weights associated with the neighbor relationship
#when points are neighbors = 1, otherwise = 0. Translated into list, numerical value put into calculation
#Here : constant criterion for all points
pond4.stan$weights
# B - In the case in which points are connected to their 8th nearest neighbors (the "queen connection") ----
#compare spatial weighted matrix with queen connection
#but different weight formulation (binary B vs standardized [line standardisation] W)
pond8.bin<- nb2listw(knn2nb(div.knear8), style="B",zero.policy=TRUE)
pond8.stan<- nb2listw(knn2nb(div.knear8), style="W", zero.policy=TRUE)
pond8.bin$weights
pond8.stan$weights
# C - In the case in which points are connected to their 4th nearest neighbors but the weight is inversely proportional to the distance (weight = 1/distance)
#To get inverse distance, we need to calculate the distances between all of the neighbours.
#for this we will use nbdists, which gives us the distances in a similar structure to our input neighbors list. To use this function we need to input the neighbours list and the coordinates.
coords <- cbind(div35$x,div35$y)
distances <- nbdists(knn2nb(div.knear4),coords)
distances[1]
distances[34]
#Calculating the inverse distances (the 1/distance gives very small values so we make an adjustement)
invd1 <- lapply(distances, function(x) (1/x))
length(invd1)
invd1[1]
invd.weights <- nb2listw(knn2nb(div.knear4),glist = invd1,style = "B")
summary(invd.weights)
invd.weights$weights[1]
```
Let's see what it looks like in the second dataset, "*mares*" :
```{r}
#Spatial series analysis : dataset mares
#Weight matrix
#For a IRREGULAR configuration of points
# 2 - Define the spatial weight matrix =========================
# Example : BINARY criterion
pondmaretri.bin <- nb2listw(mare.tri, style="B",zero.policy=TRUE)
pondmaregab.bin <- nb2listw(graph2nb(mare.gab), style="B", zero.policy=TRUE)
pondmarednear.bin <- nb2listw(mare.dnear, style="B",zero.policy=TRUE)