Modelos para Datos por areas

“All models are fiction, but some stories are better than others” - Richard Berk

“All models are wrong. We make tentative assumptions about the real world which we know are false but which we believe may be useful” - George Box

podría calcular un variograma para datos por areas?

podría emplearlo para hacer krigin?

0.1 Regresión ponderada geográficamente (Geographically weighted regression - GWR)

Los resultados de una regresión por MCe proporcionan un resumen “global” de la relación entre las dos variables. La regresión ponderada geográficamente permite explorar cómo difieren las relaciones entre regiones.

Hay varios paquetes R para GWR. Estos incluyen spgwr, GWmodel y gwrr:

• http://cran.r-project.org/web/packages/spgwr/index.html,
• http://cran.r-project.org/web/packages/GWmodel/index.html,
• http://cran.r-project.org/web/packages/gwrr/index.html

En este ejercicio, emplearemos spgwr. Para utilizar GWR, es necesario seleccionar un esquema de ponderación espacial (un kernel) que defina cuánto peso se le da a cada vecino de cada zona. El tamaño de un núcleo se indica por su ancho de banda - un ancho de banda pequeño significa que a los vecinos cercanos se les dan grandes pesos mientras que a los vecinos más distantes se les dan pesos mucho más pequeños. Con un mayor ancho de banda, los pesos están más `dispersos’ y no hay una disminución marcada de pesos con el aumento de la distancia.

0.1.1 Selección del ancho de banda (bandwidth)

El ancho de banda se puede seleccionar manualmente o se puede utilizar algún tipo de enfoque de optimización para seleccionar uno. Un posible enfoque es utilizar la validación cruzada. Esto implica eliminar cada observación por separado y estimar su valor utilizando los datos restantes. A continuación se calculan las diferencias entre los valores reales observados y los valores estimados. Este procedimiento se completa para anchos de banda múltiples y el ancho de banda que corresponde al error general más pequeño se considera óptimo. El ancho de banda identificado de esta manera se utiliza para el ajuste del modelo GWR.

## el codigo de departamento esta mal para la base 17
library(foreign)
nofetales17<-read.spss("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/nofetal2017/nofetal2017.sav",to.data.frame = TRUE)

nofetales16<-read.spss("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/nofetal2017/nofetal2017.sav",to.data.frame = TRUE)

correla667<-read.csv("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/correlativa667.csv")

names(nofetales17)<-tolower(names(nofetales17))
# names(defun_avpp13)
# table(defun_avpp13$ops_un_digito)

# names(nofetales17)[names(nofetales17)=="cod_munic"]<-"cod_mpio"
nofetales17$cod_mpio<-paste0(nofetales17$cod_dpto,nofetales17$cod_mpio)
# defun_avpp_neoplasia<-defun_avpp13[defun_avpp13$ops_un_digito==1,]
library(tidyverse)
library(readxl)

nofetales667<-left_join(nofetales17,correla667)
# defun_avpp13$ops_tres_digitos<-as.numeric(defun_avpp13$ops_desc_tres_digitos)

base<-nofetales667 %>% 
  #filter(ops_un_digito==5)%>%
  group_by(cod_mpio,ops_desc_un_digito)%>%
  #group_by(cod_mpio,ops_tres_digitos)%>%
  summarise(
    n=n())

base_rs<-reshape(data.frame(base), idvar = "cod_mpio", timevar = "ops_desc_un_digito", direction = "wide")

base_rs[is.na(base_rs)]<-0
pob<-read_xls("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/poblacion2013.xls",sheet = 1)

base<-merge(base_rs, pob, by="cod_mpio")

write.csv(x = base,"D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/667pob.csv")

library(tidyverse)
pob667<-read.csv("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/667pobrename.csv")

desempeno<-read.csv("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/Desempeno municipal mapa.csv")

base<-left_join(desempeno,pob667)

base[is.na(base)]<-0

write.csv(base,"D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/Datos 2015/Desempeno+667.csv")

mpio-dpto	mpio-dpto
Municipio	Municipio
Codigo	Codigo
Departamento	Departamento
Grupo dotaciones	GrDota
Categoria de ruralidad	Categoria de ruralidad
Educación - Cobertura media neta	ECobeMedNet
Educación - SABER 11 Matematicas	ESABEMATE
Educación - SABER 11 Lenguaje	ESABELeng
Educación - Cobertura Transición	ECoberTrans
Salud - Cobertura salud	Scobert
Salud - Vacunación Pentavalente	SVacunPenta
Salud - Mortalidad Infantil	SMortaInf
Servicios públicos - Cobertura eléctrica rural	SPCoberElectrica
Servicios públicos - Cobertura internet	SPCoberInter
Servicios públicos - Cobertura Acueducto	SPCoberAcued
Servicios públicos - Cobertura Alcantarillado	SPCoberAlcanta
Seguridad - Hurtos x 10000 hab	SegHurto
Seguridad - Homicidios x 10000 hab	SegHomici
Seguridad - Violencia Intrafamiliar x 10000 hab	SegVioIntraF
Ingresos tributarios y no tributarios per cápita	Ingresos
Desempeno Municipal	Desempeno
n.CAUSAS EXTERNAS	MortaExterna
n.CIERTAS AFECCIONES ORIGINADA EN EL PERIODO PERINATAL	MortaPeri
n.ENFERMEDADES DEL SISTEMA CIRCULATORIO	MortaCircula
n.ENFERMEDADES TRANSMISIBLES	MortaTransmiso
n.NEOPLASIAS (TUMORES)	MortaTumores
n.SINTOMAS, SIGNOS Y AFECCIONES MAL DEFINIDAS	MortaMalDefini
n.TODAS LAS DEMAS CAUSAS	MortaDemasCaus
poblacion	poblacion

# install.packages("spgwr")
library(sp)
library("spgwr")

## Loading required package: spData

## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source'))`

## NOTE: This package does not constitute approval of GWR
## as a method of spatial analysis; see example(gwr)

library(maptools)

## Checking rgeos availability: TRUE

library(lattice)

mpio <- readShapePoly("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/00. Datos/shp/mpio_desempeno_667 sin islas.shp")

## Warning: readShapePoly is deprecated; use rgdal::readOGR or sf::st_read

## Field name: 'Desempeno+' changed to: 'Desempeno.'

names(mpio)

##  [1] "AREA"       "PERIMETER"  "MUNIPIOS_"  "MUNIPIOS_I" "CODIGO"    
##  [6] "MPIOS"      "ID"         "NOMBRE"     "DPTO"       "X_COORD"   
## [11] "Y_COORD"    "cod_mpio"   "Desempeno." "Desempen_1" "Desempen_2"
## [16] "Desempen_3" "Desempen_4" "Desempen_5" "Desempen_6" "Desempen_7"
## [21] "Desempen_8" "Desempen_9" "Desempen10" "Desempen11" "Desempen12"
## [26] "Desempen13" "Desempen14" "Desempen15" "Desempen16" "Desempen17"
## [31] "Desempen18" "Desempen19" "Desempen20" "Desempen21" "Desempen22"
## [36] "Desempen23" "Desempen24" "Desempen25" "Desempen26" "Desempen27"
## [41] "Desempen28" "Desempen29" "Desempen30"

names(mpio@data)<-c("AREA", "PERIMETER", "MUNIPIOS_", "MUNIPIOS_I", "CODIGO", "MPIOS", 
"ID", "NOMBRE", "DPTO", "X_COORD", "Y_COORD", "cod_mpio","mpio.dpto", "Municipio","Departamento", "GrDota", 
"Categoria.de.ruralidad", "ECobeMedNet", "ESABEMATE", "ESABELeng", 
"ECoberTrans", "Scobert", "SVacunPenta", "SMortaInf", "SPCoberElectrica", 
"SPCoberInter", "SPCoberAcued", "SPCoberAlcanta", "SegHurto", 
"SegHomici", "SegVioIntraF", "Ingresos", "Desempeno", "MortaExterna", 
"MortaPeri", "MortaCircula", "MortaTransmiso", "MortaTumores", 
"MortaMalDefini", "MortaDemasCaus", "DPMP", "MPIO", "poblacion"
)


mpio$TMortaExterna<-mpio$MortaExterna/mpio$poblacion
mpio@data[is.na(mpio@data)]<-0

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

## Warning in `[<-.factor`(`*tmp*`, thisvar, value = 0): invalid factor level,
## NA generated

table(is.na(mpio@data))

## 
## FALSE  TRUE 
## 46369   975

summary(mpio@data)

##       AREA             PERIMETER         MUNIPIOS_        MUNIPIOS_I    
##  Min.   :1.586e+07   Min.   : 0.1407   Min.   :   2.0   Min.   :   1.0  
##  1st Qu.:1.353e+08   1st Qu.: 0.4790   1st Qu.: 274.8   1st Qu.: 267.8  
##  Median :2.917e+08   Median : 0.7474   Median : 543.5   Median : 529.5  
##  Mean   :1.064e+09   Mean   : 1.1298   Mean   : 543.5   Mean   : 529.6  
##  3rd Qu.:7.049e+08   3rd Qu.: 1.2148   3rd Qu.: 812.2   3rd Qu.: 792.2  
##  Max.   :4.970e+10   Max.   :14.9305   Max.   :1081.0   Max.   :1060.0  
##                                                                         
##             CODIGO         MPIOS            ID             NOMBRE    
##  13667000000000:   2   05001  :   1   216    :  2   BOLIVAR   :   4  
##  05001000000000:   1   05002  :   1   001    :  1   BUENAVISTA:   4  
##  05002000000000:   1   05004  :   1   002    :  1   LA UNION  :   4  
##  05004000000000:   1   05021  :   1   003    :  1   VILLANUEVA:   4  
##  05021000000000:   1   05030  :   1   004    :  1   ARGELIA   :   3  
##  05030000000000:   1   05031  :   1   (Other):209   CORDOBA   :   3  
##  (Other)       :1069   (Other):1070   NA's   :861   (Other)   :1054  
##       DPTO        X_COORD          Y_COORD          cod_mpio    
##  05     :124   Min.   :-78.83   Min.   :-4.190   Min.   : 5001  
##  15     :123   1st Qu.:-75.78   1st Qu.: 4.247   1st Qu.:15629  
##  25     :115   Median :-74.75   Median : 5.500   Median :25802  
##  68     : 87   Mean   :-74.62   Mean   : 5.557   Mean   :38705  
##  52     : 62   3rd Qu.:-73.44   3rd Qu.: 6.925   3rd Qu.:67016  
##  73     : 46   Max.   :-66.94   Max.   :11.740   Max.   :99773  
##  (Other):519                                                    
##                      mpio.dpto         Municipio          Departamento
##  ABEJORRAL-ANTIOQUIA      :   1   BUENAVISTA:   4   ANTIOQUIA   :124  
##  ABREGO-NORTE_DE_SANTANDER:   1   LA UNIÓN :   4   BOYACA      :123  
##  ABRIAQUÃ\215-ANTIOQUIA      :   1   VILLANUEVA:   4   CUNDINAMARCA:115  
##  ACACÃ\215AS-META            :   1   ARGELIA   :   3   SANTANDER   : 87  
##  ACANDÃ\215-CHOCO            :   1   BOLÃ\215VAR  :   3   NARIÑO     : 62  
##  (Other)                  :1052   (Other)   :1039   (Other)     :546  
##  NA's                     :  19   NA's      :  19   NA's        : 19  
##   GrDota                  Categoria.de.ruralidad  ECobeMedNet    
##  C   : 13   Ciudades y aglomeraciones:115        Min.   :0.0000  
##  G1  :212   Intermedios              :298        1st Qu.:0.3000  
##  G2  :211   Rural                    :357        Median :0.4000  
##  G3  :215   Rural disperso           :287        Mean   :0.3999  
##  G4  :209   NA's                     : 19        3rd Qu.:0.5000  
##  G5  :197                                        Max.   :1.0000  
##  NA's: 19                                                        
##    ESABEMATE       ESABELeng      ECoberTrans        Scobert      
##  Min.   : 0.00   Min.   : 0.00   Min.   :0.0000   Min.   :0.0000  
##  1st Qu.:44.58   1st Qu.:48.65   1st Qu.:0.4000   1st Qu.:0.6900  
##  Median :47.90   Median :50.95   Median :0.5100   Median :0.8200  
##  Mean   :46.95   Mean   :49.96   Mean   :0.5121   Mean   :0.7924  
##  3rd Qu.:50.92   3rd Qu.:53.00   3rd Qu.:0.6200   3rd Qu.:0.9500  
##  Max.   :59.60   Max.   :59.09   Max.   :1.0000   Max.   :1.0000  
##                                                                   
##   SVacunPenta       SMortaInf     SPCoberElectrica  SPCoberInter    
##  Min.   :0.0000   Min.   : 0.00   Min.   :0.0000   Min.   :0.00000  
##  1st Qu.:0.8500   1st Qu.:13.69   1st Qu.:0.8700   1st Qu.:0.00900  
##  Median :0.9300   Median :17.05   Median :0.9700   Median :0.02000  
##  Mean   :0.9227   Mean   :18.27   Mean   :0.8769   Mean   :0.03902  
##  3rd Qu.:1.0100   3rd Qu.:21.30   3rd Qu.:0.9900   3rd Qu.:0.04600  
##  Max.   :3.3900   Max.   :64.70   Max.   :1.0000   Max.   :0.68600  
##                                                                     
##   SPCoberAcued    SPCoberAlcanta      SegHurto         SegHomici     
##  Min.   :0.0000   Min.   :0.0000   Min.   :  0.000   Min.   : 0.000  
##  1st Qu.:0.3300   1st Qu.:0.1600   1st Qu.:  3.917   1st Qu.: 0.000  
##  Median :0.5600   Median :0.3300   Median :  8.955   Median : 1.505  
##  Mean   :0.5604   Mean   :0.3818   Mean   : 15.412   Mean   : 2.333  
##  3rd Qu.:0.8100   3rd Qu.:0.5700   3rd Qu.: 19.465   3rd Qu.: 3.230  
##  Max.   :1.0000   Max.   :1.0000   Max.   :157.710   Max.   :27.270  
##                                                                      
##   SegVioIntraF       Ingresos        Desempeno      MortaExterna    
##  Min.   :  0.00   Min.   :     0   Min.   : 0.00   Min.   :   0.00  
##  1st Qu.:  3.54   1st Qu.: 89592   1st Qu.:43.05   1st Qu.:   2.00  
##  Median :  8.16   Median :138929   Median :48.90   Median :   6.00  
##  Mean   : 10.48   Mean   :202347   Mean   :48.36   Mean   :  27.14  
##  3rd Qu.: 14.26   3rd Qu.:254596   3rd Qu.:54.69   3rd Qu.:  16.00  
##  Max.   :106.55   Max.   :468923   Max.   :83.24   Max.   :3048.00  
##                                                                     
##    MortaPeri        MortaCircula     MortaTransmiso     MortaTumores    
##  Min.   :  0.000   Min.   :   0.00   Min.   :   0.00   Min.   :   0.00  
##  1st Qu.:  0.000   1st Qu.:   4.00   1st Qu.:   0.00   1st Qu.:   1.00  
##  Median :  0.000   Median :  10.00   Median :   1.00   Median :   3.00  
##  Mean   :  3.566   Mean   :  57.76   Mean   :  12.91   Mean   :  36.36  
##  3rd Qu.:  0.000   3rd Qu.:  24.00   3rd Qu.:   2.00   3rd Qu.:  10.00  
##  Max.   :672.000   Max.   :9866.00   Max.   :2010.00   Max.   :8397.00  
##                                                                         
##  MortaMalDefini    MortaDemasCaus         DPMP               MPIO     
##  Min.   :  0.000   Min.   :   0.00   Min.   :    0   NA        :   6  
##  1st Qu.:  0.000   1st Qu.:   2.00   1st Qu.:15462   Buenavista:   4  
##  Median :  1.000   Median :   5.00   Median :25595   La Unión :   4  
##  Mean   :  3.494   Mean   :  46.43   Mean   :36852   Villanueva:   4  
##  3rd Qu.:  2.000   3rd Qu.:  11.00   3rd Qu.:63227   Argelia   :   3  
##  Max.   :610.000   Max.   :8519.00   Max.   :99773   (Other)   :1036  
##                                                      NA's      :  19  
##    poblacion       TMortaExterna      
##  Min.   :      0   Min.   :0.0000000  
##  1st Qu.:   6588   1st Qu.:0.0002244  
##  Median :  13302   Median :0.0004554  
##  Mean   :  43045   Mean   :0.0005428  
##  3rd Qu.:  26367   3rd Qu.:0.0007620  
##  Max.   :7674366   Max.   :0.0028307  
## 

table(is.na(mpio@data$TMortaExterna))

## 
## FALSE 
##  1076

table(is.na(mpio@data$SPCoberInter))

## 
## FALSE 
##  1076

# Identify the optimal bandwidth using cross validation; the fixed bandwidth
# Gaussian function is used
bwG <- gwr.sel(TMortaExterna ~ SPCoberInter, data = mpio, gweight = gwr.Gauss, verbose = FALSE)
# View the selected bandwidth
bwG

0.1.2 Ajuste del “modelo” (modelos)

# Fit the GWR model using the bandwidth identified in the previous step
gwrG <- gwr(TMortaExterna ~ SPCoberInter, data = mpio, bandwidth = 50000, gweight = gwr.Gauss)
# Summarise the GWR outputs
summary(gwrG$SDF)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min     max
## x 443770.49 1804280
## y  20565.69 1867933
## Is projected: NA 
## proj4string : [NA]
## Data attributes:
##      sum.w         X.Intercept.         SPCoberInter       
##  Min.   : 1.021   Min.   :-3.186e-05   Min.   :-0.0081847  
##  1st Qu.:23.228   1st Qu.: 3.291e-04   1st Qu.: 0.0003828  
##  Median :36.663   Median : 4.631e-04   Median : 0.0010650  
##  Mean   :41.397   Mean   : 4.837e-04   Mean   : 0.0051739  
##  3rd Qu.:58.467   3rd Qu.: 6.625e-04   3rd Qu.: 0.0023697  
##  Max.   :94.108   Max.   : 1.293e-03   Max.   : 0.3719067  
##      gwr.e                 pred               localR2        
##  Min.   :-7.828e-04   Min.   :-4.704e-07   Min.   :-0.14027  
##  1st Qu.:-2.041e-04   1st Qu.: 3.664e-04   1st Qu.: 0.08322  
##  Median :-3.131e-05   Median : 5.103e-04   Median : 0.11901  
##  Mean   : 1.360e-05   Mean   : 5.292e-04   Mean   : 0.18262  
##  3rd Qu.: 1.616e-04   3rd Qu.: 6.874e-04   3rd Qu.: 0.20070  
##  Max.   : 1.933e-03   Max.   : 2.016e-03   Max.   : 1.00000

Por lo tanto, la distancia de ancho de banda identificada se pasa a bwG y esto se usa con el argumento de ancho de banda usando gwr. En su lugar, podría especificar manualmente un ancho de banda (por ejemplo, ‘ancho de banda = 20000.0’) y evitar el uso de gwr.sel.

Con GWR, se proporcionan coeficientes beta de regresión para cada ubicación de entrada y éstos pueden ser mapeados para que, por ejemplo, se puedan explorar las relaciones locales entre el desempleo y las enfermedades y discapacidades de larga duración.

# Now map the coefficients (the local slope values).
spplot(gwrG$SDF, "SPCoberInter", col = "transparent")

# and the local r square values
spplot(gwrG$SDF, "localR2", col = "transparent")

El término col=“transparent” evita que los bordes se muestren como líneas sólidas distintas.

Recuerde que el cuadrado r indica la bondad del ajuste. Por lo tanto, los valores mapeados cercanos a uno indican las ubicaciones en las que el modelo encaja bien.

El esquema de ponderación (gweight) puede ser cambiado; gwr.Gauss es un núcleo de ancho de banda fijo (por ejemplo, los pesos disminuyen a casi cero después de una distancia fija). En su lugar, puede usar r gwr.bisquare que es un núcleo de ancho variable - es pequeño donde hay una densidad de observaciones alta y grande donde la densidad de observaciones es baja.

Por lo tanto, en este caso, será pequeña en las zonas con mayor densidad de municipiosy grande en el caso contrario. Para aclarar, en el caso de gwr.Gauss se especifica una distancia mientras que para gwr.bisquare es el número de vecinos más cercanos. Por ejemplo, se podrían especificar 25 vecinos más cercanos para la plaza gwr.bis y la distancia de un municipio pequeño (junto a otros pequeños) a su vecino más cercano 25 será obviamente más pequeña que la distancia de un municipio grande (rodeado por otros grandes) rural a su vecino más cercano 25. En R, los vecinos más cercanos se expresan como una fracción de los datos totales (por ejemplo, un ancho de banda de 0,000909090 = unos 1 de 1100 puntos de datos; mientras que 1 representaría un vecindario que contiene todas las demás observaciones).

# Identify the optimal bandwidth using cross validation; an adaptive
# bandwidth bi-square function is used
ad.bw <- gwr.sel(TMortaExterna ~ SPCoberInter, data = mpio, adapt = TRUE)

# Fit the GWR model using the bandwidth identified in the previous step
gwrAD <- gwr(TMortaExterna ~ SPCoberInter, data = mpio, adapt = 0.4, gweight = gwr.bisquare)
# Summarise the GWR outputs
summary(gwrAD$SDF)

## Object of class SpatialPolygonsDataFrame
## Coordinates:
##         min     max
## x 443770.49 1804280
## y  20565.69 1867933
## Is projected: NA 
## proj4string : [NA]
## Data attributes:
##      sum.w         X.Intercept.        SPCoberInter      
##  Min.   : 30.49   Min.   :0.0003047   Min.   :0.0002120  
##  1st Qu.:150.91   1st Qu.:0.0003733   1st Qu.:0.0007381  
##  Median :168.45   Median :0.0004452   Median :0.0011228  
##  Mean   :164.46   Mean   :0.0004944   Mean   :0.0013856  
##  3rd Qu.:181.65   3rd Qu.:0.0006359   3rd Qu.:0.0018568  
##  Max.   :217.84   Max.   :0.0007128   Max.   :0.0045747  
##      gwr.e                 pred              localR2       
##  Min.   :-6.981e-04   Min.   :0.0003076   Min.   :0.02654  
##  1st Qu.:-2.678e-04   1st Qu.:0.0004094   1st Qu.:0.05816  
##  Median :-6.954e-05   Median :0.0005198   Median :0.07574  
##  Mean   : 6.569e-06   Mean   :0.0005362   Mean   :0.09814  
##  3rd Qu.: 1.899e-04   3rd Qu.:0.0006594   3rd Qu.:0.12220  
##  Max.   : 2.199e-03   Max.   :0.0010617   Max.   :0.27936

spplot(gwrAD$SDF, "SPCoberInter", col = "transparent")

spplot(gwrAD$SDF, "localR2", col = "transparent")

0.1.2.1 Ejercicio

Realice una regresión geograficamente ponderada por otras variables de su interés y comente.

0.2 Modelos lineales generalizados

Los modelos lineales generalizados (GLM) son una ampliación flexible de la regresión lineal ordinaria (modelo Normal lineal) que permite variables de respuesta que tienen modelos de distribución de errores distintos de una distribución normal.

0.2.1 GLM

El modelo de regresión normal lineal es muy popular para el análisis estadístico cuando la variable respuesta es continua, simetrica y de de varianza constante. La distribución normal posee propiedades estadísticas deseables y gran cantidad de teoria desarrollada. Sin embargo, cuando la variable respuesta es discreta, o cuando la variable respuesta es continua pero los supuestos de simetria y/o varianza constante no son razonables, el modelo normal lineal puede ser inadecuado. Como alternativa para el análisis surge una extensión del modelo normal lineal, llamada modelo lineal generalizado (MLG), la cual introduce gran flexibilidad tanto en la componente aleatoria como en la sistemática, permitiendo así, abordar una variedad más amplia de situaciones prácticas.

• Componente aleatoria: Identifica la variable respuesta y su distribución de probabilidad.
• Componente sistemática: Especifica las variables explicativas (independientes o predictoras) utilizadas en la función predictora lineal

Al igual que en el modelo de regresión lineal, el supuesto de independencia entre la variable respuesta es fundamental. En el caso de suponer que los datos tienen una estructira de autocorrelación, por ejemplo temporal, es necesario recurrir a la prueba de Durbin Watson, en caso de que se rechase la hipotesis de que los datos esten autocorrelacionados, se procede a estimar el modelo normalmente, en caso contrario es necesario corregir la estructura de autocorrelación incorporandola al modelo.

Los datos con correlación espacial no son la excepción, primero deberemos diagnosticar la presencia de autocorrelaicón esparcial, si esta presente deberá incorporarse al modelo, en caso contrario, el modelo se estimará normalmente.

0.2.2 Ajustado

0.2.3 Sin ajustar (algo falta…)

0.2.4 Residuales aleatorios

0.2.5 Residuales espacialmente correlacionados

0.2.6 Un GLM Poisson

Un modelo lineal generalizado de Poisson es una forma de ajustar los datos de conteo a las variables explicativas. Se obtienen estimaciones de parámetros y errores estándar para las variables explicativas, y se pueden obtener valores ajustados y residuales.

La función glm() permite ajustar GLM con variable respuesta Poisson. Funciona igual que la función lm(), pero también se puede especificar un argumento de familia. La fórmula tiene el significado usual - respuesta a la izquierda del ~, y variables explicativas a la derecha.

Para hacer frente a los datos de conteo procedentes de poblaciones de diferentes tamaños, se especifica un argumento de desplazamiento o compensación (offset). En este caso se emplea el logaritmo de la población como término de compensación. Esto pues el número de casos depentde del tamaño de la población. En otras palabras en el offset del modelo se colocan los denominadores para modelar corregir el hecho de no estar modelando la tasa.

london<-readRDS("D:/DriveW7/2019/Sabana/Estadistica aplicada Salud Publica/Bibliografía/DataCamp Spatial Statistics in R/datasets/london_2017_2.rds", refhook = NULL)
# Fit a poisson GLM.
model_flu <- glm(
  Flu_OBS ~ HealthDeprivation, 
  offset = log(TOTAL_POP), 
  data = london, 
  family = poisson)

# Is HealthDeprivation significant?
summary(model_flu)

## 
## Call:
## glm(formula = Flu_OBS ~ HealthDeprivation, family = poisson, 
##     data = london, offset = log(TOTAL_POP))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -9.5361  -4.5285  -0.0499   2.9043   8.2194  
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)       -8.78190    0.02869 -306.043   <2e-16 ***
## HealthDeprivation  0.65689    0.06797    9.665   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 703.75  on 31  degrees of freedom
## Residual deviance: 605.03  on 30  degrees of freedom
## AIC: 762.37
## 
## Number of Fisher Scoring iterations: 5

# Put residuals into the spatial data.
london$Flu_Resid <- residuals(model_flu)

# Map the residuals using spplot
spplot(london, "Flu_Resid")

Las regiones en amarillo tenía muchos más casos de gripe de lo esperado; las que se encuentran en negro tenía muchos menos casos.

0.2.7 Residuos

Un modelo lineal debe ajustarse a los datos y dejar residuos no correlacionados. Esto se aplica a los modelos no espaciales, en los que, por ejemplo, ajustar una línea recta a través de puntos en una curva daría lugar a residuos correlacionados en serie. Un modelo sobre datos espaciales debe tener como objetivo tener residuos que no muestren una correlación espacial significativa.

Puede probar el modelo ajustado a los datos de la gripe usando moran.mc() del paquete spdep.

# Compute the neighborhood structure.
library(spdep)

## Loading required package: Matrix

borough_nb <- poly2nb(london)

# Test spatial correlation of the residuals.
moran.mc(london$Flu_Resid, listw = nb2listw(borough_nb), nsim = 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  london$Flu_Resid 
## weights: nb2listw(borough_nb)  
## number of simulations + 1: 1000 
## 
## statistic = 0.15059, observed rank = 921, p-value = 0.079
## alternative hypothesis: greater

mp <- moran.plot(london$Flu_Resid, nb2listw(borough_nb), labels = as.character(london$NAME), 
    xlab = "residuos", ylab = "residuos rezagado")

** Es deseable un valor p de menos de 0.05 para considerar que los residuos están espacialmente autocorrelacionados. En este caso, no hay suficiente evidencia en la prueba para tomar una decisión clara.**

0.3 Correlación en GLM espaciales

0.3.1 Faltan incluir variables

0.3.2 Modelo para correlación espacial

Agregar un termino espacial al modelo.

\[Y=X\beta+S(x,y)\]

0.3.3 Autocorrelación Condicional

Hay varias formas de definir \(S(x,y)\), una es suponer que los datos de una región se distribuyen normalmente con media los vecinos y varianza por ser estimada. Es un modelo sencillo, con las ventajas de la distribución normal y solo un parámetro extra.

0.3.4 Ajustar un GLM Bayesiano

Se usan los metodos bayesianos dado lo costoso que es calcular la verosimilitud del modelo convencional.

La estimación bayesiana estima la función de distribución a posteriori del parámetro dada la muestra (evidencia), con la función de distribución se calculan los intervalos de credibilidad para el parámetro.

Si el intervalo de credibilidad no contiene el cero el parámetro es significativo para el modelo (análogo a la prueba t del modelo de regresión sobre cada beta distinto de cero)

Los modelos estadísticos bayesianos devuelven muestras de los parámetros de interés (la distribución “posterior”) basadas en alguna distribución “a priori”, supuesta para el parámetro, distribución que se actualiza con los datos. El proceso de modelado bayesiano devuelve un número de muestras a partir de las cuales se puede calcular la media, o una probabilidad de superación, o cualquier otra cantidad que se pueda calcular a partir de una distribución.

Antes de ajustar un modelo con correlación espacial, primero se ajustará el mismo modelo que antes, pero utilizando la inferencia bayesiana.

Repasar Clase Bayesiana

# Use R2BayesX
# install.packages("R2BayesX")
library(R2BayesX)

## Loading required package: BayesXsrc

## Loading required package: colorspace

## Loading required package: mgcv

## Loading required package: nlme

## This is mgcv 1.8-27. For overview type 'help("mgcv-package")'.

# Fit a GLM
model_flu <- glm(Flu_OBS ~ HealthDeprivation, offset = log(TOTAL_POP),
                data = london, family = poisson)
                    
# Summarize it                    
summary(model_flu)

## 
## Call:
## glm(formula = Flu_OBS ~ HealthDeprivation, family = poisson, 
##     data = london, offset = log(TOTAL_POP))
## 
## Deviance Residuals: 
##     Min       1Q   Median       3Q      Max  
## -9.5361  -4.5285  -0.0499   2.9043   8.2194  
## 
## Coefficients:
##                   Estimate Std. Error  z value Pr(>|z|)    
## (Intercept)       -8.78190    0.02869 -306.043   <2e-16 ***
## HealthDeprivation  0.65689    0.06797    9.665   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 703.75  on 31  degrees of freedom
## Residual deviance: 605.03  on 30  degrees of freedom
## AIC: 762.37
## 
## Number of Fisher Scoring iterations: 5

# Calculate coeff confidence intervals
confint(model_flu)

## Waiting for profiling to be done...

##                       2.5 %     97.5 %
## (Intercept)       -8.838677 -8.7261843
## HealthDeprivation  0.524437  0.7908841

# Fit a Bayesian GLM
bayes_flu <- bayesx(Flu_OBS ~ HealthDeprivation, offset = log(london$TOTAL_POP), 
                    family = "poisson", data = data.frame(london), 
                    control = bayesx.control(seed = 17610407))
                    
# Summarize it                    
summary(bayes_flu)

## Call:
## bayesx(formula = Flu_OBS ~ HealthDeprivation, data = data.frame(london), 
##     offset = log(london$TOTAL_POP), control = bayesx.control(seed = 17610407), 
##     family = "poisson")
##  
## Fixed effects estimation results:
## 
## Parametric coefficients:
##                      Mean      Sd    2.5%     50%   97.5%
## (Intercept)       -8.7823  0.0278 -8.8363 -8.7828 -8.7271
## HealthDeprivation  0.6556  0.0686  0.5264  0.6578  0.7838
##  
## N = 32  burnin = 2000  method = MCMC  family = poisson  
## iterations = 12000  step = 10

# Look at the samples from the Bayesian model
plot(samples(bayes_flu))

Los coeficientes y sus intervalos de confianza son casi idénticos para los dos modelos. Esto es lo que se puede esperar, ya que en cada caso se utilizó la misma estructura de modelo. A continuación, verá cómo ampliar la versión bayesiana para incluir la correlación espacial.

0.3.5 Añadir un efecto espacialmente autocorrelacionado

Hemos ajustado un GLM no espacial usando el paquete R2BayesX. Ahora se incluirá un término correlacionado espacialmente basado en la estructura de adyacencia añadiendo un término a la fórmula que especifique un modelo correlacionado espacialmente.

# Compute adjacency objects
borough_nb <- poly2nb(london)
borough_gra <- nb2gra(borough_nb)

# Fit spatial model
flu_spatial <- bayesx(
  Flu_OBS ~ HealthDeprivation + sx(i, bs = "spatial", map = borough_gra),
  offset = log(london$TOTAL_POP),
  family = "poisson", data = data.frame(london), 
  control = bayesx.control(seed = 17610407)
)

## Note: created new output directory 'C:/Users/Daniel/AppData/Local/Temp/RtmpYh5AZ6/bayesx1'!

# Summarize the model
summary(flu_spatial)

## Call:
## bayesx(formula = Flu_OBS ~ HealthDeprivation + sx(i, bs = "spatial", 
##     map = borough_gra), data = data.frame(london), offset = log(london$TOTAL_POP), 
##     control = bayesx.control(seed = 17610407), family = "poisson")
##  
## Fixed effects estimation results:
## 
## Parametric coefficients:
##                      Mean      Sd    2.5%     50%   97.5%
## (Intercept)       -9.1773  0.1244 -9.4220 -9.1692 -8.9655
## HealthDeprivation  1.1022  0.6408 -0.1566  1.1729  2.1644
## 
## Smooth terms variances:
##             Mean     Sd   2.5%    50%  97.5%    Min    Max
## sx(i):mrf 4.9469 1.8225 2.2798 4.6614 9.3513 1.4965 16.746
##  
## N = 32  burnin = 2000  method = MCMC  family = poisson  
## iterations = 12000  step = 10

0.3.6 Mapeo de los efectos espaciales

Al igual que con glm, puede obtener los valores de ajuste y los residuos de su modelo utilizando las funciones fitted y residuals. Los modelos bayesx son un poco más complejos, ya que tiene el predictor lineal y los términos de los elementos \(sx\), como el término espacialmente correlacionado.

La función de resumen le mostrará información sobre los términos del modelo lineal y los términos de suavizado en dos tablas separadas. El término espacial se llama “sx(i):mrf”, que significa “Markov Random Field”.

El análisis bayesiano devuelve muestras de una distribución para nuestro término S(x) en cada uno de los municipios de Londres. La función de ajuste de modelos mediante bayesx proporciona un resumen estadístico de cada municipio.

# Summarize the model
summary(flu_spatial)

## Call:
## bayesx(formula = Flu_OBS ~ HealthDeprivation + sx(i, bs = "spatial", 
##     map = borough_gra), data = data.frame(london), offset = log(london$TOTAL_POP), 
##     control = bayesx.control(seed = 17610407), family = "poisson")
##  
## Fixed effects estimation results:
## 
## Parametric coefficients:
##                      Mean      Sd    2.5%     50%   97.5%
## (Intercept)       -9.1773  0.1244 -9.4220 -9.1692 -8.9655
## HealthDeprivation  1.1022  0.6408 -0.1566  1.1729  2.1644
## 
## Smooth terms variances:
##             Mean     Sd   2.5%    50%  97.5%    Min    Max
## sx(i):mrf 4.9469 1.8225 2.2798 4.6614 9.3513 1.4965 16.746
##  
## N = 32  burnin = 2000  method = MCMC  family = poisson  
## iterations = 12000  step = 10

# Map the fitted spatial term only
london$spatial <- fitted(flu_spatial, term = "sx(i):mrf")[, "Mean"]
spplot(london, zcol = "spatial")

# Map the residuals
london$spatial_resid <- residuals(flu_spatial)[, "mu"]
spplot(london, zcol = "spatial_resid")

# Test residuals for spatial correlation
moran.mc(london$spatial_resid, nb2listw(borough_nb), 999)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  london$spatial_resid 
## weights: nb2listw(borough_nb)  
## number of simulations + 1: 1000 
## 
## statistic = -0.27452, observed rank = 16, p-value = 0.984
## alternative hypothesis: greater

Con el p-valor no se rechaza la hipotesis nula de aleatoriedad espacial completa, con lo cual el modelo especifica la componente espacial.

library(spdep)
mp <- moran.plot(london$spatial_resid, nb2listw(borough_nb), labels = as.character(london$NAME), 
    xlab = "residuos", ylab = "residuos rezagado")

0.3.6.1 Ejercicio

Ajuste un GLM bayesiano con el dataset anterior y variables de su interes.

mpio@data$poblacion[mpio@data$poblacion==0]<-1
summary(mpio@data$poblacion)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       1    6588   13302   43045   26367 7674366

#str(mpio@data)
# Fit a poisson GLM.
model_flu <- glm(
  MortaExterna ~ SPCoberInter,
  offset = log(poblacion), 
  data = mpio@data, 
  family = poisson)

# Is HealthDeprivation significant?
summary(model_flu)

## 
## Call:
## glm(formula = MortaExterna ~ SPCoberInter, family = poisson, 
##     data = mpio@data, offset = log(poblacion))
## 
## Deviance Residuals: 
##      Min        1Q    Median        3Q       Max  
## -29.9584   -1.9278   -0.6338    0.7139   30.8163  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)    
## (Intercept)  -7.42020    0.01074 -691.02  < 2e-16 ***
## SPCoberInter  0.38896    0.06777    5.74 9.48e-09 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## (Dispersion parameter for poisson family taken to be 1)
## 
##     Null deviance: 8882.5  on 1075  degrees of freedom
## Residual deviance: 8849.5  on 1074  degrees of freedom
## AIC: 12640
## 
## Number of Fisher Scoring iterations: 5

# Put residuals into the spatial data.
mpio$Flu_Resid <- residuals(model_flu)

# Map the residuals using spplot
spplot(mpio, "Flu_Resid", col = "transparent")

# Compute the neighborhood structure.
library(spdep)
borough_nb <- poly2nb(mpio)

# Test spatial correlation of the residuals.
moran.mc(mpio$Flu_Resid, listw = nb2listw(borough_nb), nsim = 100)

## 
##  Monte-Carlo simulation of Moran I
## 
## data:  mpio$Flu_Resid 
## weights: nb2listw(borough_nb)  
## number of simulations + 1: 101 
## 
## statistic = 0.2974, observed rank = 101, p-value = 0.009901
## alternative hypothesis: greater

mp <- moran.plot(mpio$Flu_Resid, nb2listw(borough_nb), labels = as.character(london$NAME), 
    xlab = "residuos", ylab = "residuos rezagado")

Se rechaza la nula, por lo cual existe autocorrelación y es necesario incorporar dicha estructura al modelo.

“Statistics, the science of uncertainty, attempts to model order in disorder” - Noel Cressie