Analítica geoespacial del Censo Agropecuario 2022
Richard González
17 de mayo de 2025
1 Resumen
2 Objetivos:
2.1 General
** replicar en R (RStudio) el flujo de limpieza, agregación, matriz de pesos y modelo espacial Durbin (SDM) usando los micro‑datos y el shapefile reproyectado a UTM 21S. Todo el código es auto‑contenible y puede compilarse en Knit ➜ HTML.
2.2 Específicos
# Paquetes indispensables -------------------------------------------------
libs <- c("tidyverse", "sf", "spdep", "spatialreg", "tmap", "knitr")
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 defecto3 Materiales y métodos
3.1 Qué es el CAN 2022
enlace(www.xxxxxxxx)
3.2 Métodos estadísticos
3.2.1 Modelo espacial de Durbin (SMD)
3.2.2 Modelo alternativo
4 Gestión de dastos
4.1 Carga de datos
Ajusta
data_dira tu propia ruta de trabajo.
# Directorios --------------------------------------------------------------
data_dir <- "G:/Mi unidad/CAN2022/"
geo_dir <- file.path(data_dir, "paraguay")
# Shapefile distritos en metros (EPSG:32721) ------------------------------
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 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/CAN2022/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)
)4.2 Unipon de la base del CAN con la base geoespacial para mapear
##
## 101 102 103 104 105 106 107 108 109 110 111 112 113 201 202 203
## 1260 1158 2191 1662 60 141 1332 620 526 321 292 1648 872 2275 117 3334
## 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219
## 2206 1217 1547 810 4699 638 1969 980 1501 1114 1302 1051 3566 2837 2663 1867
## 220 221 301 302 303 304 305 306 307 308 309 310 311 312 313 314
## 2171 1616 1223 813 1812 770 1228 419 1115 560 703 1053 289 517 425 1139
## 315 316 317 318 319 320 401 402 403 404 405 406 407 408 409 410
## 831 371 702 482 732 355 2260 841 1019 525 636 1144 1673 548 614 350
## 411 412 413 414 415 416 417 418 501 502 503 504 505 506 507 508
## 573 289 170 450 180 123 1612 37 3948 1527 1795 606 985 1990 1336 581
## 509 510 511 512 513 514 515 516 517 518 519 520 521 522 601 602
## 1708 1768 2849 1924 954 1096 778 395 611 1082 933 1855 1590 514 2139 1223
## 603 604 605 606 607 608 609 610 611 701 702 703 704 705 706 707
## 537 516 889 774 1641 1568 1071 2051 1656 436 802 765 960 666 486 392
## 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723
## 746 1439 1355 503 1055 799 339 188 456 844 822 600 1947 1488 248 1562
## 724 725 726 727 728 729 730 801 802 803 804 805 806 807 808 809
## 1880 1164 161 837 93 481 956 606 390 952 311 151 795 1864 448 36
## 810 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915
## 143 677 1115 633 635 1659 623 501 754 670 1067 827 586 411 393 1429
## 916 917 918 1002 1003 1004 1005 1006 1007 1008 1009 1010 1011 1012 1013 1014
## 1091 640 262 157 667 704 635 1680 540 593 218 680 683 545 360 594
## 1015 1016 1017 1018 1019 1020 1021 1022 1101 1104 1105 1106 1111 1116 1117 1119
## 189 926 358 412 316 430 603 435 321 149 683 332 411 565 58 353
## 1201 1202 1203 1204 1205 1206 1207 1208 1209 1210 1211 1212 1213 1214 1215 1216
## 265 151 580 256 541 295 388 546 432 466 111 547 370 83 224 363
## 1301 1302 1303 1304 1305 1306 1401 1402 1403 1404 1405 1406 1407 1408 1409 1410
## 216 368 733 420 272 595 112 304 1349 1337 278 904 385 233 108 249
## 1411 1412 1413 1414 1415 1416 1502 1503 1504 1506 1507 1508 1509 1510 1602 1604
## 1715 781 1468 2186 86 75 288 423 509 124 1038 176 338 377 914 601
## 1605 1701 1702 1704 1705
## 514 205 426 162 175# 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>, …5 Resultados
5.1 Adopción de agricultura de precisión
5.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)5.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>5.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 -...5.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-graphs5.1.5 Moran I global
## Moran I statistic Expectation Variance
## 0.040163665 -0.023255814 0.0016049575.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.
5.1.7 Mapa coroplético
map <- tm_shape(shp_d) +
tm_polygons("agriprec", style = "quantile", n = 5, palette = "Blues",
title = "Prop. adopción AP") +
tm_layout(legend.outside = TRUE)## ## ── tmap v3 code detected ───────────────────────────────────────────────────────## [v3->v4] `tm_polygons()`: instead of `style = "quantile"`, use fill.scale =
## `tm_scale_intervals()`.
## ℹ Migrate the argument(s) 'style', 'n', 'palette' (rename to 'values') to
## 'tm_scale_intervals(<HERE>)'
## [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>)'## [cols4all] color palettes: use palettes from the R package cols4all. Run
## `cols4all::c4a_gui()` to explore them. The old palette name "Blues" is named
## "brewer.blues"
## Multiple palettes called "blues" found: "brewer.blues", "matplotlib.blues". The first one, "brewer.blues", is returned.