Capítulo 6 Modelos de regresión espacial

6.1 Estudio de Mercadeo

Se comparan varios tipos de modelos de regresión espacial para ver con cual se obtiene el mejor ajuste. Se consideran modelos autoregresivos y de medias móvviles así como su combinación.

6.2 Paquetes

rm(list=ls())
library(openxlsx)
library(dplyr)
library(rgdal)
library(maptools)
library(GISTools)
library(spdep)
library(readr)
library(car)
library(readxl)
library(psych)
library(rgdal)
library(FactoClass)
library(spdep)
require("GWmodel")
library("mapsRinteractive")
options(scipen = 999)

6.3 Lectura de Datos

# Lectura de Datos
BASE <- read_excel("data_3_EstudioDeMercadoEspacial/BASE.xlsx")
# Lectura del Shape de Colombia por Departamentos
Colombia = readOGR(dsn = "data_3_EstudioDeMercadoEspacial/Geodatabase Colombia", 
                   layer = "departamentos")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/jncc/Documents/Monitorias/EspacialPage/Clases-EE-UN/data_3_EstudioDeMercadoEspacial/Geodatabase Colombia", layer: "departamentos"
## with 33 features
## It has 6 fields
## Integer64 fields read as strings:  AñO_CREAC

6.3.1 Cruce de información y arreglo de coordenadas

#Cruce de información con el shape cargado
Insumo = merge(Colombia, BASE, by.x="COD_DANE", by.y="Cod")
Insumo = subset(Insumo[c(1:31,33),])
# Conversión a Coordenadas UTM
Crs.geo = CRS("+proj=tmerc +lat_0=4.599047222222222 +lon_0=-74.08091666666667 +k=1 +x_0=1000000 +y_0=1000000 +ellps=intl +towgs84=307,304,-318,0,0,0,0 +units=m +no_defs")
proj4string(Insumo) <- Crs.geo 
Insumo.utm = spTransform(Insumo, CRS("+init=epsg:3724 +units=km"))
## Warning: PROJ support is provided by the sf and terra packages among others

6.4 Matriz de vecindades

#---
# MATRIZ DE VECINDADES (W)
#---
## Centroides de las Áreas
Centros = getSpPPolygonsLabptSlots(Insumo.utm)
## Warning: use coordinates method
Centroids <- SpatialPointsDataFrame(coords = Centros, 
                                    data=Insumo.utm@data, 
                                    proj4string=CRS("+init=epsg:3724 +units=km"))
# Matriz de Distancias entre los Centriodes
Wdist = dist(Centros, up=T)
# Matriz W de vecindades
library(pgirmess)
library(HistogramTools)
library(strucchange)
library(spdep)
Insumo.nb = poly2nb(Insumo.utm, queen=T)
#n <- max(sapply(Insumo.nb, length))
#ll <- lapply(Insumo.nb, function(X) {
#  c(as.numeric(X), rep(0, times = n - length(X)))
#})
#out <- do.call(cbind, ll)
#Departamentos<-Insumo$Departamento
#MatW<-matrix(NA,32,32)
#for (i in 1:8) {
#  for (j in 1:32) {
#    if (out[i,j]!=0) {
#      MatW[out[i,j],j]<-1
#    } else{MatW[out[i,j],j]<-0}
#  }
#}
#for (i in 1:32) {
#  for (j in 1:32) {
#    if (is.na(MatW[i,j])) {
#      MatW[i,j]<-0
#    }
#  }
#}
#colnames(MatW)<-Departamentos
#rownames(MatW)<-Departamentos
#MatW1<-MatW[,1:16]
#MatW2<-MatW[,17:32]
# Martiz W (Estilos)
Insumo.lw = nb2listw(Insumo.nb)
Insumo.lwb = nb2listw(Insumo.nb, style="B")
Insumo.lwc = nb2listw(Insumo.nb, style="C")
Insumo.lwu = nb2listw(Insumo.nb, style="U")
Insumo.lww = nb2listw(Insumo.nb, style="W")

6.5 Mapa de valores observados

#  Mapa de Valores Observados
#dev.new() #windows()
choropleth(Insumo, Insumo$CAP_BAC)
shad = auto.shading(Insumo$CAP_BAC, 
                    n=5, 
                    cols=(brewer.pal(5,"Reds")), 
                    cutter = quantileCuts)

choro.legend(1555874,535165.5, 
             shad, 
             fmt="%1.1f", 
             title = "Valores Locales", 
             cex=0.7, 
             under = "Menos de", 
             between = "a", 
             over = "Mas de")

title("Valores Observados para las captaciones del banco agrario 
       en Colombia, cuarto trimestre 2020", cex.main=1)
map.scale(755874,335165.5, 250000, "km", 2, 50, sfcol='brown')

6.6 Pruebas de Autocorrelación

#----------------------------
#  PRUEBAS DE AUTOCORRELACION
#----------------------------

# Moran
moran.test(Insumo$CAP_BAC, Insumo.lw)
## 
##  Moran I test under randomisation
## 
## data:  Insumo$CAP_BAC  
## weights: Insumo.lw    
## 
## Moran I statistic standard deviate = 2.0024, p-value = 0.02262
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##       0.153081266      -0.032258065       0.008566935
# Dispersograma de Moran
#dev.new() #windows()
moran.plot(Insumo$CAP_BAC, 
           Insumo.lw, 
           labels=as.character(Insumo$Departamento), 
           xlab="Captaciones BAC", 
           ylab="Captaciones BAC rezagado", 
           las=1, 
           pch=16, 
           cex=0.5)

legend("bottomright", 
       legend=c("I de Moran: 0.1530", "Valor P:      0.02262"), 
       cex=1,
       bg='lightgreen')

title("Dispersograma de Moran para las captaciones del banco agrario en 
los Departamentos de Colombia, cuarto trimestre 2020", cex.main=1)

# Local G
nearng = dnearneigh(coordinates(Insumo.utm), 0, 550)
Insumo.lw.g = nb2listw(nearng, style="B")

localG = localG(Insumo$CAP_BAC, Insumo.lw.g); localG
##  [1]  1.66525050  0.02608278  1.33032949  1.15892050  1.85852161  0.68445519  1.49486468  0.10163662
##  [9]  1.64717068  1.32714028  2.53361281  2.21899396 -0.71311540  0.50952811  1.48899277  0.81676480
## [17]  0.98434717  2.09087184  2.08725553  1.05493906  1.32486118  2.09147517  2.16305539  1.89323276
## [25]  1.52155929  0.84992902 -1.19798594 -1.33847805  0.29701426 -1.60300117  1.67015910  1.96543367
## attr(,"internals")
##               Gi     E(Gi)       V(Gi)       Z(Gi) Pr(z != E(Gi))
##  [1,] 0.87682996 0.7096774 0.010075505  1.66525050     0.09586279
##  [2,] 0.29285361 0.2903226 0.009416471  0.02608278     0.97919132
##  [3,] 0.75219542 0.6129032 0.010963144  1.33032949     0.18340974
##  [4,] 0.82743583 0.7096774 0.010324679  1.15892050     0.24648859
##  [5,] 0.94189020 0.7741935 0.008141671  1.85852161     0.06309498
##  [6,] 0.55680866 0.4838710 0.011355702  0.68445519     0.49368778
##  [7,] 0.64933079 0.4838710 0.012251276  1.49486468     0.13494976
##  [8,] 0.49483872 0.4838710 0.011644875  0.10163662     0.91904511
##  [9,] 0.75532158 0.5806452 0.011245826  1.64717068     0.09952297
## [10,] 0.89821754 0.7741935 0.008733287  1.32714028     0.18446228
## [11,] 0.87606260 0.6129032 0.010788404  2.53361281     0.01128934
## [12,] 0.82054359 0.5806452 0.011688064  2.21899396     0.02648714
## [13,] 0.19194463 0.2580645 0.008596952 -0.71311540     0.47577435
## [14,] 0.37327079 0.3225806 0.009897162  0.50952811     0.61038210
## [15,] 0.85695474 0.7096774 0.009783327  1.48899277     0.13648927
## [16,] 0.44117384 0.3548387 0.011173289  0.81676480     0.41406285
## [17,] 0.77877751 0.6774194 0.010602806  0.98434717     0.32494485
## [18,] 0.91159284 0.7096774 0.009325758  2.09087184     0.03653955
## [19,] 0.91151422 0.7096774 0.009350815  2.08725553     0.03686504
## [20,] 0.84508265 0.7419355 0.009560043  1.05493906     0.29145320
## [21,] 0.68979418 0.5483871 0.011392042  1.32486118     0.18521720
## [22,] 0.89320523 0.6774194 0.010644875  2.09147517     0.03648549
## [23,] 0.81503466 0.5806452 0.011741972  2.16305539     0.03053692
## [24,] 0.68555602 0.4838710 0.011348525  1.89323276     0.05832692
## [25,] 0.77320627 0.6129032 0.011099561  1.52155929     0.12811954
## [26,] 0.54263403 0.4516129 0.011468829  0.84992902     0.39536455
## [27,] 0.04481324 0.1290323 0.004942162 -1.19798594     0.23092249
## [28,] 0.08047617 0.1935484 0.007136566 -1.33847805     0.18074065
## [29,] 0.64349326 0.6129032 0.010607310  0.29701426     0.76645562
## [30,] 0.08284688 0.2258065 0.007953509 -1.60300117     0.10893440
## [31,] 0.52325216 0.3548387 0.010168054  1.67015910     0.09488789
## [32,] 0.90927502 0.7741935 0.004723618  1.96543367     0.04936407
## attr(,"cluster")
##  [1] High Low  Low  High Low  Low  High Low  Low  High Low  High Low  Low  Low  High High Low  Low  High Low 
## [22] High High Low  Low  Low  Low  Low  Low  Low  Low  High
## Levels: Low High
## attr(,"gstari")
## [1] FALSE
## attr(,"call")
## localG(x = Insumo$CAP_BAC, listw = Insumo.lw.g)
## attr(,"class")
## [1] "localG"
# Simulaci?n montecarlo
sim.G = matrix(0,1000,32)
for(i in 1:1000) sim.G[i,] = localG(sample(Insumo$CAP_BAC),Insumo.lw.g)
mc.pvalor.G = (colSums(sweep(sim.G,2,localG,">="))+1)/(nrow(sim.G)+1)
mc.pvalor.G
##  [1] 0.015984016 0.453546454 0.086913087 0.083916084 0.001998002 0.278721279 0.063936064 0.465534466
##  [9] 0.048951049 0.055944056 0.000999001 0.003996004 0.740259740 0.315684316 0.026973027 0.232767233
## [17] 0.180819181 0.000999001 0.000999001 0.145854146 0.092907093 0.000999001 0.004995005 0.024975025
## [25] 0.043956044 0.236763237 0.967032967 0.963036963 0.412587413 0.995004995 0.057942058 0.000999001

6.7 Mapas

# Mapas
par(mfrow=c(1,2), mar=c(1,1,8,1)/2)
shadeg = auto.shading(localG, 
                      n=5, 
                      cols=(brewer.pal(5,"Purples")), 
                      cutter=quantileCuts)
#dev.new() #windows()
choropleth(Insumo, 
           localG, 
           shading=shadeg)

choro.legend(1555874,
             535165.5, 
             shadeg, 
             fmt="%1.2f", 
             title = "G", 
             cex=0.7, 
             under = "Menos de", 
             between = "a", 
             over = "Mas de")

title("G Getis Ord Local para las captaciones del banco agrario 
       en Colombia, cuarto trimestre 2020", cex.main=1)

map.scale(755874,335165.5, 250000, "km", 2, 50, sfcol='brown')

# Mapa de P-values
#dev.new() #windows()
shadegp = shading(c(0.01,0.05,0.1), cols = (brewer.pal(4,"Spectral")))
choropleth(Insumo, mc.pvalor.G, shading=shadegp)
choro.legend(1555874,
             535165.5, 
             shadegp, 
             fmt="%1.2f", 
             title = "P-valor de G", 
             cex=0.7, 
             under = "Menos de", 
             between = "a", 
             over = "Mas de")

title("P- Valor de G Getis Ord Local para las captaciones del banco agrario 
       en Colombia, cuarto trimestre 2020", cex.main=1)

map.scale(755874,335165.5, 250000, "km", 2, 50, sfcol='brown')

##Modelos SDEM, SDM, Manski, SARAR

####Modelos SDEM, SDM, Manski, SARAR########
#reg.eq1=CAP_BAC ~ PIB + NBI + CAP_BOG + CAP_BC + CAP_OCC + CAP_CS + Población + IPM
reg.eq1=CAP_BAC ~ PIB + NBI + CAP_BOG+CAP_BC + CAP_OCC + CAP_CS+ Población
reg1=lm(reg.eq1,data=Insumo)                                     #OLS            y=XB+e,    
reg2=lmSLX(reg.eq1,data=Insumo, Insumo.lw)                       #SLX            y=XB+WxT+e
reg3=lagsarlm(reg.eq1,data= Insumo, Insumo.lw)                   #Lag Y          y=XB+WxT+u,   u=LWu+e
## Warning in lagsarlm(reg.eq1, data = Insumo, Insumo.lw): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 5.81787e-18 - using numerical Hessian.
## Warning in sqrt(diag(fdHess)[-1]): NaNs produced
reg4=errorsarlm(reg.eq1,data=Insumo, Insumo.lw)                  #Spatial Error  y=pWy+XB+e   
## Warning in errorsarlm(reg.eq1, data = Insumo, Insumo.lw): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 9.23672e-18 - using numerical Hessian.
reg5=errorsarlm(reg.eq1, data=Insumo, Insumo.lw, etype="emixed") #SDEM Spatial Durbin Error Model y=XB+WxT+u,   u=LWu+e
## Warning in errorsarlm(reg.eq1, data = Insumo, Insumo.lw, etype = "emixed"): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 6.96475e-18 - using numerical Hessian.
reg6=lagsarlm(reg.eq1, data=Insumo,Insumo.lw, type="mixed")      #SDM Spatial Durbin Model (add lag X to SAR) y=pWy+XB+WXT+e 
## Warning in lagsarlm(reg.eq1, data = Insumo, Insumo.lw, type = "mixed"): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 4.67991e-18 - using numerical Hessian.
reg7=sacsarlm(reg.eq1,data=Insumo, Insumo.lw, type="sacmixed")   #Manski Model: y=pWy+XB+WXT+u,   u=LWu+e (no recomendado)
## Warning in sacsarlm(reg.eq1, data = Insumo, Insumo.lw, type = "sacmixed"): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 1.95939e-18 - using numerical Hessian.
reg8=sacsarlm(reg.eq1,data=Insumo,Insumo.lw, type="sac")         #SARAR o Kelejian-Prucha, Cliff-Ord, o SAC If all T=0,y=pWy+XB+u, u=LWu+e
## Warning in sacsarlm(reg.eq1, data = Insumo, Insumo.lw, type = "sac"): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 8.69628e-18 - using numerical Hessian.

6.8 Resumen de modelos

#Resumen de modelos
s=summary
s(reg1)#OLS
## 
## Call:
## lm(formula = reg.eq1, data = Insumo)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -276.51  -65.60   -7.76   46.60  396.20 
## 
## Coefficients:
##                 Estimate   Std. Error t value Pr(>|t|)  
## (Intercept) 148.21380364  79.80068638   1.857   0.0756 .
## PIB           0.00389642   0.00328986   1.184   0.2479  
## NBI          -1.28539812   1.73982368  -0.739   0.4672  
## CAP_BOG      -0.06643826   0.05411306  -1.228   0.2314  
## CAP_BC        0.00397406   0.00607852   0.654   0.5195  
## CAP_OCC      -0.04340185   0.02170799  -1.999   0.0570 .
## CAP_CS        0.47283237   0.31370238   1.507   0.1448  
## Población     0.00000137   0.00006700   0.020   0.9839  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 141.6 on 24 degrees of freedom
## Multiple R-squared:  0.8807, Adjusted R-squared:  0.8459 
## F-statistic: 25.31 on 7 and 24 DF,  p-value: 0.000000001309
s(reg2)#SLX
## 
## Call:
## lm(formula = formula(paste("y ~ ", paste(colnames(x)[-1], collapse = "+"))), 
##     data = as.data.frame(x), weights = weights)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -201.00  -74.99   -0.51   34.55  342.00 
## 
## Coefficients:
##                   Estimate   Std. Error t value Pr(>|t|)
## (Intercept)   392.97159157 284.42668099   1.382    0.185
## PIB             0.00429319   0.00383973   1.118    0.279
## NBI            -0.57603773   2.48213535  -0.232    0.819
## CAP_BOG        -0.00703635   0.07700761  -0.091    0.928
## CAP_BC         -0.00075536   0.01003615  -0.075    0.941
## CAP_OCC        -0.05372016   0.03418522  -1.571    0.135
## CAP_CS          0.12199794   0.43670734   0.279    0.783
## Población       0.00004753   0.00012149   0.391    0.701
## lag.PIB         0.00317318   0.00943185   0.336    0.741
## lag.NBI        -6.81433196   5.99552430  -1.137    0.271
## lag.CAP_BOG    -0.06868663   0.18842513  -0.365    0.720
## lag.CAP_BC      0.00684589   0.01472675   0.465    0.648
## lag.CAP_OCC    -0.00984285   0.05447122  -0.181    0.859
## lag.CAP_CS      0.33585267   1.09533771   0.307    0.763
## lag.Población  -0.00016927   0.00017459  -0.970    0.346
## 
## Residual standard error: 151.6 on 17 degrees of freedom
## Multiple R-squared:  0.9031, Adjusted R-squared:  0.8234 
## F-statistic: 11.32 on 14 and 17 DF,  p-value: 0.000005467
s(reg3)#Lag Y
## 
## Call:lagsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -213.859  -60.238  -17.811   42.960  393.389 
## 
## Type: lag 
## Coefficients: (numerical Hessian approximate standard errors) 
##                   Estimate     Std. Error z value Pr(>|z|)
## (Intercept) 77.75432883510 72.78425608641  1.0683  0.28539
## PIB          0.00330780489  0.00158809777  2.0829  0.03726
## NBI         -0.77863477697  1.47850389439 -0.5266  0.59845
## CAP_BOG     -0.05466122014  0.04390266160 -1.2451  0.21311
## CAP_BC       0.00533211182  0.00466511800  1.1430  0.25305
## CAP_OCC     -0.03421667889  0.01727641951 -1.9805  0.04764
## CAP_CS       0.40699817626  0.26143664946  1.5568  0.11952
## Población    0.00000091508            NaN     NaN      NaN
## 
## Rho: 0.22884, LR test value: 2.096, p-value: 0.14768
## Approximate (numerical Hessian) standard error: 0.15349
##     z-value: 1.4909, p-value: 0.13599
## Wald statistic: 2.2227, p-value: 0.13599
## 
## Log likelihood: -198.251 for lag model
## ML residual variance (sigma squared): 13916, (sigma: 117.97)
## Number of observations: 32 
## Number of parameters estimated: 10 
## AIC: 416.5, (AIC for lm: 416.6)
s(reg4)#Lag Error (SEM)
## 
## Call:errorsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -215.7428  -56.4299   -2.2091   46.6630  425.7867 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value Pr(>|z|)
## (Intercept) 94.500065984 73.097738889  1.2928  0.19608
## PIB          0.003656790  0.002619923  1.3958  0.16279
## NBI         -0.334674352  1.408551017 -0.2376  0.81219
## CAP_BOG     -0.032484043  0.049990128 -0.6498  0.51582
## CAP_BC       0.001202333  0.005555428  0.2164  0.82866
## CAP_OCC     -0.048739699  0.019538015 -2.4946  0.01261
## CAP_CS       0.285572487  0.287233644  0.9942  0.32012
## Población    0.000040164  0.000066648  0.6026  0.54676
## 
## Lambda: 0.50692, LR test value: 3.8031, p-value: 0.051158
## Approximate (numerical Hessian) standard error: 0.21495
##     z-value: 2.3584, p-value: 0.018356
## Wald statistic: 5.5618, p-value: 0.018356
## 
## Log likelihood: -197.3975 for error model
## ML residual variance (sigma squared): 12495, (sigma: 111.78)
## Number of observations: 32 
## Number of parameters estimated: 10 
## AIC: NA (not available for weighted model), (AIC for lm: 416.6)
s(reg5)#Durbin Error (SDEM)
## 
## Call:errorsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw,     etype = "emixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -199.7723  -69.8844   -2.0075   37.9235  362.9104 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                    Estimate    Std. Error z value Pr(>|z|)
## (Intercept)   248.934472800 220.997187634  1.1264  0.25999
## PIB             0.004499050   0.002763860  1.6278  0.10356
## NBI            -0.058410794   1.767448316 -0.0330  0.97364
## CAP_BOG        -0.026724646   0.054644416 -0.4891  0.62480
## CAP_BC         -0.000183348   0.007104434 -0.0258  0.97941
## CAP_OCC        -0.051006011   0.024157938 -2.1114  0.03474
## CAP_CS          0.227498823   0.309055945  0.7361  0.46166
## Población       0.000035177   0.000084525  0.4162  0.67729
## lag.PIB         0.003069904   0.007273446  0.4221  0.67297
## lag.NBI        -4.397300877   4.370216215 -1.0062  0.31432
## lag.CAP_BOG    -0.078935307   0.133353631 -0.5919  0.55390
## lag.CAP_BC      0.000970529   0.011031859  0.0880  0.92990
## lag.CAP_OCC    -0.013417110   0.040319111 -0.3328  0.73931
## lag.CAP_CS      0.374124930   0.777502425  0.4812  0.63038
## lag.Población  -0.000103911   0.000125213 -0.8299  0.40661
## 
## Lambda: 0.28217, LR test value: 0.34013, p-value: 0.55976
## Approximate (numerical Hessian) standard error: 0.45782
##     z-value: 0.61635, p-value: 0.53767
## Wald statistic: 0.37988, p-value: 0.53767
## 
## Log likelihood: -195.7931 for error model
## ML residual variance (sigma squared): 11856, (sigma: 108.89)
## Number of observations: 32 
## Number of parameters estimated: 17 
## AIC: NA (not available for weighted model), (AIC for lm: 423.93)
s(reg6)#Durbin (SDM)
## 
## Call:lagsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw,     type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -185.6170  -74.5630   -1.2308   34.6848  363.6597 
## 
## Type: mixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                    Estimate    Std. Error z value Pr(>|z|)
## (Intercept)   259.002885782 224.615320155  1.1531  0.24887
## PIB             0.004292096   0.002708060  1.5849  0.11298
## NBI            -0.068020238   2.132586505 -0.0319  0.97456
## CAP_BOG        -0.013498835   0.053903632 -0.2504  0.80226
## CAP_BC         -0.000674466   0.006222482 -0.1084  0.91368
## CAP_OCC        -0.052995708   0.022652248 -2.3395  0.01931
## CAP_CS          0.159111970   0.302853078  0.5254  0.59932
## Población       0.000046528   0.000082317  0.5652  0.57191
## lag.PIB         0.001006279   0.007240728  0.1390  0.88947
## lag.NBI        -5.022145755   4.587135513 -1.0948  0.27359
## lag.CAP_BOG    -0.060141055   0.131444000 -0.4575  0.64728
## lag.CAP_BC      0.003372687   0.010959544  0.3077  0.75828
## lag.CAP_OCC     0.001941497   0.053815702  0.0361  0.97122
## lag.CAP_CS      0.265435956   0.784721558  0.3383  0.73517
## lag.Población  -0.000120261   0.000146346 -0.8218  0.41122
## 
## Rho: 0.28321, LR test value: 0.99104, p-value: 0.31949
## Approximate (numerical Hessian) standard error: 0.26283
##     z-value: 1.0775, p-value: 0.28124
## Wald statistic: 1.1611, p-value: 0.28124
## 
## Log likelihood: -195.4676 for mixed model
## ML residual variance (sigma squared): 11616, (sigma: 107.78)
## Number of observations: 32 
## Number of parameters estimated: 17 
## AIC: 424.94, (AIC for lm: 423.93)
s(reg7)#Manski
## 
## Call:sacsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw,     type = "sacmixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -179.6048  -72.7882   -2.7818   33.8501  346.2488 
## 
## Type: sacmixed 
## Coefficients: (numerical Hessian approximate standard errors) 
##                    Estimate    Std. Error z value Pr(>|z|)
## (Intercept)   307.899524827 289.173507863  1.0648  0.28699
## PIB             0.004070842   0.002881048  1.4130  0.15766
## NBI            -0.195254186   1.932325107 -0.1010  0.91951
## CAP_BOG         0.002531850   0.089756232  0.0282  0.97750
## CAP_BC         -0.001838769   0.009348395 -0.1967  0.84407
## CAP_OCC        -0.056119738   0.028292541 -1.9836  0.04731
## CAP_CS          0.071595235   0.500476865  0.1431  0.88625
## Población       0.000065749   0.000125447  0.5241  0.60019
## lag.PIB         0.000817609   0.007040486  0.1161  0.90755
## lag.NBI        -6.096432791   5.918183252 -1.0301  0.30295
## lag.CAP_BOG    -0.052377434   0.139398457 -0.3757  0.70711
## lag.CAP_BC      0.006801659   0.016628203  0.4090  0.68251
## lag.CAP_OCC     0.008763485   0.049556804  0.1768  0.85964
## lag.CAP_CS      0.219189298   0.815594280  0.2687  0.78812
## lag.Población  -0.000162608   0.000233654 -0.6959  0.48647
## 
## Rho: 0.38485
## Approximate (numerical Hessian) standard error: 0.39392
##     z-value: 0.97698, p-value: 0.32858
## Lambda: -0.26343
## Approximate (numerical Hessian) standard error: 0.88846
##     z-value: -0.2965, p-value: 0.76685
## 
## LR test value: 7.7816, p-value: 0.5563
## 
## Log likelihood: -195.4082 for sacmixed model
## ML residual variance (sigma squared): 11213, (sigma: 105.89)
## Number of observations: 32 
## Number of parameters estimated: 18 
## AIC: 426.82, (AIC for lm: 416.6)
s(reg8)#SARAR lag Y and lag e (SAC)
## 
## Call:sacsarlm(formula = reg.eq1, data = Insumo, listw = Insumo.lw,     type = "sac")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -203.1991  -58.2427   -4.1109   46.5336  421.1733 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##                 Estimate   Std. Error z value Pr(>|z|)
## (Intercept) 73.149386014 83.317191361  0.8780  0.37996
## PIB          0.003487952  0.002651478  1.3155  0.18835
## NBI         -0.337466937  1.448625357 -0.2330  0.81580
## CAP_BOG     -0.038888663  0.051162705 -0.7601  0.44720
## CAP_BC       0.002799262  0.006077082  0.4606  0.64507
## CAP_OCC     -0.043084476  0.021594358 -1.9952  0.04602
## CAP_CS       0.320438586  0.293526107  1.0917  0.27497
## Población    0.000028476  0.000067569  0.4214  0.67343
## 
## Rho: 0.1104
## Approximate (numerical Hessian) standard error: 0.19009
##     z-value: 0.58081, p-value: 0.56137
## Lambda: 0.41973
## Approximate (numerical Hessian) standard error: 0.27257
##     z-value: 1.5399, p-value: 0.12359
## 
## LR test value: 4.1193, p-value: 0.1275
## 
## Log likelihood: -197.2394 for sac model
## ML residual variance (sigma squared): 12624, (sigma: 112.36)
## Number of observations: 32 
## Number of parameters estimated: 11 
## AIC: 416.48, (AIC for lm: 416.6)

6.9 Calculo de varibles significativas

#Calculo de variables signid¿ficativas
reg.eq2=CAP_BAC ~ PIB + CAP_BOG+CAP_BC + CAP_OCC + CAP_CS+ Población
reg4=errorsarlm(reg.eq2,data=Insumo, Insumo.lw)
## Warning in errorsarlm(reg.eq2, data = Insumo, Insumo.lw): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 9.29353e-18 - using numerical Hessian.
s(reg4)#Lag Error (SEM)
## 
## Call:errorsarlm(formula = reg.eq2, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -211.1287  -54.9469   -1.1316   40.5690  428.8317 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value Pr(>|z|)
## (Intercept) 82.780382772 57.200030716  1.4472  0.14784
## PIB          0.003786089  0.002550194  1.4846  0.13764
## CAP_BOG     -0.033655764  0.049317358 -0.6824  0.49497
## CAP_BC       0.001072562  0.005563108  0.1928  0.84712
## CAP_OCC     -0.049110630  0.019547548 -2.5124  0.01199
## CAP_CS       0.290079781  0.285013019  1.0178  0.30878
## Población    0.000039933  0.000066827  0.5976  0.55013
## 
## Lambda: 0.51895, LR test value: 4.4697, p-value: 0.0345
## Approximate (numerical Hessian) standard error: 0.20419
##     z-value: 2.5415, p-value: 0.011039
## Wald statistic: 6.4591, p-value: 0.011039
## 
## Log likelihood: -197.4239 for error model
## ML residual variance (sigma squared): 12470, (sigma: 111.67)
## Number of observations: 32 
## Number of parameters estimated: 9 
## AIC: NA (not available for weighted model), (AIC for lm: 415.32)
reg.eq3=CAP_BAC ~ PIB + CAP_BOG + CAP_OCC + CAP_CS+ Población
reg4=errorsarlm(reg.eq3,data=Insumo, Insumo.lw)
## Warning in errorsarlm(reg.eq3, data = Insumo, Insumo.lw): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 9.36792e-18 - using numerical Hessian.
s(reg4)#Lag Error (SEM)
## 
## Call:errorsarlm(formula = reg.eq3, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -214.6537  -56.5893   -1.9568   41.1618  430.2455 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value  Pr(>|z|)
## (Intercept) 76.659388419 48.733987996  1.5730 0.1157149
## PIB          0.003936351  0.002412095  1.6319 0.1026960
## CAP_BOG     -0.026724271  0.035542537 -0.7519 0.4521140
## CAP_OCC     -0.051550780  0.014822573 -3.4779 0.0005054
## CAP_CS       0.249952769  0.204318663  1.2233 0.2211984
## Población    0.000047471  0.000055107  0.8614 0.3890008
## 
## Lambda: 0.52489, LR test value: 4.7822, p-value: 0.028756
## Approximate (numerical Hessian) standard error: 0.19904
##     z-value: 2.6371, p-value: 0.008363
## Wald statistic: 6.9541, p-value: 0.008363
## 
## Log likelihood: -197.4421 for error model
## ML residual variance (sigma squared): 12461, (sigma: 111.63)
## Number of observations: 32 
## Number of parameters estimated: 8 
## AIC: NA (not available for weighted model), (AIC for lm: 413.67)
reg.eq4=CAP_BAC ~ PIB + CAP_OCC + CAP_CS+ Población
reg4=errorsarlm(reg.eq4,data=Insumo, Insumo.lw)
## Warning in errorsarlm(reg.eq4, data = Insumo, Insumo.lw): inversion of asymptotic covariance matrix failed for tol.solve = 0.000000000000000222044604925031 
##   reciprocal condition number = 9.46755e-18 - using numerical Hessian.
s(reg4)#Lag Error (SEM)
## 
## Call:errorsarlm(formula = reg.eq4, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -200.4867  -63.5891   -8.8979   41.1675  444.7078 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value  Pr(>|z|)
## (Intercept) 78.047294223 51.760317159  1.5079 0.1315904
## PIB          0.003727043  0.002409110  1.5471 0.1218483
## CAP_OCC     -0.050131489  0.014726792 -3.4041 0.0006638
## CAP_CS       0.099564078  0.037675929  2.6426 0.0082261
## Población    0.000057480  0.000053995  1.0646 0.2870778
## 
## Lambda: 0.56286, LR test value: 6.894, p-value: 0.0086485
## Approximate (numerical Hessian) standard error: 0.17904
##     z-value: 3.1437, p-value: 0.001668
## Wald statistic: 9.8832, p-value: 0.001668
## 
## Log likelihood: -197.7004 for error model
## ML residual variance (sigma squared): 12505, (sigma: 111.82)
## Number of observations: 32 
## Number of parameters estimated: 7 
## AIC: NA (not available for weighted model), (AIC for lm: 414.29)
reg.eq5=CAP_BAC ~ PIB + CAP_OCC + CAP_CS
reg4=errorsarlm(reg.eq5,data=Insumo, Insumo.lw)
s(reg4)#Lag Error (SEM)
## 
## Call:errorsarlm(formula = reg.eq5, data = Insumo, listw = Insumo.lw)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -189.815  -68.002  -13.434   35.002  443.595 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                 Estimate   Std. Error z value         Pr(>|z|)
## (Intercept) 100.98165936  46.44259678  2.1743         0.029680
## PIB           0.00611113   0.00092889  6.5789 0.00000000004738
## CAP_OCC      -0.04669321   0.01466331 -3.1844         0.001451
## CAP_CS        0.07261035   0.02868558  2.5312         0.011366
## 
## Lambda: 0.54225, LR test value: 6.0083, p-value: 0.014239
## Asymptotic standard error: 0.17423
##     z-value: 3.1123, p-value: 0.0018565
## Wald statistic: 9.6862, p-value: 0.0018565
## 
## Log likelihood: -198.2513 for error model
## ML residual variance (sigma squared): 13034, (sigma: 114.17)
## Number of observations: 32 
## Number of parameters estimated: 6 
## AIC: NA (not available for weighted model), (AIC for lm: 412.51)

6.10 Mapa Estimado

###Mapa estimado
fit = reg4$fitted.values
#dev.new() #windows()
shade.fit = shading(c(100,130,200,400), cols=(brewer.pal(5,"Reds")))
choropleth(Insumo, fit, shading=shade.fit)
choro.legend(1555874,
             535165.5, 
             shade.fit, 
             fmt="%1.2f", 
             title = "Estimaciones",
             cex=0.7, 
             under = "Menos de", 
             between = "a", 
             over = "Mas de")

title("Valores ajustados mediante el modelo SEM para las captaciones del banco 
      agrario en Colombia, cuarto trimestre 2020", cex.main=1)
map.scale(755874,335165.5, 250000, "km", 2, 50, sfcol='brown')

###R^2 Nagelkerke
# summary.sarlm(reg4,Nagelkerke = TRUE) TO-DO
###Test de moran residuales modelo SEM
moran.test(reg4$residuals, Insumo.lw)
## 
##  Moran I test under randomisation
## 
## data:  reg4$residuals  
## weights: Insumo.lw    
## 
## Moran I statistic standard deviate = 0.83723, p-value = 0.2012
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic       Expectation          Variance 
##        0.05266907       -0.03225806        0.01028984
#Municipal
ColombiaM = readOGR(dsn = "data_3_EstudioDeMercadoEspacial/Geodatabase Colombia", 
                    layer = "municipios")
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## Warning: OGR support is provided by the sf and terra packages among others
## OGR data source with driver: ESRI Shapefile 
## Source: "/home/jncc/Documents/Monitorias/EspacialPage/Clases-EE-UN/data_3_EstudioDeMercadoEspacial/Geodatabase Colombia", layer: "municipios"
## with 1107 features
## It has 6 fields
## Integer64 fields read as strings:  COD_MUN COD_DEPTO