CAR Regresion
Diego J. Lizcano
21 de febrero de 2019
Regresion model with CAR
Example using colombian data for municipalities. The idea is to find the best predictor for IPM.
The response variable in this case is Gausian (Normal), but can be binomial, Gaussian, multinomial, Poisson or zero-inflated Poisson (ZIP). Spatial autocorrelation is acounted and modelled by a set of random effects that are assigned as a conditional autoregressive (CAR) prior distribution.
See more examples in: https://r-spatial.github.io/spdep/articles/CO69.html
library (sf)
library (readxl) # read excel data
library (tidyverse) # Data handling
library (tmap)
library (spdep)
library (MASS)
library (ggiraphExtra)
library (sjPlot)
library (sjmisc)
library (sjlabelled)
library (rsq)
library (ggeffects)Read Data
Read a shp file for colombian municipalities and a excel file with the economic data for each municipalitie.
## Reading layer `municipios' from data source `C:\Users\diego.lizcano\Box Sync\CodigoR\Andres\shp\municipios.shp' using driver `ESRI Shapefile'
## Simple feature collection with 1128 features and 8 fields
## geometry type: POLYGON
## dimension: XY
## bbox: xmin: -233755.1 ymin: -469226.2 xmax: 1410054 ymax: 1490018
## epsg (SRID): 32618
## proj4string: +proj=utm +zone=18 +datum=WGS84 +units=m +no_defs
Merge shp and excel, then plot
# merge both data sets
Full_data <- merge(col_map ,
regre_data,
by="ID")
# see is IMP is Normal ?
hist(Full_data$IPM, main="Rural_MPI") # seems not but ignore
nombres <- c("ID",
"ID_ESPACIA" ,
"AREA_OFICI",
"ENTIDAD_TE" ,
"NOM_DEPART" ,
"NOM_MUNICI",
"COD_DEPTO" ,
"Area_ha" ,
"Region",
"Rural_MPI",
"Intervened_land_area_grazed",
"Total_land_area_grazed",
"Area_excluded_from_grazing",
"Grazed_area_to_manage",
"Grazing_Suitable",
"Cattle_ranching_GDP",
"Productivity_per_capita",
"Productivity_per_hectare",
"Cattle_ranching_jobs",
"Property_taxes_per_capita",
"Land_GINI_index",
"Municipal_fiscal_performance",
"Homicides",
"Kidnappings",
"Ilicit_crops",
"Thefts",
"Refugiees",
"geometry" )
names(Full_data) <- nombres
# fix desplazados Na to 0
index_d <- which(is.na(Full_data$Refugiees))
Full_data$Refugiees[index_d] <- 0
# fix Land_GINI_index Na to mean
index_g <- which(is.na(Full_data$Land_GINI_index))
Full_data$Land_GINI_index[index_g] <- mean (Full_data$Land_GINI_index, na.rm = T)
# fix Municipal_fiscal_performance Na to mean
index_g <- which(is.na(Full_data$Municipal_fiscal_performance))
Full_data$Municipal_fiscal_performance[index_g] <- mean (Full_data$Municipal_fiscal_performance, na.rm = T)
#### plot Rural_MPI
brks <- c(0, 10, 20, 30, 40, 50, 60, 70, 80, 90, 100)
tm_shape(Full_data) + tm_fill("Rural_MPI", breaks=brks, midpoint=50, palette="YlOrRd") + tm_layout(legend.outside=TRUE)
Stepwise multiple regresion, using AIC 2
glm1 starts with all covariates and ends with statistically good predictors.
glm1 <- glm(Rural_MPI ~ Intervened_land_area_grazed+
Total_land_area_grazed+
Area_excluded_from_grazing+
Grazed_area_to_manage+
Grazing_Suitable+
Cattle_ranching_GDP+
Productivity_per_capita+
Productivity_per_hectare+
Cattle_ranching_jobs+
Property_taxes_per_capita+
Land_GINI_index+
Municipal_fiscal_performance+
Homicides+
Kidnappings+
Ilicit_crops+
Thefts+
Refugiees+
Region,
data= Full_data )
step <- stepAIC(glm1, direction="backward", k = 2)## Start: AIC=7714.52
## Rural_MPI ~ Intervened_land_area_grazed + Total_land_area_grazed +
## Area_excluded_from_grazing + Grazed_area_to_manage + Grazing_Suitable +
## Cattle_ranching_GDP + Productivity_per_capita + Productivity_per_hectare +
## Cattle_ranching_jobs + Property_taxes_per_capita + Land_GINI_index +
## Municipal_fiscal_performance + Homicides + Kidnappings +
## Ilicit_crops + Thefts + Refugiees + Region
##
## Df Deviance AIC
## - Productivity_per_capita 1 76391 7712.5
## - Kidnappings 1 76397 7712.6
## - Cattle_ranching_GDP 1 76408 7712.8
## - Total_land_area_grazed 1 76442 7713.2
## - Refugiees 1 76462 7713.5
## - Ilicit_crops 1 76463 7713.5
## <none> 76391 7714.5
## - Cattle_ranching_jobs 1 76627 7715.9
## - Homicides 1 76670 7716.5
## - Grazing_Suitable 1 76856 7719.1
## - Area_excluded_from_grazing 1 77075 7722.2
## - Grazed_area_to_manage 1 77117 7722.7
## - Intervened_land_area_grazed 1 77650 7730.2
## - Land_GINI_index 1 77990 7734.9
## - Productivity_per_hectare 1 78427 7740.9
## - Thefts 1 78884 7747.2
## - Municipal_fiscal_performance 1 79323 7753.2
## - Region 6 89387 7872.2
## - Property_taxes_per_capita 1 93406 7929.7
##
## Step: AIC=7712.52
## Rural_MPI ~ Intervened_land_area_grazed + Total_land_area_grazed +
## Area_excluded_from_grazing + Grazed_area_to_manage + Grazing_Suitable +
## Cattle_ranching_GDP + Productivity_per_hectare + Cattle_ranching_jobs +
## Property_taxes_per_capita + Land_GINI_index + Municipal_fiscal_performance +
## Homicides + Kidnappings + Ilicit_crops + Thefts + Refugiees +
## Region
##
## Df Deviance AIC
## - Kidnappings 1 76397 7710.6
## - Total_land_area_grazed 1 76442 7711.2
## - Refugiees 1 76462 7711.5
## - Ilicit_crops 1 76463 7711.5
## <none> 76391 7712.5
## - Cattle_ranching_jobs 1 76628 7713.9
## - Homicides 1 76670 7714.5
## - Grazing_Suitable 1 76856 7717.1
## - Cattle_ranching_GDP 1 76988 7718.9
## - Area_excluded_from_grazing 1 77075 7720.2
## - Grazed_area_to_manage 1 77119 7720.8
## - Intervened_land_area_grazed 1 77658 7728.3
## - Land_GINI_index 1 78005 7733.1
## - Productivity_per_hectare 1 78429 7739.0
## - Thefts 1 78885 7745.2
## - Municipal_fiscal_performance 1 79328 7751.3
## - Region 6 89443 7870.9
## - Property_taxes_per_capita 1 93411 7927.8
##
## Step: AIC=7710.61
## Rural_MPI ~ Intervened_land_area_grazed + Total_land_area_grazed +
## Area_excluded_from_grazing + Grazed_area_to_manage + Grazing_Suitable +
## Cattle_ranching_GDP + Productivity_per_hectare + Cattle_ranching_jobs +
## Property_taxes_per_capita + Land_GINI_index + Municipal_fiscal_performance +
## Homicides + Ilicit_crops + Thefts + Refugiees + Region
##
## Df Deviance AIC
## - Total_land_area_grazed 1 76447 7709.3
## - Refugiees 1 76469 7709.6
## - Ilicit_crops 1 76471 7709.7
## <none> 76397 7710.6
## - Cattle_ranching_jobs 1 76640 7712.0
## - Homicides 1 76672 7712.5
## - Grazing_Suitable 1 76858 7715.1
## - Cattle_ranching_GDP 1 77002 7717.1
## - Area_excluded_from_grazing 1 77076 7718.2
## - Grazed_area_to_manage 1 77120 7718.8
## - Intervened_land_area_grazed 1 77658 7726.3
## - Land_GINI_index 1 78009 7731.2
## - Productivity_per_hectare 1 78429 7737.0
## - Thefts 1 78902 7743.5
## - Municipal_fiscal_performance 1 79332 7749.3
## - Region 6 89469 7869.2
## - Property_taxes_per_capita 1 93412 7925.8
##
## Step: AIC=7709.32
## Rural_MPI ~ Intervened_land_area_grazed + Area_excluded_from_grazing +
## Grazed_area_to_manage + Grazing_Suitable + Cattle_ranching_GDP +
## Productivity_per_hectare + Cattle_ranching_jobs + Property_taxes_per_capita +
## Land_GINI_index + Municipal_fiscal_performance + Homicides +
## Ilicit_crops + Thefts + Refugiees + Region
##
## Df Deviance AIC
## - Ilicit_crops 1 76511 7708.2
## - Refugiees 1 76520 7708.3
## <none> 76447 7709.3
## - Cattle_ranching_jobs 1 76674 7710.5
## - Homicides 1 76739 7711.4
## - Cattle_ranching_GDP 1 77070 7716.1
## - Grazing_Suitable 1 77545 7722.7
## - Intervened_land_area_grazed 1 77703 7724.9
## - Land_GINI_index 1 78009 7729.2
## - Grazed_area_to_manage 1 78015 7729.2
## - Area_excluded_from_grazing 1 78063 7729.9
## - Productivity_per_hectare 1 78594 7737.2
## - Thefts 1 78911 7741.6
## - Municipal_fiscal_performance 1 79366 7747.8
## - Region 6 92034 7897.7
## - Property_taxes_per_capita 1 93444 7924.1
##
## Step: AIC=7708.22
## Rural_MPI ~ Intervened_land_area_grazed + Area_excluded_from_grazing +
## Grazed_area_to_manage + Grazing_Suitable + Cattle_ranching_GDP +
## Productivity_per_hectare + Cattle_ranching_jobs + Property_taxes_per_capita +
## Land_GINI_index + Municipal_fiscal_performance + Homicides +
## Thefts + Refugiees + Region
##
## Df Deviance AIC
## - Refugiees 1 76600 7707.5
## <none> 76511 7708.2
## - Cattle_ranching_jobs 1 76740 7709.4
## - Homicides 1 76872 7711.3
## - Cattle_ranching_GDP 1 77130 7714.9
## - Grazing_Suitable 1 77652 7722.2
## - Intervened_land_area_grazed 1 77773 7723.9
## - Land_GINI_index 1 78059 7727.8
## - Grazed_area_to_manage 1 78128 7728.8
## - Area_excluded_from_grazing 1 78171 7729.4
## - Productivity_per_hectare 1 78691 7736.6
## - Thefts 1 79062 7741.6
## - Municipal_fiscal_performance 1 79379 7746.0
## - Region 6 92817 7904.9
## - Property_taxes_per_capita 1 93690 7925.0
##
## Step: AIC=7707.48
## Rural_MPI ~ Intervened_land_area_grazed + Area_excluded_from_grazing +
## Grazed_area_to_manage + Grazing_Suitable + Cattle_ranching_GDP +
## Productivity_per_hectare + Cattle_ranching_jobs + Property_taxes_per_capita +
## Land_GINI_index + Municipal_fiscal_performance + Homicides +
## Thefts + Region
##
## Df Deviance AIC
## <none> 76600 7707.5
## - Cattle_ranching_jobs 1 76824 7708.6
## - Homicides 1 76976 7710.8
## - Cattle_ranching_GDP 1 77203 7713.9
## - Grazing_Suitable 1 77790 7722.1
## - Intervened_land_area_grazed 1 77887 7723.5
## - Land_GINI_index 1 78174 7727.4
## - Grazed_area_to_manage 1 78239 7728.3
## - Area_excluded_from_grazing 1 78281 7728.9
## - Productivity_per_hectare 1 78830 7736.5
## - Thefts 1 79069 7739.7
## - Municipal_fiscal_performance 1 79414 7744.4
## - Region 6 93067 7905.8
## - Property_taxes_per_capita 1 93692 7923.0GLM (regresion) con las variables filtradas por el stepwise regresion con AIC 2
glm2 is assembled with just the statistically significant variables.
glm2 <- glm (Rural_MPI ~ Intervened_land_area_grazed +
Area_excluded_from_grazing +
Grazed_area_to_manage +
Grazing_Suitable +
Cattle_ranching_GDP +
Productivity_per_hectare +
Cattle_ranching_jobs +
Property_taxes_per_capita +
Land_GINI_index +
Municipal_fiscal_performance +
Homicides +
Thefts +
Region,
data = Full_data)
######## coeficientes para el mejor modelo
summary(glm2)##
## Call:
## glm(formula = Rural_MPI ~ Intervened_land_area_grazed + Area_excluded_from_grazing +
## Grazed_area_to_manage + Grazing_Suitable + Cattle_ranching_GDP +
## Productivity_per_hectare + Cattle_ranching_jobs + Property_taxes_per_capita +
## Land_GINI_index + Municipal_fiscal_performance + Homicides +
## Thefts + Region, data = Full_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -35.217 -4.888 0.516 5.371 42.354
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.078e+02 3.709e+00 29.076 < 2e-16 ***
## Intervened_land_area_grazed 8.600e+00 2.037e+00 4.221 2.64e-05 ***
## Area_excluded_from_grazing -8.608e+00 1.784e+00 -4.825 1.61e-06 ***
## Grazed_area_to_manage -9.227e+00 1.937e+00 -4.764 2.16e-06 ***
## Grazing_Suitable -6.312e+00 1.555e+00 -4.059 5.29e-05 ***
## Cattle_ranching_GDP 3.212e-06 1.112e-06 2.888 0.003950 **
## Productivity_per_hectare -3.083e+00 5.547e-01 -5.558 3.45e-08 ***
## Cattle_ranching_jobs -1.142e-02 6.490e-03 -1.760 0.078776 .
## Property_taxes_per_capita -1.551e+02 1.008e+01 -15.387 < 2e-16 ***
## Land_GINI_index -1.359e+01 2.910e+00 -4.669 3.41e-06 ***
## Municipal_fiscal_performance -2.251e-01 3.606e-02 -6.242 6.22e-10 ***
## Homicides 1.583e-02 6.934e-03 2.283 0.022640 *
## Thefts -1.835e-02 3.138e-03 -5.847 6.64e-09 ***
## RegionANDINA -4.073e+00 1.809e+00 -2.252 0.024528 *
## RegionCARIBE 7.250e+00 1.876e+00 3.864 0.000118 ***
## RegionCATATUMBO 4.528e+00 3.635e+00 1.246 0.213089
## RegionMAGDALENA 2.631e+00 2.118e+00 1.242 0.214411
## RegionORINOQUIA 6.247e-01 2.486e+00 0.251 0.801629
## RegionPACIFICO 7.193e+00 2.251e+00 3.196 0.001436 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 72.19629)
##
## Null deviance: 193614 on 1079 degrees of freedom
## Residual deviance: 76600 on 1061 degrees of freedom
## (44 observations deleted due to missingness)
## AIC: 7707.5
##
## Number of Fisher Scoring iterations: 2tab_model(glm2, show.aic=TRUE)| Rural MPI | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 107.85 | 100.58 – 115.12 | <0.001 |
|
Intervened land area grazed |
8.60 | 4.61 – 12.59 | <0.001 |
|
Area excluded from grazing |
-8.61 | -12.11 – -5.11 | <0.001 |
| Grazed area to manage | -9.23 | -13.02 – -5.43 | <0.001 |
| Grazing Suitable | -6.31 | -9.36 – -3.26 | <0.001 |
| Cattle ranching GDP | 0.00 | 0.00 – 0.00 | 0.004 |
| Productivity per hectare | -3.08 | -4.17 – -2.00 | <0.001 |
| Cattle ranching jobs | -0.01 | -0.02 – 0.00 | 0.079 |
| Property taxes per capita | -155.09 | -174.85 – -135.34 | <0.001 |
| Land GINI index | -13.59 | -19.29 – -7.88 | <0.001 |
|
Municipal fiscal performance |
-0.23 | -0.30 – -0.15 | <0.001 |
| Homicides | 0.02 | 0.00 – 0.03 | 0.023 |
| Thefts | -0.02 | -0.02 – -0.01 | <0.001 |
| RegionANDINA | -4.07 | -7.62 – -0.53 | 0.025 |
| RegionCARIBE | 7.25 | 3.57 – 10.93 | <0.001 |
| RegionCATATUMBO | 4.53 | -2.60 – 11.65 | 0.213 |
| RegionMAGDALENA | 2.63 | -1.52 – 6.78 | 0.214 |
| RegionORINOQUIA | 0.62 | -4.25 – 5.50 | 0.802 |
| RegionPACIFICO | 7.19 | 2.78 – 11.60 | 0.001 |
| Observations | 1080 | ||
| Cox & Snell’s R2 / Nagelkerke’s R2 | 1.000 / 1.000 | ||
| AIC | 7707.478 | ||
######################################
####### calcular el R cuadrado #######
######################################
rsq(glm2)## [1] 0.6043655###### Graficas
###### By default, estimates are sorted in descending order, with the highest effect at the top.
plot_model(glm2, type="std",
title = "MPI (standardized beta values)",
show.values = TRUE, value.offset = .4) # standardized beta values.
plot_model(glm2, type="slope") 
plot_model(glm2, type="pred") ## $Intervened_land_area_grazed
##
## $Area_excluded_from_grazing
##
## $Grazed_area_to_manage
##
## $Grazing_Suitable
##
## $Cattle_ranching_GDP
##
## $Productivity_per_hectare
##
## $Cattle_ranching_jobs
##
## $Property_taxes_per_capita
##
## $Land_GINI_index
##
## $Municipal_fiscal_performance
##
## $Homicides
##
## $Thefts
##
## $Region
p <- ggpredict(glm2, c("Area_excluded_from_grazing", "Region"))
# plot(p, facet = TRUE)glm3. Regresion con solo las variables ganaderas
glm3 <- glm (Rural_MPI ~ Intervened_land_area_grazed +
Area_excluded_from_grazing +
Grazed_area_to_manage + Grazing_Suitable +
Cattle_ranching_GDP +
Productivity_per_hectare +
Cattle_ranching_jobs,
data = Full_data)
######## coeficientes para el mejor modelo
summary(glm3)##
## Call:
## glm(formula = Rural_MPI ~ Intervened_land_area_grazed + Area_excluded_from_grazing +
## Grazed_area_to_manage + Grazing_Suitable + Cattle_ranching_GDP +
## Productivity_per_hectare + Cattle_ranching_jobs, data = Full_data)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -57.361 -6.103 1.631 8.159 59.857
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 6.918e+01 2.478e+00 27.915 < 2e-16 ***
## Intervened_land_area_grazed 2.569e+01 2.714e+00 9.464 < 2e-16 ***
## Area_excluded_from_grazing -1.816e+01 2.173e+00 -8.359 < 2e-16 ***
## Grazed_area_to_manage -1.090e+01 2.502e+00 -4.355 1.46e-05 ***
## Grazing_Suitable -1.403e+01 1.937e+00 -7.244 8.31e-13 ***
## Cattle_ranching_GDP 1.760e-06 1.394e-06 1.263 0.207
## Productivity_per_hectare -4.701e+00 7.622e-01 -6.167 9.83e-10 ***
## Cattle_ranching_jobs 1.236e-04 8.744e-03 0.014 0.989
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 145.175)
##
## Null deviance: 193614 on 1079 degrees of freedom
## Residual deviance: 155628 on 1072 degrees of freedom
## (44 observations deleted due to missingness)
## AIC: 8451.1
##
## Number of Fisher Scoring iterations: 2tab_model(glm3, show.aic=TRUE)| Rural MPI | |||
|---|---|---|---|
| Predictors | Estimates | CI | p |
| (Intercept) | 69.18 | 64.32 – 74.04 | <0.001 |
|
Intervened land area grazed |
25.69 | 20.37 – 31.01 | <0.001 |
|
Area excluded from grazing |
-18.16 | -22.42 – -13.90 | <0.001 |
| Grazed area to manage | -10.90 | -15.80 – -5.99 | <0.001 |
| Grazing Suitable | -14.03 | -17.82 – -10.23 | <0.001 |
| Cattle ranching GDP | 0.00 | -0.00 – 0.00 | 0.207 |
| Productivity per hectare | -4.70 | -6.19 – -3.21 | <0.001 |
| Cattle ranching jobs | 0.00 | -0.02 – 0.02 | 0.989 |
| Observations | 1080 | ||
| Cox & Snell’s R2 / Nagelkerke’s R2 | 1.000 / 1.000 | ||
| AIC | 8451.053 | ||
######################################
####### calcular el R cuadrado #######
######################################
rsq(glm3)## [1] 0.1961953###### Graficas
plot_model(glm3, type="std",
title = "MPI (standardized beta values)",
show.values = TRUE, value.offset = .4) # standardized beta values.
plot_model(glm3, type="slope") 
plot_model(glm3, type="pred") ## $Intervened_land_area_grazed
##
## $Area_excluded_from_grazing
##
## $Grazed_area_to_manage
##
## $Grazing_Suitable
##
## $Cattle_ranching_GDP
##
## $Productivity_per_hectare
##
## $Cattle_ranching_jobs
CAR model
Spatial Autoregression Model Estimation by Maximum Likelihood. The implementation in the package spautolm use Genenarlised Least Squares (GLS), as a single spatial coefficient value, here termed lambda.
#### drop datos con NA
# drop Municipios sin IPM
outl <- which(is.na( Full_data$Rural_MPI))
as.character(Full_data$NOM_MUNICI[outl])## [1] "NOROSI"
## [2] "GUACHENÉ"
## [3] "SAN JOSÉ DE URÉ"
## [4] "TUCHÍN"
## [5] "MAPIRIPÁN"
## [6] "PUERTO RONDÓN"
## [7] "EL ENCANTO (Cor. Departamental)"
## [8] "LA CHORRERA (Cor. Departamental)"
## [9] "LA PEDRERA (Cor. Departamental)"
## [10] "LA VICTORIA (Pacoa) (Cor. Departamental)"
## [11] "MIRITÍ-PARANÁ (Campoamor) (Cor. Departamental)"
## [12] "PUERTO ALEGRÍA (Cor. Departamental)"
## [13] "PUERTO ARICA (Cor. Departamental)"
## [14] "SANTANDER (Araracuara) (Cor. Departamental)"
## [15] "TARAPACÁ (Cor. Departamental)"
## [16] "BARRANCO MINA (Cor. Departamental)"
## [17] "MAPIRIPANA (Cor. Departamental)"
## [18] "SAN FELIPE (Cor. Departamental)"
## [19] "PUERTO COLOMBIA (Cor. Departamental)"
## [20] "LA GUADALUPE (Cor. Departamental)"
## [21] "CACAHUAL (Cor. Departamental)"
## [22] "PANÁ PANÁ (Campo Alegre) (Cor. Departamental)"
## [23] "MORICHAL (Morichal Nuevo) (Cor. Departamental)"
## [24] "CARURÚ"
## [25] "PACOA (Cor. Departamental)"
## [26] "PAPUNAUA (Cor. Departamental)"
## [27] "YAVARATÉ (Cor. Departamental)"data_sin_na1 <- Full_data[-outl, -outl]
# drop Municipios sin Empleos
out2 <- which(is.na( data_sin_na1$Cattle_ranching_jobs))
as.character(data_sin_na1$NOM_MUNICI[out2])## [1] "ATRATO (Yuto)" "BOJAYÁ (Bellavista)"
## [3] "CARMEN DEL DARIÉN (Curbaradó)" "CERTEGUI"
## [5] "MEDIO ATRATO (Beté)" "MEDIO BAUDÓ (Boca de Pepé)"
## [7] "MEDIO SAN JUAN (Andagoya)" "RÍO IRÓ (Santa Rita)"
## [9] "RÍO QUITO (Paimadó)" "SAN JOSÉ DEL PALMAR"
## [11] "UNIÓN PANAMERICANA (Ánimas)" "NARIÑO"
## [13] "TARAIRA"data_sin_na2 <- data_sin_na1[-out2, -out2]
# drop Municipios sin PIB_Ganadero
out3 <- which(is.na( data_sin_na2$Cattle_ranching_GDP))
as.character(data_sin_na2$NOM_MUNICI[out3])## [1] "EL LITORAL DEL SAN JUAN"data_sin_na3 <- data_sin_na2[-out3, -out3]
# drop Municipios sin Productividad_ha
out4 <- which(is.na( data_sin_na3$Productivity_per_hectare))
as.character(data_sin_na3$NOM_MUNICI[out4])## [1] "LÓPEZ" "TIMBIQUÍ"
## [3] "SANTA BÁRBARA (Iscuandé)"data_sin_na4 <- data_sin_na3[-out4, -out4]
out_total <- c(outl, out2, out3, out4, 1118, 1115, 578)
out5 <- c(1118, 1115, 730, 578)
data_sin_na5 <- data_sin_na4[!data_sin_na4$NOM_MUNICI %in%
c("ACANDÍ", "MITÚ", "EL CHARCO"),]
nb_Full_data_2 <- poly2nb(data_sin_na5)
nb_B_2 <- nb2listw(nb_Full_data_2, style="B", zero.policy=TRUE)
Full_data_2 <- Full_data[!(1:length(Full_data$ID) %in% out_total),]
########################################
#################### Moran test
########################################
form1 <- Rural_MPI ~ Intervened_land_area_grazed +
Area_excluded_from_grazing +
Grazed_area_to_manage + Grazing_Suitable +
Cattle_ranching_GDP +
Productivity_per_hectare +
Cattle_ranching_jobs +
Property_taxes_per_capita +
Land_GINI_index +
Municipal_fiscal_performance +
Homicides +
Thefts +
Refugiees +
Region
model1 <- lm(formula=form1, data=data.frame(data_sin_na5))
W.list <- nb_B_2
resid.model <- residuals(model1)
# moran.mc(x=resid.model, listw=W.list, nsim=1000,
# zero.policy=TRUE, na.action=na.omit, adjust.n=TRUE)
# The Moran's I test has a p-value much less than 0.05, which suggests that the residuals contain
# substantial positive spatial autocorrelation.
# Convert to vecinos y lista
nb_Full_data <- poly2nb(Full_data)
nb_B_Full <- nb2listw(nb_Full_data, style="B", zero.policy=TRUE)
#B <- as(nb_B, "CsparseMatrix")
#B <- as(nb_B, "CsparseMatrix")
# plot vecinos
plot(st_geometry(data_sin_na5), border="grey")
plot(nb_B_2, st_centroid(st_geometry(data_sin_na5), of_largest_polygon), add=TRUE, col="blue")
final_data <- data.frame(data_sin_na5)
############## Names
########### Model 1 test
# CAR_mod1 <- spautolm(IPM ~ PIB_Ganadero +
# Proporcion_municipio +
# Productividad_ha, data=final_data, listw=nb_B_2, # family="CAR")
# summary(CAR_mod1)
# data_sin_na5$IPM_3covs <- fitted.values(CAR_mod1)
# tm_shape(data_sin_na5) + tm_fill("IPM_3covs")
########### Model 2 Full model
CAR_mod2 <- spautolm(Rural_MPI ~ Intervened_land_area_grazed +
Area_excluded_from_grazing +
Grazed_area_to_manage + Grazing_Suitable +
Cattle_ranching_GDP +
Productivity_per_hectare +
Cattle_ranching_jobs +
Property_taxes_per_capita +
Land_GINI_index +
Municipal_fiscal_performance +
Homicides +
Thefts +
Refugiees +
Region,
data=final_data,
listw=nb_B_2,
na.action="na.omit",
family="CAR")
summary(CAR_mod2, adj.se=TRUE, Nagelkerke=TRUE)##
## Call: spautolm(formula = Rural_MPI ~ Intervened_land_area_grazed +
## Area_excluded_from_grazing + Grazed_area_to_manage + Grazing_Suitable +
## Cattle_ranching_GDP + Productivity_per_hectare + Cattle_ranching_jobs +
## Property_taxes_per_capita + Land_GINI_index + Municipal_fiscal_performance +
## Homicides + Thefts + Refugiees + Region, data = final_data,
## listw = nb_B_2, na.action = "na.omit", family = "CAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.87263 -4.21067 0.36398 4.47418 30.94798
##
## Coefficients:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 9.9318e+01 3.7281e+00 26.6405 < 2.2e-16
## Intervened_land_area_grazed 7.9551e+00 2.0940e+00 3.7990 0.0001453
## Area_excluded_from_grazing -4.5851e+00 1.8667e+00 -2.4562 0.0140406
## Grazed_area_to_manage -5.3557e+00 1.9758e+00 -2.7107 0.0067142
## Grazing_Suitable -6.4894e+00 1.6713e+00 -3.8828 0.0001032
## Cattle_ranching_GDP 2.6754e-06 1.1162e-06 2.3970 0.0165303
## Productivity_per_hectare -1.7110e+00 4.9588e-01 -3.4505 0.0005596
## Cattle_ranching_jobs -8.3518e-03 6.5428e-03 -1.2765 0.2017862
## Property_taxes_per_capita -9.0222e+01 9.4539e+00 -9.5434 < 2.2e-16
## Land_GINI_index -1.3875e+01 2.8118e+00 -4.9346 8.031e-07
## Municipal_fiscal_performance -1.1407e-01 3.2803e-02 -3.4776 0.0005060
## Homicides 1.1471e-02 6.8429e-03 1.6763 0.0936847
## Thefts -2.1815e-02 3.1440e-03 -6.9387 3.959e-12
## Refugiees 3.7725e-04 1.0193e-04 3.7009 0.0002148
## RegionANDINA -4.8891e+00 2.0351e+00 -2.4023 0.0162912
## RegionCARIBE 5.0744e+00 2.2957e+00 2.2104 0.0270774
## RegionCATATUMBO 1.6782e+00 3.9374e+00 0.4262 0.6699397
## RegionMAGDALENA -1.6440e+00 2.4986e+00 -0.6580 0.5105569
## RegionORINOQUIA -1.9774e+00 2.6947e+00 -0.7338 0.4630759
## RegionPACIFICO 5.9503e+00 2.5869e+00 2.3002 0.0214359
##
## Lambda: 0.14871 LR test value: 265.78 p-value: < 2.22e-16
## Numerical Hessian standard error of lambda: 0.001162
##
## Log likelihood: -3690.16
## Residual variance (sigma squared): 50.298, (sigma: 7.0921)
## Number of observations: 1077
## Number of parameters estimated: 22
## AIC: 7424.3
## Nagelkerke pseudo-R-squared: 0.69056# add fited values to dataframe
data_sin_na5$MPI_predicted <- fitted.values(CAR_mod2)
# plot tm
tm_shape(data_sin_na5) + tm_fill("MPI_predicted", breaks=brks, midpoint=50, palette="YlOrRd") + tm_layout(legend.outside=TRUE)
CAR Solo variables ganaderas
CAR_mod3 <- spautolm(Rural_MPI ~ Intervened_land_area_grazed +
Area_excluded_from_grazing +
Grazed_area_to_manage + Grazing_Suitable +
Cattle_ranching_GDP +
Productivity_per_hectare +
Cattle_ranching_jobs,
data=final_data,
listw=nb_B_2,
na.action="na.omit",
family="CAR")
summary(CAR_mod3, adj.se=TRUE, Nagelkerke=TRUE)##
## Call: spautolm(formula = Rural_MPI ~ Intervened_land_area_grazed +
## Area_excluded_from_grazing + Grazed_area_to_manage + Grazing_Suitable +
## Cattle_ranching_GDP + Productivity_per_hectare + Cattle_ranching_jobs,
## data = final_data, listw = nb_B_2, na.action = "na.omit",
## family = "CAR")
##
## Residuals:
## Min 1Q Median 3Q Max
## -37.40347 -4.30243 0.56039 4.89731 33.80127
##
## Coefficients:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 7.7129e+01 2.2108e+00 34.8867 < 2.2e-16
## Intervened_land_area_grazed 1.2713e+01 2.4071e+00 5.2815 1.281e-07
## Area_excluded_from_grazing -6.9912e+00 2.0419e+00 -3.4238 0.0006175
## Grazed_area_to_manage -6.7548e+00 2.2132e+00 -3.0520 0.0022730
## Grazing_Suitable -9.2529e+00 1.8188e+00 -5.0875 3.629e-07
## Cattle_ranching_GDP 1.7271e-06 1.2552e-06 1.3760 0.1688310
## Productivity_per_hectare -2.4525e+00 5.7064e-01 -4.2977 1.725e-05
## Cattle_ranching_jobs -1.3936e-03 7.4600e-03 -0.1868 0.8518102
##
## Lambda: 0.14945 LR test value: 677.04 p-value: < 2.22e-16
## Numerical Hessian standard error of lambda: 0.00024547
##
## Log likelihood: -3867.266
## Residual variance (sigma squared): 68.822, (sigma: 8.2959)
## Number of observations: 1077
## Number of parameters estimated: 10
## AIC: 7754.5
## Nagelkerke pseudo-R-squared: 0.57006# add fited values to dataframe
data_sin_na5$MPI_predicted_cattle <- fitted.values(CAR_mod3)
# plot tm
tm_shape(data_sin_na5) + tm_fill("MPI_predicted_cattle", breaks=brks, midpoint=50, palette="YlOrRd") + tm_layout(legend.outside=TRUE)
R sesion info
print(sessionInfo(), locale = FALSE)## R version 3.5.2 (2018-12-20)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 17134)
##
## Matrix products: default
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] ggeffects_0.8.0 rsq_1.1 sjlabelled_1.0.16
## [4] sjmisc_2.7.7.9000 sjPlot_2.6.2.9000 ggiraphExtra_0.2.9
## [7] MASS_7.3-51.1 spdep_1.0-2 spData_0.3.0
## [10] Matrix_1.2-15 sp_1.3-1 tmap_2.2
## [13] forcats_0.4.0 stringr_1.4.0 dplyr_0.8.0.1
## [16] purrr_0.3.0 readr_1.3.1 tidyr_0.8.2
## [19] tibble_2.0.1 ggplot2_3.1.0 tidyverse_1.2.1
## [22] readxl_1.2.0 sf_0.7-2
##
## loaded via a namespace (and not attached):
## [1] backports_1.1.3 lwgeom_0.1-5 plyr_1.8.4
## [4] lazyeval_0.2.1 TMB_1.7.15 splines_3.5.2
## [7] mycor_0.1.1 crosstalk_1.0.0 leaflet_2.0.2
## [10] TH.data_1.0-10 digest_0.6.18 htmltools_0.3.6
## [13] gdata_2.18.0 magrittr_1.5 modelr_0.1.4
## [16] gmodels_2.18.1 sandwich_2.5-0 colorspace_1.4-0
## [19] rvest_0.3.2 haven_2.1.0 xfun_0.5
## [22] rgdal_1.3-6 crayon_1.3.4 jsonlite_1.6
## [25] lme4_1.1-20 survival_2.43-3 zoo_1.8-4
## [28] glue_1.3.0 gtable_0.2.0 ppcor_1.1
## [31] emmeans_1.3.2 webshot_0.5.1 sjstats_0.17.3
## [34] scales_1.0.0 mvtnorm_1.0-8 DBI_1.0.0
## [37] Rcpp_1.0.0 viridisLite_0.3.0 xtable_1.8-3
## [40] units_0.6-2 foreign_0.8-71 stats4_3.5.2
## [43] prediction_0.3.6.2 htmlwidgets_1.3 httr_1.4.0
## [46] RColorBrewer_1.1-2 modeltools_0.2-22 pkgconfig_2.0.2
## [49] XML_3.98-1.16 deldir_0.1-16 labeling_0.3
## [52] tidyselect_0.2.5 rlang_0.3.1 reshape2_1.4.3
## [55] later_0.7.5 tmaptools_2.0-1 munsell_0.5.0
## [58] cellranger_1.1.0 tools_3.5.2 cli_1.0.1
## [61] generics_0.0.2 broom_0.5.1 ggridges_0.5.1
## [64] evaluate_0.13 yaml_2.2.0 knitr_1.21
## [67] satellite_1.0.1 coin_1.2-2 nlme_3.1-137
## [70] mime_0.6 ggiraph_0.6.0 xml2_1.2.0
## [73] compiler_3.5.2 bayesplot_1.6.0 rstudioapi_0.9.0
## [76] png_0.1-7 e1071_1.7-0.1 stringi_1.3.1
## [79] gdtools_0.1.7 rgeos_0.4-2 lattice_0.20-38
## [82] classInt_0.3-1 psych_1.8.12 nloptr_1.2.1
## [85] pillar_1.3.1 LearnBayes_2.15.1 pwr_1.2-2
## [88] estimability_1.3 data.table_1.12.0 raster_2.8-19
## [91] mapview_2.6.3 httpuv_1.4.5.1 R6_2.4.0
## [94] promises_1.0.1 KernSmooth_2.23-15 codetools_0.2-16
## [97] dichromat_2.0-0 boot_1.3-20 gtools_3.8.1
## [100] assertthat_0.2.0 withr_2.1.2 mnormt_1.5-5
## [103] multcomp_1.4-8 parallel_3.5.2 mgcv_1.8-26
## [106] expm_0.999-3 hms_0.4.2 grid_3.5.2
## [109] coda_0.19-2 glmmTMB_0.2.3 class_7.3-15
## [112] minqa_1.2.4 rmarkdown_1.11 snakecase_0.9.2
## [115] shiny_1.2.0 lubridate_1.7.4 base64enc_0.1-3