Estad´ıstica Espacial 7 30' W MAR CARIBE Boca de la Barra

C. Clarín

R. Sevilla C. Grande

R.

Ar ac at ac a

10 45' N

R. Fundacion

NOTAS DE CLASE Ram´on Giraldo Henao Profesor Asociado Departamento de Estad´ıstica Universidad Nacional de Colombia

Bogot´a, Diciembre de 2011

Prefacio La necesidad de acudir a herramientas estad´ısticas para el an´alisis de datos en todas las ´areas del conocimiento, ha hecho que aparezcan con el correr de los a˜ nos nuevas metodolog´ıas que, no obstante se centran en fundamentos probabil´ısticos comunes, son espec´ıficas para cada una de las diversas disciplinas del saber. Algunos ejemplos son econometr´ıa, psicometr´ıa o bioestad´ıstica. Desde hace unas dos d´ecadas se ha dedicado mucho inter´es a la modelaci´on de datos espaciales, entendidos como aquellos en los que adem´as de los valores de las variables de inter´es se registran tambi´en lo sitios donde estos ocurren. El an´alisis de datos espaciales se ha hecho frecuente en muchas ´areas del conocimiento, incluyendo ente otras geograf´ıa, geolog´ıa, ciencias ambientales, econom´ıa epidemiolog´ıa o medicina. Las t´ecnicas estad´ısticas cl´asicas suponen que al estudiar un fen´omeno se realizan observaciones bajo circunstancias id´enticas y adem´as que las observaciones son independientes y por ello no son convenientes para analizar fen´omenos que var´ıan en tiempo o espacio. Cuando se tienen datos espaciales, intuitivamente se tiene la noci´on de que observaciones cercanas est´an correlacionadas y por ello es necesario acudir a herramientas de an´alisis que contemplen dicha estructura. El presente documento ha sido preparado como gu´ıa de clase para estudiantes del curso de estad´ıstica espacial del pregrado en estad´ıstica de la Universidad Nacional de Colombia. Este cubre los elementos esenciales de cada una de las t´ecnicas de an´alisis de datos espaciales, es decir, geoestad´ıstica, datos de ´areas (regionales) y patrones puntuales. Estas notas tambi´en tienen el prop´osito servir de consulta a ge´ologos, bi´ologos, agr´onomos, ingenieros, meteor´ologos, ec´ologos, epidemi´ologos y todos aquellos profesionales que se encargan del estudio de informaci´on georreferenciada. El documento tiene un enfoque te´orico-pr´actico. i

ii Para el seguimiento completo de la teor´ıa descrita se requiere tener conocimientos b´asicos de ´algebra lineal, estad´ıstica matem´atica (univariable y multivariable) y de modelos lineales. Un resumen no exhaustivo de estos temas es hecho al final en el ap´endice. Aquellos lectores que est´en poco familiarizados con los m´etodos estad´ısticos, podr´an obviar las secciones en las que se hacen desarrollos te´oricos y centrar su atenci´on en la filosof´ıa de los m´etodos presentados y en las aplicaciones mostradas en cada uno de los cap´ıtulos del documento. No obstante en el escrito se cubren diversos temas de la estad´ıstica espacial y se hacen aplicaciones de m´etodos recientes, es necesario acudir a la lectura de art´ıculos cient´ıficos y textos avanzados para lograr un buen dominio de estas metodolog´ıas. En el texto no se hace una redacci´on en castellano riguroso, en el sentido de que en muchas ocasiones se prefiere evitar la traducci´on de palabras t´ecnicas. Por ello en muchas partes del documento se encuentran expresiones como kriging, test, nugget o kernel, en vez de traducciones como krigeage, prueba, pepita o n´ ucleo, respectivamente..

´Indice general 1. Datos Espaciales y An´ alisis Exploratorio

1

1.1. Conceptos b´asicos de probabilidad y procesos estoc´asticos . . . . . . . . . .

1

1.2. Datos espaciales y ´areas de la estad´ıstica espacial . . . . . . . . . . . . . .

3

1.3. Medidas de dependencia espacial . . . . . . . . . . . . . . . . . . . . . . .

8

1.3.1. Test de Mantel . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

9

1.3.2. Test de Moran . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

10

1.3.3. Variograma . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

12

1.4. Efectos de la correlaci´on en inferencia estad´ıstica . . . . . . . . . . . . . . .

14

1.4.1. Efecto en la estimaci´on . . . . . . . . . . . . . . . . . . . . . . . . .

14

1.4.2. Efecto en la predicci´on . . . . . . . . . . . . . . . . . . . . . . . . .

15

1.5. Implementaci´on en R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

18

2. Patrones Puntuales

22

2.1. Tipos de Patrones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.2. Estimaci´on de la intensidad . . . . . . . . . . . . . . . . . . . . . . . . . .

23

2.2.1. Definici´on de Intensidad . . . . . . . . . . . . . . . . . . . . . . . .

25

2.2.2. Estimaci´on Kernel . . . . . . . . . . . . . . . . . . . . . . . . . . .

27

2.2.3. Estimaci´on Kernel de la intensidad λ(s) . . . . . . . . . . . . . . .

29

2.3. M´etodos basados en cuadrantes . . . . . . . . . . . . . . . . . . . . . . . .

30

2.3.1. Test de bondad de ajuste χ2 . . . . . . . . . . . . . . . . . . . . . .

31

2.4. M´etodos basados en distancias . . . . . . . . . . . . . . . . . . . . . . . . .

32

2.4.1. Funci´on G(h) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

32

iii

´INDICE GENERAL

iv

2.4.2. Funci´on F(h) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

2.4.3. Funci´on K de Ripley . . . . . . . . . . . . . . . . . . . . . . . . . .

36

2.4.4. Funci´on L . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.5. Modelos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

2.5.1. Proceso Poisson no Homog´eneo . . . . . . . . . . . . . . . . . . . .

38

2.5.2. Proceso Poisson Homog´eneo . . . . . . . . . . . . . . . . . . . . . .

41

2.5.3. Proceso de Cox . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

42

2.5.4. Proceso de Poisson Cluster . . . . . . . . . . . . . . . . . . . . . . .

43

2.5.5. Implementaci´on en R . . . . . . . . . . . . . . . . . . . . . . . . . .

44

3. Datos de Areas

58

3.1. Visualizaci´on de datos de ´areas . . . . . . . . . . . . . . . . . . . . . . . .

58

3.1.1. Mapas de datos originales y tasas . . . . . . . . . . . . . . . . . . .

58

3.1.2. Mapas de riesgos relativos . . . . . . . . . . . . . . . . . . . . . . .

59

3.1.3. Mapas Probabil´ısticos y Suavizados con Estimaci´on Bayesiana . . .

62

3.2. Exploraci´on de datos de ´areas . . . . . . . . . . . . . . . . . . . . . . . . .

64

3.2.1. Medidas de proximidad espacial . . . . . . . . . . . . . . . . . . . .

64

3.2.2. Medias m´oviles . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

3.2.3. Estimaci´on kernel . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

3.3. Correlaci´on espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

64

3.3.1. I de Moran y C de Geary . . . . . . . . . . . . . . . . . . . . . . . .

67

3.3.2. Indice Bayesiano Emp`ırico (IBE) . . . . . . . . . . . . . . . . . . .

68

3.3.3. Funci´on de semivariograma

. . . . . . . . . . . . . . . . . . . . . .

69

3.4. Modelos de regresi´on espacial . . . . . . . . . . . . . . . . . . . . . . . . .

69

3.4.1. Modelo de regresi´on polinomial local y ponderado geogr´aficamente .

69

3.4.2. Modelos lineales con errores correlacionados . . . . . . . . . . . . .

72

3.4.3. Modelos autorregresivos espaciales

. . . . . . . . . . . . . . . . . .

74

3.4.4. Modelo autorregresivo simult´aneo (SAR) . . . . . . . . . . . . . . .

74

3.4.5. Modelo autorregresivo condicional (CAR) . . . . . . . . . . . . . .

77

´INDICE GENERAL

v

4. Geoestad´ıstica

1

4.1. Definici´on de Geoestad´ıstica . . . . . . . . . . . . . . . . . . . . . . . . . .

1

4.2. Variable Regionalizada . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2

4.2.1. Momentos de una Variable Regionalizada . . . . . . . . . . . . . . .

2

4.2.2. Estacionariedad . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

3

4.3. Correlaci´on espacial muestral y ajuste de modelos . . . . . . . . . . . . . .

6

4.3.1. Funciones de correlaci´on espacial . . . . . . . . . . . . . . . . . . .

6

4.3.2. Modelos te´oricos de semivarianza . . . . . . . . . . . . . . . . . . .

8

4.4. Estimaci´on por m´ınimos cuadrados . . . . . . . . . . . . . . . . . . . . . .

14

4.4.1. Estimaci´on por M´ınimos cuadrados . . . . . . . . . . . . . . . . . .

14

4.4.2. M´ınimos cuadrados ponderados . . . . . . . . . . . . . . . . . . . .

15

4.4.3. Cuasi-verosimilitud . . . . . . . . . . . . . . . . . . . . . . . . . . .

16

4.5. Estimaci´on asumiendo normalidad . . . . . . . . . . . . . . . . . . . . . . .

18

4.5.1. Estimaci´on por M´axima Verosimilitud . . . . . . . . . . . . . . . .

18

4.6. Predicci´on espacial . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21

4.6.1. Predicci´on espacial optima . . . . . . . . . . . . . . . . . . . . . . .

21

4.6.2. Definici´on de Kriging . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.6.3. Kriging Ordinario . . . . . . . . . . . . . . . . . . . . . . . . . . . .

22

4.6.4. Otros M´etodos Kriging . . . . . . . . . . . . . . . . . . . . . . . . .

30

4.6.5. Kriging en Bloques. . . . . . . . . . . . . . . . . . . . . . . . . . . .

33

4.6.6. Kriging Universal . . . . . . . . . . . . . . . . . . . . . . . . . . . .

34

4.6.7. Kriging Residual . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.6.8. Kriging Indicador . . . . . . . . . . . . . . . . . . . . . . . . . . . .

37

4.6.9. Kriging Log-Normal y Multi-Gaussiano . . . . . . . . . . . . . . . .

39

4.6.10. Cokriging Ordinario . . . . . . . . . . . . . . . . . . . . . . . . . .

40

´Indice de figuras 1.1. Distribuci´on espacial de clorofila en la Ci´enaga Grande de Santa Marta (costa norte de Colombia). Datos medidos en un jornada de muestreo realizada en marzo de 1997.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

1.2. Mapa cloropl´etico de la tasa de delitos en Colombia en el a˜ no 2003. . . . .

7

1.3. Ubicaci´on de deslizamientos en el corredor Ca˜ no Lim´on-Cove˜ nas en 2008 (panel izquierdo) y ubicaci´on de sismos de baja magnitud en Colombia en el periodo Julio a Diciembre de 2008 (panel derecho). . . . . . . . . . . . .

9

1.4. Comportamiento t´ıpico de un semivariograma acotado con una representaci´on de los par´ametros b´asicos. SEMEXP corresponde al semivariograma experimental y MODELO al ajuste de un modelo te´orico.

. . . . . . . . . . . .

13

2.1. Ejemplos simulados de patrones puntuales con distribuci´on aleatoria (izquierda), regular (derecha) y agregada (centro).

. . . . . . . . . . . . . . . . .

24

2.2. a.)Ejemplos de patrones que incluyen covariables y b.) efecto de borde (derecha). En la figura de la izquierda la escala de valores corresponde a la elevaci´on en metros

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3. Marcas asociadas a los datos de sismos en Colombia

25

. . . . . . . . . . . .

26

2.4. Estimaci´on de la intensidad de sismos en Colombia (Julio-Diciembre 2008).

30

2.5. Gr´aficas del la funci´on G para cada tipo de proceso . . . . . . . . . . . . .

34

2.6. Gr´aficas del la funci´on F para cada tipo de proceso . . . . . . . . . . . . .

35

2.7. Gr´aficas del la funci´on K para cada tipo de proceso . . . . . . . . . . . . .

38

2.8. Gr´aficas del la funci´on L para cada tipo de proceso . . . . . . . . . . . . .

39

vi

´INDICE DE FIGURAS

vii

3.1. Mapas tem´aticos de mortalidad infantil en Colombia en 2003. . . . . . . .

59

3.2. Mapas tem´aticos de tasas de mortalidad infantil en Colombia en 2003.

60

. .

3.3. Distribuci´on de las tasa de violencia familiar (arriba) y de hurto calificado (abajo) seg´ un unidades seccionales de la Fiscal´ıa General de la Naci´on, en el a˜ no 2003. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

61

3.4. Dispersograma entre tasas de violencia familiar de 2003 y tasas suavizadas por el m´etodo bayesiano emp´ırico. ). . . . . . . . . . . . . . . . . . . . . .

63

3.5. Las ´areas sombreadas corresponden a los departamentos de la costa Caribe de Colombia (sin incluir San Andr´es y Providencia). . . . . . . . . . . . . .

65

4.1. Representaci´on de una superficie interpolada para una variable regionalizada estacionaria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

4

4.2. Representaci´on de una superficie interpolada para una variable regionalizada no estacionaria. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

6

4.3. Comportamiento t´ıpico de un semivariograma acotado con una representaci´on de los par´ ametros b´ asicos. SEMEXP corresponde al semivariograma experimental y MODELO al ajuste de un modelo te´ orico. . . . . . . . . . . . . . . . . . . . .

10

4.4. Comparaci´on de los modelos exponencial, esf´erico y Gaussiano. La l´ınea punteada vertical representa el rango en el caso del modelo esf´erico y el rango efectivo en el de los modelos exponencial y gaussiano. Este tiene un valor de 210, respecto a una escala simulada entre 0 y 300. El valor de la meseta es 30 y el de la pepita 0. El 95 % de la meseta es igual a 28.5. . . . . . . . . . . . . . . . . . . . . . .

11

4.5. CComportamiento t´ıpico de los modelos de semivarianza mon´omicos. . . . . . .

13

4.6. Modelo de semivarianza te´orico para variables sin correlaci´on espacial. . . . . .

13

4.7. Representaci´on del procedimiento kriging residual. La superficie interpolada (arriba) es igual a la suma de la tendencia lineal (centro) m´ as la predicci´ on de los errores (abajo). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

38

4.8. Representaci´on de la transformaci´on a scores normales. . . . . . . . . . . . . .

40

Cap´ıtulo 1 Datos Espaciales y An´ alisis Exploratorio En este cap´ıtulo se presentan conceptos b´asicos de probabilidad que permiten posteriormente enmarcar las ´areas de la estad´ıstica espacial dentro del contexto de los procesos estoc´asticos. Se definen algunas medidas de autocorrelaci´on espacial y se dan dos ejemplos de como la dependencia espacial afecta la inferencia estad´ıstica.

1.1.

Conceptos b´ asicos de probabilidad y procesos estoc´ asticos

Definition 1.1. Sea Ω 6= ∅. Un sistema F de subconjuntos de Ω se llama σ-´algebra si satisface las siguientes condiciones 1. Ω ∈ F 2. Si A ∈ F ⇒ Ac ∈ F 3. Si A1 , · · · , An ∈ F ⇒

Sn

i=1

Ai ∈ F .

Definition 1.2. Sea Ω 6= ∅, F una σ-´ algebra de subconjuntos de Ω. La pareja (Ω, F ) se llama espacio medible.

1

´ ´ 21.1. CONCEPTOS BASICOS DE PROBABILIDAD Y PROCESOS ESTOCASTICOS Definition 1.3. Sea (Ω, F ) un espacio medible. P : F → [0, 1] se llama medida de prob-

abilidad si satisface

1. P (A) ≥ 0, ∀A ∈ F 2. P (Ω) = 1 3. Si A1 , · · · , An ∈ F con Ai ∩ Aj = ∅, ∀i 6= j ⇒ P ( Propiedades de P

Sn

i=1

Ai ) =

Pn

i=1

P (Ai ).

1. P (∅) = 0 2. P (A ∪ B) = P (A) + P (B) − P (A ∩ B) 3. Si A ⊆ B ⇒ P (A) ≤ P (B). 4. P (A) = 1 − P (Ac ). Definition 1.4. La tripla (Ω, F , P ), donde Ω 6= ∅, F σ-´ algebra sobre Ω y P es una medida

de probabilidad sobre (Ω, F ), se denomina espacio de probabilidad.

Definition 1.5. Sea (Ω, F , P ), un espacio de probabilidad. X : Ω → R se llama variable aleatoria.

Definition 1.6 (Proceso Estoc´ astico). Es una familia de variables aleatorias {Z(s) :

s ∈ D ⊂ RP } definida sobre un espacio de probabilidad (Ω, F , P ). El conjunto D de´ındices

del procesos se denomina espacio de par´ ametros. Los valores que toma Z(s) se llaman estados y el conjunto de todos los posibles valores de Z(s) se llama espacio de estados. Los procesos estoc´asticos son clasificados de acuerdo con el espacio de par´ametros (discreto y continuo) y el espacio de estados (discreto y continuo). Algunos ejemplos de procesos estoc´asticos no espaciales son los siguientes 1. Espacio de par´ametros discreto y espacio de estados discreto Z(n): Preferencia del consumidor en el n-´esimo mes, con n ∈ N.

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

3

Z(n): N´ umero de individuos de la n-´esima generaci´on de una poblaci´on, con n ∈ N. 2. Espacio de par´ametros continuo y espacio de estados discreto Z(t): N´ umero de part´ıculas de una sustancia acuosa de volumen t, con t ∈ T ⊂

R.

Z(t): N´ umero de individuos que esperan el bus por periodo de tiempo t, con t ∈ T ⊂ R. 3. Espacio de par´ametros discreto y espacio de estados continuo Z(n): Tiempo de espera hasta que el n-´esimo estudiante arribe a la parada de bus, con n ∈ N. Z(n): Utilidad en pesos de un jugador despu´es del n-´esimo lanzamiento de una moneda, con n ∈ N. 4. Espacio de par´ametros continuo y espacio de estados continuo Z(t): Contenido de un embalse sobre un periodo de tiempo t, con t ∈ T ⊂ R. Z(t): Temperatura en el instante t, con t∈ T ⊂ R.

1.2.

Datos espaciales y ´ areas de la estad´ıstica espacial

Estad´ıstica espacial es la reuni´on de un conjunto de metodolog´ıas apropiadas para el an´alisis de datos que corresponden a la medici´on de variables aleatorias en diversos sitios (puntos del espacio o agregaciones espaciales) de una regi´on. De manera m´as formal se puede decir que la estad´ıstica espacial trata con el an´alisis de realizaciones de un proceso estoc´astico {Z(s) : s ∈ D ⊂ RP }, en el que s es la ubicaci´on en el espacio Euclidiano P-dimensional y Z(s) es una variable aleatoria en la ubicaci´on s. Observaciones

4

´ 1.2. DATOS ESPACIALES Y AREAS DE LA ESTAD´ISTICA ESPACIAL El proceso estoc´astico {Z(s) : s ∈ D ⊂ RP }, en el que s es sitio del espacio, tambi´en

se denomina proceso aleatorio o campo aleatorio.

El proceso estoc´astico {Z(s) : s ∈ D ⊂ RP }, en el que Z(s) es un vector aleatorio se

denomina en el contexto espacial proceso aleatorio multivariable o campo aleatorio multivariable

El proceso estoc´astico {Z(s) : s ∈ D ⊂ RP } en el que tanto el espacio de estados

como el espacio de par´ametros es continuo (es decir que las variables aleatorias Z(s) son continuas y D ⊂ RP es un conjunto continuo) se denomina variable regionalizada.

Este t´ermino es particularmente usado en aplicaciones de la estad´ıstica espacial en ingenier´ıa y geolog´ıa. Cuando se tiene una observaci´on del proceso estoc´astico {Z(s) : s ∈ D ⊂ RP } se

dispone de una muestra de tama˜ no {Z(s) = (Z(s1 ), Z(s2 ), · · · , Z(sn ))} (con n el

n´ umero de sitios donde se hace la medici´on de la variable aleatoria Z(s)) y no de una muestra de tama˜ no n de una variable aleatoria. Por ello puede ser carente de sentido pr´actico hacer inferencia estad´ıstica cl´asica (intervalos de confianza, pruebas de normalidad, etc) con los datos obtenidos. Desconocer esto hace que se cometan errores intentando validar los supuestos necesarios para la aplicaci´on de m´etodos estad´ısticos espaciales. En general en estad´ıstica espacial, como en el caso cl´asico, es deseable tener normalidad para hacer inferencia. Sin embargo lo que se asume en este contexto es que la muestra corresponde a la observaci´on de vector aleatorio con distribuci´on normal multivaluada y no que se tiene una muestra n-variada de una variable aleatoria con distribuci´on normal. Usar una prueba de normalidad univariada (por ejemplo la de Shapiro-Wilk) para comprobar si los datos siguen una distribuci´on normal es ciertamente equivocado en el contexto espacial, puesto que adem´as de desconocer que no se tiene una muestra iid (puesto que hay dependencia espacial), lo que en realidad habr´ıa que probar es normalidad multivariada. La estad´ıstica espacial se subdivide en tres grandes ´areas. La pertinencia de cada una de ellas est´a asociada a las caracter´ısticas del conjunto D de ´ındices del proceso estoc´astico

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

5

de inter´es. A continuaci´on se mencionan dichas ´areas y se describen las propiedades de D en cada una de ´estas. Geoestad´ıstica: Estudia datos de procesos estoc´asticos en los que el espacio de par´ametros D ⊂ RP es continuo. Algunos ejemplos de datos espaciales que son tratados con

m´etodos geoestad´ısticos son

{Z(s) : s ∈ D ⊂ RP }, donde Z(s) mide el contenido de nitr´ogeno en sitios de un parcela experimental. En este caso los sitios pertenecen a D ⊂ R2 .

{Z(s) : s ∈ D ⊂ RP }, donde Z(s) corresponde a la precipitaci´on en sitios de Colombia.

En los dos ejemplos anteriores hay infinitos sitios donde medir y por ello el conjunto de par´ametros es continuo. Sin embargo en la practica es potestad del investigador seleccionar en que sitios de la regi´on de inter´es hace la medici´on de las variables, es decir, el investigador puede hacer selecci´on de puntos del espacio a conveniencia o puede seleccionar los sitios bajo alg´ un esquema de muestreo probabil´ıstico. En este sentido se dice que el conjunto D ⊂ RP es fijo. Un ejemplo de un conjunto de datos analizado con m´etodos geoestad´ısticos es presentado en la Figura 1.1. Es importante resaltar que en geoestad´ısti-

ca el prop´osito esencial es la interpolaci´on y si no hay continuidad espacial pueden hacerse predicciones carentes de sentido. Datos de ´ areas o regionales: En este caso el proceso estoc´astico tiene espacio de par´ametros D ⊂ RP discreto y la selecci´on de los sitios de medici´on depende del investigador (D fijo). Las ubicaciones de muestreo pueden estar regular o irregularmente espaciadas. Dos ejemplos de datos regionales son {Z(s) : s ∈ D ⊂ RP }, donde Z(s) es la variable aleatoria correspondiente a la tasa de mortalidad y los sitios son departamentos de Colombia, es decir D es el conjunto discreto formado por los departamentos del pa´ıs.

6

´ 1.2. DATOS ESPACIALES Y AREAS DE LA ESTAD´ISTICA ESPACIAL

7 30' W MAR CARIBE Boca de la Barra

C. Clarín

R. Sevilla C. Grande

R.

Ar ac at ac a

10 45' N

R. Fundacion

Figura 1.1: Distribuci´on espacial de clorofila en la Ci´enaga Grande de Santa Marta (costa norte de Colombia). Datos medidos en un jornada de muestreo realizada en marzo de 1997.

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

7

Figura 1.2: Mapa cloropl´etico de la tasa de delitos en Colombia en el a˜ no 2003. {Z(s) : s ∈ D ⊂ RP }, donde Z(s) corresponde a la producci´on cafetera (en kilogramos) y D es el conjunto de todas las fincas productoras de caf´e del pa´ıs.

Nuevamente el investigador puede decidir donde (en que departamentos o en que fincas en los ejemplos) hace la medici´on de las variables de inter´es, es decir en datos de ´areas tambi´en el conjunto D ⊂ RP es fijo. En la Figura 1.2 se presenta un ejemplo de un conjunto

de datos que corresponde a la observaci´on de un proceso aleatorio de datos regionales. Un ejemplo de datos de ´area con sitios regularmente espaciados es el de colores de pixeles en interpretaci´on de im´agenes de sat´elite. En ese caso el conjunto de ubicaciones de inter´es es discreto y estas corresponden a agregaciones espaciales m´as que a un conjunto de puntos del espacio. Es obvio que la interpolaci´on espacial puede ser carente de sentido con este tipo de datos. Sus principales aplicaciones se encuentran en el campo epidemiol´ogico. Patrones Puntuales: La diferencia central del an´alisis de patrones puntuales con las

8

1.3. MEDIDAS DE DEPENDENCIA ESPACIAL

t´ecnicas geoestad´ısticas y el an´alisis de datos de ´areas radica en el hecho de que el conjunto D ⊂ RP es aleatorio, es decir que la decisi´on al respecto de donde se hace la medici´on no

depende del investigador. Dicho conjunto puede ser discreto o continuo, pero la ubicaci´on

de los sitios donde ocurre el fen´omeno a estudiar es dada por la naturaleza. En general el prop´osito de an´alisis en estos casos es el de determinar si la distribuci´on de los individuos dentro de la regi´on es aleatoria, agregada o uniforme. Algunos ejemplos de datos correspondientes a patrones puntuales son dados a continuaci´on Ubicaci´on de nidos de p´ajaros en una regi´on dada. Localizaci´on de imperfectos en una placa met´alica Sitios de terremotos en Colombia Municipios de Colombia con mayor´ıas negras En los tres primeros ejemplos D ⊂ RP es continuo y en el u ´ ltimo es discreto. Cuando

en cada sitio se hace medici´on de alguna variable (por ejemplo del n´ umero de huevos en los nidos, de la forma del imperfecto en la placa o de la tasa de analfabetismo de los municipios de mayor´ıas negras) se dice que ese tiene un patr´ on espacial marcado. Dos ejemplos de datos correspondientes a patrones espaciales son dados en la Figura 1.3.

1.3.

Medidas de dependencia espacial

La dependencia espacial hace referencia a la estructura de correlaci´on de las variables aleatorias del proceso {Z(s) : s ∈ D ⊂ RP }. Cuando hay dependencia espacial los sitios

cercanos tienen valores m´as similares que los distantes. Por el contrario la ausencia de correlaci´on espacial se refleja en el hecho de que la distancia entre los sitios no tiene influencia en la relaci´on de sus valores. A continuaci´on se presentan algunos test y funciones que permiten establecer estad´ısticamente o de manera emp´ırica si existe dependencia (correlaci´on) espacial.

400000

1000000

1500000

9

0

500000

Norte

1000000 0

500000

Norte

1500000

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

600000

800000

1000000

1400000

1800000

400000

750000

1100000

Este

1450000

1800000

Este

Figura 1.3: Ubicaci´on de deslizamientos en el corredor Ca˜ no Lim´on-Cove˜ nas en 2008 (panel izquierdo) y ubicaci´on de sismos de baja magnitud en Colombia en el periodo Julio a Diciembre de 2008 (panel derecho).

1.3.1.

Test de Mantel

Permite comprobar estad´ısticamente si las observaciones provienen de un proceso estoc´astico en el que las variables son correlacionadas espacialmente. Hip´otesis H0 : Hay aleatoriedad espacial Ha : Hay correlaci´on espacial Estad´ıstica de prueba M=

n X n X

Wij Uij ,

i=1 i=1

donde W ij = ksi − sj k y Uij = (Z(si ) − Z(sj ))2 . La estad´ıstica de mantel est´a rela-

cionada con la pendiente del modelo de regresi´on simple Uij = βWij + eij a trav´es de P P β = M/ ni=1 ni=1 Wij2 , es decir que intuitivamente se tiene que a mayor M, mayor de-

pendencia espacial positiva. La significancia de la prueba puede establecerse por varios caminos. Puede emplearse un test de permutaciones en el que asumiendo aleatoriedad

10

1.3. MEDIDAS DE DEPENDENCIA ESPACIAL

se encuentran las n! posibles asignaciones de sitios a valores y con cada una de ellas se calcula M, obteni´endose por consiguiente su distribuci´on bajo H0 . Tambi´en en el caso de n grande puede usarse un test de Monte Carlo en el que solo se toman k de las asignaciones aleatorias de sitios a valores de la variable. En ambos casos (permutaciones, Monte Carlo) podr´ıa usarse una aproximaci´on a la normal estimando E(M) y V (M) a trav´es de P P ¯ = 1/n n Mi y s2 = 1/(n − 1) n (Mi − M) ¯ 2 . En el caso de asumir normalidad y M M i=1 i=1 aleatoriedad, es decir si Z(s1 ), · · · .Z(sn ) son iid, con Z(si) ∼ N(µ, σ 2 ), pueden obtenerse

expresiones para E(M) y V (M) y establecer el nivel de significancia bas´andose en un test

normal.

1.3.2.

Test de Moran

Este test es especialmente usado en datos de ´areas. Sean Z(s1 ), · · · , Z(sn ), las variables

medidas en las n ´areas. La noci´on de autocorrelaci´on espacial de estas variables est´a asociada con la idea de que valores observados en ´areas geogr´aficas adyacentes ser´an m´as similares que los esperados bajo el supuesto de independencia espacial. El ´ındice de autocorrelaci´on de Moran considerando la informaci´on de los vecinos m´as cercanos es definida como n I= P n P n

i=1 j=1

n P n P

i=1 j=1

Wij

¯ ¯ Wij (Z(si ) − Z)(Z(s j ) − Z) n P

¯ 2 (Z(si ) − Z)

(1.1)

i=1

Valores positivos (entre 0 y 1) indican autocorrelaci´on directa (similitud entre valores cercanos) y valores negativos (entre -1 y 0) indican autocorrelaci´on inversa (disimilitud entre las ´areas cercanas). Valores del coeficiente cercanos a cero apoyan la hip´otesis de aleatoriedad espacial. Para el c´alculo del ´ındice de Moran es necesario definir la proximidad entre las ´areas. Lo anterior se lleva a cabo por medio del c´alculo de una matriz de proximidad espacial. Dado un conjunto de n ´areas (A1 , · · · An ) se construye una matriz W (1) de orden (n × n)

donde cada uno de los elementos Wij representa una medida de proximidad entre Ai y Aj j. Dicha medida puede ser calculada con alguno de los siguientes criterios:

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

11

Wij = 1 si el centro de Ai se encuentra a una distancia determinada de Aj o Wij = 0 en caso contrario. Wij = 1 si Ai comparte frontera con Aj y en caso contrario Wij = 0. Wij = Iij /Ii , donde Iij es la longitud de la frontera entre Ai y Aj y Ii es el per´ımetro de Ai . Wij = dij , con dij la distancia entre los centros de las dos ´areas. En todos los casos anteriores Wii = 0. La idea de la matriz de proximidad espacial puede ser generalizada a vecinos de mayor orden (vecinos de vecinos) construy´endose as´ı las matrices W (2) , · · · , W (n) . Se acostumbra a normalizar las filas de la matriz, es decir que la suma por fila de los Wij sea igual a uno.

Una vez obtenido el valor del coeficiente es necesario evaluar su significancia estad´ıstica. En otras palabras se requiere probar la hip´otesis de aleatoriedad espacial con base en el valor observado. Para llevar a cabo esto es necesario establecer la correspondiente distribuci´on de probabilidad de la estad´ıstica de prueba I. Bajo normalidad, es decir asumiendo que Z(s1 ), · · · , Z(sn ) son iid con Z( si ) ∼ N(µ, σ 2 ), la estad´ıstica

I − E(I) Z= p V(I) sigue una distribuci´on normal est´andar, en la que el valor esperado y la varianza est´an dados por E(I) = − donde S0 = Wi0 =

1 n2 S1 − n2 S2 + 3S02 1 , V(I) = − , 2 2 (n + 1) (n − 1)S0 (n − 1)2

n X

i6=j n X j=1

Wij , S1 =

n X

2

(Wij + Wji ) , S2 =

i6=j n X

Wij , W0i =

n X

(Wi0 + W0i )2 ,

i=1

Wji .

j=1

Otra posibilidad para establecer la significancia estad´ıstica, con menos supuestos, es llevando a cabo un test de permutaci´on o de Monte Carlo como los descritos para la estad´ıstica de Mantel.

12

1.3. MEDIDAS DE DEPENDENCIA ESPACIAL

1.3.3.

Variograma

El variograma, denotado por 2γ(h), se define como la varianza de la diferencia entre variables separadas por una distancia h = ksi − sj k. Asumiendo que E(Z(s)) = µ se tiene 2γ(h) = V(Z(s + h) − Z(s)) = E(Z(s + h) − Z(s))2 .

(1.2)

La mitad del variograma se llama semivariograma y caracteriza las propiedades de dependencia espacial de un fen´omeno espacial. Esta funci´on es usualmente empleada para tratar datos de un fen´omeno con continuidad espacial (datos geoestad´ısticos). Usando el m´etodo de momentos se tiene que un estimador del semivariograma es n(h)

1 X (Z(s + h) − Z(s))2 , γ¯ (h) = n(h)

(1.3)

donde n(h) representa el n´ umero de parejas de sitios (si , sj ) que se encuentran separados por una distancia h. En la pr´actica, debido a irregularidad en el muestreo y por ende en las distancias entre los sitios, se toman intervalos de distancia {[0, h], (h, 2h], (2h, 3h], · · · } y el

semivariograma experimental corresponde a una distancia promedio entre parejas de sitios

dentro de cada intervalo y no a una distancia h espec´ıfica. Obviamente el n´ umero de parejas de puntos n dentro de los intervalos no es constante. Para interpretar el semivariograma experimental se parte del criterio de que a menor distancia entre los sitios mayor similitud o correlaci´on espacial entre las observaciones. Por ello en presencia de autocorrelaci´on se espera que para valores de h peque˜ nos el semivariograma experimental tenga magnitudes menores a las que este toma cuando las distancias se incrementan. Como se ver´a en el cap´ıtulo 4 la soluci´on del problema de predicci´on espacial requiere del conocimiento de la estructura de autocorrelaci´on para cualquier posible distancia entre sitios dentro del ´area de estudio. De la ecuaci´on (1.3) es claro que el semivariograma muestral es calculado s´olo para algunas distancias promedios particulares. Por ello se hace necesario el ajuste de modelos que generalicen la dependencia espacial para cualquier distancia (Figura 1.3). Existen diversos modelos te´oricos de semivarianza que pueden ajustarse al semivariograma muestral. En Cressie (1993) se presenta una discusi´on respecto

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

Semivarianza

2,0

13

2

Sill ( σ )

1,6 1,2

SEMEXP MODELO

0,8 0,4

Rango ( φ )

Nugget ( τ )

0,0

0

10000

20000

30000

Distancia

Figura 1.4: Comportamiento t´ıpico de un semivariograma acotado con una representaci´on de los par´ametros b´asicos. SEMEXP corresponde al semivariograma experimental y MODELO al ajuste de un modelo te´orico. a las caracter´ısticas y condiciones que ´estos deben cumplir. En general dichos modelos pueden dividirse en no acotados (lineal, logar´ıtmico, potencial) y acotados (esf´erico, exponencial, Gaussiano) (?). Los del segundo grupo garantizan que la covarianza de los incrementos es finita, por lo cual son ampliamente usados cuando hay evidencia de que presentan buen ajuste. La mayor´ıa de modelos empleados para ajustar el semivariograma muestral, tienen tres par´ametros en com´ un (Figura 1.4) que son descritos a continuaci´on: Nugget (τ ): Representa una discontinuidad puntual del semivariograma en el origen (Figura 1.3). Puede ser debido a errores de medici´on en la variable o a la escala de la misma. En algunas ocasiones puede ser indicativo de que parte de la estructura espacial se concentra a distancias inferiores a las observadas. Sill (σ 2 ): Es un estimador de la varianza de las variables del proceso. Tambi´en puede definirse como el limite del semivariograma cuando la distancia h tiende a infinito. Rango(φ). En t´erminos pr´acticos corresponde a la distancia a partir de la cual dos observaciones son independientes. El rango se interpreta como la zona de influencia. Existen algunos modelos de semivariograma en los que no existe una distancia finita para la cual dos observaciones sean independientes; por ello se llama rango efectivo a la distancia para la cual el semivariograma alcanza el 95 % de la meseta (sill ).

´ EN INFERENCIA ESTAD´ISTICA 1.4. EFECTOS DE LA CORRELACION

14

1.4.

Efectos de la correlaci´ on en inferencia estad´ıstica

Muchas m´etodos estad´ısticos est´an basados en el supuesto de que las variables aleatorias involucradas en la muestra son independientes. La violaci´on de dicho supuesto tiene consecuencias en todos los procesos inferenciales. En esta secci´on se ilustra como la correlaci´on entre las variables (por consiguiente la no independencia entre las mismas) afecta la estimaci´on y la predicci´on en el modelo de regresi´on simple (sin covariables).

1.4.1.

Efecto en la estimaci´ on

Sea Y1 , · · · , Yn una muestra aleatoria de Y ∼ N(µ, σ 2 ). El estimador de µ es Y¯ = Pn 1 2 i=1 Yi . El valor esperado y la varianza de este estimador son µ y σ /n, respectivamente. n

Ahora suponga que las variables Y1 , · · · , Yn son correlacionadas y que Cov(Yi , Yj ) = σ 2 ρ. P En este caso nuevamente el estimador de µ es Y¯ = n1 ni=1 Yi y su valor esperado es µ, sin

embargo la correlaci´on aumenta (en este caso) la varianza del estimador. Veamos n

V (Y¯ ) = V ( 1 = 2 n = = = =

1X Yi ) n i=1

n X n X i=1 j=1

Cov(Yi , Yj )

!

 1  2 2 2 2 2 2 (σ + σ ρ, · · · , +σ ρ), · · · , (σ + σ ρ, · · · , +σ ρ) n2
1 2 2 2 nσ + (n − 1)σ ρ, · · · , (n − 1)σ ρ n2
1 2 2 nσ + n(n − 1)σ ρ n2 σ2 (1 + (n − 1)ρ) . n

(1.4)

Si ρ > 0 en (1.4), V (Y¯ ) > σ 2 /n, es decir la varianza del estimador de µ cuando hay correlaci´on es mayor que la de este mismo cuando las variables son independientes.

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

1.4.2.

15

Efecto en la predicci´ on

Sean Y1 , · · · , Yn variables aleatorias tales que Yi ∼ N(µ, σ 2 ) y Cov(Yi, Yj ) = σ 2 ρ. Un

modelo lineal para representar este escenario es 

  Y1  .   .   Y=  . = Yn

  µ  ..   + .    µ

= µ1 + ǫ,

 ǫ1 ..  .   ǫn

(1.5)

donde 

  2 V (ǫ) = Σ = σ   

1 ρ ··· ρ



 1 ··· ρ   . .. . . ..  . .  .  ρ ρ ··· 1 ρ .. .

Suponga que se quiere predecir una nueva observaci´on Y0 . Definiendo el predictor por Y0∗

=

n X

λi Y i ,

(1.6)

i=1

los pesos λi se obtienen de tal forma que se minimice la esperanza de una funci´on de p´erdida. Bajo p´erdida cuadr´atica, el mejor predictor (lineal en este caso), ser´a el que minimiza la funci´on m´ın E(Y0∗ − Y0 )2 , sujeto a E(Y0∗ ) = E(Y0 ).

λ1 ,...,λn

De acuerdo con lo anterior, la funci´on objetivo es m´ınV(Y0∗ − Y0 ), sujeto a λ,m

n X i=1

λi = 1.

16

´ EN INFERENCIA ESTAD´ISTICA 1.4. EFECTOS DE LA CORRELACION

Se tiene que Cov(Y, Y0) = σ 2 ρ1. Desarrollando la varianza e incluyendo un multiplicador de Lagrange para la condici´on de insesgadez la funci´on a optimizar es m´ınV(Y0∗ ) λ,m

+ V(Y0 ) −

2Cov(Y0∗ , Y0)

− 2m(

n X i=1

λi − 1)

m´ınV(λT Y) + σ 2 − 2Cov(λT Y, Y0) − 2m(λT 1 − 1) λ,m  σ2ρ  .. m´ınλT Σλ + σ 2 − 2λT c − 2m(λT 1 − 1), c =  .  λ,m σ2ρ



  = σ 2 ρ1. 

Tomando derivadas respecto a λ y m se obtiene el siguiente sistema Σλ − c − m1 = 0

λT 1 − 1 = 0.

(1.7)

Despejando λ en la primera ecuaci´on del sistema (1.7), se obtiene λ = Σ−1 (c + m1).

(1.8)

Reemplazando esta expresi´on en la segunda ecuaci´on del sistema (1.7) se encuentra (Σ−1 (c + m1))T 1 = 1 (Σ−1 c + Σ−1 m1)T 1 = 1 1T (Σ−1 c) + 1T (Σ−1 m1) = 1 1T (Σ−1 m1) = 1 − 1T (Σ−1 c)
m = 1 − 1T (Σ−1 c) (1T Σ−1 1)−1 m=

1 − 1T (Σ−1 c) 1T Σ−1 1

Sustituyendo (1.9) en la ecuaci´on (1.8) se obtiene

(1.9)

´ CAP´ITULO 1. DATOS ESPACIALES Y ANALISIS EXPLORATORIO

  1 − 1T (Σ−1 c) λ=Σ c+ 1 1T Σ−1 1 T  1 − 1T (Σ−1 c) T λ = c+1 (Σ−1 )T T −1 1 Σ 1  T 1 − 1T (Σ−1 c) T λ = c+1 Σ−1 . 1T Σ−1 1

17

−1

(1.10)

De acuerdo con la soluci´on obtenida en (1.10), el predictor en (1.6) es definido por Y0∗

=

n X

λi Y i

i=1

= λT Y " # T 1 − 1T (Σ−1 c) = c+1 Σ−1 Y 1T Σ−1 1 Haciendo algunas manipulaiones de ´algebra se obtiene que Y0∗ = µ ˆ + cT Σ−1 (Y − 1ˆ µ) ,

(1.11)

donde µ ˆ es el estimador de m´ınimos cuadrados generalizados de µ en la ecuaci´on (1.5). La varianza del predictor en (1.11) est´a dada por σp2 = σ 2 − cT Σ1 c +

(1 − 1T Σ−1 c)2 . (1T Σ−1 1)

(1.12)

Observaci´ on Del modelo lineal general Y = Xβ+ǫ se tiene que el estimador de m´ınimos cuadrados
−1
generalizados del vector de par´ametros es β = XT Σ−1 X XT Σ−1 Y . Definiendo
−1 T −1
X = 1 y β = µ, se obtiene que µ ˆ = 1T Σ−1 1 1 Σ Y .

Ahora consid´erese el caso de predicci´on teniendo una muestra aleatoria. Sean Y1 , · · · , Yn

variables aleatorias independientes e id´enticamente distribuidas, con Yi ∼ N(µ, σ 2 ). Planteando el mismo predictor dado en (1.6) y reemplazando Σ = σ 2 I y c = 0 en las ecuaciones

´ EN R 1.5. IMPLEMENTACION

18

(1.8) y (1.9), se encuentra que m = (1T (σ 2 I)−1 1)−1 y que λ = (σ 2 I)−1 (1T (σ 2 I)−1 1)−1 1   1/σ 2 · · · 0 σ 2 /n  .  . ..  ..  . . = . .   .  . 0 · · · 1/σ 2 σ 2 /n

Al sustituir (1.13) en (1.6) se obtiene Y0∗

=

n X





 1/n   .   =  ..  .    1/n

(1.13)

λi Y i

i=1

n

1X Yi = Y¯ . = n i=1

(1.14)

Tomando Σ = σ 2 I y c = 0 en (1.12) se obtiene que σp2 = σ 2 (1 + 1/n), es decir la varianza de predicci´on del modelo bajo independencia.

1.5.

Implementaci´ on en R

######################################################### # Ejemplos del test de Mantel para aleatoriedad versus # asociaci´ on espacial #########################################################

library(ape)

# Para el test de Mantel

library(spdep) library(geoR) library (gstat) library(MASS) # Para simulaci´ on de normales # Dos ejemplos bajo aleatoriedad q1