Comparación de modelos espaciales globales y locales para la toma de decisiones agropecuarias en Paraguay: SDM versus GWR con evidencias del Censo Agropecuario 2022
Richard González
16 de junio de 2025
1 Resumen
Este anteproyecto plantea un estudio comparativo entre dos enfoques de análisis espacial —el Modelo Espacial de Durbin (SDM) y la Regresión Geográficamente Ponderada (GWR)— aplicado a indicadores agropecuarios disgregados al nivel de distrito en Paraguay. El objetivo central es evaluar cómo varía la intensidad de mecanización y la diversificación ganadera entre distritos, identificando diferencias en los efectos globales (SDM) frente a efectos locales (GWR), y determinar cuál herramienta ofrece mejor ajuste y explicación de la variabilidad espacial.
2 Introducción
3 Objetivos:
3.1 General
Comparar, a nivel de distrito, la capacidad explicativa y predictiva de los Modelos Espaciales de Durbin (SDM) versus la Regresión Geográficamente Ponderada (GWR) sobre índices de mecanización y diversificación ganadera derivados del CAN 2022.
3.2 Específicos
Construir índices distritales de mecanización (proporción de tecnologías usadas) y diversificación ganadera (conteo de especies).
Estimar un SDM para cada índice, evaluando impactos directos e indirectos globales.
Ajustar una GWR para cada índice, mapeando coeficientes locales de escala (ln(sup agrícola+1)).
Contrastar resultados: bondad de ajuste, significancia espacial y utilidad práctica para políticas focalizadas.
# Paquetes indispensables -------------------------------------------------
#libs <- c("tidyverse", "sf", "spdep", "spatialreg", "tmap", "knitr", "GWmodel")
#inst <- libs[!libs %in% installed.packages()]
#if(length(inst)) install.packages(inst)
#invisible(lapply(libs, library, character.only = TRUE))
#options(tmap.mode = "plot") # mapas estáticos por defecto## Warning: package 'tidyverse' was built under R version 4.3.3## Warning: package 'ggplot2' was built under R version 4.3.3## Warning: package 'tibble' was built under R version 4.3.3## Warning: package 'tidyr' was built under R version 4.3.3## Warning: package 'readr' was built under R version 4.3.3## Warning: package 'purrr' was built under R version 4.3.3## Warning: package 'dplyr' was built under R version 4.3.3## Warning: package 'stringr' was built under R version 4.3.3## Warning: package 'forcats' was built under R version 4.3.3## Warning: package 'lubridate' was built under R version 4.3.3## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.3 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors## Warning: package 'sf' was built under R version 4.3.3## Linking to GEOS 3.11.2, GDAL 3.8.2, PROJ 9.3.1; sf_use_s2() is TRUE## Warning: package 'spdep' was built under R version 4.3.3## Loading required package: spData## Warning: package 'spData' was built under R version 4.3.3## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`## Warning: package 'spatialreg' was built under R version 4.3.3## Loading required package: Matrix
##
## Attaching package: 'Matrix'
##
## The following objects are masked from 'package:tidyr':
##
## expand, pack, unpack
##
##
## Attaching package: 'spatialreg'
##
## The following objects are masked from 'package:spdep':
##
## get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption## Warning: package 'sp' was built under R version 4.3.3## Warning: package 'GWmodel' was built under R version 4.3.3## Loading required package: robustbase## Warning: package 'robustbase' was built under R version 4.3.3## Loading required package: Rcpp## Warning: package 'Rcpp' was built under R version 4.3.3## Welcome to GWmodel version 2.4-2.## Warning: package 'tmap' was built under R version 4.3.34 Materiales y métodos
4.1 Qué es el CAN 2022
El Censo Agropecuario Nacional (CAN) 2022 es la operación estadística más completa realizada en Paraguay para caracterizar la estructura, el uso de la tierra y las condiciones de los productores y sus hogares en el ámbito rural. Tras 14 años sin actualización (el anterior fue en 2008), el MAG (Ministerio de Agricultura y Ganadería) —con la coordinación técnica de la DGEEC— definió una “finca agropecuaria” como aquella que cumple al menos uno de estos criterios:
- ≥ 0,1 ha con cultivos para venta
- ≥ 3 vacunos mayores
- ≥ 5 ovinos, caprinos o porcinos
- ≥ 100 aves de corral
- ≥ 10 colmenas
- ≥ 500 m² de espejo de agua para piscicultura
- ≥ 2 ha de plantaciones forestales comerciales ([henoi.org.py][1])
Cobertura y despliegue El levantamiento de campo comenzó el 2 de agosto de 2022 y, tras sucesivas etapas, culminó a fines de ese año. Se visitaron 454 950 unidades de observación, de las cuales 336 742 resultaron parcelas con actividad agropecuaria y 291 497 se clasificaron como fincas agropecuarias propiamente dichas. El 97 % de esas fincas se encuentra en la Región Oriental (44 % de la superficie total) y el 3 % en la Región Occidental (56 % de la superficie) ([opsaa.iica.int][2], [henoi.org.py][1]).
Principales resultados
- Estructura de tenencia: 95,4 % de las fincas son manejadas por un solo productor (96,7 % paraguayos), 55 % cuentan con título definitivo y 71,4 % cultivos temporales.
- Composición productiva: 69 % de las fincas poseen ganado vacuno y 72 % aves de corral; cultivos temporales dominados por soja y maíz, y dentro de los permanentes, yerba mate y banano.
- Asistencia y financiamiento: Solo el 15,2 % recibió asistencia técnica y el 14,9 % accedió a créditos agropecuarios ([opsaa.iica.int][2], [henoi.org.py][1]).
Alcance del microdato La base resultante incluye variables detalladas sobre:
- Territorio y equipamiento: superficies (agrícola, total, regadío), tipo de manejo (propia, arrendada), infraestructura (riego, silos, cercos).
- Tecnificación: uso de tractores con/sin GPS, drones, sembradoras con sensor, pulverizadores, etc.
- Insumos y prácticas agronómicas: semillas híbridas/transgénicas, fertilizantes químicos y orgánicos, rotación de cultivos, sistemas silvopastoriles.
- Hogar y capital humano: nivel educativo (CLASIFICANEW), tamaño y composición familiar, seguridad alimentaria (p30000), mano de obra familiar y asalariada.
- Ganadería: presencia y totales de vacunos, búfalos, ovinos, caprinos, porcinos y aves.
- Servicios: asistencia técnica, crédito y rechazo de crédito; mercados de comercialización y contratos de venta. ([datos.gov.py][3])
4.2 Métodos estadísticos
4.2.1 Modelo espacial de Durbin (SMD)
4.2.2 Modelo alternativo
5 Gestión de dastos
5.1 Carga de datos
# Directorios --------------------------------------------------------------
# Shapefile distritos en metros (EPSG:32721) ------------------------------
shp <- st_read(file.path("G:/Mi unidad/TESIS_RICHARD_GONZALEZ/paraguay/gadm40_PRY_2.shp")) %>%
st_make_valid()## Reading layer `gadm40_PRY_2' from data source
## `G:\Mi unidad\TESIS_RICHARD_GONZALEZ\paraguay\gadm40_PRY_2.shp'
## using driver `ESRI Shapefile'
## Simple feature collection with 218 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -62.64652 ymin: -27.60586 xmax: -54.25863 ymax: -19.29137
## Geodetic CRS: WGS 84## [1] "ID_0" "COUNTRY" "NAME_1" "NL_NAME_1" "ID_2" "NAME_2"
## [7] "VARNAME_2" "NL_NAME_2" "TYPE_2" "ENGTYPE_2" "CC_2" "HASC_2"
## [13] "geometry"can <- read_delim("G:/Mi unidad/TESIS_RICHARD_GONZALEZ/BASE_CAN2022_DATOS_GENERALES_P2.csv",
delim = ",",
show_col_types = FALSE
)## Warning: One or more parsing issues, call `problems()` on your data frame for details,
## e.g.:
## dat <- vroom(...)
## problems(dat)## [1] "ID_FCA_ALEATORIO" "CLASIFICANEW" "DPTO_PROC"
## [4] "DSTO_PROC" "tipoi" "tipofinca"
## [7] "p5000" "menores_total" "mayores_total"
## [10] "p7001_sn" "p7021p" "p7021p_nosabe"
## [13] "p7033p" "p7033p_nosabe" "p7221t"
## [16] "p7221t_nosabe" "p7223t" "p7223t_nosabe"
## [19] "p8000_1" "p8000_3" "p9001"
## [22] "p9000_total" "p10001" "p10002"
## [25] "p10003" "p10004" "p10005"
## [28] "p10006" "p10007" "p10008"
## [31] "p10009" "p10000_total" "p12000_sn"
## [34] "p12000_total" "p12000b_total" "p13300_total"
## [37] "p14000_sn" "p14000_total" "p15000_sn"
## [40] "p15000_total" "p16000_sn" "p16000_total"
## [43] "p18000_sn" "p18002" "XP18003_1"
## [46] "XP18003_2" "XP18003_3" "p18003_esp"
## [49] "p18003_0" "p18003_1" "p18003_2"
## [52] "p18003_3" "p19100_sn" "p19101"
## [55] "p19102" "p19103" "p19104"
## [58] "p19105" "p19106" "p19107"
## [61] "p19108" "p19109" "p19110"
## [64] "p19111" "p19112" "p19119"
## [67] "p19121" "p19122" "p19123"
## [70] "p19123_nosabe" "xp19124_1" "xp19124_2"
## [73] "xp19124_3" "xp19124_4" "p19125_1"
## [76] "p19131" "p19132" "p19133"
## [79] "p19135_1_sn" "p19136_sn" "p19139_sn"
## [82] "p19140_sn" "p19300_sn" "p19303_total"
## [85] "raza_oveja_total" "p19400_sn" "p19405_total"
## [88] "raza_cabra_total" "p19500_sn" "p19509_total"
## [91] "p19700_sn" "p19701" "p19710"
## [94] "p19800_sn" "p19802" "p19802_nosabe"
## [97] "p210051_sn" "p21010_sn" "p21011_sn"
## [100] "p21013_sn" "p21014_sn" "p21016_sn"
## [103] "p21017_sn" "p21018_sn" "p21018_1"
## [106] "XP21018_2_1" "XP21018_2_2" "XP21018_2_3"
## [109] "XP21018_2_4" "p22000_sn" "p23100_sn"
## [112] "XP23101_1" "XP23101_2" "XP23101_3"
## [115] "XP23101_4" "XP23101_5" "XP23101_6"
## [118] "XP23101_7" "XP23102_1_1" "XP23102_1_2"
## [121] "XP23102_1_3" "XP23102_1_4" "XP23102_1_5"
## [124] "XP23102_1_6" "XP23102_1_7" "XP23102_1_8"
## [127] "XP23102_1_9" "XP23102_1_10" "XP23102_1_11"
## [130] "XP23102_1_12" "XP23102_3_1" "XP23102_3_2"
## [133] "XP23102_3_3" "XP23102_3_4" "XP23102_3_5"
## [136] "XP23102_3_6" "XP23102_3_7" "XP23102_3_8"
## [139] "XP23102_3_9" "XP23102_3_10" "XP23102_3_11"
## [142] "XP23102_3_12" "p23150_sn" "p23200_sn"
## [145] "XP23201_1" "XP23201_2" "XP23201_3"
## [148] "XP23201_4" "XP23201_5" "XP23201_6"
## [151] "XP23202_1_1" "XP23202_1_2" "XP23202_1_3"
## [154] "XP23202_1_4" "XP23202_1_5" "XP23202_1_6"
## [157] "XP23202_1_7" "XP23202_31" "XP23202_32"
## [160] "XP23202_33" "XP23202_34" "XP23202_35"
## [163] "XP23202_36" "XP23202_37" "p23250_sn"
## [166] "XP23251_1" "XP23251_2" "XP23251_3"
## [169] "XP23251_4" "XP23251_5" "XP23251_6"
## [172] "XP23252_11" "XP23252_12" "XP23252_13"
## [175] "XP23252_14" "XP23252_15" "XP23252_16"
## [178] "XP23252_17" "XP23252_31" "XP23252_32"
## [181] "XP23252_33" "XP23252_34" "XP23252_35"
## [184] "XP23252_36" "p23300_sn" "p24000_1_sn"
## [187] "p24000_2_sn" "p24000_3_sn" "p24000_4_sn"
## [190] "p24000_5_sn" "p24000_13_sn" "p26002_sn"
## [193] "p26101_sn" "p26102_sn" "p26103_sn"
## [196] "p27000_sn" "p28100_sn" "XP28100_1"
## [199] "XP28100_2" "XP28100_3" "XP28100_4"
## [202] "XP28100_5" "XP28100_6" "XP28100_7"
## [205] "XP28100_9" "p28200_sn" "XP29000_1"
## [208] "XP29000_2" "XP29000_3" "XP29000_4"
## [211] "XP29000_5" "XP29000_6" "XP29000_9"
## [214] "p29000_esp" "p30000_sn" "p30000_1"
## [217] "SEMBRAMANUAL" "PULVEMOCHILA" "PULVEMOTO"
## [220] "ARADOS" "MOTOCULTOR" "MOTOSIERRA"
## [223] "CARRETA" "PICADORA" "TRACTOR"
## [226] "SEMBRADIRECTA" "PULVEARRASTRE" "COSEAUTOMOTRIZ"
## [229] "ACOPLADOS" "CAMIONES" "AERONAVES"
## [232] "TRACTORGPS" "SEMBRASENSOR" "PULVECONTROL"
## [235] "COSECHASENSOR" "DRONES"# Vector key → nombre de distrito para los 218 distritos (códigos DSTO_PROC)
dept_map <- c(
`0` = "Asunción",
`1` = "Concepción",
`2` = "San Pedro",
`3` = "Cordillera",
`4` = "Guairá",
`5` = "Caaguazú",
`6` = "Caazapá",
`7` = "Itapúa",
`8` = "Misiones",
`9` = "Paraguarí",
`10` = "Alto Paraná",
`11` = "Central",
`12` = "Ñeembucú",
`13` = "Amambay",
`14` = "Canindeyú",
`15` = "Presidente Hayes",
`16` = "Boquerón",
`17` = "Alto Paraguay"
)
distrito_labels <- c(
# ── Capital ────────────────────────────────────────────
`0` = "Asunción",
# ── Concepción (10x) ───────────────────────────────────
`101` = "Concepción",
`102` = "Belén",
`103` = "Horqueta",
`104` = "Loreto",
`105` = "San Carlos del Apa",
`106` = "San Lázaro",
`107` = "Yby Ya'u",
`108` = "Azote'y",
`109` = "Sargento José Félix López",
`110` = "San Alfredo",
`111` = "Paso Barreto",
`112` = "Arroyito",
`113` = "Paso Horqueta",
# ── San Pedro (20x) ───────────────────────────────────
`201` = "San Pedro del Ycuamandyyú",
`202` = "Antequera",
`203` = "Choré",
`204` = "General Elizardo Aquino",
`205` = "Itacurubí del Rosario",
`206` = "Lima",
`207` = "Nueva Germania",
`208` = "San Estanislao",
`209` = "San Pablo",
`210` = "Tacuatí",
`211` = "Unión",
`212` = "25 de Diciembre",
`213` = "Villa del Rosario",
`214` = "General Francisco Isidoro Resquín",
`215` = "Yataity del Norte",
`216` = "Guajayví",
`217` = "Capiibary",
`218` = "Santa Rosa del Aguaray",
`219` = "Yrybucuá",
`220` = "Liberación",
`221` = "San Vicente Pancholo",
`222` = "San José del Rosario",
# ── Cordillera (30x) ──────────────────────────────────
`301` = "Caacupé",
`302` = "Altos",
`303` = "Arroyos y Esteros",
`304` = "Atyrá",
`305` = "Caraguatay",
`306` = "Emboscada",
`307` = "Eusebio Ayala",
`308` = "Isla Pucú",
`309` = "Itacurubí de la Cordillera",
`310` = "Juan de Mena",
`311` = "Loma Grande",
`312` = "Mbocayaty del Yhaguy",
`313` = "Nueva Colombia",
`314` = "Piribebuy",
`315` = "Primero de Marzo",
`316` = "San Bernardino",
`317` = "Santa Elena",
`318` = "Tobatí",
`319` = "Valenzuela",
`320` = "San José Obrero",
# ── Guairá (40x) ──────────────────────────────────────
`401` = "Villarrica",
`402` = "Borja",
`403` = "Capitán Mauricio José Troche",
`404` = "Coronel Martínez",
`405` = "Félix Pérez Cardozo",
`406` = "General Eugenio A. Garay",
`407` = "Independencia",
`408` = "Itapé",
`409` = "Iturbe",
`410` = "José Fassardi",
`411` = "Mbocayaty",
`412` = "Natalicio Talavera",
`413` = "Ñumí",
`414` = "San Salvador",
`415` = "Yataity",
`416` = "Doctor Botrell",
`417` = "Paso Yobai",
`418` = "Tebicuary",
# ── Caaguazú (50x) ────────────────────────────────────
`501` = "Coronel Oviedo",
`502` = "Caaguazú",
`503` = "Carayaó",
`504` = "Doctor Cecilio Báez",
`505` = "Santa Rosa del Mbutuy",
`506` = "Doctor Juan Manuel Frutos",
`507` = "Repatriación",
`508` = "Nueva Londres",
`509` = "San Joaquín",
`510` = "San José de los Arroyos",
`511` = "Yhú",
`512` = "Doctor J. Eulogio Estigarribia",
`513` = "R.I. 3 Corrales",
`514` = "Raúl Arsenio Oviedo",
`515` = "José Domingo Ocampos",
`516` = "Mariscal Francisco Solano López",
`517` = "La Pastora",
`518` = "3 de Febrero",
`519` = "Simón Bolívar",
`520` = "Vaquería",
`521` = "Tembiaporá",
`522` = "Nueva Toledo",
# ── Caazapá (60x) ─────────────────────────────────────
`601` = "Caazapá",
`602` = "Abai",
`603` = "Buena Vista",
`604` = "Doctor Moisés S. Bertoni",
`605` = "General Higinio Morínigo",
`606` = "Maciel",
`607` = "San Juan Nepomuceno",
`608` = "Tavaí",
`609` = "Yegros",
`610` = "Yuty",
`611` = "3 de Mayo",
# ── Itapúa (70x) ──────────────────────────────────────
`701` = "Encarnación",
`702` = "Bella Vista",
`703` = "Cambyretá",
`704` = "Capitán Meza",
`705` = "Capitán Miranda",
`706` = "Nueva Alborada",
`707` = "Carmen del Paraná",
`708` = "Coronel Bogado",
`709` = "Carlos Antonio López",
`710` = "Natalio",
`711` = "Fram",
`712` = "General Artigas",
`713` = "General Delgado",
`714` = "Hohenau",
`715` = "Jesús",
`716` = "José Leandro Oviedo",
`717` = "Obligado",
`718` = "Mayor Julio Dionisio Otaño",
`719` = "San Cosme y Damián",
`720` = "San Pedro del Paraná",
`721` = "San Rafael del Paraná",
`722` = "Trinidad",
`723` = "Edelira",
`724` = "Tomás Romero Pereira",
`725` = "Alto Vera",
`726` = "La Paz",
`727` = "Yatytay",
`728` = "San Juan del Paraná",
`729` = "Pirapó",
`730` = "Itapúa Poty",
# ── Misiones (80x) ────────────────────────────────────
`801` = "San Juan Bautista de las Misiones",
`802` = "Ayolas",
`803` = "San Ignacio",
`804` = "San Miguel",
`805` = "San Patricio",
`806` = "Santa María",
`807` = "Santa Rosa",
`808` = "Santiago",
`809` = "Villa Florida",
`810` = "Yabebyry",
# ── Paraguarí (90x) ───────────────────────────────────
`901` = "Paraguarí",
`902` = "Acahay",
`903` = "Caapucú",
`904` = "Caballero",
`905` = "Carapeguá",
`906` = "Escobar",
`907` = "La Colmena",
`908` = "Mbuyapey",
`909` = "Pirayú",
`910` = "Quiindy",
`911` = "Quyquyhó",
`912` = "Roque González de Santa Cruz",
`913` = "Sapucái",
`914` = "Tebicuary-mi",
`915` = "Yaguarón",
`916` = "Ybycuí",
`917` = "Ybytymi",
`918` = "María Antonia",
# ── Alto Paraná (100x) ────────────────────────────────
`1001` = "Ciudad del Este",
`1002` = "Presidente Franco",
`1003` = "Domingo Martínez de Irala",
`1004` = "Doctor Juan León Mallorquín",
`1005` = "Hernandarias",
`1006` = "Itakyry",
`1007` = "Juan E. O'Leary",
`1008` = "Ñacunday",
`1009` = "Yguazú",
`1010` = "Los Cedrales",
`1011` = "Minga Guazú",
`1012` = "San Cristóbal",
`1013` = "Santa Rita",
`1014` = "Naranjal",
`1015` = "Santa Rosa del Monday",
`1016` = "Minga Porã",
`1017` = "Mbaracayú",
`1018` = "San Alberto",
`1019` = "Iruña",
`1020` = "Santa Fe del Paraná",
`1021` = "Tavapy",
`1022` = "Doctor Raúl Peña",
# ── Central (110x) ────────────────────────────────────
`1101` = "Areguá",
`1102` = "Capiatá",
`1103` = "Fernando de la Mora",
`1104` = "Guarambaré",
`1105` = "Itá",
`1106` = "Itauguá",
`1107` = "Lambaré",
`1108` = "Limpio",
`1109` = "Luque",
`1110` = "Mariano Roque Alonso",
`1111` = "Nueva Italia",
`1112` = "Ñemby",
`1113` = "San Antonio",
`1114` = "San Lorenzo",
`1115` = "Villa Elisa",
`1116` = "Villeta",
`1117` = "Ypacaraí",
`1118` = "Ypané",
`1119` = "J. Augusto Saldívar",
# ── Ñeembucú (120x) ───────────────────────────────────
`1201` = "Pilar",
`1202` = "Alberdi",
`1203` = "Cerrito",
`1204` = "Desmochados",
`1205` = "General José Eduvigis Díaz",
`1206` = "Guazú-Cuá",
`1207` = "Humaitá",
`1208` = "Isla Umbú",
`1209` = "Laureles",
`1210` = "Mayor José de Jesús Martínez",
`1211` = "Paso de la Patria",
`1212` = "San Juan Bautista de Ñeembucú",
`1213` = "Tacuaras",
`1214` = "Villa Franca",
`1215` = "Villa Oliva",
`1216` = "Villalbín",
# ── Amambay (130x) ───────────────────────────────────
`1301` = "Pedro Juan Caballero",
`1302` = "Bella Vista",
`1303` = "Capitán Bado",
`1304` = "Zanja Pytá",
`1305` = "Karapaí",
`1306` = "Cerro Corá",
# ── Canindeyú (140x) ─────────────────────────────────
`1401` = "Salto del Guairá",
`1402` = "Corpus Christi",
`1403` = "Curuguaty",
`1404` = "Ygatimí",
`1405` = "Itanará",
`1406` = "Ypejhú",
`1407` = "Francisco Caballero Álvarez",
`1408` = "Katueté",
`1409` = "La Paloma del Espíritu Santo",
`1410` = "Nueva Esperanza",
`1411` = "Yasy Cañy",
`1412` = "Ybyrarobana",
`1413` = "Yby Pytá",
`1414` = "Maracaná",
`1415` = "Puerto Adela",
`1416` = "Laurel",
# ── Presidente Hayes (150x) ───────────────────────────
`1502` = "Benjamín Aceval",
`1503` = "Puerto Pinasco",
`1504` = "Villa Hayes",
`1505` = "Nanawa",
`1506` = "José Falcón",
`1507` = "Teniente 1.º Manuel Irala Fernández",
`1508` = "Teniente Esteban Martínez",
`1509` = "General José María Bruguez",
`1510` = "Campo Aceval",
# ── Boquerón (160x) ──────────────────────────────────
`1602` = "Mariscal José Félix Estigarribia",
`1604` = "Filadelfia",
`1605` = "Loma Plata",
`1606` = "Boquerón",
# ── Alto Paraguay (170x) ─────────────────────────────
`1701` = "Fuerte Olimpo",
`1702` = "Puerto Casado",
`1704` = "Bahía Negra",
`1705` = "Carmelo Peralta"
)
# 5. Asignar etiqueta a cada fila según DSTO_PROC
can <- can %>%
mutate(
DSTO_PROC = as.integer(DSTO_PROC), # asegurarse de tipo entero
distrito_nombre = recode(DSTO_PROC, !!!distrito_labels)
)
can <- can %>%
mutate(
DPTO_PROC = as.integer(DPTO_PROC),
depto_nombre = recode(as.character(DPTO_PROC), !!!dept_map)
)5.2 Unipon de la base del CAN con la base geoespacial para mapear
# 1. Librerías (dplyr debe cargarse después de paquetes que exporten select S4)
library(sf)
library(stringr)
library(dplyr)
# 4. Preparar tu df 'can' para el join
can <- can %>%
mutate(dstocode = as.integer(DSTO_PROC))# 5. Hacer el join usando dplyr::select() para evitar conflictos
shp <- shp %>%
mutate(
raw_id2 = as.character(ID_2), # o ID_2 si así se llama
dstocode = raw_id2 %>%
str_remove("^PRY\\.") %>% # quita "PRY."
str_remove("_.*$") %>% # quita el "_" y lo que sigue
str_replace_all("\\.", "") %>% # elimina cualquier punto
as.integer() # convierte a integer (ej. 1613)
)
can <- can %>%
left_join(
shp %>%
sf::st_drop_geometry() %>%
dplyr::select(dstocode, NAME_1, NAME_2),
by = "dstocode"
) %>%
dplyr::rename(
departamento = NAME_1,
distrito = NAME_2
)
# 7. Vista rápida del resultado
head(can)## # A tibble: 6 × 241
## ID_FCA_ALEATORIO CLASIFICANEW DPTO_PROC DSTO_PROC tipoi tipofinca p5000
## <dbl> <dbl> <int> <int> <dbl> <dbl> <dbl>
## 1 1 11 15 1503 NA 3 NA
## 2 2 11 4 402 NA 1 1
## 3 3 11 7 723 NA 1 1
## 4 4 11 2 208 NA 1 1
## 5 5 11 2 221 NA 1 1
## 6 7 11 5 519 NA 1 1
## # ℹ 234 more variables: menores_total <dbl>, mayores_total <dbl>,
## # p7001_sn <dbl>, p7021p <dbl>, p7021p_nosabe <dbl>, p7033p <dbl>,
## # p7033p_nosabe <dbl>, p7221t <dbl>, p7221t_nosabe <dbl>, p7223t <dbl>,
## # p7223t_nosabe <dbl>, p8000_1 <dbl>, p8000_3 <dbl>, p9001 <dbl>,
## # p9000_total <dbl>, p10001 <dbl>, p10002 <dbl>, p10003 <dbl>, p10004 <dbl>,
## # p10005 <dbl>, p10006 <dbl>, p10007 <dbl>, p10008 <dbl>, p10009 <dbl>,
## # p10000_total <dbl>, p12000_sn <dbl>, p12000_total <dbl>, …6 Resultados
6.1 Adopción de agricultura de precisión
6.1.1 Variable dependiente y covariables
Las columnas del RDS se llaman según el listado entregado. Creamos las variables clave:
can <- can %>%
# Dependiente: adopción de agricultura de precisión ----------------------
mutate(agriprec = ifelse(p21016_sn == 1, 1, 0)) %>%
# Escala: superficie agrícola (p10001 ≈ tierra agrícola) -----------------
mutate(ln_supag = log(p10001 + 1)) %>%
# Capital humano: CLASIFICANEW (proxy, factor de clasificación nueva) ----
mutate(educ_cat = factor(CLASIFICANEW)) %>%
# Liquidez / crédito -----------------------------------------------------
mutate(credito = ifelse(p23200_sn == 1, 1, 0)) %>%
# Asistencia técnica -----------------------------------------------------
mutate(asistencia = ifelse(p23100_sn == 1, 1, 0)) %>%
# Uso de semilla híbrida -------------------------------------------------
mutate(sem_hibrida = ifelse(p24000_1_sn == 1, 1, 0)) %>%
# Comercialización -------------------------------------------------------
mutate(vende = ifelse(p28100_sn == 1, 1, 0)) %>%
# Riego ------------------------------------------------------------------
mutate(riego = ifelse(p21018_sn == 1, 1, 0))
# Diversificación cultivos: columnas p12000* y p14000* ---------------------
cult_cols <- grep("^(p12000|p14000).*_total$", names(can), value = TRUE)
can$cultivos <- rowSums(can[cult_cols] > 0, na.rm = TRUE)6.1.2 Agregación distrital
distrital <- can %>%
group_by(dstocode, DSTO_PROC) %>% # DSTO_PROC ≈ distrito procesado
summarise(across(c(agriprec, ln_supag, credito, asistencia,
sem_hibrida, vende, riego, cultivos),
mean, na.rm = TRUE),
educ_mode = as.numeric(names(sort(table(educ_cat), TRUE)[1]))) %>%
ungroup()## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(...)`.
## ℹ In group 1: `dstocode = 101` and `DSTO_PROC = 101`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
##
## # Previously
## across(a:b, mean, na.rm = TRUE)
##
## # Now
## across(a:b, \(x) mean(x, na.rm = TRUE))## `summarise()` has grouped output by 'dstocode'. You can override using the
## `.groups` argument.## # A tibble: 245 × 11
## dstocode DSTO_PROC agriprec ln_supag credito asistencia sem_hibrida vende
## <int> <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 101 101 0.00159 0.0344 0.162 0.169 0.265 0.934
## 2 102 102 0.0320 0.0210 0.226 0.262 0.0777 0.964
## 3 103 103 0.00183 0.101 0.189 0.170 0.0858 0.905
## 4 104 104 0.00361 0.0717 0.452 0.468 0.314 0.830
## 5 105 105 0 0.0400 0.1 0.0833 0 0.983
## 6 106 106 0 0 0.0922 0.0142 0 0.482
## 7 107 107 0.0420 0.0513 0.155 0.172 0.0736 0.856
## 8 108 108 0.0274 0.0393 0.145 0.187 0.292 0.952
## 9 109 109 0.00380 0.0702 0.152 0.120 0.0304 0.601
## 10 110 110 0.00312 0.0115 0.112 0.103 0.0374 0.333
## # ℹ 235 more rows
## # ℹ 3 more variables: riego <dbl>, cultivos <dbl>, educ_mode <dbl>6.1.3 Unir con shapefile
## Simple feature collection with 44 features and 24 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -58.65899 ymin: -27.38545 xmax: -54.35558 ymax: -22.07954
## Geodetic CRS: WGS 84
## First 10 features:
## ID_0 COUNTRY NAME_1 NL_NAME_1 ID_2
## 1 PRY Paraguay Alto Paraná <NA> PRY.2.10_1
## 2 PRY Paraguay Alto Paraná <NA> PRY.2.11_1
## 3 PRY Paraguay Alto Paraná <NA> PRY.2.12_1
## 4 PRY Paraguay Alto Paraná <NA> PRY.2.13_1
## 5 PRY Paraguay Alto Paraná <NA> PRY.2.14_1
## 6 PRY Paraguay Alto Paraná <NA> PRY.2.15_1
## 7 PRY Paraguay Alto Paraná <NA> PRY.2.16_1
## 8 PRY Paraguay Alto Paraná <NA> PRY.2.17_1
## 9 PRY Paraguay Alto Paraná <NA> PRY.2.18_1
## 10 PRY Paraguay Caaguazú <NA> PRY.6.10_1
## NAME_2 VARNAME_2 NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2
## 1 Minga Porá <NA> <NA> Distrito District <NA>
## 2 Ñacunday <NA> <NA> Distrito District <NA>
## 3 Naranjal <NA> <NA> Distrito District <NA>
## 4 Presidente Franco <NA> <NA> Distrito District <NA>
## 5 San Alberto <NA> <NA> Distrito District <NA>
## 6 San Cristóbal <NA> <NA> Distrito District <NA>
## 7 Santa Rita <NA> <NA> Distrito District <NA>
## 8 Santa Rosa del Monday <NA> <NA> Distrito District <NA>
## 9 Yguazú <NA> <NA> Distrito District <NA>
## 10 Mariscal Francisco Solano López <NA> <NA> Distrito District <NA>
## HASC_2 raw_id2 dstocode DSTO_PROC agriprec ln_supag credito
## 1 PY.AA.MP PRY.2.10_1 210 210 0.0553580498 0.084298839 0.08684611
## 2 PY.AA.NC PRY.2.11_1 211 211 0.0010204082 0.001534773 0.01020408
## 3 PY.AA.NR PRY.2.12_1 212 212 0.0079946702 0.103158619 0.03131246
## 4 <NA> PRY.2.13_1 213 213 0.0000000000 0.023697256 0.21992819
## 5 PY.AA.SA PRY.2.14_1 214 214 0.0122887865 0.103408848 0.45238095
## 6 PY.AA.SC PRY.2.15_1 215 215 0.0009514748 0.041745015 0.04947669
## 7 PY.AA.SR PRY.2.16_1 216 216 0.0100953449 0.340068725 0.16993831
## 8 PY.AA.SM PRY.2.17_1 217 217 0.0035248502 0.037502102 0.17729996
## 9 PY.AA.YG PRY.2.18_1 218 218 0.0157716861 0.067034407 0.15921893
## 10 PY.CG.ML PRY.6.10_1 610 610 0.0029254022 0.046777553 0.04729400
## asistencia sem_hibrida vende riego cultivos educ_mode
## 1 0.056373794 0.38649060 0.8928390 0.0020314881 0.9740985 11
## 2 0.009183673 0.04183673 0.4479592 0.0000000000 0.5051020 11
## 3 0.033977348 0.11125916 0.9600266 0.0059960027 0.8500999 11
## 4 0.269299820 0.17414722 0.8285458 0.0000000000 0.5421903 11
## 5 0.421658986 0.23809524 0.9301075 0.0069124424 1.0721966 11
## 6 0.026641294 0.03996194 0.5889629 0.0009514748 0.6841104 11
## 7 0.173022995 0.21508693 0.9610208 0.0033651150 1.1365676 11
## 8 0.156150864 0.27634826 0.8854424 0.0088121255 0.8741628 11
## 9 0.179872324 0.25647766 0.8460383 0.0067592940 0.8434097 11
## 10 0.032179425 0.27888835 0.5738664 0.0019502682 0.9819600 11
## geometry
## 1 MULTIPOLYGON (((-55.02682 -...
## 2 MULTIPOLYGON (((-54.65066 -...
## 3 MULTIPOLYGON (((-54.99545 -...
## 4 MULTIPOLYGON (((-54.59611 -...
## 5 MULTIPOLYGON (((-54.37466 -...
## 6 MULTIPOLYGON (((-55.34438 -...
## 7 MULTIPOLYGON (((-55.01954 -...
## 8 MULTIPOLYGON (((-55.22297 -...
## 9 MULTIPOLYGON (((-55.0987 -2...
## 10 MULTIPOLYGON (((-55.02353 -...6.1.4 Matriz de pesos (k = 5 vecinos más cercanos)
# Cargar los paquetes necesarios
library(sf)
library(spdep)
library(class)
# Suponiendo que 'shp_d' es un objeto sf (un shapefile cargado con sf)
coords <- st_coordinates(st_centroid(shp_d))## Warning: st_centroid assumes attributes are constant over geometries# Calcular los vecinos más cercanos con k = 5
knn5 <- knearneigh(coords, k = 5)
# Convertir los índices de vecinos a un objeto de clase 'nb'
nb_list <- knn2nb(knn5)## Warning in knn2nb(knn5): neighbour object has 2 sub-graphs6.1.5 Moran I global
## Moran I statistic Expectation Variance
## 0.040163665 -0.023255814 0.0016049576.1.6 Modelo espacial Durbin (SDM)
form_sdm <- agriprec ~ ln_supag + educ_mode + credito + asistencia +
sem_hibrida + vende + riego + cultivos
sdm <- lagsarlm(form_sdm, data = shp_d, listw = listw, Durbin = TRUE, type = "mixed")## Warning in lagsarlm(form_sdm, data = shp_d, listw = listw, Durbin = TRUE, :
## Aliased variables found: educ_mode lag.educ_mode##
## Call:lagsarlm(formula = form_sdm, data = shp_d, listw = listw, Durbin = TRUE,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0603298 -0.0129453 0.0020318 0.0124543 0.1480843
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## (2 not defined because of singularities)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.128846 0.108978 -1.1823 0.2371
## ln_supag -0.017593 0.084263 -0.2088 0.8346
## educ_mode NA NA NA NA
## credito 0.092002 0.122902 0.7486 0.4541
## asistencia -0.014909 0.118984 -0.1253 0.9003
## sem_hibrida 0.042909 0.056460 0.7600 0.4473
## vende 0.017665 0.042370 0.4169 0.6767
## riego -0.040033 0.078799 -0.5080 0.6114
## cultivos 0.019985 0.022526 0.8872 0.3750
## lag.ln_supag 0.343545 0.281177 1.2218 0.2218
## lag.educ_mode NA NA NA NA
## lag.credito -0.331792 0.304409 -1.0900 0.2757
## lag.asistencia 0.319760 0.338647 0.9442 0.3451
## lag.sem_hibrida -0.015740 0.135238 -0.1164 0.9073
## lag.vende 0.133811 0.095645 1.3990 0.1618
## lag.riego -0.271691 0.264056 -1.0289 0.3035
## lag.cultivos -0.040151 0.046311 -0.8670 0.3860
##
## Rho: 0.06393, LR test value: 0.035859, p-value: 0.84981
## Asymptotic standard error: 0.23735
## z-value: 0.26935, p-value: 0.78766
## Wald statistic: 0.072549, p-value: 0.78766
##
## Log likelihood: 91.90396 for mixed model
## ML residual variance (sigma squared): 0.00089746, (sigma: 0.029958)
## Number of observations: 44
## Number of parameters estimated: 17
## AIC: -149.81, (AIC for lm: -151.77)
## LM test for residual autocorrelation
## test value: 1.3103, p-value: 0.25235## Impact measures (mixed, exact):
## Direct Indirect Total
## ln_supag -0.01425875 0.36247143 0.34821268
## credito 0.08882854 -0.34499444 -0.25616590
## asistencia -0.01180449 0.33747508 0.32567059
## sem_hibrida 0.04278263 -0.01375842 0.02902421
## vende 0.01897957 0.14284245 0.16182201
## riego -0.04270323 -0.29031002 -0.33301325
## cultivos 0.01960674 -0.04114926 -0.02154252Discusión:
El análisis sugiere que las variables incluidas en el modelo no tienen un impacto estadísticamente significativo sobre la variable dependiente, ya que la mayoría de los valores p son altos. El modelo presenta problemas con variables como educ_mode que tienen colinealidad. La autocorrelación espacial en el conjunto de datos es leve, y los impactos directos e indirectos también son pequeños.
Es posible que necesites revisar los datos o incluir variables adicionales que podrían mejorar el ajuste del modelo o considerar la posibilidad de usar otro tipo de modelo si hay más complejidades en los datos que no se están capturando adecuadamente.
6.1.7 Mapa coroplético
library(dplyr)
library(tmap)
library(RColorBrewer)
# 1) Reconstruye el sf con TODOS los distritos y NA→0 solo para el plot
shp_all <- shp %>%
left_join(distrital, by = "dstocode") %>%
mutate(agriprec0 = coalesce(agriprec, 0))
# 2) Define cortes idénticos a antes
max_val <- max(shp_all$agriprec0, na.rm = TRUE)
breaks <- c(-1e-6, 0.0000001, 0.002, 0.004, 0.010, max_val)
# 3) Prepara paleta: primera clase blanca + 4 tonos de “Blues”
pal <- c(
"white",
brewer.pal(4, "Blues")
)
# 4) Dibuja
tmap_mode("plot")## ℹ tmap mode set to "plot".tm_shape(shp_all) +
tm_polygons(
"agriprec0",
style = "fixed",
breaks = breaks,
palette = pal,
border.col= "grey50",
lwd = 0.3,
title = "Prop. adopción AP"
) +
tm_layout(legend.outside = TRUE)##
## ── tmap v3 code detected ───────────────────────────────────────────────────────
## [v3->v4] `tm_polygons()`: instead of `style = "fixed"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'breaks', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'[v3->v4] `tm_polygons()`: use 'fill' for the fill color of polygons/symbols
## (instead of 'col'), and 'col' for the outlines (instead of 'border.col').[v3->v4] `tm_polygons()`: migrate the argument(s) related to the legend of the
## visual variable `fill` namely 'title' to 'fill.legend = tm_legend(<HERE>)'
6.1.8 Modelo alternativo
library(sf)
library(sp)
library(GWmodel)
# 1) Partimos de tu shp_d (un sf con datos ya filtrados)
shp_sp <- as(shp_d, "Spatial") # convierte a SpatialPointsDataFrame
coords_sp<- coordinates(shp_sp) # extrae la matriz n×2
# 2) Matriz de distancias (longlat=FALSE porque ya está en UTM)
dMat <- gw.dist(dp.locat = coords_sp, longlat = FALSE)
# 3) Selección de bandwidth (CV, kernel bisquare, fijo)
bw_gwr <- bw.gwr(
formula = form_sdm, # tu fórmula: agriprec ~ ...
data = shp_sp,
approach = "CV",
kernel = "bisquare",
adaptive = FALSE,
dMat = dMat
)## Fixed bandwidth: 3.122368 CV score: 1.733227e+25
## Fixed bandwidth: 1.930116 CV score: 2.45474e+25
## Fixed bandwidth: 3.859221 CV score: 3.084886e+25
## Fixed bandwidth: 2.666968 CV score: 3.438473e+23
## Fixed bandwidth: 2.385515 CV score: 2.020429e+26
## Fixed bandwidth: 2.840915 CV score: 3.104956e+24
## Fixed bandwidth: 2.559463 CV score: 2.068696e+24
## Fixed bandwidth: 2.73341 CV score: 3.271124e+24
## Fixed bandwidth: 2.625905 CV score: 2.318367e+25
## Fixed bandwidth: 2.692347 CV score: 5.339226e+23
## Fixed bandwidth: 2.651283 CV score: 1.783002e+26
## Fixed bandwidth: 2.676662 CV score: 5.039946e+24
## Fixed bandwidth: 2.660977 CV score: 3.485697e+23
## Fixed bandwidth: 2.670671 CV score: 4.975707e+27
## Fixed bandwidth: 2.66468 CV score: 5.898951e+25
## Fixed bandwidth: 2.668382 CV score: 1.022687e+24
## Fixed bandwidth: 2.666094 CV score: 4.459832e+25
## Fixed bandwidth: 2.667508 CV score: 4.352202e+27
## Fixed bandwidth: 2.666634 CV score: 3.344616e+24
## Fixed bandwidth: 2.667174 CV score: 3.969739e+22
## Fixed bandwidth: 2.667302 CV score: 3.253216e+25
## Fixed bandwidth: 2.667096 CV score: 2.76917e+24
## Fixed bandwidth: 2.667223 CV score: 2.080793e+23
## Fixed bandwidth: 2.667144 CV score: 7.219613e+24# 4) Ajuste del GWR
gwr_res <- gwr.basic(
formula = form_sdm,
data = shp_sp,
bw = bw_gwr,
kernel = "bisquare",
adaptive = FALSE,
dMat = dMat
)library(sf)
library(dplyr)
library(tmap)
library(RColorBrewer)
# — 1) Extraer los coeficientes del objeto GWR (SpatialDataFrame) —
gwr_df <- data.frame(
dstocode = shp_d$dstocode,
beta_ln_supag = gwr_res$SDF@data$ln_supag
)
gwr_df## dstocode beta_ln_supag
## 1 210 0.0066222392
## 2 211 0.0470489951
## 3 212 0.0283602590
## 4 213 0.0346216118
## 5 214 0.0098244101
## 6 215 0.0313762152
## 7 216 0.0401496048
## 8 217 0.0398877015
## 9 218 0.0340983229
## 10 610 0.0306019613
## 11 611 0.0042891066
## 12 710 -0.1114501622
## 13 910 0.0092686530
## 14 911 0.0171014531
## 15 912 -0.0176901426
## 16 913 -0.0798846884
## 17 914 -0.0216907909
## 18 915 -0.0088804868
## 19 916 -0.0144418716
## 20 917 -0.0742270288
## 21 918 0.0005098001
## 22 101 -0.0617374697
## 23 102 -0.2923641133
## 24 103 -0.1581290939
## 25 104 -0.1786297262
## 26 105 -0.5156250000
## 27 106 -0.4140625000
## 28 107 -0.2120330155
## 29 111 0.0099912931
## 30 112 0.0422048765
## 31 113 0.0081616524
## 32 1111 0.0264522027
## 33 1116 0.0047479507
## 34 1117 0.0320015326
## 35 1119 0.0167874192
## 36 1210 -0.0089936380
## 37 1211 -0.0238795483
## 38 1212 -0.0142509358
## 39 1213 -0.0439856875
## 40 1214 -0.0566064327
## 41 1215 -0.0407452101
## 42 1216 -0.0257248907
## 43 1410 -0.2584617783
## 44 1510 -0.2781749158# — 2) Reconstruir un sf con todos los distritos y unir el GWR —
shp_all_gwr <- shp %>%
# asegurar dstocode en el shapefile original
mutate(
raw_id2 = as.character(ID_2),
dstocode = raw_id2 %>%
str_remove("^PRY\\.") %>%
str_remove("_.*$") %>%
str_replace_all("\\.", "") %>%
as.integer()
) %>%
# left join para traer beta_ln_supag; los ausentes quedan NA
left_join(gwr_df, by = "dstocode")
shp_all_gwr## Simple feature collection with 218 features and 15 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -62.64652 ymin: -27.60586 xmax: -54.25863 ymax: -19.29137
## Geodetic CRS: WGS 84
## First 10 features:
## ID_0 COUNTRY NAME_1 NL_NAME_1 ID_2 NAME_2
## 1 PRY Paraguay Alto Paraguay <NA> PRY.1.1_1 Fuerte Olimpo
## 2 PRY Paraguay Alto Paraguay <NA> PRY.1.2_1 La Victoria
## 3 PRY Paraguay Alto Paraguay <NA> PRY.1.3_1 Mayor Pablo Lagerenza
## 4 PRY Paraguay Alto Paraná <NA> PRY.2.1_1 Ciudad del Este
## 5 PRY Paraguay Alto Paraná <NA> PRY.2.2_1 Doctor Juan León Mallorquín
## 6 PRY Paraguay Alto Paraná <NA> PRY.2.3_1 Domingo Martínez de Irala
## 7 PRY Paraguay Alto Paraná <NA> PRY.2.4_1 Hernandarias
## 8 PRY Paraguay Alto Paraná <NA> PRY.2.5_1 Itakyry
## 9 PRY Paraguay Alto Paraná <NA> PRY.2.6_1 Juan Emilio O'Leary
## 10 PRY Paraguay Alto Paraná <NA> PRY.2.7_1 Los Cedrales
## VARNAME_2 NL_NAME_2 TYPE_2 ENGTYPE_2 CC_2 HASC_2 raw_id2
## 1 <NA> <NA> Distrito District <NA> PY.AG.FO PRY.1.1_1
## 2 Puerto Casado <NA> Distrito District <NA> PY.AG.LV PRY.1.2_1
## 3 <NA> <NA> Distrito District <NA> PY.AG.ML PRY.1.3_1
## 4 <NA> <NA> Distrito District <NA> PY.AA.CE PRY.2.1_1
## 5 <NA> <NA> Distrito District <NA> PY.AA.DM PRY.2.2_1
## 6 <NA> <NA> Distrito District <NA> <NA> PRY.2.3_1
## 7 <NA> <NA> Distrito District <NA> <NA> PRY.2.4_1
## 8 Itaquyry|Ytakyry <NA> Distrito District <NA> PY.AA.IT PRY.2.5_1
## 9 <NA> <NA> Distrito District <NA> PY.AA.JO PRY.2.6_1
## 10 <NA> <NA> Distrito District <NA> <NA> PRY.2.7_1
## dstocode beta_ln_supag geometry
## 1 11 NA MULTIPOLYGON (((-59.14466 -...
## 2 12 NA MULTIPOLYGON (((-57.87265 -...
## 3 13 NA MULTIPOLYGON (((-59.1192 -1...
## 4 21 NA MULTIPOLYGON (((-54.60418 -...
## 5 22 NA MULTIPOLYGON (((-55.34259 -...
## 6 23 NA MULTIPOLYGON (((-54.62878 -...
## 7 24 NA MULTIPOLYGON (((-54.46827 -...
## 8 25 NA MULTIPOLYGON (((-54.97109 -...
## 9 26 NA MULTIPOLYGON (((-55.29816 -...
## 10 27 NA MULTIPOLYGON (((-54.65395 -...# — 3) Mapa de coeficiente local, con NAs en gris claro —
tmap_mode("plot")
tm_shape(shp_all_gwr) +
tm_polygons(
"beta_ln_supag",
style = "quantile",
n = 5,
palette = brewer.pal(5, "RdBu"),
title = expression(beta[ln~supag]),
colorNA = "lightgrey", # relleno de los sin estimación
textNA = "No estimado", # etiqueta en la leyenda
border.col= "grey50", # contorno de cada distrito
lwd = 0.3
) +
tm_layout(
main.title = "Coeficiente local de ln_supag (GWR)",
legend.outside = TRUE
)
6.2 Mecanización
# Paquetes --------------------------------------------------------------
# 1) Carga de datos --------------------------------------------------------
data_dir <- "G:/Mi unidad/CAN2022/"
geo_dir <- file.path(data_dir, "paraguay")
shp <- st_read(file.path(geo_dir, "gadm40_PRY_2.shp")) %>% st_make_valid()## Reading layer `gadm40_PRY_2' from data source
## `G:\Mi unidad\CAN2022\paraguay\gadm40_PRY_2.shp' using driver `ESRI Shapefile'
## Simple feature collection with 218 features and 12 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: -62.64652 ymin: -27.60586 xmax: -54.25863 ymax: -19.29137
## Geodetic CRS: WGS 84can <- read_delim(
file.path(data_dir, "BASE_CAN2022_DATOS_GENERALES_P2.csv"),
delim = ",", show_col_types = FALSE
)
# 2) Índice de mecanización por finca -------------------------------------
mech_vars <- c(
"TRACTORGPS", "TRACTOR", "MOTOSIERRA", "MOTOCULTOR",
"PICADORA", "CARRETA", "CAMIONES", "AERONAVES",
"SEMBRASENSOR", "SEMBRAMANUAL", "SEMBRADIRECTA",
"PULVECONTROL", "PULVEMOCHILA", "PULVEARRASTRE", "DRONES"
)can <- can %>%
# Asegurar que los binaries estén en 0/1
mutate(across(all_of(mech_vars), ~ as.numeric(. == 1))) %>%
# Índice de mecanización: proporción de prácticas usadas
rowwise() %>%
mutate(
mech_index = mean(c_across(all_of(mech_vars)), na.rm = TRUE),
# Covariables de control
ln_supag = log(p10001 + 1),
credito = p23200_sn,
asistencia = p23100_sn,
semilla_hib = p24000_1_sn,
abono_org = p24000_13_sn,
riego = p21018_sn
) %>%
ungroup() %>%
mutate(DSTO_PROC = as.integer(DSTO_PROC))# 3) Agregación distrital --------------------------------------------------
distr_mech <- can %>%
group_by(DSTO_PROC) %>%
summarise(
mech_prop = mean(mech_index, na.rm = TRUE),
ln_supag = mean(ln_supag, na.rm = TRUE),
credito = mean(credito, na.rm = TRUE),
asistencia = mean(asistencia, na.rm = TRUE),
semilla_hib = mean(semilla_hib, na.rm = TRUE),
abono_org = mean(abono_org, na.rm = TRUE),
riego = mean(riego, na.rm = TRUE),
.groups = "drop"
)
distr_mech## # A tibble: 245 × 8
## DSTO_PROC mech_prop ln_supag credito asistencia semilla_hib abono_org riego
## <int> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 101 0.0433 0.0344 5.21 5.16 4.68 5.00 5.91
## 2 102 0.0538 0.0210 4.87 4.69 5.61 5.70 5.96
## 3 103 0.0410 0.101 5.08 5.17 5.59 5.65 5.98
## 4 104 0.0773 0.0717 3.76 3.68 4.45 5.77 5.89
## 5 105 0.0211 0.0400 5.55 5.63 6.05 5.97 6.05
## 6 106 0.0180 0 5.54 5.95 6 6 6
## 7 107 0.0349 0.0513 5.24 5.14 5.65 5.80 5.99
## 8 108 0.0526 0.0393 5.29 5.07 4.55 5.16 5.97
## 9 109 0.0393 0.0702 5.26 5.41 5.85 5.83 6
## 10 110 0.0114 0.0115 5.46 5.50 5.84 5.86 6
## # ℹ 235 more rows# 4) Unir con shapefile (todos los distritos) ------------------------------
shp_all <- shp %>%
mutate(
raw_id2 = as.character(ID_2),
dstocode = raw_id2 %>%
str_remove("^PRY\\.") %>%
str_remove("_.*$") %>%
str_replace_all("\\.", "") %>%
as.integer()
) %>%
left_join(distr_mech, by = c("dstocode" = "DSTO_PROC"))# 5) Datos para análisis (sin NA en la dependiente) -----------------------
shp_ana <- shp_all %>% filter(!is.na(mech_prop))
# 6) Matriz de pesos (k = 5 vecinos más cercanos) ------------------------
coords <- st_coordinates(st_centroid(shp_ana))
knn5 <- knearneigh(coords, k = 5)
nb <- knn2nb(knn5)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# 7) Moran’s I -------------------------------------------------------------
print(moran.test(shp_ana$mech_prop, listw, zero.policy = TRUE))##
## Moran I test under randomisation
##
## data: shp_ana$mech_prop
## weights: listw
##
## Moran I statistic standard deviate = 0.73215, p-value = 0.232
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.03554115 -0.02325581 0.00644933# 8) SDM (Spatial Durbin) -------------------------------------------------
form_mech <- mech_prop ~ ln_supag + credito + asistencia +
semilla_hib + abono_org + riego
sdm_mech <- lagsarlm(
formula = form_mech,
data = shp_ana,
listw = listw,
Durbin = TRUE,
type = "mixed",
zero.policy = TRUE
)
summary(sdm_mech)##
## Call:lagsarlm(formula = form_mech, data = shp_ana, listw = listw,
## Durbin = TRUE, type = "mixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0337824 -0.0089612 0.0021560 0.0083216 0.0335333
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.1295640 0.1485876 -0.8720 0.38322
## ln_supag 0.0115541 0.0355172 0.3253 0.74495
## credito 0.0173381 0.0114552 1.5136 0.13014
## asistencia -0.0182110 0.0117561 -1.5491 0.12137
## semilla_hib -0.0304586 0.0068798 -4.4272 9.545e-06
## abono_org 0.0092040 0.0044721 2.0581 0.03958
## riego 0.0092883 0.0075168 1.2357 0.21658
## lag.ln_supag -0.2826114 0.1178457 -2.3981 0.01648
## lag.credito -0.0535892 0.0361012 -1.4844 0.13770
## lag.asistencia 0.0624600 0.0394722 1.5824 0.11356
## lag.semilla_hib -0.0187537 0.0157971 -1.1872 0.23517
## lag.abono_org -0.0013221 0.0104832 -0.1261 0.89964
## lag.riego 0.0522749 0.0247377 2.1132 0.03459
##
## Rho: 0.032716, LR test value: 0.012194, p-value: 0.91207
## Asymptotic standard error: 0.23883
## z-value: 0.13698, p-value: 0.89104
## Wald statistic: 0.018765, p-value: 0.89104
##
## Log likelihood: 120.9908 for mixed model
## ML residual variance (sigma squared): 0.00023934, (sigma: 0.015471)
## Number of observations: 44
## Number of parameters estimated: 15
## AIC: -211.98, (AIC for lm: -213.97)
## LM test for residual autocorrelation
## test value: 0.77101, p-value: 0.3799## Impact measures (mixed, exact):
## Direct Indirect Total
## ln_supag 0.010179603 -0.290404784 -0.280225182
## credito 0.017079917 -0.054557144 -0.037477228
## asistencia -0.017909726 0.063655351 0.045745624
## semilla_hib -0.030554795 -0.020321975 -0.050876770
## abono_org 0.009199014 -0.001050587 0.008148427
## riego 0.009544352 0.054101096 0.063645448# 9) GWR -------------------------------------------------------------------
shp_sp <- as(shp_ana, "Spatial")
coords_sp<- coordinates(shp_sp)
dMat <- gw.dist(dp.locat = coords_sp, longlat = FALSE)
bw_mech <- bw.gwr(
formula = form_mech,
data = shp_sp,
dMat = dMat,
kernel = "bisquare",
adaptive = FALSE
)## Fixed bandwidth: 3.122368 CV score: 0.0537551
## Fixed bandwidth: 1.930116 CV score: 0.200666
## Fixed bandwidth: 3.859221 CV score: 0.05227186
## Fixed bandwidth: 4.314621 CV score: 0.05039002
## Fixed bandwidth: 4.596073 CV score: 0.04938544
## Fixed bandwidth: 4.770021 CV score: 0.04881354
## Fixed bandwidth: 4.877526 CV score: 0.04849189
## Fixed bandwidth: 4.943968 CV score: 0.04830765
## Fixed bandwidth: 4.985031 CV score: 0.04819923
## Fixed bandwidth: 5.01041 CV score: 0.04813427gwr_mech <- gwr.basic(
formula = form_mech,
data = shp_sp,
bw = bw_mech,
dMat = dMat,
kernel = "bisquare",
adaptive = FALSE
)
# Extraer coeficiente local de ln_supag
gwr_coef <- data.frame(
dstocode = shp_ana$dstocode,
beta_ln_supag = gwr_mech$SDF@data$ln_supag
)
shp_all <- shp_all %>% left_join(gwr_coef, by = "dstocode")# 10) Mapas ---------------------------------------------------------------
tmap_mode("plot")
# 10.1 Coroplético de mecanización
tm_shape(shp_all) +
tm_polygons(
"mech_prop",
style = "quantile",
n = 5,
palette = brewer.pal(5, "Purples"),
title = "Prop. mecanización",
colorNA = "lightgrey",
textNA = "Sin dato",
border.col= "grey50", lwd = 0.3
) +
tm_layout(legend.outside = TRUE,
main.title = "Intensidad de mecanización por distrito")
# 10.2 Coeficiente local ln_supag (GWR)
tm_shape(shp_all) +
tm_polygons(
"beta_ln_supag",
style = "quantile",
n = 5,
palette = brewer.pal(5, "RdBu"),
title = expression(beta[ln~supag]),
colorNA = "lightgrey",
textNA = "No estimado",
border.col= "grey50", lwd = 0.3
) +
tm_layout(legend.outside = TRUE,
main.title = "Coef. local de escala (mecanización GWR)")
6.3 Diversificación ganadera
# Paquetes ----------------------------------------------------------------
library(tidyverse)
library(sf)
library(spdep)
library(spatialreg)
library(sp)
library(GWmodel)
library(tmap)
library(RColorBrewer)
library(stringr)
# 2) Crear índice de diversificación ganado -------------------------------
can <- can %>%
mutate(DSTO_PROC = as.integer(DSTO_PROC)) %>%
mutate(
vacunos = p19100_sn, # 1/0
ovinos = p19300_sn,
caprinos = p19400_sn,
porcinos = p19500_sn,
aves = p19700_sn,
diversidad = vacunos + ovinos + caprinos + porcinos + aves,
# covariables de control
ln_supag = log(p10001 + 1),
credito = p23200_sn,
asistencia = p23100_sn,
riego = p21018_sn,
semilla_hib = p24000_1_sn
)
# 3) Agregación distrital --------------------------------------------------
distr_div <- can %>%
group_by(DSTO_PROC) %>%
summarise(
div_med = mean(diversidad, na.rm = TRUE),
ln_supag = mean(ln_supag, na.rm = TRUE),
credito = mean(credito, na.rm = TRUE),
asistencia = mean(asistencia, na.rm = TRUE),
riego = mean(riego, na.rm = TRUE),
semilla_hib= mean(semilla_hib, na.rm = TRUE),
.groups = "drop"
)# Paquetes necesarios
library(sf)
library(dplyr)
library(stringr)
library(spdep)
library(spatialreg)
# 4) Unir con shapefile (todos los distritos)
shp_all <- shp %>%
mutate(
raw_id2 = as.character(ID_2),
dstocode = raw_id2 %>%
str_remove("^PRY\\.") %>%
str_remove("_.*$") %>%
str_replace_all("\\.", "") %>%
as.integer()
) %>%
left_join(
distr_div,
by = c("dstocode" = "DSTO_PROC") # ojo: 'DSTO_PROC' viene de distr_div
)
# 5) Filtrar para análisis ------------------------------------------------
shp_ana <- shp_all %>%
filter(!is.na(div_med))
# 6) Matriz de pesos (k = 5 vecinos más cercanos) --------------------------
coords <- st_coordinates(st_centroid(shp_ana))
knn5 <- knearneigh(coords, k = 5)
nb <- knn2nb(knn5)
listw <- nb2listw(nb, style = "W", zero.policy = TRUE)
# 7) Moran’s I -------------------------------------------------------------
print(moran.test(shp_ana$div_med, listw, zero.policy = TRUE))##
## Moran I test under randomisation
##
## data: shp_ana$div_med
## weights: listw
##
## Moran I statistic standard deviate = 3.107, p-value = 0.0009449
## alternative hypothesis: greater
## sample estimates:
## Moran I statistic Expectation Variance
## 0.21822032 -0.02325581 0.00604028# 8) SDM (Spatial Durbin) -------------------------------------------------
form_div <- div_med ~ ln_supag + credito + asistencia + riego + semilla_hib
sdm_div <- lagsarlm(
formula = form_div,
data = shp_ana,
listw = listw,
Durbin = TRUE,
type = "mixed",
zero.policy = TRUE
)
summary(sdm_div)##
## Call:
## lagsarlm(formula = form_div, data = shp_ana, listw = listw, Durbin = TRUE,
## type = "mixed", zero.policy = TRUE)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.88128 -0.89636 -0.13858 1.10593 4.16946
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 43.91890 16.32852 2.6897 0.007152
## ln_supag 1.37517 3.52494 0.3901 0.696444
## credito -1.66320 1.02094 -1.6291 0.103294
## asistencia 1.05565 1.07758 0.9797 0.327258
## riego -2.31834 0.71763 -3.2305 0.001236
## semilla_hib -1.00378 0.52730 -1.9036 0.056962
## lag.ln_supag 3.54498 11.33591 0.3127 0.754493
## lag.credito 3.52600 2.60200 1.3551 0.175382
## lag.asistencia -2.38280 2.97479 -0.8010 0.423134
## lag.riego -3.77975 2.46287 -1.5347 0.124859
## lag.semilla_hib 1.21697 1.28732 0.9453 0.344481
##
## Rho: 0.36356, LR test value: 3.556, p-value: 0.059331
## Asymptotic standard error: 0.17375
## z-value: 2.0925, p-value: 0.036398
## Wald statistic: 4.3784, p-value: 0.036398
##
## Log likelihood: -82.23469 for mixed model
## ML residual variance (sigma squared): 2.4021, (sigma: 1.5499)
## Number of observations: 44
## Number of parameters estimated: 13
## AIC: 190.47, (AIC for lm: 192.03)
## LM test for residual autocorrelation
## test value: 6.0947, p-value: 0.013559## Impact measures (mixed, exact):
## Direct Indirect Total
## ln_supag 1.6740494 6.056707 7.7307569
## credito -1.4473440 4.374256 2.9269121
## asistencia 0.9079451 -2.993217 -2.0852718
## riego -2.6599023 -6.921699 -9.5816011
## semilla_hib -0.9408203 1.275792 0.3349718# 9) GWR -------------------------------------------------------------------
shp_sp <- as(shp_ana, "Spatial")
coords_sp <- coordinates(shp_sp)
dMat <- gw.dist(dp.locat = coords_sp, longlat = FALSE)
bw_div <- bw.gwr(
formula = form_div,
data = shp_sp,
dMat = dMat,
kernel = "bisquare",
adaptive = FALSE
)## Fixed bandwidth: 3.122368 CV score: 603.0396
## Fixed bandwidth: 1.930116 CV score: 3822.033
## Fixed bandwidth: 3.859221 CV score: 576.5295
## Fixed bandwidth: 4.314621 CV score: 564.8775
## Fixed bandwidth: 4.596073 CV score: 556.9044
## Fixed bandwidth: 4.770021 CV score: 550.6573
## Fixed bandwidth: 4.877526 CV score: 547.2613
## Fixed bandwidth: 4.943968 CV score: 545.4041
## Fixed bandwidth: 4.985031 CV score: 544.339
## Fixed bandwidth: 5.01041 CV score: 543.7094
## Fixed bandwidth: 5.026095 CV score: 543.3305
## Fixed bandwidth: 5.035789 CV score: 543.1
## Fixed bandwidth: 5.04178 CV score: 542.959
## Fixed bandwidth: 5.045482 CV score: 542.8723
## Fixed bandwidth: 5.047771 CV score: 542.819
## Fixed bandwidth: 5.049185 CV score: 542.7861
## Fixed bandwidth: 5.050059 CV score: 542.7658
## Fixed bandwidth: 5.050599 CV score: 542.7532
## Fixed bandwidth: 5.050933 CV score: 542.7455
## Fixed bandwidth: 5.051139 CV score: 542.7407
## Fixed bandwidth: 5.051267 CV score: 542.7378
## Fixed bandwidth: 5.051346 CV score: 542.7359
## Fixed bandwidth: 5.051395 CV score: 542.7348
## Fixed bandwidth: 5.051425 CV score: 542.7341gwr_div <- gwr.basic(
formula = form_div,
data = shp_sp,
bw = bw_div,
dMat = dMat,
kernel = "bisquare",
adaptive = FALSE
)
# Extraer coeficiente local de ln_supag
gwr_coef <- data.frame(
dstocode = shp_ana$dstocode,
beta_ln_supag = gwr_div$SDF@data$ln_supag
)
shp_all <- shp_all %>% left_join(gwr_coef, by = "dstocode")# 10) Mapas ---------------------------------------------------------------
tmap_mode("plot")
# 10.1 Coroplético diversificación ganadera
tm_shape(shp_all) +
tm_polygons(
"div_med",
style = "quantile",
n = 5,
palette = brewer.pal(5, "Oranges"),
title = "Índice diversificación",
colorNA = "lightgrey", textNA = "Sin dato",
border.col= "grey50", lwd = 0.3
) +
tm_layout(legend.outside = TRUE,
main.title = "Diversificación ganadera media por distrito")
# 10.2 Coeficiente local ln_supag (GWR)
tm_shape(shp_all) +
tm_polygons(
"beta_ln_supag",
style = "quantile",
n = 5,
palette = brewer.pal(5, "RdBu"),
title = expression(beta[ln~supag]),
colorNA = "lightgrey", textNA = "No estimado",
border.col= "grey50", lwd = 0.3
) +
tm_layout(legend.outside = TRUE,
main.title = "Coef. local de escala (GWR sobre diversificación)")