Modelos de regresión para datos espaciales

FACULTAD DE MATEM´ ATICAS DEPARTAMENTO DE ESTAD´ ISTICA E INVESTIGACI´ ON OPERATIVA Trabajo Fin de Grado Modelos de Regresi´on para Datos Espaciales Jos´e ´ Angel Borrego S´anchez Dirigido por: Juan Manuel Mu˜noz Pichardo 2018 2 ´ Indice general Resumen 5 Abstract 7 1. Introducci´on 9 1.1. Modelado del Espacio Geogr´afico . . . . . . . . . . . . . . . . . . . 10 1.2. Tipos de Datos Espaciales . . . . . . . . . . . . . . . . . . . . . . . 11 1.3. Fiabilidad de los Datos Espaciales . . . . . . . . . . . . . . . . . . . 12 2. An´alisis Exploratorio Datos Espaciales 15 2.1. Visualizaci´on de Datos . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.2. Matriz de Datos Espaciales . . . . . . . . . . . . . . . . . . . . . . . 18 2.3. Autocorrelaci´on Espacial . . . . . . . . . . . . . . . . . . . . . . . . 19 2.4. Estad´ıstico General de Productos Cruzados . . . . . . . . . . . . . . 22 3. Medidas de Autocorrelaci´on Espacial 25 3.1. Medidasglobales ............................ 25 3.2. Medidaslocales............................. 31 4. Modelado de los datos de ´area 35 4.1. Tipos de modelos de regresi´on espacial . . . . . . . . . . . . . . . . 35 4.1.1. Modelos de retardo espacial . . . . . . . . . . . . . . . . . . 37 4.1.2. Modelos de error espacial . . . . . . . . . . . . . . . . . . . . 38 4.1.3. Modelos de orden superior . . . . . . . . . . . . . . . . . . . 38 4.1.4. Modelo espacial de Durbin . . . . . . . . . . . . . . . . . . . 39 4.2. Pruebas para Dependencia Espacial . . . . . . . . . . . . . . . . . . 40 4.3. Estimaci´on de modelos espaciales . . . . . . . . . . . . . . . . . . . 44 5. Implementaci´on de datos espaciales en R 49 5.1. Paquetesprevios ............................ 49 5.2. Paquete spdep .............................. 51 5.2.1. Vecindad de datos de ´area . . . . . . . . . . . . . . . . . . . 51 3 5.2.2. Pruebas para Dependencia Espacial . . . . . . . . . . . . . . 53 5.2.3. Modelos de Regresi´on Espacial en R . . . . . . . . . . . . . 55 5.3. Ilustraci´on de los Modelos en R . . . . . . . . . . . . . . . . . . . . 57 Bibliograf´ıa 71 4 Resumen En este trabajo nuestro objeto de estudio ser´an los datos espaciales, de los cuales nos centraremos con detalle en los datos de tipo ´area. Nuestro objetivo final ser´a proponer diferentes tipos de modelos para los datos de ´area y aplicarlos a un conjunto de datos real. Para ello necesitaremos modelar el espacio donde se ubican los datos y modificarlo para conseguir un modelado discreto en el cual podamos distinguir diferentes ´areas. A continuaci´on, procederemos a realizar un an´alisis exploratorio de los datos que consiste en crear una matriz de datos con las diferentes propiedades de los datos y su ubicaci´on en el ´area de estudio. A partir de esta matriz, y con las t´ecnicas que se expondr´an en el cap´ıtulo 2, podremos analizar si existe dependencia espacial en los datos, es decir, si los valores de la variable bajo estudio en un punto muestral est´an relacionados con los valores en puntos muestrales cercanos. Para ver este an´alisis mostraremos medidas de dependencia espacial y contraste de hip´otesis con algunos estad´ısticos como el ´ındice Ide Moran o el cde Geary (cap´ıtulo 3). Para finalizar el an´alisis te´orico de los datos, detallaremos diferentes modelos de regresi´on espacial y la estimaci´on de sus par´ametros (cap´ıtulo 4). Por ´ultimo, el cap´ıtulo final est´a reservado para analizar los datos de ´area con el programa estad´ıstico R. En ´el veremos los diferentes paquetes, librer´ıas y funciones que necesitaremos para realizar este an´alisis y expondremos un ejemplo pr´actico con datos reales para ver el uso de algunas funciones definidas en este conjunto de datos. 5 6 Abstract In this work our study objetive will be the spatial datas in which we will be centered with detail in the type area data. Our final aim will be to propose differents types of models for area datas and apply them to a real data set. For this will be need to model the space where datas are located acd modify it to get a discrete modelling in which we can distinguish differents areas. Bellow we will proceed to make an exploratory analysis about datas that consist on creating a datas matrix with the differents datas properties and their location in the study area. From this matrix on, and with techniques that we will see in chapter 2, we could analyze if there is a spatial dependence in the datas, that is, if the variable values under study in a sample point are related to the values in nearby sampling points. To see this analysis we will show measures of spatial dependence and contrast tests with some statistics like Moran’s Ior Geary’s c(chapter 3). To finish the datas theoric analysis, we will detail differents spatial regression models and the estimation of its parameters (chapter 4). Finally, the last chapter is reserved to analyze the area datas with the statistic programme R. In this programme we will see the different packages, libraries and functions that we are going to be needed to make this analysis and we will expose a practical example with real datas to see the use of some functions which are defined in this datas set. 7 8 Cap´ıtulo 1 Introducci´on En este cap´ıtulo comenzaremos introduciendo una definici´on de datos espaciales proporcionado por Fischer y Wang [5], en el cual nos apoyaremos bastante para la redacci´on de este trabajo. Dentro de este cap´ıtulo veremos la importancia de la ubicaci´on geogr´afica por lo que necesitaremos modelar nuestra regi´on de estudio y tratarla como un espacio continuo o discreto. Debido a este modelado del espacio podremos distinguir varios tipos de datos que veremos en la secci´on 1.2. Para finalizar este cap´ıtulo de introducci´on trataremos de alcanzar una fiabilidad y calidad de los datos y veremos algunos de los problemas m´as comunes que pueden surgir en el an´alisis a partir de unos datos defectuosos o de baja calidad. Seg´un M. Fisher y J. Wang, los datos consisten en n´umeros o s´ımbolos que, en cierto sentido, son neutrales y, en contraste con la informaci´on, casi libres de contexto. Los datos geogr´aficos brutos, como la temperatura en un momento y lugar espec´ıficos, son ejemplos de datos. En Longley et al. (2001, p. 64) podemos ver los datos espaciales como construidos a partir de elementos o hechos sobre el mundo geogr´afico. Un dato espacial vincula una ubicaci´on geogr´afica (espacio geogr´afico), a menudo un atributo de tiempo, y un atributo (o propiedad descriptiva) de la entidad entre s´ı. Estos atributos pueden ser de naturaleza ambiental (por ejemplo; la temperatura), econ´omicos (capital de una regi´on concreta), otros pueden identificar una ubicaci´on (direcciones postales). En particular, en el an´alisis de datos espaciales el tiempo es opcional, sin embargo la ubicaci´on geogr´afica es esencial ya que es la que distingue este an´alisis de otros con datos no espaciales. Para llevar a cabo el an´alisis de datos espaciales requerimos, al menos, la informaci´on de la ubicaci´on y de los atributos, independientemente de c´omo se midan los atributos. Este tipo de an´alisis requiere un marco espacial sobre el cual ubicaremos los fen´omenos espaciales que vamos a estudiar. 9 16 CAP´ ITULO 2. AN ´ ALISIS EXPLORATORIO DATOS ESPACIALES Divisiones naturales. Las clases se definen de acuerdo con algunas agrupaciones naturales de los valores de los datos. Los diferentes cortes pueden imponerse a partir de puntos de ruptura que se sabe que son relevantes (asignaci´on deductiva) o mediante el uso de herramientas del sistema GIS1que buscan saltos importantes en los valores de datos. Divisiones por cuantiles. Cada clase contiene un n´umero igual de observaciones. Las clasificaciones de cuartiles (cuatro categor´ıas) y quintiles (cinco categor´ıas) son las m´as com´unmente utilizadas en la pr´actica. Divisiones de intervalos iguales. Tienen validez cuando las observaciones se distribuyen de manera razonablemente uniforme en su rango. Pero si los datos son marcadamente sesgados, dar´an un gran n´umero de observaciones en unas pocas clases. Divisi´on seg´un desviaci´on est´andar. Se basan en intervalos distribuidos alrededor de la media en unidades de desviaci´on est´andar. El uso de mapa de coropletas tiene algunos problemas que deben ser tenidos en cuenta antes de su elaboraci´on: 1. Los mapas de coropletas traen la implicaci´on visual de la uniformidad dentro del ´area de los valores variables. Adem´as, el mapa de coropletas convencional permite que cualquier ´area f´ısicamente grande domine la pantalla, de una manera que puede ser bastante inapropiada para el tipo de datos que se mapean. Por ejemplo, al mapear datos socioecon´omicos, las ´areas rurales grandes y escasamente pobladas pueden dominar el mapa pero el inter´es real puede estar en ´areas f´ısicamente m´as peque˜nas, como las ´areas urbanas m´as densamente pobladas. 2. La variable de inter´es ha surgido de la agregaci´on de datos individuales a las ´areas de forma arbitraria y cualquier patr´on que observemos, tanto en los l´ımites de ´area elegidos como en la distribuci´on espacial, proviene de dicha arbitrariedad de los valores de las variables, no teniendo ning´un rigor ni formalidad a la hora de estudiar esos datos con la hip´otesis de seguir dicho patr´on encontrado. Este es el problema de unidad de ´area modificable (MAUP) que surge a menudo en el mapa de coropletas. 1Las siglas GIS provienen de sistema de informaci´on geogr´afica que es un conjunto de herramientas que integra y relaciona diversos componentes que permiten la organizaci´on, almacenamiento, an´alisis y modelizaci´on de grandes cantidades de datos procedentes del mundo real que est´an vinculados a una referencia espacial. 2.1. VISUALIZACI ´ ON DE DATOS 17 Figura 2.1: Tipos de intervalos de clases para un mapa de coropletas 3. Es importante darse cuenta de que los resultados estad´ısticos de cualquier an´alisis de patrones y relaciones depender´an inevitablemente de la configuraci´on de ´area particular que se utilice, como hemos visto en el p´arrafo anterior. En general, para conseguir un buen an´alisis y, por tanto unos resultados fiables de lo que estudiemos, los datos deben analizarse sobre la base de las unidades de ´area m´as peque˜nas para las cuales est´an disponibles 18 CAP´ ITULO 2. AN ´ ALISIS EXPLORATORIO DATOS ESPACIALES y debe evitarse la agregaci´on a ´areas m´as grandes arbitrarias, a menos que haya buenas razones para hacerlo. Tambi´en es importante verificar cualquier inferencia extra´ıda de los datos mediante el uso de diferentes configuraciones de ´area de los mismos datos. 2.2. Matriz de Datos Espaciales Todas las t´ecnicas anal´ıticas que utilizamos en este trabajo usan una matriz de datos que recoge los datos espaciales necesarios para la realizaci´on del an´alisis. Sean Z1, Z2, . . . , ZKvariables aleatorias y sea Sla ubicaci´on del punto o ´area. La matriz de datos espaciales viene generada por:      z1(1) z2(1) ··· zK(1) s(1) z1(2) z2(2) ··· zK(2) s(2) . . .. . .. . .. . . z1(n)z2(n)··· zK(n)s(n)      la cual puede ser expresada de la siguente forma z1(i), z2(i),··· zK(i)s(i)i=1...n (2.1) donde zkdenota una realizaci´on muestral (valor de datos real) de la variable Zk(k= 1, . . . , K) mientras que el s´ımbolo identro de los par´entesis hace referencia al caso u observaci´on particular (punto, ´area, pol´ıgono o zona espacial). Para cada caso, i= 1, . . . , n,s(i) representa la ubicaci´on del objeto espacial. La referencia implicar´a dos coordenadas geogr´aficas. Generalmente nos referiremos a un espacio bidimensional, por tanto s(i) = (s1(i), s2(i))t. En el caso de datos referentes a objetos de puntos en un espacio bidimensional, la ubicaci´on del punto i-´esimo puede estar dada por un par de coordenadas cartesianas (ortogonales). Los ejes del sistema de coordenadas generalmente se han construido para el conjunto de datos particular, pero se puede usar cualquier sistema de referencia. En el caso de datos referentes a objetos de ´area de forma irregular, una opci´on para representar la ubicaci´on es seleccionar un punto representativo de cada ´area denominado centroide y luego usar el mismo procedimiento que para un objeto de punto para identificar s(i). Hay situaciones en las que la informaci´on de georreferencia proporcionada por {s(i)}en la expresi´on (2.1) tiene que completarse con informaci´on de vecindad que define no solo qu´e pares de ´areas son adyacentes entre s´ı, sino que tambi´en puede cuantificar la cercan´ıa de esa adyacencia. Esta informaci´on es necesaria para la especificaci´on de muchos modelos 2.3. AUTOCORRELACI ´ ON ESPACIAL 19 estad´ısticos espaciales, como los modelos de regresi´on espacial. Dicha informaci´on est´a recogida en el concepto de autocorrelaci´on espacial que veremos en la siguiente secci´on. 2.3. Autocorrelaci´on Espacial Para comenzar con esta secci´on veamos algunas definiciones que han dado diferentes autores de la autocorrelaci´on espacial que quedan recogidas en el estudio del laboratorio de paleontolog´ıa y climatolog´ıa de la universidad de Ottawa [12]. Para Cliff y Ord (1973) es ´util ver si la distribuci´on de cierta cualidad o cantidad en los condados o estados de un pa´ıs hace que su presencia en los condados vecinos sea m´as o menos probable. En tal caso, exponen que el fen´omeno exhibe autocorrelaci´on espacial. Sokal y Oden (1978) argumentaron que el an´alisis de autocorrelaci´on espacial prueba si el valor observado de una variable nominal, ordinal o de intervalo en una localidad es independiente de los valores de esa misma variable en las localidades vecinas. Upton y Fingleton (1985) la definen como una propiedad que los datos poseen cuando muestran un patr´on de comportamiento. Estos autores exponen que la autocorrelaci´on espacial existe siempre que haya una variaci´on espacial sistem´atica en los valores a lo largo de un mapa, o patrones en los valores registrados en las localizaciones. Goodchild (1987) y Griffith (1991) destacaron que la autocorrelaci´on espacial trata simult´aneamente con la informaci´on de ubicaci´on y de atributos. Goodchild (1987) dice que en su sentido m´as general la autocorrelaci´on espacial se refiere al grado en que los objetos o actividades en alg´un lugar de la superficie son similares a otros objetos o actividades ubicados cerca y refleja la primera ley de geograf´ıa de Tobler [13] “todo est´a relacionado con todo lo dem´as, pero las cosas cercanas est´an m´as relacionadas que las cosas distantes”. As´ı, seg´un las diferentes aportaciones sobre el concepto de autocorrelaci´on espacial, se puede llegar a la conclusi´on de que si el valor de una o varias variables en una ubicaci´on son similares a los valores de dichas variables en ubicaciones cercanas, entonces se dice que el patr´on en conjunto exhibe una autocorrelaci´on espacial positiva (autocorrelaci´on). Por el contrario, se dice que existe autocorrelaci´on espacial negativa cuando las observaciones que est´an cerca en el espacio tienden a ser m´as diferentes en los valores de las variables que las observaciones que est´an m´as separadas (en contradicci´on con la ley de Tobler). La autocorrelaci´on cero se produce cuando los valores variables son independientes de la ubicaci´on (en esta situaci´on los datos no ser´ıan espaciales y por tanto no aplicaremos las t´ecnicas que estudiaremos en este trabajo). En el an´alisis de autocorrelaci´on espacial se necesita una medida de contig¨uidad que podemos definirla de manera general como una relaci´on de vecindad. Estas 20 CAP´ ITULO 2. AN ´ ALISIS EXPLORATORIO DATOS ESPACIALES relaciones pueden ser de tres tipos principalmente, caso de torre (Rook’s case), caso de alfil (Bishop’s case) y caso de la reina (Queen’s case) que reciben este nombre por analog´ıa con los movimientos de las figuras de ajedrez. La contig¨uidad de tipo torre es por un vecindario de 4 ubicaciones adyacentes a cada celda (posiciones N, S, E, W), la de tipo alfil solo considera las diagonales de cada celda (posiciones NE, SE, NW, SW) mientras que el ´ultimo tipo considera un vecindario de ocho celdas (posiciones N, S, E, W, NE, NW, SW, SE). En la siguiente figura podemos observar m´as intuitivamente estos tipos de contig¨uidad. Figura 2.2: Tipos de relaciones de vecindad Un aspecto crucial de la definici´on de la autocorrelaci´on espacial es la determinaci´on de ubicaciones cercanas, es decir, aquellas ubicaciones que rodean un punto de datos dado que podr´ıa considerarse que influyen en la observaci´on en ese punto de datos. Sin embargo, la determinaci´on de este vecindario tiene un cierto grado de arbitrariedad. El n´umero de observaciones en el vecindario establecido para cada ubicaci´on puede expresarse mediante una matriz de ponderaciones W: W=     W11 W12 ··· W1n W21 W22 ··· W2n . . .. . .. . . Wn1Wn2··· Wnn      donde nrepresenta el n´umero de ubicaciones (observaciones). La entrada en la fila i(i= 1, . . . , n) y columna j(j= 1, . . . , n), denotado como Wij, corresponde al par (i, j) de ubicaciones. Los elementos diagonales de la matriz son cero, por convenio, mientras que los elementos no diagonales Wij(i6=j) toman valores distintos de cero (uno, para una matriz binaria) cuando las ubicaciones iyjse consideran vecinas. En concreto para los datos de ´areas, nrepresentar´ıa el n´umero de ´areas donde cada ´area se identifica con un punto (centroide) del que se conocen sus coordenadas 2.3. AUTOCORRELACI ´ ON ESPACIAL 21 cartesianas y cada elemento de la matriz Wcorresponde a la relaci´on de dos ´areas. Dicha matriz, W, a menudo est´a estandarizada por filas, es decir, cada suma de filas en la matriz se hace igual a uno, as´ı, los valores individuales Wij est´an proporcionalmente representados. Esto se hace para que cada vecino de un ´area tenga el mismo peso y la suma de todos los Wij (sobre j) sea igual a uno n X j=1 Wij = 1!. Una forma de representar las relaciones espaciales con datos de ´area es a trav´es del concepto de contig¨uidad. Los vecinos contiguos de primer orden se definen como ´areas que tienen un l´ımite com´un. Formalmente: Wij =1 si el ´area jcomparte un l´ımite com´un con el ´area i 0 caso contrario Alternativamente, dos ´areas iyjpueden definirse como vecinas cuando la distancia dij entre sus centroides es menor que un valor cr´ıtico dado, pongamos d, donde las distancias se calculan a partir de la informaci´on sobre latitud y longitud, s(i),de las ubicaciones del centroide: Wij =1 si dij < d, (d < 0) 0 caso contrario (2.2) La especificaci´on basada en la distancia (2.2) de la matriz de ponderaciones depende de un valor de distancia cr´ıtica dado, d. Sin embargo, cuando hay un alto grado de heterogeneidad en el tama˜no de las unidades de ´area, puede ser dif´ıcil encontrar una distancia cr´ıtica satisfactoria. En tales circunstancias, una peque˜na distancia tender´a a conducir a muchas “islas”, mientras que una distancia elegida para garantizar que cada unidad de ´area tenga al menos un vecino puede producir un tama˜no inaceptablemente grande de n´umero de vecinos para las unidades de ´area m´as peque˜nas. Una soluci´on com´un a este problema es restringir la estructura contigua a los k-vecinos m´as cercanos, y por lo tanto excluir las “islas”(´areas que por no estar a una distancia dde otras ´areas se podr´ıa decir que no tiene vecinos) y forzar a cada unidad de ´area a tener el mismo n´umero kde vecinos. Formalmente: Wij =1 si el centroide de jes uno de los kcentroides m´as cercanos al de i 0 caso contrario En este caso el n´umero de vecinos, k, es el par´ametro de este esquema de ponderaci´on. Tambi´en podemos cambiar la ponderaci´on para que los vecinos m´as distantes obtengan menos peso introduciendo un par´ametro θque permita indicar la tasa de disminuci´on de los pesos. 22 CAP´ ITULO 2. AN ´ ALISIS EXPLORATORIO DATOS ESPACIALES Un esquema de ponderaci´on continua com´unmente utilizado se basa en la funci´on de distancia inversa, de modo que los pesos est´an inversamente relacionados con el ´area de separaci´on de distancia iy el ´area jdonde el par´ametro θse estima o se establece a priori: Wij =d−θ ij si la distancia entre centroides dij < d (d > 0, θ > 0) 0 caso contrario Otro esquema de ponderaci´on continua se deriva de la funci´on exponencial negativa que viene dado por: Wij =exp(−θdij ) si la distancia entre centroides dij < d (d > 0, θ > 0) 0 caso contrario donde θes un par´ametro que puede ser estimado, pero generalmente es elegido a priori por el investigador. Una elecci´on muy com´un es θ= 2. Evidentemente, existe una gran cantidad de matrices de ponderaciones espaciales para el mismo modelo espacial. Es importante tener siempre en cuenta que los resultados de cualquier an´alisis estad´ıstico espacial est´an condicionados a la matriz de pesos espaciales elegida. Con frecuencia, es una buena pr´actica verificar la sensibilidad de las conclusiones con la elecci´on de diferentes matrices de ponderaciones espaciales y ver si dichas conclusiones se siguen cumpliendo para la mayor´ıa de las matrices de ponderaciones elegidas, a menos que exista una raz´on convincente sobre bases te´oricas para considerar solo una. 2.4. Estad´ıstico General de Productos Cruzados La mayor´ıa de medidas de autocorrelaci´on espacial utilizan el producto cruzado de matrices definido de la siguiente forma: n X i=1 n X j=1 MijWij donde Wij son las entradas de la matriz de pesos espaciales definidas en la secci´on anterior y Mij es una medida de proximidad entre las ´areas iyjen alguna otra dimensi´on, por ejemplo, la distancia eucl´ıdea, distancia esf´erica o la distancia Manhattan. Veamos un ejemplo simple y pr´actico que servir´a para visualizar todos estos conceptos m´as f´acilmente. 2.4. ESTAD´ ISTICO GENERAL DE PRODUCTOS CRUZADOS 23 Ubicaci´on A B C D Valores 1 4 5 2 En este ejemplo estamos utilizando la contig¨uidad de tipo torre y utilizaremos la distancia eucl´ıdea para el c´alculo de la matriz de medida de proximidad entre las ´areas. El primer cuadro nos da informaci´on de c´omo est´an repartidas cuatro ´areas (A,B,C,D), por ejemplo el ´area A es contig¨ua a las ´areas B y C, mientras que el segundo nos indica los valores de una variable de cada ´area. La matriz de pesos Wla formar´ıamos creando una matriz binaria 4 ×4 donde las columnas y las filas est´an etiquetadas por las ubicaciones de las celdas. Cuando dos celdas son adyacentes en los datos originales, ingresamos un 1 en la entrada (i, j) de la matriz W, si dos celdas no son adyacentes ingresamos un 0. W=    0 1 1 0 1 0 0 1 1 0 0 1 0 1 1 0     El siguiente paso antes de determinar el producto cruzado es crear la matriz M, que tiene la misma dimensi´on que la matriz W, donde las entradas se basan en alguna otra medida de distancia entre todas las celdas. Tomemos la distancia eucl´ıdea al cuadrado d(xi, xj) = (xi−xj)2. Para ello escogemos la celda A y restamos su valor al de cada una de las otras celdas y lo elevamos al cuadrado. Por tanto la matriz Mnos queda: M=    0 9 16 1 9014 16 1 0 9 1490     Veremos ahora la multiplicaci´on de Wij ∗Mij para despu´es sumar las entradas de esta nueva matriz 4×4 y as´ı obtener el valor del estad´ıstico general de productos cruzados. Wij ∗Mij =    0 9 16 0 9004 16 0 0 9 0490     24 CAP´ ITULO 2. AN ´ ALISIS EXPLORATORIO DATOS ESPACIALES Tomando la suma de las entradas de la ´ultima matriz llega a un producto cruzado de 76 n X i=1 n X j=1 WijMij = 76!. Sin embargo, este valor no significa nada porque ser´ıa necesario un intervalo de referencia. El estad´ıstico de productos cruzados no proporciona informaci´on sobre la estructura del mapa original porque no sabemos en qu´e contexto poner el valor de 76. El producto cruzado general es un estad´ıstico en el sentido de que la matriz resulta de una muestra de varias matrices posibles. El valor del producto cruzado se puede comparar con el rango de valores que se pueden producir si se crean varios mapas con el mismo conjunto de valores mediante la asignaci´on aleatoria de esos valores a las ´areas que tenemos. En nuestro ejemplo tenemos cuatro valores (1,4,5,2) que asignamos a las ´areas (A,B,C,D), respectivamente. Sin embargo este mapa original lo podemos considerar como una posible permutaci´on de los valores en esas ´areas. Tenemos 4! mapas distintos que podr´ıan producirse aleatoriamente si cada uno de los cuatro valores se asignan al azar, con reemplazo, a un ´area individual. Esto nos dar´ıa 4! ´ındices de productos cruzados con los que podemos crear una distribuci´on a este estad´ıstico. Obviamente en un mapa con mayor n´umero de ´areas no ser´ıa computacionalmente posible obtener el valor utilizando todos los ´ındices de productos cruzados que tendr´ıamos con todas las permutaciones posibles de los valores en las ´areas del mapa. Para ello se utilizan unas t´ecnicas de muestreo como el enfoque de Monte Carlo donde se toma una muestra aleatoria de todas las permutaciones posibles del mapa. Esta distribuci´on o valores del producto cruzado se utiliza mayoritariamente en las medidas de autocorrelaci´on espacial como por ejemplo el estad´ıstico Ide Moran o cde Geary que veremos en profundidad en el siguiente cap´ıtulo. Cap´ıtulo 3 Medidas de Autocorrelaci´on Espacial Las medidas de autocorrelaci´on espacial se distinguen entre globales y locales. Global implica que todos los elementos de la matriz Wse aplican a una evaluaci´on de la autocorrelaci´on, es decir, todas las asociaciones espaciales de ´areas se incluyen en el c´alculo de la autocorrelaci´on espacial. En otras palabras, el an´alisis espacial global es un estad´ıstico para resumir toda la zona de estudio, por tanto asume homogeneidad, y esto permite que se produzca un valor para cualquier matriz de pesos espaciales. Por el contrario, las medidas locales eval´uan la autocorrelaci´on espacial asociada con una o algunas unidades de ´area particulares. 3.1. Medidas globales Las dos medidas que m´as se han utilizado para el caso de las unidades de ´area son los estad´ısticos Ide Moran y cde Geary. Moran [10] introdujo la primera medida de autocorrelaci´on espacial para estudiar fen´omenos estoc´asticos que se distribuyen en el espacio en dos o m´as dimensiones. La medida Ide Moran viene dada por la siguiente expresi´on: I=n n X i=1 n X j=1 Wij n X i=1 n X j=1 Wij(zi−¯z)(zj−¯z) n X i=1 (zi−¯z)2 (3.1) donde nes el n´umero de ´areas que consideramos y Wij es una matriz de pesos definida anteriormente. En esta medida concreta, la distancia que utiliza para 25 32 CAP´ ITULO 3. MEDIDAS DE AUTOCORRELACI ´ ON ESPACIAL expresar formalmente en una matriz de pesos binarios sim´etrica W(d) (v´ease 2.3), con los elementos Wij(d) indexados por la distancia d. Los estad´ısticos GiyG∗ imiden el grado de asociaci´on local para cada observaci´on ien un conjunto de datos que contiene nobservaciones. Formalmente, estas medidas para cada observaci´on (´area) ise pueden expresar como: Gi(d) = n X j6=i Wij(d)zj n X j6=i zj G∗ i(d) = n X j=1 Wij(d)zj n X j=1 zj El estad´ıstico Gipuede interpretarse como una medida de la agrupaci´on de valores similares entorno a una observaci´on particular i, independientemente del valor en esa ´area, mientras que la estad´ıstica G∗ iincluye el valor dentro de la medida de la agrupaci´on. Un valor positivo indica agrupamiento de valores altos y un valor negativo indica un grupo de valores bajos. La distribuci´on de estos estad´ısticos es normal si la distribuci´on subyacente de las observaciones es normal. Pero si la distribuci´on es asim´etrica, la prueba solo se aproxima a la normalidad a medida que aumenta la distancia cr´ıtica d. Por tanto, bajo estas circunstancias, la normalidad del estad´ıstico de contraste solo puede garantizarse cuando el n´umero de ´areas vecinas es grande. El estad´ıstico local Iide Moran para cada ´area ise define: Ii= (zi−¯z) n X j∈Ji Wij(zj−¯z)2 donde Jidenota el conjunto de vecinos del ´area iy ¯zdenota la media de las obseraciones que son vecinas al ´area i. Podemos comprobar que la suma de los Iipara todas las observaciones ies proporcional al estad´ıstico Ide Moran dado por la ecuaci´on (3.1) n X i=1 Ii= n X i=1 (zi−¯z) n X j∈Ji Wij(zj−¯z) 3.2. MEDIDAS LOCALES 33 El estad´ıstico local cide Geary para cada ´area ise define: ci= n X j∈Ji Wij(zi−zj)2 De la misma forma que para el estad´ıstico Iide Moran es f´acil comprobar que la suma de los cipara todas las observaciones ies proporcional al estad´ıstico cde Geary dado por la ecuaci´on (3.2) n X i=1 ci= n X i=1 n X j∈Ji Wij(zi−zj)2 Estos estad´ısticos de LISA1,Iiyci, tienen dos prop´ositos. Por un lado, pueden ser vistos como indicadores de puntos calientes, similares a GiyG∗ i. Por otro lado, se pueden usar para evaluar la influencia de ubicaciones individuales (observaciones) en la magnitud del estad´ıstico de autocorrelaci´on espacial global correspondiente, Ide Moran y cde Geary. 1LISA: (Local Indicator of Spatial Association),son estad´ısticas que descompone el ´ındice global de autocorrelaci´on y verifica en cu´anto contribuye cada unidad espacial a la formaci´on del valor general, permitiendo capturar de forma simult´anea el grado de asociaci´on espacial y la heterogeneidad resultante del aporte de cada unidad espacial. 34 CAP´ ITULO 3. MEDIDAS DE AUTOCORRELACI ´ ON ESPACIAL Cap´ıtulo 4 Modelado de los datos de ´area El an´alisis exploratorio visto en el cap´ıtulo 2 es un paso preliminar para abordar el modelo que busca establecer relaciones entre las observaciones de una variable y las observaciones de otras variables, registradas para cada unidad de ´area. Para estos modelos de regresi´on espacial se utilizan datos recogidos a partir de un dise˜no transversal simple, sin tener en cuenta datos de panel1. Este dise˜no toma una muesta de datos de la poblaci´on objetivo y se obtiene informaci´on de esta muestra una ´unica vez. Tambi´en consideraremos que los datos en cuesti´on se distribuyen aproximadamente como una distribuci´on normal. Para los datos en los que la variable de inter´es sea un recuento o una proporci´on esperamos que los modelos para dichos datos tengan una distribuci´on de Poisson o Binomial ya que a medida que el n´umero de ´areas crezca, estas distribuciones tienden asint´oticamente a una distribuci´on normal. En este cap´ıtulo, a partir del modelo de regresi´on lineal (m´ultiple) veremos diferentes modelos de regresi´on espacial, tales como: modelo de retardo espacial, modelo de error espacial, modelo de orden superior y el modelo espacial de Durbin. En la secci´on 4.2 podremos ver pruebas para la dependencia espacial que nos ayudar´a a saber si los datos que estamos analizando poseen autocorrelaci´on espacial. Tambi´en estudiaremos las estimaciones de los par´ametros de algunos de los modelos nombrados. 4.1. Tipos de modelos de regresi´on espacial Para comenzar esta secci´on, recordaremos el modelo de regresi´on lineal que ser´a nuestro punto de partida para estudiar los principales modelos de regresi´on espacial. 1Un conjunto de datos de panel recoge observaciones sobre m´ultiples fen´omenos a lo largo de determinados per´ıodos. 35 36 CAP´ ITULO 4. MODELADO DE LOS DATOS DE ´ AREA El modelo de regresi´on lineal (m´ultiple) tiene como forma funcional: Y= Q X q=1 Xqβq+ε(4.1) y su versi´on muestral yi= Q X q=1 xiqβq+εii= 1,···n(4.2) donde yies una observaci´on de la variable dependiente o de inter´es, xiq es una observaci´on en una variable explicativa con q= 1, . . . , Q,βqes el coeficiente de regresi´on que mide la influencia por s´ı sola de la q-variable explicativa en la variable dependiente, es decir, mide el cambio en Ypor cada cambio unitario en Xqmanteniendo las restantes variables explicativas constantes. El t´ermino εies el error aleatorio, que puede ser debido a variables no controladas y a la variabilidad muestral. Para estos t´erminos de error asumimos que εison variables independientes e id´enticamente distribuidas a una variable normal con media cero y varianza una constante σ2, esto es, εi∼ N(0, σ2) donde E(εi·εj) = E(εi)·E(εj) = 0. En notaci´on matricial este modelo de regresi´on lineal puede ser escrito de la siguiente forma: Y=Xβ+ε(4.3) donde Yes un vector n×1 en el que se encuentran las nobservaciones de la variable dependiente, Xes una matriz n×Qque muestra las observaciones de las variables explicativas, βes el vector Q×1 de par´ametros de regresi´on asociados a dichas variables explicativas y el vector εde dimensi´on n×1 de t´erminos de error. Este modelo es v´alido para entender la relaci´on funcional entre la variable dependiente y las variables explicativas y estudiar cu´ales pueden ser las causas de la variaci´on de Y. Para ello, a partir de la informaci´on muestral, obtenemos un estimador de β(ˆ β) para as´ı obtener una predicci´on de Ya partir de las variables explicativas X1··· , XQ. Podemos ver que bajo la suposici´on de que las observaciones son independientes, el modelo se simplifica notablemente. Sin embargo a la hora de estudiar los datos espaciales, esta simplificaci´on puede desembocar a unos resultados parciales inconsistentes debido a la dependencia espacial. Esta dependencia espacial puede estar presente en variables explicativas, variable dependiente o en los residuos (t´erminos de error). Cuando la dependecia espacial se encuentra en la varible dependiente los modelos se denominan modelos de retardo espacial mientras que si est´a en los residuos se denominan modelos de error espacial. Cuando 4.1. TIPOS DE MODELOS DE REGRESI ´ ON ESPACIAL 37 est´a presente en las variables explicativas se llaman modelos de regresi´on cruzada omodelos X espacialmente retardados pero, en contraste con los otros dos modelos, no precisan de procedimientos especiales para la estimaci´on. 4.1.1. Modelos de retardo espacial Estos modelos presentan la correlaci´on espacial (dependencia) en la variable dependiente (Y). Estos modelos son extensiones de modelos de regresi´on de tipo (4.2). Permiten a las observaciones de la variable dependiente Yen el ´area i(i= 1, . . . , n) depender de observaciones en ´areas vecinas. El modelo de retardo espacial b´asico, llamado modelo autorregresivo espacial de primer orden (SAR), toma la forma muestral: yi=ρ n X j=1 Wijyj+ Q X q=1 xiqβq+εii= 1,···n(4.4) donde el t´ermino error εies independiente e id´enticamente distribuido. Wij es el elemento (i, j)th de la matriz de pesos n×nW(v´ease 2.2). El escalar ρes un par´ametro a estimar que determinar´a el nivel de relaci´on autorregresiva espacial entre yiyX j Wijyj(combinaci´on lineal de observaciones espacialmente relacionadas basadas en elementos distintos de cero en la i-´esima fila de W). En notaci´on matricial se escribe: Y=ρWY+Xβ+ε(4.5) WYse conoce como una variable dependiente espacialmente retardada. El dominio de ρviene definido por (w−1 min, w−1 max) donde wmin ywmax representan el m´aximo y m´ınimo autovalor de la matriz W. Si la matriz de pesos Wes estoc´astica, se puede probar que los autovalores de dicha matriz est´an en el intervalo [−1,1]. Podemos ver la prueba de este lema en el cap´ıtulo 2 de Baris y Mete [8]. Cuando Wes estoc´astica por filas, una entrada (i, j) distinta de cero nos indica que la j-´esima observaci´on se usar´a para ajustar la predicci´on de la fila i. Este es un modelo estructural, que puede ser expresado en forma reducida como sigue: Y= (I−ρW)−1(Xβ+ε).(4.6) En consecuencia su valor esperado viene dado por la expresi´on: E[Y] = (I−ρW)−1(Xβ) 38 CAP´ ITULO 4. MODELADO DE LOS DATOS DE ´ AREA El t´ermino (I−ρW)−1se denomina multiplicador espacial e indica que el valor esperado de cada observaci´on yidepender´a de una combinaci´on lineal de valores Xtomados por observaciones vecinas, escalado por el par´ametro de dependencia ρ. 4.1.2. Modelos de error espacial Estos modelos explican la dependencia espacial en el t´ermino de error. La dependencia del error espacial puede surgir de variables latentes no observables que est´an correlacionadas espacialmente. Tambi´en puede surgir de los l´ımites del ´area que no reflejan con precisi´on la vecindad que dan lugar a las variables recopiladas para el an´alisis. La especificaci´on m´as com´un es un proceso autorregresivo espacial de primer orden, dado por: εi=λ n X j=1 Wijεj+uii= 1,···n donde λes el par´ametro autorregresivo, y uiun t´ermino aleatorio de error, que asumimos que es independiente e id´enticamente distribuido. Si |λ|<1, entonces nos queda: ε= (I−λW)−1u Insertando esta ecuaci´on en el modelo de regresi´on est´andar obtenemos: Y=Xβ+ (I−λW)−1u(4.7) con E[uu0] = σ2I. En consecuencia la matriz de varianza-covarianza del vector de errores aleatorio es: E[εεt] = σ2(I−λW)−1(I−λWt)−1 El modelo (4.7), denominado modelo de error espacial (SEM) puede verse como una combinaci´on de un modelo de regresi´on est´andar con un modelo autorregresivo espacial en el t´ermino de error ε. 4.1.3. Modelos de orden superior Estos modelos usan m´ultiples matrices de pesos proporcionando una generalizaci´on sencilla de los modelos SAR y SEM. Por ejemplo, si usamos dos matrices espaciales W1yW2(no necesariamente una distinta de la otra) para combinar los modelos de retardo espacial y de error b´asicos da lugar: 4.1. TIPOS DE MODELOS DE REGRESI ´ ON ESPACIAL 39 Y=ρW1Y+Xβ+ε ε=λW2ε+u u∼ N(0, σ2 uI) En este ejemplo con dos matrices de pesos, si establecemos ρ= 0, queda eliminado la influencia de las observaciones de la variable objetivo en los vecinos, resultando el modelo de error espacial dado por (4.7). Por otro lado, si λ= 0 se elimina el t´ermino de perturbaci´on de retardo espacial generando el modelo de retardo espacial dado por (4.4). 4.1.4. Modelo espacial de Durbin En algunos estudios puede sospecharse que la dependencia influye sobre la variable objetivo o dependiente a trav´es de la propia variable objetivo y de las variables predictoras o explicativas. Se puede construir un modelo tipo SAR, visto antes dado por la ecuaci´on (4.4), aumentado por variables explicativas espacialmente retardadas: Y=ρWY+Xβ+W¯ Xγ+ε(4.8) donde ¯ Xes la matriz n×(Q−1) no constante de variables explicativas. El modelo puede reducirse como: Y= (I−ρW)−1(Xβ+W¯ Xγ+ε) con ε∼ N(0, σ2I) donde γes un vector (Q−1)×1 de par´ametros que miden el impacto marginal de las variables explicativas de observaciones vecinas (´areas) sobre la variable dependiente Y.W¯ Xproduce variables explicativas retardadas espacialmente que reflejan un promedio de observaciones vecinas. Una motivaci´on para el uso de SDM se basa en dos circunstancias que suelen surgir cuando aplicamos un modelo de regresi´on espacial a datos de ´area. Una de esas circunstancia es la dependencia espacial en las perturbaciones de un modelo de regresi´on MCO (m´ınimos cuadrados ordinarios). La segunda circunstancia es la existencia de una variable explicativa omitida que muestra una covarianza no nula con una variable incluida en el modelo, y las variables omitidas son probables 40 CAP´ ITULO 4. MODELADO DE LOS DATOS DE ´ AREA cuando se trata de muestras de datos de ´area. Tambi´en es importante este modelo de Durbin porque generaliza muchos de los modelos que se utilizan: Si imponemos γ= 0 nos lleva a un modelo tipo SAR (4.5). Si imponemos γ=−ρ¯ βnos lleva a un modelo tipo SEM (4.7). La restricci´on ρ= 0 da como resultado un modelo de regresi´on Xespacialmente retardados que asume la independencia entre las observaciones de la variable dependiente, pero incluye caracter´ısticas de las ´areas vecinas, en forma de variables explicativas espacialmente retardadas. Por ´ultimo, si ρ= 0 y γ= 0 obtenemos el modelo de regresi´on est´andar dado por (4.3). Si definimos Z= [X W ¯ X] y δ=β γse puede observar de forma m´as evidente que es un modelo tipo SAR (4.4) quedando el modelo escrito de la siguiente forma: Y=ρWY+Zδ+ε(4.9) que si agrupamos nos queda: Y= (I−ρW)−1(Zδ+ε) Finalmente, podemos generalizar m´as el modelo espacial de Durbin de la siguiente manera: Y=ρW1Y+Xβ+W1¯ Xγ+ε ε=λW2ε+u u∼ N(0, σ2 uI) donde W1yW2son matrices n×nde pesos espaciales que pueden ser iguales o distintas. 4.2. Pruebas para Dependencia Espacial La presencia de dependencia espacial en un modelo de regresi´on se puede detectar mediante pruebas de diagn´ostico. A continuaci´on se incluyen tres pruebas para detectar la dependencia espacial en t´erminos de error: el test de Moran y los test de los multiplicadores de Lagrange para los modelos espaciales de retardo (4.4) y de error (4.7) visto en la secci´on anterior. 4.2. PRUEBAS PARA DEPENDENCIA ESPACIAL 41 El estad´ıstico m´as conocido para hacer el contraste de hip´otesis de existencia de autocorrelaci´on espacial es el estad´ıstico Ide Moran, definido anteriormente por la ecuaci´on (3.1), aplicada a los residuos de regresi´on. Si denotamos aecomo el vector n×1 de m´ınimos cuadrados ordinarios (MCO) definido como e=Y−Xb β, el estad´ıstico Ide Moran queda redefinido de la siguiente forma: I=n W0 e0We e0e con W0el vector normalizado que aparece en el estad´ıstico Ide Moran, W0= n X i=1 n X j=1 Wij ye0ela suma de los cuadrados de los residuos. El estad´ıstico Ise interpreta como el coeficiente de una regresi´on de MCO de Wesobre e. En la pr´actica, la inferencia mediante la prueba Ide Moran, se basa en una aproximaci´on normal utilizando un valor estandarizado obtenido restando la media, bajo la hip´otesis nula de no dependencia espacial, y dividiendo por la ra´ız cuadrada de la varianza. El test estad´ıstico queda: Z(I) = I−E(I) pV ar(I)∼ N(0,1) donde el valor esperado y la varianza del estad´ıstico Ide Moran est´an expresados con m´as detalle en la secci´on (3.1) La causa de la dependencia espacial bajo la hip´otesis alternativa no est´a especificada, por lo tanto, la prueba de Moran es una prueba general para detectar autocorrelaci´on espacial. Al aplicar esta prueba a los residuos, debemos tener cuidado, ya que si estimamos los Qcoeficientes de regresi´on, los residuos observados est´an sujetos a Qrestricciones lineales, lo cual quiere decir que dichos residuos estar´an correlacionados y el procedimiento de prueba mediante Ide Moran no ser´a v´alido. Sin embargo, en caso de que el n´umero de coeficientes Qfuese peque˜no en comparaci´on con n(n´umero de ´areas) podr´ıamos ignorar esta correlaci´on y de este modo validar la prueba. Una prueba alternativa se basa en el principio del multiplicador de Lagrange (LM), sugerido por Burridge [4], que se calcula a trav´es de los residuos MCO y que es diferente seg´un qu´e tipo de dependencia espacial queramos probar. Para una dependencia de tipo error dado por la ecuaci´on (4.7), el estad´ıstico LM(error) viene dado por: 48 CAP´ ITULO 4. MODELADO DE LOS DATOS DE ´ AREA Cap´ıtulo 5 Implementaci´on de datos espaciales en R En este ´ultimo cap´ıtulo nos centraremos en ver los modelos de regresi´on espacial, estudiados en el cap´ıtulo anterior, en el programa estad´ıstico R. Analizaremos los paquetes que deben instalarse y las funciones m´as importantes que utilizaremos, as´ı como alg´un ejemplo con datos reales con la ayuda de Bivand y otros [1] y la p´agina web The Geospatial and Farming Systems Research Consortium [15]. Los paquetes m´as esenciales y usados para los datos espaciales de tipo ´area son rgdal, raster, sp yspdep. Los tres primeros paquetes son necesarios para utilizar las funciones del paquete spdep, el cual, como veremos a continuaci´on, es el que visualmente m´as se usa para el an´alisis espacial. Es decir, para utilizar las funciones que analizar´an la autocorrelaci´on espacial y los distintos modelos vistos en el cap´ıtulo 3 y 4, respectivamente, debemos cargar los dem´as paquetes y librerias. 5.1. Paquetes previos Paquete rgdal Este paquete proporciona enlaces a la librer´ıa de abstracci´on de datos geoespaciales (’GDAL’) y da acceso a las operaciones de proyecci´on/transformaci´on desde la librer´ıa ’PROJ.4’. Su uso m´as general es leer los datos para devolver un objeto espacial y proporcionar la proyecci´on de coordenadas o matrices de coordenadas. Para m´as detalles sobre el uso de este paquete en datos espaciales ver Bivand [2]. Paquete r´aster El paquete r´aster proporciona clases y funciones para manipular datos geogr´aficos (espaciales) en formato ’raster’. Los datos r´aster dividen el espacio en 49 50 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R celdas (rect´angulos, p´ıxeles) de igual tama˜no (en unidades del sistema de referencia de coordenadas). Dichos datos espaciales continuos tambi´en se conocen como datos de “cuadr´ıcula” y se contrastan con datos espaciales discretos (basados en objetos como puntos, l´ıneas o pol´ıgonos). El paquete es ´util cuando se usan conjuntos de datos muy grandes que no se pueden cargar en la memoria de la computadora. Las funciones se ejecutar´an correctamente, ya que procesan archivos grandes en trozos, es decir, leen, computan y escriben bloques de datos, sin cargar todos los valores en la memoria a la vez. El paquete implementa clases para datos r´aster y su uso es muy variado: •Creaci´on de objetos r´aster desde cero o desde archivo. •Manejo de archivos r´aster computacionalmente grandes. •´ Algebra r´aster y funciones de superposici´on. •Funciones de distancia y vecindad. •Conversi´on de pol´ıgonos, l´ıneas y puntos a datos r´aster. •Modelo de predicciones. •Resumir valores r´aster. •F´acil acceso a valores de celda rasterizados. •Dibujar (hacer mapas). •C´alculo de los n´umeros de fila, columna y celda a coordenadas y viceversa. •Leer y escribir varios tipos de archivos r´aster. Con ayuda de Hijmans [7] vemos las diferentes funciones y sus usos con mayor detalle. Paquete sp Clases y m´etodos para datos espaciales. El documento de clases donde reside la informaci´on de ubicaci´on espacial, para datos 2D o 3D. Se proporcionan funciones de utilidad, por ejemplo, para trazar datos como mapas, selecci´on espacial, as´ı como m´etodos para recuperar coordenadas, subconjuntos, impresi´on, resumen, etc. Este paquete no proporciona muchas funciones para modificar o analizar datos espaciales, pero las clases que se definen son muy utilizadas por otros muchos paquetes de R. El paquete sp introduce una serie de clases: Para datos vectoriales, los tipos b´asicos son SpatialPoints, SpatialLines y SpatialPolygons. Para almacenar tambi´en los atributos, 5.2. PAQUETE SPDEP 51 las clases est´an disponibles con estos nombres m´as DataFrame, por ejemplo, SpatialPolygonsDataFrame y SpatialPointsDataFrame. Cuando se hace referencia a cualquier objeto con un nombre que comienza con Spatial, es com´un escribir Spatial*. Cuando se hace referencia a un objeto SpatialPolygons o SpatialPolygonsDataFrame, es com´un escribir SpatialPolygons*. Para m´as informaci´on sobre este paquete y la utilidad de las funciones mencionadas se recomienda ver Pebesma y otros [11] y Bivand y otros [1]. 5.2. Paquete spdep Es uno de los paquetes m´as importantes a la hora del an´alisis exploratorio y la inferencia estad´ıstica de los datos espaciales. Posee una colecci´on de funciones para crear matrices de pesos espaciales a partir de la contig¨uidad de los pol´ıgonos, para el estad´ıstico general de productos cruzados, para pruebas para el contraste de existencia de autocorrelaci´on espacial como el estad´ıstico Ide Moran o cde Geary y sus versiones locales (ver cap´ıtulos 2 y 3). Tambi´en contiene funciones para estimar modelos del tipo SAR o SEM (ver cap´ıtulo 4). A continuaci´on veremos las funciones en R m´as relevantes para llevar a cabo este an´alisis e inferencia de datos de ´area. Para mayor detalle de las funciones, a continuci´on descritas y otras funciones que incluye el paquete, se recomienda Bivand y Wong [3]. 5.2.1. Vecindad de datos de ´area Como hemos visto en teor´ıa, antes de poder hacer un contraste de hip´otesis de autocorrelaci´on espacial o modelar los datos, debemos disponer de una organizaci´on de dichos datos, es decir, debemos tener en cuenta la contig¨uidad que se utiliza y a partir de ello, y de los datos, crear una matriz de pesos para su posterior uso en la inferencia estad´ıstica de los datos de ´area. Para conseguir estas matrices de peso, dentro del paquete spdep usamos la clase nb. Veamos las funciones m´as frecuentemente usadas en R de dicha clase para crear matrices de pesos a partir de unos datos: poly2nb: La funci´on crea una lista de vecinos basada en regiones con l´ımites contiguos, que comparte uno o m´as puntos de l´ımite. Uso poly2nb(pl, row.names = NULL, snap=sqrt(.Machine$double.eps), queen=TRUE,...) 52 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R Argumentos pl: Lista de pol´ıgonos de la clase SpatialPolygons. row.names: Vector de identidad de las regiones para ser a˜nadidas a la lista de los vecinos. snap: Los puntos de l´ımites inferiores a la distancia de separaci´on para indicar la contig¨uidad. queen: Si TRUE el tipo de contig¨uidad que se utiliza es la de tipo reina. Si FALSE, la contig¨uidad utlizada es la de tipo torre. Valor Nos devuelve una lista de vecinos de la clase nb. knearneigh: La funci´on devuelve una matriz con los ´ındices de puntos que pertenecen al conjunto de los kvecinos m´as cercanos entre s´ı. Uso knearneigh(x, k=1, longlat = NULL,...) Argumentos x: Matriz de puntos coordenados o un objeto de SpatialPoints. k: N´umero de vecinos m´as cercanos para cada ´area. longlat: TRUE si las coordenadas del punto son grados decimales de longitudlatitud, en cuyo caso las distancias se miden en kil´ometros; si x es un objeto SpatialPoints, el valor se toma del propio objeto. Valor Una lista con los siguientes elementos: nn: Matriz entera de identificadores de n´umero de regi´on. np: N´umero de puntos de entrada. k: Valor de entrada requerido. dimension: N´umero de columnas de x x: Coordenadas de entrada knn2nb: La funci´on convierte un objeto knn devuelto por knearneigh en una lista de vecinos de clase nb con una lista de vectores de n´umeros enteros que contienen identificadores de n´umero de regi´on vecina. Uso knn2nb(knn, row.names = NULL, sym = FALSE) 5.2. PAQUETE SPDEP 53 Argumentos knn: Un objeto knn devuelto po la funci´on knearneigh. sym: Fuerza a la lista de regiones vecinas de salida a ser sim´etrica. nbdists: Dada una lista de enlaces vecinos adyacentes (una lista de vecinos de tipo de objeto nb), la funci´on devuelve las distancias euclideas a lo largo de los enlaces en una lista de la misma forma que la lista de vecinos. Si longlat = TRUE, se usan las distancias geod´esicas (great-circle distance). Uso nbdists(nb, coords, longlat = NULL) nb2listw: Esta funci´on toma un objeto de lista de vecinos y lo convierte en un objeto o matriz de pesos. El estilo de conversi´on predeterminado es W, donde los pesos para cada entidad de ´area se estandarizan para sumar por filas igual a la unidad, es decir, est´a estandarizado por filas. Uso nb2listw(neighbours, glist=NULL, style="W",...) Argumentos neighbours: Un objeto de clase nb. glist: Lista de pesos generales correspondientes a vecinos. style: A partir de una lista binaria de vecinos, en la que las regiones est´an listadas como vecinas o est´an ausentes (por lo tanto, no est´an en el conjunto de vecinos para alguna definici´on), la funci´on agrega una lista de ponderaciones con valores dados por el estilo de esquema de codificaci´on elegido. Bes la codificaci´on binaria b´asica mientras que West´a estandarizada por filas. Hay m´as opciones para el esquema de codificaci´on que no utilizaremos puesto que no lo hemos visto en teor´ıa. 5.2.2. Pruebas para Dependencia Espacial En la parte de teor´ıa del trabajo hemos detallado las pruebas mediante los estad´ısticos globales Ide Moran y cde Geary, as´ı como sus formas locales y los estad´ısticos locales de Getis y Ord. En la pr´actica y, en concreto, en el paquete de R que estamos utilizando, el m´as usado es el estad´ıstico global Ide Moran ya que es efectivo en la mayor´ıa de las situaciones que se pueden plantear. Por este motivo, analizaremos detalladamente las funciones del paquete spdep que utilizan este estad´ıstico para la prueba de dependencia espacial. Tambi´en indicaremos algunas funciones de otros estad´ısticos pero sin entrar en mucho detalle ya que su uso ser´a similar al del estad´ıstico I. 54 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R moran.test: La prueba de Moran para la autocorrelaci´on espacial utilizando una matriz de ponderaciones espaciales en forma de lista de ponderaciones. Las hip´otesis que subyacen a la prueba son sensibles a la forma de la gr´afica de relaciones entre vecinos y otros factores, y los resultados pueden compararse con la funci´on de “moran.mc”basada en permutaciones. Uso moran.test(x, listw, randomisation=TRUE, alternative="greater", rank = FALSE,...) Argumentos x: Vector num´erico de la misma longitud que la lista de vecinas de listw. litw: Un objeto tipo listw creado por ejemplo por nb2listw. randomisation: Varianza de Icalculada bajo la hip´otesis de aleatoriedad. Si FALSE, la hip´otesis utilizada ser´a normalidad. alternative: Una cadena de caracteres que especifique la hip´otesis alternativa, debe ser una de mayor (por defecto), menor o doble. rank: Valor l´ogico. Si FALSE utiliza variables continuas. Si TRUE utiliza la adaptaci´on de Ide Moran para los rangos sugeridos por Cliff y Ord. Valor La funci´on nos devuelve el nombre del m´etodo y el de los datos usados, asi como una descripci´on de la hip´otesis alternativa. Nos devuelve el valor del estad´ıstico Iy el p-valor del test para saber cu´al es el resultado del contraste. lm.morantest: La prueba Ide Moran para autocorrelaci´on espacial en residuos de un modelo lineal estimado. Uso lm.morantest(model, listw, alternative = "greater", resfun=weighted.residuals,.... Argumentos model: Un objeto de la clase lm devuelto por lm. Los pesos pueden especificarse en el ajuste de lm, pero no se deben usar las compensaciones. resfun: Por defecto: weighted.residuals. La funci´on se utiliza para extraer los residuos del objeto lm. Valor Al igual que la funci´on anterior, esta funci´on nos devolver´a los nombres de los m´etodos y datos utilizados, el p-valor, y el estad´ıstico I, su esperanza y varianza. 5.2. PAQUETE SPDEP 55 moran.mc: Una prueba de permutaci´on para el estad´ıstico Ide Moran calculado usando nsim permutaciones aleatorias de x para el esquema de ponderaci´on espacial dado, para establecer el rango de la estad´ıstica observada en relaci´on con los valores simulados nsim. Uso moran.mc(x, listw, nsim, alternative="greater",...) Argumentos Son los mismos argumentos de la funci´on moran.test. localmoran: El estad´ıstico espacial local Ide Moran se calcula para cada zona seg´un el objeto de ponderaciones espaciales utilizado. Uso localmoran(x, listw, alternative = "greater", p.adjust.method="none", mlvar=TRUE) Argumentos mlvar: Predeterminado TRUE: los valores de la Ide Moran local se informan utilizando la varianza de la variable de inter´es (suma de desviaciones al cuadrado sobre n), pero se pueden informar como la varianza de la muestra, dividiendo por (n−1) en su lugar. adjust.x: Si TRUE, los valores de las observaciones de x que no tengan vecinos se omiten para calcular la media de x. A la hora de utilizar cualquiera de estas funciones, no es necesario usar todos los argumentos, es decir, en cualquiera de ellos los fundamentales son los datos que vamos a utilizar para hacer la inferencia y el tipo de matriz de pesos que vamos a usar, xylistw respectivamente. Para el estad´ıstico cde Geary, el paquete spdep tambi´en proporciona unas funciones para su utilizaci´on. Estas son geary.test ygeary.mc que, como en el caso de I, la primera funci´on es para realizar un test para la autocorrelaci´on espacial y la segunda a partir de un muestreo aleatorio mediantes el enfoque de Monte Carlo para tomar varias matrices de pesos y realizar su posterior an´alisis. 5.2.3. Modelos de Regresi´on Espacial en R En este ´ultimo apartado del cap´ıtulo vamos a aplicar los modelos de regresi´on, que hemos visto en la parte de teor´ıa, a los datos en el programa R. Dentro del paquete spdep podemos encontrar la funci´on lagsarlm que se utiliza tanto para el modelo tipo SAR dado por la ecuaci´on (4.4) como para el modelo espacial de 56 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R Durbin, ecuaci´on (4.8), ya que el modelo SDM se puede expresar como un modelo tipo SAR. Tambi´en estudiaremos con detalle la funci´on errorsarlm que la usaremos para el modelo tipo SEM dado por la ecuaci´on (4.7). lagsarlm: Estimaci´on del modelo de regresi´on espacial de retardo. Uso lagsarlm(formula, data = list(), listw, Durbin, method= " ", interval=NULL, tol.solve=1.0e-10,...) Argumentos formula: Una descripci´on simb´olica del modelo a ser ajustado. Los detalles de la especificaci´on del modelo se dan para lm(). data: Un marco de datos opcional que contiene las variables en el modelo. Por defecto, las variables se toman del entorno al que se llama la funci´on. listw: Una lista de datos creados por ejemplo por nb2listw. Durbin: Si FALSE, el modelo utilizado ser´a el de retraso. Si TRUE, el modelo que utiliza es el de Durbin. method: El m´etodo predeterminado de method = “eigen“ usa valores propios y, por lo tanto, tambi´en puede establecer los l´ımites inferior y superior para la b´usqueda de l´ıneas para ρpero es no es factible para grandes n. El otro m´etodo utilizado es method = “Matrix“ para calcular directamente el determinante de la matriz I−ρW. interval: Intervalo de b´usqueda para el par´ametro autorregresivo ρ. tol.solve: La tolerancia para detectar dependencias lineales en las columnas de matrices que se van a invertir (por defecto = 1.0e-10). errorsarlm: Estimaci´on del modelo de regresi´on espacial de error. Uso errorsarlm(formula, data=list(), listw, method=" ", interval = NULL, tol.solve=1.0e-10,...) Argumentos Al ser los mismos argumentos que para el caso de retardo espacial, omitimos su detalle. Como hemos especificado anteriormente, no es necesario utilizar todos los argumentos descritos en ambos modelos. En el ejemplo que daremos a continuaci´on veremos que a la hora de estudiar los modelos y estimar sus coeficientes utilizaremos los siguientes argumentos: formula, data, listw ytol.solve. 5.3. ILUSTRACI ´ ON DE LOS MODELOS EN R 57 5.3. Ilustraci´on de los Modelos en R El conjunto de datos que vamos a utilizar para nuestro ejemplo es houses20001 que son datos de las viviendas de California del Censo 2000. En este conjunto de datos hay 7049 casos con 29 variables: TRACT: Identificador del distrito censal GEOID: Identificador geogr´afico label: Etiqueta de identificaci´on HouseValue: Valor de la vivienda NhousingUn: N´umero de viviendas RecHouses: N´umero de casas para uso recreativo NMobileHom: N´umero de casas m´oviles NBadPlumbi: N´umero de casas con plomer´ıa incompleta NBadKitche: N´umero de viviendas con cocinas incompletas Population: Poblaci´on total Males: N´umero de hombres Females: N´umero de mujeres Under5: N´umero de personas menores de cinco a˜nos White: N´umero de personas que se identifican como blancas (s´olo) Black: N´umero de personas que se identifican como afro-ameriano (s´olo) AmericanIn: N´umero de personas que se identifican como indio americano (s´olo) Asian: N´umero de personas que se identifican como asi´atico americano (s´olo) Hispan: N´umero de personas que se identifican como hispanas (solo) PopInHouse: N´umero de personas que viven en hogares nHousehold: N´umero de hogares 1Datos house2000. Fuente: package spdep 64 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R Para ilustrar los modelos veamos un caso sencillo para este conjunto de datos, veremos el valor de las viviendas a partir de la antig¨uedad de construcci´on de la vivienda y del n´umero de habitaciones de ella. Para ello crearemos una nueva variable llamada age que ser´an los a˜nos en media que tienen las viviendas en los distintos condados (houseV alue ∼age +nBedrooms). Realizaremos este simple modelo lineal (m´ultiple) y as´ı obtendremos los residuos que analizaremos posteriormente: > hh$age <- 2000 - hh$ > f1 <- houseValue ~ age + nBedrooms > m1 <- lm(f1, data=hh) > summary(m1) Call: lm(formula = f1, data = hh) Residuals: Min 1Q Median 3Q Max -222541 -67489 -6128 60509 217655 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -628578 233217 -2.695 0.00931 ** age 12695 2480 5.119 4.05e-06 *** nBedrooms 191889 76756 2.500 0.01543 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 94740 on 55 degrees of freedom Multiple R-squared: 0.3235, Adjusted R-squared: 0.2989 F-statistic: 13.15 on 2 and 55 DF, p-value: 2.147e-05 El modelo lineal (m´ultiple) queda de la siguiente forma: houseV aluei=β0+β1∗agei+β2∗nBedroomsi+εi donde los par´ametros βqestimados son β0=−628578, β1= 12695 y β2= 191889. Por tanto, como hab´ıamos mencionado los a˜nos que lleve construida la casa y el n´umero de habitaciones influyen en el valor de ella. Podemos ver incluso que el n´umero de habitaciones hace que aumente cerca de 200.000$ el precio de la vivienda. Veamos si los residuos poseen autocorrelaci´on espacial a partir del ´ındice 5.3. ILUSTRACI ´ ON DE LOS MODELOS EN R 65 Ide Moran. Para ello, primero tendremos que crear la matriz de ponderaciones espaciales, W: > nb <- poly2nb(hh) > nb Neighbour list object: Number of regions: 58 Number of nonzero links: 276 Percentage nonzero weights: 8.204518 Average number of links: 4.758621 > par(mai=c(0,0,0,0)) > plot(hh) > plot(nb, coordinates(hh), col=’red’, lwd=2, add=TRUE) Con la orden anterior se obtiene la figura 5.1 que muestra los enlaces de los condados vecinos. A continuaci´on obtenemos la inferencia sobre la dependencia espacial a trav´es del estad´ıstico Ide Moran: > lw<- nb2listw(nb) > moran.mc(hh$residuals, lw, 999) Monte-Carlo simulation of Moran I data: hh$residuals weights: lw number of simulations + 1: 1000 statistic = 0.40431, observed rank = 1000, p-value = 0.001 alternative hypothesis: greater 66 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R Figura 5.1: Mapa de los enlaces de los condados vecinos Como tenemos un p-valor de 0,001, se verifica la hip´otesis alternativa y, por tanto, los residuos poseen autocorrelaci´on espacial. Lo ´ultimo que veremos son los tipos de modelos que hemos visto aplicados a estos datos, el modelo tipo SAR (4.4), el modelo de error (4.7) o el modelo de Durbin (4.8): —–Modelo de Retardo Espacial—– houseV aluei=ρ∗lag.housesV aluei+β0+β1∗agei+β2∗nBedroomsi+εi donde εi∼ N(0, σ2) lag.housesV aluei= n X j=1 Wij ∗housesV aluej > m1s = lagsarlm(f1, data=hh, lw, tol.solve=1.0e-30) > summary(m1s) 5.3. ILUSTRACI ´ ON DE LOS MODELOS EN R 67 Call:lagsarlm(formula = f1, data = hh, listw = lw, tol.solve = 1e-30) Residuals: Min 1Q Median 3Q Max -122914.4 -57887.7 -3330.4 41620.1 210803.7 Type: lag Coefficients: (asymptotic standard errors) Estimate Std. Error z value Pr(>|z|) (Intercept) -488761.6 165471.9 -2.9537 0.0031394 age 6796.5 1831.9 3.7101 0.0002072 nBedrooms 146301.3 54732.7 2.6730 0.0075174 Rho: 0.7291, LR test value: 27.797, p-value: 1.3471e-07 Asymptotic standard error: 0.090506 z-value: 8.0559, p-value: 8.8818e-16 Wald statistic: 64.897, p-value: 7.7716e-16 Log likelihood: -731.4781 for lag model ML residual variance (sigma squared): 4494700000, (sigma: 67042) Number of observations: 58 Number of parameters estimated: 5 AIC: 1473, (AIC for lm: 1498.8) LM test for residual autocorrelation test value: 0.46008, p-value: 0.49759 Observamos que nos devuelve 5 par´ametros estimados que son los que en la parte de teor´ıa vimos que necesit´abamos para conseguir el modelo. Estos son: el par´ametro ρque determina el nivel de relaci´on autorregresiva espacial es estimado por 0.7291, los coeficientes βqdados por β0=−488761,6, β1= 6796,5 yβ2= 146301,3 y por ´ultimo la desviaci´on t´ıpica de los residuos generados por este m´etodo, es decir, σ= 67042. Tambi´en nos devuelve el p-valor de un test realizado a los residuos que generan y como dicho p-valor es 0.497 no podemos rechazar la hip´otesis nula de no autocorrelaci´on espacial, es decir, los residuos se pueden considerar independientes y por tanto, tras la aplicaci´on del modelo SAR, la dependencia de los residuos desaparece. —–Modelo de Durbin —– houseV aluei=ρ∗lag.housesV aluei+β0+β1∗agei+β2∗nBedroomsi+ +γ1∗lag.agei+γ2∗lag.nBedroomsi+εi 68 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R donde lag.agei= n X j=1 Wij ∗agej lag.nBedroomsi= n X j=1 Wij ∗nBedroomsj y las variables lag.housesV alueiyεison las mismas que para el modelo SAR. > m1d = lagsarlm(f1, data=hh, lw, Durbin=TRUE,tol.solve=1.0e-30) >summary(m1d) Call:lagsarlm(formula = f1, data = hh, listw = lw, Durbin = TRUE, tol.solve = 1e-30) Residuals: Min 1Q Median 3Q Max -119926.32 -54437.79 -900.84 38916.33 214614.24 Type: mixed Coefficients: (asymptotic standard errors) Estimate Std. Error z value Pr(>|z|) (Intercept) -770005.2 389556.9 -1.9766 0.048085 age 5929.2 2263.2 2.6198 0.008799 nBedrooms 135170.0 56760.2 2.3814 0.017246 lag.age 3096.9 4221.2 0.7337 0.463151 lag.nBedrooms 102897.1 135275.1 0.7607 0.446866 Rho: 0.70085, LR test value: 22.59, p-value: 2.005e-06 Asymptotic standard error: 0.099823 z-value: 7.021, p-value: 2.2036e-12 Wald statistic: 49.294, p-value: 2.2036e-12 Log likelihood: -731.1205 for mixed model ML residual variance (sigma squared): 4510300000, (sigma: 67159) Number of observations: 58 Number of parameters estimated: 7 AIC: 1476.2, (AIC for lm: 1496.8) LM test for residual autocorrelation test value: 1.2435, p-value: 0.26479 5.3. ILUSTRACI ´ ON DE LOS MODELOS EN R 69 En este modelo se estiman 7 par´ametros que son, ρ= 0,7, los coeficientes βq dados por β0=−770005,2, β1= 5929,2 y β2= 135170, los coeficientes γ1= 3096,9 yγ2= 102897,1 y por ´ultimo la desviaci´on t´ıpica de los residuos, σ= 67159. Al igual que el modelo SAR nos devuelve un p-valor para un test de dependencia espacial para los residuos generados, este es, 0,264, el cual nos indica que los nuevos residuos son independientes desapareciendo as´ı la dependencia espacial de ellos. —–Modelo de Error Espacial—– houseV aluei=β0+β1∗agei+β2∗nBedroomsi+λ∗lag.residi+ui donde ui∼ N(0, σ2) lag.residi= n X j=1 Wij ∗εj > m1e = errorsarlm(f1, data=hh, lw, tol.solve=1.0e-30) > summary(m1e) Call:errorsarlm(formula = f1, data = hh, listw = lw, tol.solve = 1e-30) Residuals: Min 1Q Median 3Q Max -121933.170 -58532.837 47.423 40304.197 216867.129 Type: error Coefficients: (asymptotic standard errors) Estimate Std. Error z value Pr(>|z|) (Intercept) -302291.4 188876.7 -1.6005 0.109494 age 6989.1 2234.6 3.1277 0.001762 nBedrooms 132526.7 57138.4 2.3194 0.020374 Lambda: 0.75099, LR test value: 23.731, p-value: 1.1078e-06 Asymptotic standard error: 0.088919 z-value: 8.4457, p-value: < 2.22e-16 Wald statistic: 71.331, p-value: < 2.22e-16 70 CAP´ ITULO 5. IMPLEMENTACI ´ ON DE DATOS ESPACIALES EN R Log likelihood: -733.5113 for error model ML residual variance (sigma squared): 4757300000, (sigma: 68973) Number of observations: 58 Number of parameters estimated: 5 AIC: 1477, (AIC for lm: 1498.8) En este modelo, de nuevo, el n´umero de par´ametros estimados son 5 donde el par´ametro autorregresivo λes estimado por 0.75. Los coeficientes son β0= −302291,4, β1= 6989,1 y β2= 132526,7 y la desviaci´on t´ıpica de los residuos ui es σ= 68973. Tambi´en nos devuelve unos residuos y aunque no comprobamos que sean independientes, por lo estudiado en teor´ıa han de serlo, de lo contrario habr´ıa que modelizar de nuevo los datos teniendo en cuenta estos nuevos residuos. Como conclusi´on, los modelos ajustados se resumen en la siguiente tabla y se incluye la medida AIC (Criterio de Informaci´on de Akaike)2que nos puede permitir seleccionar el modelo m´as adecuado: Modelo Ecuaci´on AIC Lineal houseV alue =−628578 + 1498.8 (m´ultiple) + 12695 ∗age + 191889 ∗nBedrooms Retardo houseV alue = 0,729 ∗lag.houseV alue −488761,6 + 1473 Espacial + 6796,5∗age + 146301,3∗nBedrooms Espacial houseV alue = 0,7∗lag.houseV alue −770005,1 + 1476.2 de + 5929,2∗age + 135170 ∗nBedrooms + Durbin + 3096,9∗lag.age + 102897,1∗lag.nBedrooms Error houseV alue = 0,75 ∗lag.resid −302291,4 + 1477 Espacial + 6989,1∗age + 132526,7∗nBedrooms A partir de la tabla anterior podemos comprobar que, siguiendo el criterio AIC, ambos modelos de regresi´on espacial se ajustan mejor a los datos que el modelo de regresi´on lineal. Para una elecci´on de los modelos espaciales, utilizar´ıamos el modelo de retardo espacial ya que el ´ındice AIC es menor. 2AIC es una medida de la calidad relativa de un modelo estad´ıstico. Su f´ormula viene dado por AIC = 2k−2ln(L) donde kes el n´umero de par´ametros del modelo y Les el m´aximo de la funci´on de verosimilitud. Bibliograf´ıa [1] Bivand, R.A.; Pebesma, E. y G´omez-Rubio, V. (2013). Applied Spatial Data. Analysis with R. Springer. [2] Bivand, R; Keitt, T. y Rowlingson, B. (2018). rgdal: Bindings for the ’Geospatial’ Data Abstraction Library. R package version1.3-6. Recuperado de https://CRAN.R-project.org/package=rgdal. [3] Bivand, R.S y Wong, D. (2018) Comparing implementations of global and local indicators of spatial association TEST, 27(3), 716-748. Recuperado de https://doi.org/10.1007/s11749-018-0599-x. [4] Burridge, P. (1980). On the Cliff-Ord Test for Spatial Correlation. Journal of the Royal Statistical Society. Series B (Methodological), Vol. 42, No. 1. [5] Fischer, M.M. y Wang, J. (2011). Spatial Data Analysis. Models, Methods and Techniques. Springer. [6] Gangodagam, C., Zhou, X. y Lin, H. (2008). Autocorrelation Spatial. Encyclopedia of GIS. Eds. Shekthar, S. y Xiong, H. Springer. [7] Hijmans, R.J. (2017). raster: Geographic Data Analysis and Modeling. R package version 2.6-7. https://CRAN.Rproject.org/package=raster. [8] Kazar, B.M. y Celik, M. (2012). Spatial Autoregression (SAR) Model Parameter Estimation Techniques. Springer. [9] LeSage, J. y Pace, R.K., (2009). Introduction to Spatial Econometrics Statistics. Chapman and Hall/ CRC. [10] Moran, P.A.P. (1950). Notes on Continuous Stochastic Phenomena.Biometrika, Vol. 37, No. 1/2, pp. 17-23. [11] Pebesma, E.J. y Bivand, R.S. (2005). Classes and methods for spatial data in R. R News 5 (2). Recuperado de https://cran.rproject.org/doc/Rnews/. 71 72 BIBLIOGRAF´ IA [12] Sawada M. Global Spatial Autocorrelation Indices - Moran’s I, Geary’s C and the General Cross-Product Statistic. Laboratory for Paleoclimatology and Climatology, Department of Geography University of Ottawa. Recuperado de http://www.lpc.uottawa.ca/publications/moransi/moran.htm. [13] Tobler, W.R. (1970). A computer model simulation of urban growth in the Detroit region. Economic Geaography 46 (2), pp. 234-240. [14] Vitturini M. ; Fillottrani P. y Castro S. (2003) Modelo de datos para datos espaciales. V Workshop de Investigadores en Ciencias de la Computaci´on (WICC 2003). [15] The Geospatial and Farming Systems Research Consortium. Recuperado de http://rspatial.org/.