Proyecto análisis estadístico sobre la compensación por aprovechamiento de fauna silvestre en Colombia
María Paula Gómez Silva - Angie Maria Moreno Mantilla
Facultad de Ciencias Naturales e Ingeniería
San Gil, Santander - 2026
1 Problema de investigación
A nivel nacional, cuantificar el aprovechamiento de los recursos naturales renovables es un problema técnico recurrente que involucra procesos de liquidación económica, los cuales dependen de múltiples variables técnicas y ecológicas medibles. En Colombia, uno de estos procesos es el cálculo de la tasa compensatoria por fauna silvestre nativa, reglamentada en el artículo 42 de la Ley 99 de 1993 y administrada por el Ministerio de Ambiente y Desarrollo Sostenible - Minambiente.
El problema que motiva este estudio surge en un contexto de ingeniería ambiental en donde, aunque la fórmula de liquidación se encuentra definida normativamente, en la práctica, el sistema de información puede generar registros con alta variabilidad en los coeficientes, diferencias regionales y montos que oscilan entre valores muy bajos y sumamente elevados, los cuales impiden identificar con claridad qué variables tienen mayor influencia en el monto final.
Pregunta principal: ¿Cuál es el comportamiento estadístico de las variables cuantitativas que determinan el monto a pagar en la tasa compensatoria por caza de fauna silvestre en Colombia, y en qué medida el Factor Regional (FR) y el número de especímenes aprobados (Es) influyen linealmente sobre dicho monto?
1.1 Identificación de variables
Variable dependiente (Y):
- Monto a pagar (MP): variable cuantitativa continua expresada en pesos colombianos (COP). Representa el valor económico que debe cancelar cada usuario autorizado para ejercer la caza de fauna silvestre nativa.
Variables independientes (X):
- X1 – Factor Regional (FR): pondera la presión de uso y el estado de conservación del hábitat en la región donde se ejerce la caza. A mayor intervención del ecosistema, mayor es su valor.
- X2 – Coeficiente Biótico (Cb): refleja el estado de conservación de la especie autorizada y de su hábitat. Especies más vulnerables reciben valores más altos.
- X3 – Grupo Trófico (Gt): coeficiente numérico que pondera el impacto ecológico según la posición de la especie en la red trófica (0 < Gt ≤ 1).
- X4 – Coeficiente de Valoración (V): diferencia el impacto esperado según la finalidad de la caza: científica, comercial, de fomento o deportiva.
- X5 – Número de especímenes autorizados (Es): cantidad de individuos o muestras aprobados en cada permiso. A mayor número, mayor extracción sobre el recurso.
2 Cargar librerías
Para este proyecto, se cargarán herramientas especializadas en el ecosistema de R que cubren las fases críticas del análisis de datos: desde la importación y manipulación eficiente de registros (con el núcleo de dplyr y readr), pasando por la exploración estadística avanzada y visual (usando psych y skimr), hasta la creación de gráficas de alta calidad y diagnósticos de modelado estadístico.
library(psych) # Estadísticos descriptivos avanzados
library(corrplot) # Visualización de matrices de correlación
library(readr) # Importación eficiente de datos
library(GGally) # Gráficos multivariados (ggpairs)
library(car) # Diagnóstico de modelos de regresión
library(MASS) # Selección de modelos (stepAIC)
library(ggplot2) # Visualización gráfica de datos
library(skimr) # Resumen más visual
library(dplyr) # Manipulación y transformación de datos3 Carga de datos
Se importa el archivo CSV con los registros de la Tasa Compensatoria por Caza de Fauna Silvestre entre los años 2016 y 2022. Se especifican explícitamente los tipos de cada columna para garantizar una correcta interpretación de las variables numéricas y de texto, evitando errores de conversión automática por parte de R.
Datos_Tasa_compensatoria_por_Caza_de_Fauna <- read_csv(
"Datos_de_la_Tasa_compensatoria_por_Caza_de_Fauna_Silvestre_20260509.csv",
col_types = cols(
Año = col_character(),
`FECHA DE OTORGAMIENTO` = col_character(),
`VIGENCIA (meses)` = col_number(),
`VALOR TIPO DE CAZA` = col_number(),
`COEFICIENTE BIÓTICO (Cb)` = col_number(),
`NACIONALIDAD (N)` = col_number(),
`GRUPO TRÓFICO (Gt)` = col_number(),
`COEFICIENTE DE VALORACIÓN (V)` = col_number(),
`FACTOR REGIONAL (FR)` = col_number(),
`No. ESPECÍMENES O MUESTRAS (Es)` = col_number(),
`TARIFA MÍNIMA (TM)` = col_number(),
`COSTO DE IMPLEMENTACIÓN (CI)` = col_number(),
`MONTO A PAGAR (MP)` = col_number(),
`MP DESPUÉS DE RECLAMACIÓN` = col_number(),
RECAUDO = col_number()
)
)4 Exploración inicial
Antes de realizar cualquier análisis, se explora la base de datos
para comprender su estructura. Las funciones head(),
summary() y dim() permiten obtener una visión
general del contenido, el tamaño y la distribución de las variables.
Adicionalmente, colSums(is.na()) y sapply()
identifican la cantidad y porcentaje de valores faltantes por columna, y
names() lista los nombres de todas las variables
disponibles.
4.1 Preparación de datos
# Asignamos nombre corto para facilitar el trabajo
datostasaCF <- Datos_Tasa_compensatoria_por_Caza_de_Fauna
head(datostasaCF)Tipo de datos:
## spc_tbl_ [1,526 × 34] (S3: spec_tbl_df/tbl_df/tbl/data.frame)
## $ Año : chr [1:1526] "2,021" "2,021" "2,021" "2,021" ...
## $ Autoridad Ambiental : chr [1:1526] "CORTOLIMA" "CORTOLIMA" "CORTOLIMA" "CORTOLIMA" ...
## $ NOMBRE/RAZÓN SOCIAL : chr [1:1526] "UNVERSIDAD DEL TOLIMA" "UNVERSIDAD DEL TOLIMA" "UNVERSIDAD DEL TOLIMA" "UNVERSIDAD DEL TOLIMA" ...
## $ PERSONA NATURAL O JURÍDICA : chr [1:1526] "Persona jurídica" "Persona jurídica" "Persona jurídica" "Persona jurídica" ...
## $ REPRESENTANTE LEGAL : chr [1:1526] "Omar A. Mejía Patiño" "Omar A. Mejía Patiño" "Omar A. Mejía Patiño" "Omar A. Mejía Patiño" ...
## $ TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA : chr [1:1526] "Permiso / Autoridad Ambiental diferente a ANLA" "Permiso / Autoridad Ambiental diferente a ANLA" "Permiso / Autoridad Ambiental diferente a ANLA" "Permiso / Autoridad Ambiental diferente a ANLA" ...
## $ No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO: chr [1:1526] "Resolución 3758 de 2016" "Resolución 3758 de 2016" "Resolución 3758 de 2016" "Resolución 3758 de 2016" ...
## $ FECHA DE OTORGAMIENTO : chr [1:1526] "2016-11-16T00:00:00.000" "2016-11-16T00:00:00.000" "2016-11-16T00:00:00.000" "2016-11-16T00:00:00.000" ...
## $ VIGENCIA (meses) : num [1:1526] 120 120 120 120 120 120 120 6 120 120 ...
## $ TIPO DE CAZA : chr [1:1526] "Científica estudios ambientales" "Científica estudios ambientales" "Científica estudios ambientales" "Científica estudios ambientales" ...
## $ VALOR TIPO DE CAZA : num [1:1526] 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 0.6 ...
## $ ESPECIE O GRUPO TAXONÓMICO : chr [1:1526] "Panthiades bitias" "Actinote anteas" "Altinote dicaeus" "Catasticta tomyris tomyris" ...
## $ ORDEN : chr [1:1526] "Lepidoptera" "Lepidoptera" "Lepidoptera" "Lepidoptera" ...
## $ CLASE : chr [1:1526] "Insecta" "Insecta" "Insecta" "Insecta" ...
## $ No. INDIVIDUOS y/o MUESTRAS APROBADOS : chr [1:1526] "No registra" "No registra" "No registra" "No registra" ...
## $ DEPARTAMENTO : chr [1:1526] "Tolima" "Tolima" "Tolima" "Tolima" ...
## $ MUNICIPIO : chr [1:1526] "Dolores, Natagaima, Prado, Purificación, Cunday,Cunday, Venadillo, Alvarado, Piedras, Ortega, Espinal, Guamo, A"| __truncated__ "Dolores, Natagaima, Prado, Purificación, Cunday,Cunday, Venadillo, Alvarado, Piedras, Ortega, Espinal, Guamo, A"| __truncated__ "Dolores, Natagaima, Prado, Purificación, Cunday,Cunday, Venadillo, Alvarado, Piedras, Ortega, Espinal, Guamo, A"| __truncated__ "Dolores, Natagaima, Prado, Purificación, Cunday,Cunday, Venadillo, Alvarado, Piedras, Ortega, Espinal, Guamo, A"| __truncated__ ...
## $ No. CERTIFICADO SiB : chr [1:1526] "1634DA1EAAC -16341C01996 -6321985BD7- 164F610CBDC- 16DBB9126D7- 16321A6F91C- 16321C92741- 16364C22A54- 1632181A"| __truncated__ "1634DA1EAAC -16341C01996 -6321985BD7- 164F610CBDC- 16DBB9126D7- 16321A6F91C- 16321C92741- 16364C22A54- 1632181A"| __truncated__ "1634DA1EAAC -16341C01996 -6321985BD7- 164F610CBDC- 16DBB9126D7- 16321A6F91C- 16321C92741- 16364C22A54- 1632181A"| __truncated__ "1634DA1EAAC -16341C01996 -6321985BD7- 164F610CBDC- 16DBB9126D7- 16321A6F91C- 16321C92741- 16364C22A54- 1632181A"| __truncated__ ...
## $ ESTADO DE CONSERVACIÓN DE LA ESPECIE : chr [1:1526] "Datos insuficientes (DD) o especies No evaluadas (NE)" "Datos insuficientes (DD) o especies No evaluadas (NE)" "Datos insuficientes (DD) o especies No evaluadas (NE)" "Datos insuficientes (DD) o especies No evaluadas (NE)" ...
## $ ESTADO DE CONSERVACIÓN DEL HÁBITAT : chr [1:1526] "Pobremente conservado" "Pobremente conservado" "Pobremente conservado" "Pobremente conservado" ...
## $ PRESIÓN POR USO : chr [1:1526] NA NA NA NA ...
## $ COEFICIENTE BIÓTICO (Cb) : num [1:1526] 1 1 1 1 1 1 1 1 1 1 ...
## $ NACIONALIDAD (N) : num [1:1526] 0 0 0 0 0 0 0 0 0 0 ...
## $ GRUPO TRÓFICO (Gt) : num [1:1526] 0.15 0.15 0.15 0.15 0.15 0.15 0.15 0.8 0.15 0.15 ...
## $ COEFICIENTE DE VALORACIÓN
## (V) : num [1:1526] 1 1 1 1 1 1 1 1 1 1 ...
## $ FACTOR REGIONAL (FR) : num [1:1526] 0.09 0.09 0.09 0.09 0.09 0.09 0.09 0.48 0.09 0.09 ...
## $ No. ESPECÍMENES O MUESTRAS (Es) : num [1:1526] 1 1 1 1 1 1 1 1 1 1 ...
## $ UNIDAD DE Es : chr [1:1526] "Individuo" "Individuo" "Individuo" "Individuo" ...
## $ TARIFA MÍNIMA (TM) : num [1:1526] 10903 10903 10903 10903 10903 ...
## $ COSTO DE IMPLEMENTACIÓN (CI) : num [1:1526] 29530 29530 29530 29530 29530 ...
## $ MONTO A PAGAR (MP) : num [1:1526] 981 981 981 981 981 ...
## $ RECLAMACIONES : chr [1:1526] "No" "No" "No" "No" ...
## $ MP DESPUÉS DE RECLAMACIÓN : num [1:1526] NA NA NA NA NA NA NA NA NA NA ...
## $ RECAUDO : num [1:1526] 981 981 981 981 981 ...
## - attr(*, "spec")=
## .. cols(
## .. Año = col_character(),
## .. `Autoridad Ambiental` = col_character(),
## .. `NOMBRE/RAZÓN SOCIAL` = col_character(),
## .. `PERSONA NATURAL O JURÍDICA` = col_character(),
## .. `REPRESENTANTE LEGAL` = col_character(),
## .. `TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA` = col_character(),
## .. `No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO` = col_character(),
## .. `FECHA DE OTORGAMIENTO` = col_character(),
## .. `VIGENCIA (meses)` = col_number(),
## .. `TIPO DE CAZA` = col_character(),
## .. `VALOR TIPO DE CAZA` = col_number(),
## .. `ESPECIE O GRUPO TAXONÓMICO` = col_character(),
## .. ORDEN = col_character(),
## .. CLASE = col_character(),
## .. `No. INDIVIDUOS y/o MUESTRAS APROBADOS` = col_character(),
## .. DEPARTAMENTO = col_character(),
## .. MUNICIPIO = col_character(),
## .. `No. CERTIFICADO SiB` = col_character(),
## .. `ESTADO DE CONSERVACIÓN DE LA ESPECIE` = col_character(),
## .. `ESTADO DE CONSERVACIÓN DEL HÁBITAT` = col_character(),
## .. `PRESIÓN POR USO` = col_character(),
## .. `COEFICIENTE BIÓTICO (Cb)` = col_number(),
## .. `NACIONALIDAD (N)` = col_number(),
## .. `GRUPO TRÓFICO (Gt)` = col_number(),
## .. `COEFICIENTE DE VALORACIÓN
## .. (V)` = col_double(),
## .. `FACTOR REGIONAL (FR)` = col_number(),
## .. `No. ESPECÍMENES O MUESTRAS (Es)` = col_number(),
## .. `UNIDAD DE Es` = col_character(),
## .. `TARIFA MÍNIMA (TM)` = col_number(),
## .. `COSTO DE IMPLEMENTACIÓN (CI)` = col_number(),
## .. `MONTO A PAGAR (MP)` = col_number(),
## .. RECLAMACIONES = col_character(),
## .. `MP DESPUÉS DE RECLAMACIÓN` = col_number(),
## .. RECAUDO = col_number()
## .. )
## - attr(*, "problems")=<externalptr>Tamaño (número de filas y número de columnas):
## [1] 1526 34Nombres de las variables:
## [1] "Año"
## [2] "Autoridad Ambiental"
## [3] "NOMBRE/RAZÓN SOCIAL"
## [4] "PERSONA NATURAL O JURÍDICA"
## [5] "REPRESENTANTE LEGAL"
## [6] "TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA"
## [7] "No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO"
## [8] "FECHA DE OTORGAMIENTO"
## [9] "VIGENCIA (meses)"
## [10] "TIPO DE CAZA"
## [11] "VALOR TIPO DE CAZA"
## [12] "ESPECIE O GRUPO TAXONÓMICO"
## [13] "ORDEN"
## [14] "CLASE"
## [15] "No. INDIVIDUOS y/o MUESTRAS APROBADOS"
## [16] "DEPARTAMENTO"
## [17] "MUNICIPIO"
## [18] "No. CERTIFICADO SiB"
## [19] "ESTADO DE CONSERVACIÓN DE LA ESPECIE"
## [20] "ESTADO DE CONSERVACIÓN DEL HÁBITAT"
## [21] "PRESIÓN POR USO"
## [22] "COEFICIENTE BIÓTICO (Cb)"
## [23] "NACIONALIDAD (N)"
## [24] "GRUPO TRÓFICO (Gt)"
## [25] "COEFICIENTE DE VALORACIÓN \n(V)"
## [26] "FACTOR REGIONAL (FR)"
## [27] "No. ESPECÍMENES O MUESTRAS (Es)"
## [28] "UNIDAD DE Es"
## [29] "TARIFA MÍNIMA (TM)"
## [30] "COSTO DE IMPLEMENTACIÓN (CI)"
## [31] "MONTO A PAGAR (MP)"
## [32] "RECLAMACIONES"
## [33] "MP DESPUÉS DE RECLAMACIÓN"
## [34] "RECAUDO"Resumen de las variables:
## Año Autoridad Ambiental NOMBRE/RAZÓN SOCIAL
## Length:1526 Length:1526 Length:1526
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## PERSONA NATURAL O JURÍDICA REPRESENTANTE LEGAL
## Length:1526 Length:1526
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA
## Length:1526
## Class :character
## Mode :character
##
##
##
##
## No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO FECHA DE OTORGAMIENTO
## Length:1526 Length:1526
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## VIGENCIA (meses) TIPO DE CAZA VALOR TIPO DE CAZA
## Min. : 0.00 Length:1526 Min. :0.1000
## 1st Qu.: 6.00 Class :character 1st Qu.:0.6000
## Median : 24.00 Mode :character Median :0.6000
## Mean : 36.48 Mean :0.5965
## 3rd Qu.: 24.00 3rd Qu.:0.6000
## Max. :120.00 Max. :1.2000
##
## ESPECIE O GRUPO TAXONÓMICO ORDEN CLASE
## Length:1526 Length:1526 Length:1526
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## No. INDIVIDUOS y/o MUESTRAS APROBADOS DEPARTAMENTO MUNICIPIO
## Length:1526 Length:1526 Length:1526
## Class :character Class :character Class :character
## Mode :character Mode :character Mode :character
##
##
##
##
## No. CERTIFICADO SiB ESTADO DE CONSERVACIÓN DE LA ESPECIE
## Length:1526 Length:1526
## Class :character Class :character
## Mode :character Mode :character
##
##
##
##
## ESTADO DE CONSERVACIÓN DEL HÁBITAT PRESIÓN POR USO COEFICIENTE BIÓTICO (Cb)
## Length:1526 Length:1526 Min. :0.000
## Class :character Class :character 1st Qu.:1.000
## Mode :character Mode :character Median :1.000
## Mean :1.209
## 3rd Qu.:1.000
## Max. :5.000
##
## NACIONALIDAD (N) GRUPO TRÓFICO (Gt) COEFICIENTE DE VALORACIÓN \n(V)
## Min. :0 Min. :0.0000 Min. : 0.0000
## 1st Qu.:0 1st Qu.:0.8000 1st Qu.: 1.0000
## Median :0 Median :0.8000 Median : 1.0000
## Mean :0 Mean :0.7282 Mean : 0.9784
## 3rd Qu.:0 3rd Qu.:0.8000 3rd Qu.: 1.0000
## Max. :0 Max. :1.0000 Max. :20.0000
##
## FACTOR REGIONAL (FR) No. ESPECÍMENES O MUESTRAS (Es) UNIDAD DE Es
## Min. : 0.0000 Min. : 0.00 Length:1526
## 1st Qu.: 0.4800 1st Qu.: 1.00 Class :character
## Median : 0.4800 Median : 2.00 Mode :character
## Mean : 0.6805 Mean : 14.92
## 3rd Qu.: 0.4800 3rd Qu.: 5.00
## Max. :76.8000 Max. :7261.00
##
## TARIFA MÍNIMA (TM) COSTO DE IMPLEMENTACIÓN (CI) MONTO A PAGAR (MP)
## Min. : 9600 Min. :26000 Min. : 966
## 1st Qu.:10567 1st Qu.:28620 1st Qu.: 5433
## Median :10903 Median :29530 Median : 11278
## Mean :10961 Mean :29686 Mean : 59518
## 3rd Qu.:11318 3rd Qu.:30652 3rd Qu.: 40577
## Max. :11748 Max. :31817 Max. :7334339
## NA's :50 NA's :50 NA's :385
## RECLAMACIONES MP DESPUÉS DE RECLAMACIÓN RECAUDO
## Length:1526 Min. : 0 Min. : 981
## Class :character 1st Qu.: 6791 1st Qu.: 5433
## Mode :character Median : 16298 Median : 10865
## Mean : 41668 Mean : 37914
## 3rd Qu.: 38028 3rd Qu.: 31817
## Max. :1249507 Max. :1962540
## NA's :1012 NA's :565Interpretación de la exploración inicial:
head(): Los primeros 6 registros
confirman que la carga fue exitosa. Se observan variables de tipo
carácter como Autoridad Ambiental, Tipo de Caza y Especie, y variables
numéricas como los coeficientes de valoración y el Monto a Pagar, lo que
es consistente con la especificación de tipos definida en la carga.
summary(): La variable MP (Monto a
Pagar) presenta una media muy superior a su mediana, lo que evidencia
una distribución fuertemente asimétrica hacia la derecha, con presencia
de valores extremos que elevan el promedio. La variable RECAUDO muestra
un comportamiento similar. La variable VIGENCIA (meses) oscila entre 0 y
120 meses, con una media de 36.48 meses, lo que indica una amplia
variedad en la duración de los permisos otorgados. La variable MP
DESPUÉS DE RECLAMACIÓN presenta 1.012 valores faltantes de 1.526,
indicando que aproximadamente el 66% de los permisos no tuvieron
reclamaciones.
dim(): El conjunto de datos tiene 1.526
filas y 34 columnas. Este tamaño muestral es adecuado para la aplicación
de técnicas estadísticas como el análisis de correlación y la regresión
lineal múltiple, permitiendo estimaciones estables de los coeficientes
del modelo.
names(): Las 34 variables cubren
información administrativa (autoridad, permiso, razón social), biológica
(especie, orden, clase, estado de conservación) y económica
(coeficientes, monto, recaudo), lo que refleja la naturaleza
multidimensional del sistema de liquidación de la tasa
compensatoria.
4.2 Datos faltantes
Conteo de NA por columna:
## Año
## 0
## Autoridad Ambiental
## 0
## NOMBRE/RAZÓN SOCIAL
## 0
## PERSONA NATURAL O JURÍDICA
## 0
## REPRESENTANTE LEGAL
## 42
## TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA
## 0
## No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO
## 0
## FECHA DE OTORGAMIENTO
## 0
## VIGENCIA (meses)
## 0
## TIPO DE CAZA
## 0
## VALOR TIPO DE CAZA
## 0
## ESPECIE O GRUPO TAXONÓMICO
## 0
## ORDEN
## 4
## CLASE
## 4
## No. INDIVIDUOS y/o MUESTRAS APROBADOS
## 405
## DEPARTAMENTO
## 0
## MUNICIPIO
## 0
## No. CERTIFICADO SiB
## 0
## ESTADO DE CONSERVACIÓN DE LA ESPECIE
## 383
## ESTADO DE CONSERVACIÓN DEL HÁBITAT
## 352
## PRESIÓN POR USO
## 866
## COEFICIENTE BIÓTICO (Cb)
## 0
## NACIONALIDAD (N)
## 0
## GRUPO TRÓFICO (Gt)
## 0
## COEFICIENTE DE VALORACIÓN \n(V)
## 0
## FACTOR REGIONAL (FR)
## 0
## No. ESPECÍMENES O MUESTRAS (Es)
## 0
## UNIDAD DE Es
## 0
## TARIFA MÍNIMA (TM)
## 50
## COSTO DE IMPLEMENTACIÓN (CI)
## 50
## MONTO A PAGAR (MP)
## 385
## RECLAMACIONES
## 0
## MP DESPUÉS DE RECLAMACIÓN
## 1012
## RECAUDO
## 565Porcentaje de NA por columna:
## Año
## 0.0000000
## Autoridad Ambiental
## 0.0000000
## NOMBRE/RAZÓN SOCIAL
## 0.0000000
## PERSONA NATURAL O JURÍDICA
## 0.0000000
## REPRESENTANTE LEGAL
## 2.7522936
## TIPO DE AUTORIZACIÓN / ENTIDAD QUE OTORGA
## 0.0000000
## No. PERMISO/ LICENCIA AMBIENTAL/ SANCIONATORIO
## 0.0000000
## FECHA DE OTORGAMIENTO
## 0.0000000
## VIGENCIA (meses)
## 0.0000000
## TIPO DE CAZA
## 0.0000000
## VALOR TIPO DE CAZA
## 0.0000000
## ESPECIE O GRUPO TAXONÓMICO
## 0.0000000
## ORDEN
## 0.2621232
## CLASE
## 0.2621232
## No. INDIVIDUOS y/o MUESTRAS APROBADOS
## 26.5399738
## DEPARTAMENTO
## 0.0000000
## MUNICIPIO
## 0.0000000
## No. CERTIFICADO SiB
## 0.0000000
## ESTADO DE CONSERVACIÓN DE LA ESPECIE
## 25.0982962
## ESTADO DE CONSERVACIÓN DEL HÁBITAT
## 23.0668414
## PRESIÓN POR USO
## 56.7496723
## COEFICIENTE BIÓTICO (Cb)
## 0.0000000
## NACIONALIDAD (N)
## 0.0000000
## GRUPO TRÓFICO (Gt)
## 0.0000000
## COEFICIENTE DE VALORACIÓN \n(V)
## 0.0000000
## FACTOR REGIONAL (FR)
## 0.0000000
## No. ESPECÍMENES O MUESTRAS (Es)
## 0.0000000
## UNIDAD DE Es
## 0.0000000
## TARIFA MÍNIMA (TM)
## 3.2765400
## COSTO DE IMPLEMENTACIÓN (CI)
## 3.2765400
## MONTO A PAGAR (MP)
## 25.2293578
## RECLAMACIONES
## 0.0000000
## MP DESPUÉS DE RECLAMACIÓN
## 66.3171691
## RECAUDO
## 37.0249017Interpretación de los datos faltantes:
Se detectaron valores faltantes principalmente en TARIFA MÍNIMA (TM), COSTO DE IMPLEMENTACIÓN (CI), MONTO A PAGAR (MP) (~25%) y MP DESPUÉS DE RECLAMACIÓN (~66%).
- Si < 5%: impacto bajo.
- Si > 5%: considerar imputación o eliminación.
Los faltantes en MP son los más críticos para el análisis, ya que corresponden a la variable dependiente del estudio. Los faltantes en las demás variables tienen un impacto bajo dado que están por debajo del 5%.
5 Selección de variables
De las 34 variables disponibles en el conjunto de datos, se seleccionan únicamente las 6 que participan directamente en la fórmula de cálculo de la tasa compensatoria. Se renombran con nombres cortos para facilitar el análisis posterior. La variable MP actúa como variable dependiente (Y) y las demás como variables independientes (X).
tabla_analisis <- dplyr::select(datostasaCF,
"MONTO A PAGAR (MP)", # Variable depediente
"FACTOR REGIONAL (FR)",
"COEFICIENTE BIÓTICO (Cb)",
"GRUPO TRÓFICO (Gt)",
"COEFICIENTE DE VALORACIÓN \n(V)",
"No. ESPECÍMENES O MUESTRAS (Es)"
)
#Renombrar columnas
colnames(tabla_analisis) <- c("MP", "FR", "Cb", "Gt", "V", "Es")
#Vista previa
head(tabla_analisis)Resumen de variables seleccionadas:
## MP FR Cb Gt
## Min. : 966 Min. : 0.0000 Min. :0.000 Min. :0.0000
## 1st Qu.: 5433 1st Qu.: 0.4800 1st Qu.:1.000 1st Qu.:0.8000
## Median : 11278 Median : 0.4800 Median :1.000 Median :0.8000
## Mean : 59518 Mean : 0.6805 Mean :1.209 Mean :0.7282
## 3rd Qu.: 40577 3rd Qu.: 0.4800 3rd Qu.:1.000 3rd Qu.:0.8000
## Max. :7334339 Max. :76.8000 Max. :5.000 Max. :1.0000
## NA's :385
## V Es
## Min. : 0.0000 Min. : 0.00
## 1st Qu.: 1.0000 1st Qu.: 1.00
## Median : 1.0000 Median : 2.00
## Mean : 0.9784 Mean : 14.92
## 3rd Qu.: 1.0000 3rd Qu.: 5.00
## Max. :20.0000 Max. :7261.00
## ## [1] 1526 6Datos faltantes en variables seleccionadas:
## MP FR Cb Gt V Es
## 385 0 0 0 0 0## MP FR Cb Gt V Es
## 25.22936 0.00000 0.00000 0.00000 0.00000 0.00000Interpretación:
Se seleccionaron 6 variables de las 34 disponibles, todas directamente vinculadas a la fórmula normativa de cálculo de la tasa compensatoria. La variable MP actúa como variable dependiente (Y) y representa el resultado económico del proceso de liquidación. Las variables FR, Cb, Gt, V y Es actúan como variables independientes (X) y corresponden a los factores técnicos y ecológicos que determinan el monto. Tras renombrar las columnas, el subconjunto conserva las 1.526 filas originales, con valores faltantes únicamente en MP (~25%), los cuales serán eliminados antes del análisis para garantizar la integridad de los resultados.
5.1 Visualización de filas vacías
Se identifican y visualizan las filas donde la variable dependiente MP no tiene valor registrado. Estas observaciones no pueden usarse en el análisis de correlación ni en el modelo de regresión, por lo que serán eliminadas en el paso siguiente.
5.2 Limpieza de datos
Se eliminan todas las filas con valores faltantes mediante
na.omit(), obteniendo un conjunto de datos completo listo
para el análisis. La función describe() del paquete
psych entrega estadísticos descriptivos avanzados
incluyendo media, desviación estándar, asimetría y curtosis.
#Eliminar filas con valores faltantes
tabla_analisis <- na.omit(tabla_analisis)
# Configurar para ver números normales y no científicos
options(scipen = 999)
resumen_estadistico <- describe(tabla_analisis)
round(resumen_estadistico, 2)Interpretación:
Tras aplicar na.omit(), el conjunto de análisis quedó con 1.141 observaciones completas, eliminando 385 filas que contenían al menos un valor faltante en alguna de las variables seleccionadas. La función describe() revela que la variable MP tiene una media de $59.518 COP con una desviación estándar muy elevada, confirmando la alta variabilidad y la presencia de valores extremos. La variable Es presenta la mayor dispersión relativa de todas las variables, con un rango que va desde 0 hasta 7.261 especímenes, lo que la convierte en el predictor con mayor variabilidad y, posiblemente, con mayor capacidad explicativa sobre el monto final..
6 Estadísticos descriptivos
En esta sección se presentan dos tablas complementarias. La primera,
generada con la función describe() del paquete
psych, ofrece un resumen estadístico completo que incluye
número de observaciones, media, desviación estándar, mediana, mínimo,
máximo, rango, error estándar, asimetría y curtosis. La segunda tabla
consolida de forma más organizada las medidas de tendencia central,
dispersión y variabilidad relativa, añadiendo la moda y el coeficiente
de variación que describe() no calcula directamente.
6.1 Tabla 1 — Resumen estadístico completo
6.2 Tabla 2 — Medidas de tendencia central, dispersión y variabilidad
# 1. Calcular la varianza
varianzas <- apply(tabla_analisis, 2, var, na.rm = TRUE)
# 2. Calcular el coeficiente de variación (CV%)
cv_porcentaje <- (resumen_estadistico$sd / resumen_estadistico$mean) * 100
# 3. Función para calcular la moda
calcular_moda <- function(x) {
ux <- unique(x)
ux[which.max(tabulate(match(x, ux)))]
}
# Calcular las modas de cada variable
modas <- apply(tabla_analisis, 2, calcular_moda)
# 4. Consolidar en una sola tabla
tabla_final <- data.frame(
Media = resumen_estadistico$mean,
Mediana = resumen_estadistico$median,
Moda = modas,
Minimo = resumen_estadistico$min,
Maximo = resumen_estadistico$max,
Varianza = varianzas,
Sd = resumen_estadistico$sd,
CV_porcentaje = cv_porcentaje
)
# Sin notación científica y redondeado a 2 decimales
options(scipen = 999)
tabla_final <- round(tabla_final, 2)
tabla_finalInterpretación de los estadísticos descriptivos:
La variable MP (Monto a pagar) presenta una media de $59.518 COP con una desviación estándar muy elevada, indicando una variabilidad extrema entre permisos. Su asimetría positiva pronunciada sugiere que la mayoría de los permisos generan montos bajos, mientras que un grupo reducido de permisos con muchos especímenes genera montos muy altos que elevan considerablemente el promedio. El coeficiente de variación (CV) confirma esta alta dispersión relativa respecto a su media.
La variable Es (Número de especímenes) muestra una media de 14.92 con una desviación estándar altísima, reflejo de permisos que van desde 1 espécimen hasta más de 7.000. Esta variabilidad hace de Es el predictor con mayor potencial explicativo del monto a pagar.
La variable FR (Factor regional) tiene una media de 0.68 con valores entre 0 y 76.8, lo que indica que la mayoría de los permisos se otorgan en regiones con presión ecológica moderada, aunque existen casos extremos en regiones de alta intervención.
Las variables Cb (Coeficiente biótico) y V (Coeficiente de valoración) presentan medias cercanas a 1 con poca dispersión, dado que sus valores están definidos normativamente en rangos acotados, lo que se refleja en coeficientes de variación bajos. Gt (Grupo trófico) tiene una media de 0.73, con la mayoría de los valores concentrados en 0.80, correspondiente a mamíferos y aves que ocupan niveles tróficos superiores.
En conjunto, la alta variabilidad de MP y Es, combinada con la menor dispersión de los coeficientes normativos, sugiere que el número de especímenes es el factor con mayor influencia práctica sobre el monto final. Los valores de asimetría positiva en MP y Es, visibles en la Tabla 1, anticipan la necesidad de considerar transformaciones logarítmicas para el modelamiento de regresión.
7 Visualización de datos
La representación gráfica permite identificar visualmente la distribución, tendencia, dispersión y posibles relaciones entre variables. Se utilizan histogramas para variables cuantitativas continuas, diagramas de caja para detectar valores atípicos, diagramas de barras para variables categóricas y tablas cruzadas para analizar relaciones entre categorías.
7.1 Histogramas
El histograma divide el rango de los datos en intervalos y representa la frecuencia de cada uno mediante barras. Permite identificar la forma de la distribución: si es simétrica, sesgada, unimodal o multimodal. Se genera un histograma para cada variable cuantitativa del análisis.
ggplot(tabla_analisis, aes(x = MP)) +
geom_histogram(fill = "#21618C", color = "black", bins = 30) +
scale_x_log10(labels = scales::comma) + # Se transforma el eje X con log10
theme_minimal() +
labs(title = "Monto a pagar (MP) (Escala Logarítmica)")
ggplot(tabla_analisis, aes(x = Es)) +
geom_histogram(fill = "#2E86C1", color = "black", bins = 30) +
scale_x_log10(labels = scales::comma) + # Se transforma el eje X con log10
theme_minimal() +
labs(title = "Número de especímenes (Es) (Escala Logarítmica)")
ggplot(tabla_analisis, aes(x = FR)) +
geom_histogram(fill = "#F1C40F", color = "black", bins = 30) +
scale_x_log10(labels = scales::comma) + # Se transforma el eje X con log10
theme_minimal() +
labs(title = "Factor Regional (FR) (Escala Logarítmica)")
ggplot(tabla_analisis, aes(x = Cb)) +
geom_histogram(fill = "#5499C7", color = "black", bins = 30) +
scale_x_log10(labels = scales::comma) + # Se transforma el eje X con log10
theme_minimal() +
labs(title = "Coeficiente Biótico (Cb) (Escala Logarítmica)")
ggplot(tabla_analisis, aes(x = V)) +
geom_histogram(fill = "#154360", color = "black", bins = 30) +
scale_x_log10(labels = scales::comma) + # Se transforma el eje X con log10
theme_minimal() +
labs(title = "Coeficiente de Valoración (V) (Escala Logarítmica)")
Interpretación de los histogramas:
Se aplicó transformación logarítmica (scale_x_log10) en todos los histogramas dado que las variables presentan rangos muy amplios y distribuciones asimétricas que en escala lineal comprimían la mayoría de los datos en un extremo, impidiendo una visualización adecuada.
MP (Monto a pagar): En escala logarítmica se observa una distribución con concentración entre $5.000 y $10.000 COP, con una cola hacia valores superiores al millón. Esto confirma que la mayoría de los permisos generan montos relativamente bajos, mientras que un grupo pequeño de permisos con grandes cantidades de especímenes genera montos extremadamente altos. La distribución es claramente asimétrica positiva incluso en escala logarítmica.
Es (Número de especímenes): La gran mayoría de los permisos autorizan entre 1 y 10 especímenes, con frecuencias decrecientes a medida que aumenta el número. Los casos con más de 100 especímenes son excepcionales pero existen, y son los responsables de los montos más altos del conjunto de datos.
FR (Factor regional): Los valores se concentran en dos picos principales alrededor de 0.48 y 1.0, que corresponden a las categorías normativas más frecuentes. La escala logarítmica revela también la presencia de valores en el rango 0.01–100, lo que refleja la diversidad de condiciones ecológicas regionales representadas en los datos.
Cb (Coeficiente biótico): La distribución muestra que la inmensa mayoría de los registros (más del 85%) tienen Cb = 1, con una minoría en valores 3 y 5. Esto indica que la mayoría de las especies autorizadas para caza no presentan alto grado de vulnerabilidad según los criterios normativos.
V (Coeficiente de valoración): Casi la totalidad de los valores se concentra en V = 1.0, con muy pocos registros en valores distintos. Esto refleja que la caza científica, a la que le corresponde este coeficiente, es la modalidad ampliamente dominante en el conjunto de datos analizados.
7.2 Diagrama de caja y bigotes
El diagrama de caja y bigotes resume cinco estadísticos clave: valor mínimo, primer cuartil (Q1), mediana (Q2), tercer cuartil (Q3) y valor máximo. Los puntos por fuera de los bigotes representan valores atípicos (outliers). Es especialmente útil para detectar la presencia de datos extremos y evaluar la simetría de la distribución.
ggplot(tabla_analisis, aes(x = "", y = MP)) +
geom_boxplot(fill = "#2E86C1", alpha = 0.7, outlier.color = "#154360", outlier.shape = 1, outlier.size = 2, outlier.alpha = 0.5) +
# Transformación logarítmica y formato de moneda/comas
scale_y_log10(labels = scales::label_comma()) +
theme_minimal() +
labs(title = "Distribución del Monto a Pagar (MP)",
subtitle = "Escala logarítmica para visualizar la dispersión de valores atípicos",
x = "",
y = "Monto a Pagar (COP)")
Interpretación del diagrama de caja y bigotes:
Se aplicó escala logarítmica en el eje Y (scale_y_log10) dado que en escala lineal los valores extremos de MP comprimían completamente la caja, haciendo ilegible el gráfico. En escala logarítmica se aprecia que el rango intercuartílico (caja) se ubica aproximadamente entre $10.000 y $30.000 COP, con la mediana cercana a $11.000 COP. Los bigotes se extienden desde aproximadamente $1.000 hasta $100.000 COP, cubriendo el rango no atípico de los datos. Los puntos por encima del bigote superior corresponden a outliers que superan el millón de pesos, asociados a permisos con un número muy elevado de especímenes autorizados. La asimetría positiva sigue siendo evidente incluso en escala logarítmica, lo que refuerza la necesidad de considerar transformaciones de variables antes del modelamiento de regresión.
7.3 Tablas de frecuencias — Variables categóricas
Las tablas de frecuencias permiten conocer cuántos registros corresponden a cada categoría de las variables cualitativas del conjunto de datos original. Se analizan el tipo de caza, la clase taxonómica de las especies y los departamentos con mayor número de permisos otorgados.
Tipo de caza:
##
## Científica estudios ambientales Científica no comercial
## 1370 90
## Comercial Fomento
## 65 1Permisos por clase taxonómica:
##
## Aves No registra Mammalia
## 550 426 145
## Insecta Amphibia Sauropsida
## 142 94 56
## Reptilia Actinopterygii GB
## 44 37 8
## Anfibia Teleostei Aves, Reptiles
## 5 5 2
## insecta Anphibia aves
## 2 1 1
## Aves, Herpetos aves, herpetos, mamiferos herpetos, mamiferos
## 1 1 1
## mamalia
## 1Permisos por departamento (top 10):
##
## Tolima Risaralda Caldas Valle del Cauca Quindío
## 603 495 347 16 14
## Boyacá Cauca No registra La Guajira Cundinamarca
## 12 10 9 7 47.4 Diagrama de barras — Tipo de caza
El diagrama de barras representa la frecuencia de cada categoría de la variable Tipo de Caza. Cada barra corresponde a una modalidad y su altura es proporcional al número de permisos otorgados, facilitando la comparación visual entre categorías.
datostasaCF %>%
count(`TIPO DE CAZA`) %>%
ggplot(aes(x = reorder(`TIPO DE CAZA`, n), y = n, fill = `TIPO DE CAZA`)) +
geom_col( show.legend = FALSE) +
scale_fill_manual(values = c("#AED6F1", "#A9DFBF", "#F9E79F", "#F5B7B1" )) +
coord_flip() +
theme_minimal() +
labs(title = "Número de permisos por tipo de caza",
x = "", y = "Cantidad de permisos")
Interpretación de tablas de frecuencias y barras:
Las tablas de frecuencias revelan que la caza científica es la modalidad más frecuente. En cuanto a la clase taxonómica, los mamíferos e insectos concentran la mayor cantidad de permisos. A nivel geográfico, los departamentos con mayor número de autorizaciones corresponden a regiones con mayor actividad investigativa y presencia institucional ambiental.
8 Análisis de correlación
El análisis de correlación mide la fuerza y dirección de la relación
lineal entre las variables cuantitativas seleccionadas. Se calcula la
matriz de correlación de Pearson con cor(), se redondea a
dos decimales para facilitar su lectura, y se visualiza mediante
corrplot(). Posteriormente, pairs() y
ggpairs() permiten explorar gráficamente las relaciones
entre todas las variables simultáneamente.
8.1 Matriz de correlaciones
## MP FR Cb Gt V Es
## MP 1.00 0.18 0.10 -0.02 0.17 0.85
## FR 0.18 1.00 0.41 0.11 0.93 -0.02
## Cb 0.10 0.41 1.00 0.14 0.08 -0.02
## Gt -0.02 0.11 0.14 1.00 0.01 -0.12
## V 0.17 0.93 0.08 0.01 1.00 0.00
## Es 0.85 -0.02 -0.02 -0.12 0.00 1.00#Visualización de correlaciones
corrplot(matriz_cor,method = "color", addCoef.col = "black",main= "Matriz de correlaciones (Taza caza)", mar = c(0, 0, 1, 0))
8.2 Diagrama de dispersión
El diagrama de dispersión representa visualmente la relación entre cada par de variables mediante puntos en un plano cartesiano. Es una herramienta fundamental antes de aplicar el modelo de regresión, ya que permite detectar la dirección, forma e intensidad de las asociaciones entre variables.
8.3 Matriz de correlación multivariada
La función ggpairs() del paquete GGally
genera una matriz gráfica completa que combina en un solo panel: las
distribuciones individuales de cada variable (diagonal), los diagramas
de dispersión entre pares (triángulo inferior) y los coeficientes de
correlación exactos (triángulo superior). Es la visualización más
completa para el análisis multivariado.
Interpretación de la matriz de correlación:
La matriz de correlación permite analizar la relación lineal entre las variables cuantitativas del estudio. En particular, se observa que la variable MP presenta la correlación positiva más fuerte, siendo casi perfecta con Es (número de especímenes), lo que indica que a mayor número de individuos o muestras autorizados, mayor es el monto a pagar. Esta relación es la más relevante del conjunto y es consistente con la fórmula normativa, donde Es actúa como factor multiplicador directo del monto.
Contrario a la hipótesis normativa, las variables FR (-0.482), Cb (-0.470) y Gt (-0.600) presentan correlaciones negativas moderadas con el monto. Aunque la norma sugiere que a mayores coeficientes el cobro debería aumentar, los datos reflejan un comportamiento inverso. Es fundamental señalar que, aunque la tarifa mínima no es una variable independiente analizada en este modelo, su existencia en el marco normativo explica estas correlaciones negativas. El sistema de información aplica un “piso” económico o tarifa plana que neutraliza el efecto de los coeficientes FR y Cb, especialmente en la caza científica (modalidad dominante según el análisis de frecuencias). Al aplicarse cobros mínimos fijos a registros con coeficientes ecológicos altos, se rompe la proporcionalidad lineal esperada, generando la tendencia negativa observada.
Se identifica una correlación extremadamente alta de 0.955 entre el Factor Regional (FR) y el Coeficiente de Valoración (V). Desde la perspectiva de ingeniería, esto advierte sobre una redundancia de información; ambas variables aportan casi el mismo peso estadístico, lo que sugiere que el modelo de regresión podría simplificarse sin perder capacidad explicativa.
En conjunto, estos resultados indican que Es es la variable con mayor potencial predictivo del monto a pagar (MP). La baja correlación entre la mayoría de las variables independientes favorece la construcción de un modelo de regresión estable, aunque la fuerte asociación entre FR y V requiere especial atención para evitar la multicolinealidad.
Rangos de referencia:
- Correlaciones muy fuertes: |r| ≥ 0.80
- Correlaciones fuertes: 0.60 ≤ |r| < 0.80
- Correlaciones moderadas: 0.40 ≤ |r| < 0.60
- Correlaciones débiles: 0.20 ≤ |r| < 0.40
- Correlaciones muy débiles o nulas: |r| < 0.20
8.4 Formulación de hipótesis
La prueba cor.test() evalúa formalmente si la
correlación entre dos variables es estadísticamente significativa. Se
plantean una hipótesis nula (H₀) y una hipótesis alternativa (H₁). Si el
p-value obtenido es menor a 0.05, se rechaza H₀ y se concluye que existe
evidencia de relación lineal entre las variables. Se evalúa la
correlación entre Es y MP por ser la
relación más relevante del estudio.
##
## Pearson's product-moment correlation
##
## data: tabla_analisis$Es and tabla_analisis$MP
## t = 54.559, df = 1139, p-value < 0.00000000000000022
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.8335486 0.8657489
## sample estimates:
## cor
## 0.8504434Hipótesis:
- H₀: No existe relación lineal entre Es y MP (ρ = 0)
- H₁: Existe relación lineal entre Es y MP (ρ ≠ 0)
Interpretación: Si p-value < 0.05, se rechaza H₀ y se concluye que existe evidencia estadística de una relación lineal significativa entre el número de especímenes autorizados y el monto a pagar por concepto de la tasa compensatoria.
9 Evaluación de normalidad
La normalidad se evalúa sobre la variable respuesta (Monto a Pagar) como uno de los supuestos fundamentales del modelo de regresión lineal. Verificar si la variable dependiente sigue una distribución normal permite determinar la validez de las inferencias estadísticas del modelo y anticipar la necesidad de transformaciones sobre los datos.
9.1 Prueba de Kolmogorov-Smirnov
La prueba de Kolmogorov-Smirnov (K-S) es una prueba no paramétrica
que evalúa si una variable cuantitativa sigue una distribución teórica
específica, en este caso la distribución normal. A diferencia de la
prueba de Shapiro-Wilk, diseñada para muestras pequeñas, la prueba K-S
es más adecuada para conjuntos de datos con un mayor número de
observaciones como el presente estudio (n = 1.141). La función
ks.test() compara la distribución empírica acumulada de los
datos con la distribución normal teórica construida a partir de la media
y desviación estándar de la propia variable.
El resultado se interpreta mediante el estadístico D y el valor p asociado.
9.1.1 Hipótesis
-H0: La variable sigue una distribución normal. -H1: La variable no sigue una distribución normal.
##
## Asymptotic one-sample Kolmogorov-Smirnov test
##
## data: tabla_analisis$MP
## D = 0.41574, p-value < 0.00000000000000022
## alternative hypothesis: two-sidedInterpretación de la Prueba de Kolmogorov-Smirnov:
El estadístico D = 0.41574 indica la máxima diferencia entre la función de distribución acumulada empírica de los datos y la función de distribución normal teórica. Un valor de 0.416 es considerablemente alto, evidenciando que los datos se alejan notablemente de la distribución normal. El p-value < 2.2e-16 es prácticamente cero, muy por debajo del nivel de significancia α = 0.05, por lo que se rechaza H₀ y se concluye que la variable MP no sigue una distribución normal. Este resultado es consistente con lo observado en el histograma y el diagrama de caja, donde se evidenció una fuerte asimetría positiva y la presencia de valores atípicos extremos.
Este hallazgo es relevante para el modelamiento, ya que sugiere la conveniencia de evaluar transformaciones logarítmicas sobre la variable dependiente en etapas posteriores del análisis.
10 División en entrenamiento y prueba
Antes de ajustar el modelo de regresión, el conjunto de datos se
divide en dos subconjuntos: datos de entrenamiento
(70%) utilizados para estimar los coeficientes del modelo, y
datos de prueba (30%) reservados para evaluar su
capacidad predictiva sobre observaciones no vistas. La función
set.seed(123) garantiza que la división aleatoria sea
reproducible en cualquier ejecución del código.
set.seed(123)
muestra <- sample(
1:nrow(tabla_analisis),
size = 0.7*nrow(tabla_analisis)
)
datos_entrena <- tabla_analisis[muestra, ]
datos_test <- tabla_analisis[-muestra, ]
dim(datos_entrena)## [1] 798 6## [1] 343 6Interpretación de la división:
El conjunto de entrenamiento quedó conformado por 798 observaciones y 6 variables, correspondiente al 70% del total. El conjunto de prueba contiene 343 observaciones, equivalente al 30% restante. Esta proporción 70/30 es ampliamente utilizada en modelamiento estadístico y garantiza que el modelo tenga suficientes datos para estimar los coeficientes de forma estable, a la vez que se reserva una muestra representativa para evaluar su desempeño predictivo.
11 Ajuste de modelos de regresión
Se ajustan dos modelos de regresión lineal múltiple sobre los datos de entrenamiento. La comparación entre ambos modelos permite identificar el modelo con mejor ajuste.
11.1 Modelo de regresión múltiple
Modelo 1 — todas las variables:
##
## Call:
## lm(formula = MP ~ FR + Cb + Gt + V + Es, data = datos_entrena)
##
## Residuals:
## Min 1Q Median 3Q Max
## -525099 -29775 -25156 10201 1715719
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -263703.26 80259.91 -3.286 0.001062 **
## FR -26737.00 16553.54 -1.615 0.106669
## Cb 56740.34 15593.81 3.639 0.000292 ***
## Gt 89777.36 21553.23 4.165 0.0000345 ***
## V 181969.32 64630.70 2.816 0.004991 **
## Es 1013.13 19.45 52.082 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 142100 on 792 degrees of freedom
## Multiple R-squared: 0.7839, Adjusted R-squared: 0.7826
## F-statistic: 574.7 on 5 and 792 DF, p-value: < 0.00000000000000022Interpretación Modelo 1:
El modelo con las 5 variables independientes presenta un R² = 0.7839, explicando el 78.39% de la variabilidad del Monto a Pagar. La prueba F global es significativa (F = 574.7, p < 2.2e-16), confirmando que el modelo en conjunto tiene poder explicativo. Sin embargo, la variable FR no es estadísticamente significativa (p = 0.107), lo que sugiere que su aporte al modelo es marginal cuando las demás variables están incluidas. Las variables Es, Gt, Cb y V sí resultan significativas, siendo Es la de mayor relevancia con un coeficiente de 1.013 (p < 2.2e-16).
Modelo 2 — sin la variable FR:
##
## Call:
## lm(formula = MP ~ Cb + Gt + V + Es, data = datos_entrena)
##
## Residuals:
## Min 1Q Median 3Q Max
## -524683 -30021 -25403 9950 1737919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -137241.27 17665.34 -7.769 0.0000000000000245 ***
## Cb 33114.79 5410.38 6.121 0.0000000014646584 ***
## Gt 75098.71 19563.40 3.839 0.000134 ***
## V 78288.35 7529.97 10.397 < 0.0000000000000002 ***
## Es 1013.55 19.47 52.055 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 142300 on 793 degrees of freedom
## Multiple R-squared: 0.7832, Adjusted R-squared: 0.7821
## F-statistic: 716.3 on 4 and 793 DF, p-value: < 0.00000000000000022Interpretación Modelo 2:
Al excluir FR, el modelo mantiene un R² = 0.7832, prácticamente idéntico al Modelo 1 (diferencia de apenas 0.0007), lo que confirma que FR no aportaba información relevante. El R² ajustado = 0.7821 penaliza correctamente la reducción de una variable, siendo ligeramente inferior al del Modelo 1 pero con la ventaja de mayor simplicidad. Todas las variables incluidas son altamente significativas (p < 0.001): Es (coef = 1.013, p < 2.2e-16), V (coef = 78.288, p < 2.2e-16), Cb (coef = 33.115, p < 0.0001) y Gt (coef = 75.099, p = 0.000134). El Modelo 2 es preferible al Modelo 1 por ser más simple sin pérdida significativa de ajuste.
12 Interpretación del modelo
La ecuación estimada del modelo seleccionado (Modelo 2) se obtiene
con la función coef(), que extrae los coeficientes
estimados por mínimos cuadrados ordinarios.
## (Intercept) Cb Gt V Es
## -137241.27 33114.79 75098.71 78288.35 1013.55Ecuación del modelo estimado:
MP = -137241.27 + 33114.79·Cb + 75098.71·Gt + 78288.35·V + 1013.55·Es
Interpretación de coeficientes:
- Intercepto (-137241.27): valor estimado de MP cuando todos los predictores son cero. No tiene interpretación práctica directa en el contexto del sistema de liquidación.
- Cb (33114.79): por cada unidad adicional en el Coeficiente Biótico, el monto a pagar aumenta en promedio $33.115 COP, lo que refleja que las especies más vulnerables generan tasas más altas.
- Gt (75098.71): por cada unidad adicional en el Grupo Trófico, el monto aumenta en promedio $75.099 COP, evidenciando que las especies en niveles tróficos superiores (depredadores) conllevan tasas más elevadas.
- V (78288.35): por cada unidad adicional en el Coeficiente de Valoración, el monto aumenta en promedio $78.288 COP, reflejando el impacto de la finalidad de la caza sobre el monto liquidado.
- Es (1013.55): por cada espécimen adicional autorizado, el monto a pagar aumenta en promedio $1.013 COP. Es la variable con mayor impacto acumulado dado el amplio rango de valores que puede tomar.
13 Evaluación del ajuste del modelo
Una vez ajustado el modelo, se evalúa su capacidad explicativa mediante el coeficiente de determinación R² y el R² ajustado. Estas métricas permiten cuantificar qué proporción de la variabilidad total del Monto a Pagar es explicada por las variables independientes incluidas en el modelo.
13.1 Coeficiente de determinación
El R² representa el porcentaje de variabilidad de MP explicado por el modelo. Toma valores entre 0 y 1, donde valores más cercanos a 1 indican mejor ajuste.
## [1] 0.7832305Interpretación: El R² = 0.7832 indica que el modelo explica el 78.32% de la variabilidad total del Monto a Pagar, lo cual representa un ajuste alto y satisfactorio considerando la complejidad y variabilidad del fenómeno estudiado.
13.2 R cuadrado ajustado
El R² ajustado penaliza la inclusión innecesaria de variables en el modelo, corrigiendo el R² según el número de predictores y el tamaño muestral. Es la métrica más adecuada para comparar modelos con diferente número de variables.
## [1] 0.7821371Interpretación: El R² ajustado = 0.7821 es muy cercano al R², lo que confirma que las 4 variables incluidas en el Modelo 2 son todas relevantes y que el modelo no presenta sobreajuste. La diferencia mínima entre R² y R² ajustado (0.0011) indica que ninguna variable está inflando artificialmente el ajuste del modelo.
14 Evaluación de hipótesis del modelo
14.1 Prueba F global
Evalúa si el modelo es significativo.
14.1.1 Hipótesis
- \(H_0\): Todos los coeficientes son iguales a cero.
- \(H_1\): Al menos un coeficiente es diferente de cero.
##
## Call:
## lm(formula = MP ~ Cb + Gt + V + Es, data = datos_entrena)
##
## Residuals:
## Min 1Q Median 3Q Max
## -524683 -30021 -25403 9950 1737919
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -137241.27 17665.34 -7.769 0.0000000000000245 ***
## Cb 33114.79 5410.38 6.121 0.0000000014646584 ***
## Gt 75098.71 19563.40 3.839 0.000134 ***
## V 78288.35 7529.97 10.397 < 0.0000000000000002 ***
## Es 1013.55 19.47 52.055 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 142300 on 793 degrees of freedom
## Multiple R-squared: 0.7832, Adjusted R-squared: 0.7821
## F-statistic: 716.3 on 4 and 793 DF, p-value: < 0.00000000000000022Interpretación: El estadístico F = 716.3 con 4 y 793 grados de libertad y un p-value < 2.2e-16 indica que el modelo en su conjunto es altamente significativo. Se rechaza \(H_0\) y se concluye que al menos una de las variables independientes (Cb, Gt, V, Es) tiene un efecto lineal real y significativo sobre el Monto a Pagar. Este resultado valida la estructura del modelo y justifica el uso de la regresión lineal múltiple para explicar el comportamiento de MP.
15 Evaluación individual de coeficientes
Cada predictor tiene una prueba t.
15.1 Hipótesis
- \(H_0 : \beta_i = 0\)
- \(H_1 : \beta_i \neq 0\)
Interpretación: Todas las variables incluidas en el Modelo 2 son estadísticamente significativas:
- Es (coef = 1013.55, p < 2.2e-16 ***): es la variable con mayor significancia. Por cada espécimen adicional autorizado, el monto a pagar aumenta en promedio $1.013 COP. Su altísimo estadístico t (52.055) confirma que es el predictor dominante del modelo.
- V (coef = 78288.35, p < 2.2e-16 ***): por cada unidad adicional en el Coeficiente de Valoración, el monto aumenta en promedio $78.288 COP, reflejando el fuerte impacto de la finalidad de la caza.
- Cb (coef = 33114.79, p < 0.0001 ***): por cada unidad adicional en el Coeficiente Biótico, el monto aumenta en $33.115 COP, consistente con que especies más vulnerables generan tasas más altas.
- Gt (coef = 75098.71, p = 0.000134 ***): por cada unidad adicional en el Grupo Trófico, el monto aumenta en $75.099 COP, indicando que especies en niveles tróficos superiores conllevan tasas más elevadas.
16 Predicciones
La función predict() aplica los coeficientes estimados
del modelo a los datos de prueba para generar valores predichos de MP.
Se construye una tabla comparativa que incluye el valor real, el valor
predicho, el error simple y el error absoluto para evaluar la calidad de
las predicciones de forma individual.
predicciones <- predict(modelo2, newdata = datos_test)
comparativa <- data.frame(
Real = datos_test$MP,
Predicho = predicciones,
Error = datos_test$MP - predicciones,
Error_Absoluto = abs(datos_test$MP - predicciones)
)
set.seed(123)
muestra_reporte <- comparativa[sample(1:nrow(comparativa), 6), ]
rownames(muestra_reporte) <- NULL
print("Muestra de predicciones - Modelo 2:")## [1] "Muestra de predicciones - Modelo 2:"Interpretación:Del análisis de las predicciones se observa que el modelo 2 presenta un desempeño aceptable en algunos casos, principalmente cuando el monto real es alto. Sin embargo, se identfica que el modelo genera predicciones negativas en algunos registros, lo cual carece de sentido práctico para un monto económico. Adicionalmente, existen errores absolutos elevados en varias situaciones; estos resultados son consistentes con la distribución fuertemente asimétrica y no normal de la variable MP, observada previamente.
17 Evaluación predictiva
La evaluación predictiva cuantifica el desempeño del modelo sobre datos no vistos (conjunto de prueba), permitiendo detectar problemas de sobreajuste y medir la precisión real de las estimaciones.
17.1 Error Cuadrático Medio (MSE)
El MSE mide el promedio de los errores al cuadrado entre los valores reales y los predichos. Penaliza con mayor fuerza los errores grandes. Se define como: \[MSE = \frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2\]
## [1] 1194982808517.2 Raíz del Error Cuadrático Medio (RMSE)
El RMSE es la raíz cuadrada del MSE y tiene la ventaja de expresarse en las mismas unidades que la variable dependiente (COP), lo que facilita su interpretación práctica: \[RMSE = \sqrt{MSE} = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(y_i - \hat{y}_i)^2}\]
## [1] 109315.317.3 Error Absoluto Medio (MAE)
El MAE mide el promedio de los errores absolutos, siendo menos sensible a valores extremos que el MSE. Se interpreta directamente en COP: \[MAE = \frac{1}{n}\sum_{i=1}^{n}|y_i - \hat{y}_i|\]
## [1] 43177.1417.4 Coeficiente R2 en prueba
El \(R^2\) en el conjunto de prueba indica qué proporción de la variabilidad de MP es explicada por el modelo sobre datos nuevos. Se calcula como: \[R^2 = 1 - \frac{SSE}{SST} = 1 - \frac{\sum(y_i - \hat{y}_i)^2}{\sum(y_i - \bar{y})^2}\]
sst <- sum((datos_test$MP -
mean(datos_test$MP))^2)
sse <- sum((datos_test$MP -
predicciones)^2)
r_squared <- 1 - (sse/sst)
r_squared## [1] 0.6650201Interpretación de la evaluación predictiva: El modelo presenta un R² de entrenamiento = 0.7832 que cae a 0.665 en el conjunto de prueba, evidenciando una pérdida de capacidad explicativa del 11.8% al generalizar a datos no vistos. Esto indica un leve sobreajuste del modelo a los datos de entrenamiento.
El RMSE = $109.315 COP representa el error típico de predicción en la misma escala que MP. Considerando que la mediana de MP es aproximadamente $11.000 COP, este error equivale a casi 10 veces la mediana, lo que confirma que el modelo tiene dificultades para predecir con precisión en la escala original de la variable.
El MAE = $43.177 COP es más representativo del error típico real, ya que no penaliza los errores extremos al cuadrado como el RMSE. La diferencia considerable entre RMSE ($109.315) y MAE ($43.177) confirma que existen algunos errores muy grandes que inflan el RMSE, mientras que la mayoría de las predicciones tienen errores menores. Esto es consistente con las predicciones negativas y los errores absolutos elevados detectados en la tabla comparativa.
En conjunto, estos resultados confirman que el modelo lineal en escala original tiene limitaciones importantes derivadas de la distribución asimétrica de MP y la presencia de observaciones influyentes, lo que justifica buscar alternativas de solución para obtener un modelo más estable y confiable.
18 Validación de supuestos
Los supuestos del modelo de regresión lineal múltiple se evalúan mediante pruebas formales y gráficos diagnósticos. El cumplimiento de estos supuestos es necesario para que los estimadores de mínimos cuadrados ordinarios sean los mejores estimadores lineales insesgados (teorema de Gauss-Markov).
18.1 Gráficos diagnósticos
Los cuatro gráficos diagnósticos estándar de R permiten evaluar visualmente los supuestos de linealidad, normalidad, homocedasticidad e identificación de observaciones influyentes.
Interpretación de gráficos: - Residuals vs Fitted: Se observa una clara concentración de residuos en valores bajos de predicción, con varios puntos extremos. La línea roja horizontal muestra que no hay una relación perfectamente lineal y existe cierta heteroscedasticidad (la dispersión no es constante).
Normal Q-Q: Se evidencia que los residuos se desvían notablemente de la línea diagonal, especialmente en la cola derecha (valores altos). Esto confirma la falta de normalidad detectada previamente.
Scale-Location: El gráfico nos muestra que la línea roja presenta una clara pendiente descendente, indicando que la varianza de los residuos disminuye a medida que aumentan los valores ajustados. Existe un problema importante de heterocedasticidad.
Residuals vs Leverage: Se identifican varios puntos con alto leverage y algunos con alta influencia (Cook’s distance). Estos puntos están afectando fuertemente el ajuste del modelo.
Gracias a los resultados mencionados, se determina que el modelo presenta violaciones importantes a los supuestos de normalidad y homocedasticidad, además de observaciones influyentes que justifican su tratamiento.
18.2 Normalidad de residuos
La prueba de Shapiro-Wilk evalúa formalmente si los residuos del modelo siguen una distribución normal. Es la prueba más potente para muestras de tamaño moderado. El estadístico W se define como:
18.2.1 Hipótesis
- \(H_0:\) Los residuos siguen una distribución normal
- \(H_1:\) Los residuos no siguen una distribución normal
Criterio de decisión:
- Si valor-\(p > 0.05\): No se rechaza \(H_0\). Existe evidencia de normalidad.
- Si valor-\(p < 0.05\): Se rechaza \(H_0\). Los residuos no son normales.
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo2)
## W = 0.33032, p-value < 0.00000000000000022Interpretación:
El estadístico W = 0.33032 con p-value < 2.2e-16 indica que se rechaza \(H_0\) con un nivel de confianza del 95%. Los residuos del modelo no siguen una distribución normal, lo cual es consistente con lo observado en el gráfico Q-Q y con la distribución asimétrica de la variable MP. Este resultado confirma que el supuesto de normalidad de residuos está violado en el Modelo 2.
18.3 Homocedasticidad
Prueba de no constancia de varianza (Non-Constant Variance Test), utilizada para evaluar el supuesto de homocedasticidad en modelos de regresión lineal.
18.3.1 Hipótesis
- \(H_0\): La varianza de los residuos es constante (homocedasticidad)
- \(H_1\): La varianza de los residuos no es constante (heterocedasticidad)
Criterio de decisión:
- Si valor-\(p > 0.05\): No se rechaza \(H_0\). Existe homocedasticidad.
- Si valor-\(p < 0.05\): Se rechaza \(H_0\). Existe heterocedasticidad.
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 33.91343, Df = 1, p = 0.0000000057619Interpretación:
El estadístico Chisquare = 33.913 con p = 0.0000000057619 indica que se rechaza \(H_0\). Existe evidencia estadística significativa de heterocedasticidad en los residuos del modelo. La varianza de los errores no es constante, lo que viola uno de los supuestos fundamentales de la regresión lineal y puede afectar la validez de las inferencias sobre los coeficientes estimados.
18.4 Independencia de residuos
Se avalua con el test de Durvin Watson.
18.4.1 Hipótesis
- \(H_0\): Los residuos son independientes (\(\rho = 0\))
- \(H_1\): Existe autocorrelación (\(\rho \neq 0\))
Criterio de decisión:
- Si valor-\(p > 0.05\): No se rechaza \(H_0\). Los residuos son independientes.
- Si valor-\(p < 0.05\): Se rechaza \(H_0\). Existe autocorrelación.
## lag Autocorrelation D-W Statistic p-value
## 1 -0.02201962 2.044008 0.408
## Alternative hypothesis: rho != 0Interpretación:
El estadístico D-W = 2.044 con p-value = 0.408 indica que no se rechaza \(H_0\). No existe evidencia de autocorrelación en los residuos del modelo. Este es el único supuesto que se cumple satisfactoriamente en el Modelo 2, lo que indica que las observaciones son independientes entre sí, como era esperado dado que los registros corresponden a permisos distintos sin relación temporal secuencial.
18.5 Multicolinealidad
La multicolinealidad se evalúa mediante el Factor de Inflación de la Varianza (VIF).
Criterio de interpretación:
- \(VIF < 5\): baja multicolinealidad, aceptable.
- \(5 \leq VIF \leq 10\): multicolinealidad moderada, requiere atención.
- \(VIF > 10\): multicolinealidad problemática, compromete la estabilidad de los coeficientes.
## Cb Gt V Es
## 1.030521 1.030853 1.010044 1.010654Interpretación:
Los valores VIF obtenidos son: Cb = 1.031, Gt = 1.031, V = 1.010, Es = 1.011, todos muy cercanos a 1.0 y muy por debajo del umbral de 5. Esto indica una ausencia total de multicolinealidad entre las variables independientes del modelo. Cada predictor aporta información estadísticamente independiente de los demás, lo que garantiza la estabilidad en la estimación de los coeficientes y facilita su interpretación individual. Este resultado es el supuesto mejor cumplido del Modelo 2.
19 Observaciones influyentes
El leverage (apalancamiento) mide cuánto influye cada observación sobre los valores ajustados del modelo en función de sus valores en las variables independientes. Una observación tiene alto leverage cuando sus valores en X son inusuales respecto al resto del conjunto de datos.
Umbral de leverage alto (regla práctica):
\[h_{ii} > \frac{2(p+1)}{n}\]
donde \(p\) es el número de predictores y \(n\) el número de observaciones.
Interpretación:
- Valores pequeños de leverage: observaciones cercanas al centro de los datos, poca influencia.
- Valores grandes de leverage: observaciones alejadas del patrón general, pueden distorsionar los coeficientes.
leverage <- hatvalues(modelo2)
umbral <- 2*(length(coef(modelo2))/nrow(tabla_analisis))
puntos_alto_leverage <- which(leverage > umbral)
puntos_alto_leverage## 12 15 25 45 87 89 95 104 113 145 147 155 170 183 184 185 188 192 197 202
## 12 15 25 45 87 89 95 104 113 145 147 155 170 183 184 185 188 192 197 202
## 204 216 230 231 240 241 246 257 261 263 266 275 280 296 299 331 335 350 357 363
## 204 216 230 231 240 241 246 257 261 263 266 275 280 296 299 331 335 350 357 363
## 377 383 384 396 416 421 423 440 443 471 480 490 512 517 520 523 524 533 542 550
## 377 383 384 396 416 421 423 440 443 471 480 490 512 517 520 523 524 533 542 550
## 563 575 592 611 620 626 629 654 658 664 673 693 701 708 732 742 746 757 766 774
## 563 575 592 611 620 626 629 654 658 664 673 693 701 708 732 742 746 757 766 774
## 782 788
## 782 78819.1 Gráfico de leverage
plot(leverage,
main = "Valores de leverage",
ylab = "Leverage",
xlab = "Observaciones")
abline(h = umbral,
col = "red",
lty = 2)
Interpretación del leverage:
El gráfico revela que la gran mayoría de las observaciones tienen leverage muy bajo (cercano a 0), pero existen dos observaciones extremas con leverage cercano a 1.0 (aproximadamente las observaciones 592 y 673 del conjunto de entrenamiento). Un leverage de 1.0 indica que esa observación tiene una influencia absoluta sobre su propio valor ajustado, lo que representa una situación problemática para la estabilidad del modelo. Estas observaciones corresponden muy probablemente a permisos con un número de especímenes extremadamente alto, alejado del patrón general del conjunto de datos.
19.2 Distancia de Cook
La distancia de Cook combina el leverage y el tamaño del residuo para medir la influencia global de cada observación sobre todos los coeficientes estimados del modelo simultáneamente.
Umbral de referencia:
\[D_i > \frac{4}{n}\]
Interpretación:
- \(D_i \approx 0\): observación sin influencia importante sobre el modelo.
- \(D_i > \frac{4}{n}\): posible observación influyente que merece análisis.
- \(D_i > 1\): observación altamente influyente.
cooks_dist <- cooks.distance(modelo2)
puntos_influyentes <- which(
cooks_dist > 4/nrow(tabla_analisis)
)
puntos_influyentes## 25 37 76 104 145 155 192 221 241 296 312 335 422 423 443 457 475 490 518 557
## 25 37 76 104 145 155 192 221 241 296 312 335 422 423 443 457 475 490 518 557
## 592 673 796
## 592 673 796Interpretación de la distancia de Cook:
Se identificaron 23 observaciones influyentes con distancia de Cook superior al umbral \(\frac{4}{n}\), incluyendo las observaciones 25, 37, 76, 104, 145, 155, 192, 221, 241, 296, 312, 335, 422, 423, 443, 457, 475, 490, 518, 557, 592, 673 y 796. Estas observaciones modifican considerablemente los coeficientes estimados, las predicciones y el comportamiento general del modelo cuando son incluidas. Las observaciones 592 y 673 resultan especialmente críticas al combinar alto leverage con residuos elevados, cruzando las líneas de distancia de Cook en el gráfico Residuals vs Leverage. Su presencia en el conjunto de entrenamiento justifica plenamente el tratamiento realizado en la siguiente sección mediante la eliminación de puntos influyentes para obtener un modelo más estable y confiable.
20 Transformación logarítmica
Debido a los resultados obtenidos anteriormente y a la no normalidad de la variable MP, se busca probar la aplicación de una transformación logarítmica (log(MP)) sobre la variable dependiente con el fin de mejorar el ajuste y la capacidad predictiva del modelo seleccionado.
# Crear variable transformada
datos_entrena$logMP <- log(datos_entrena$MP + 1) # +1 por si hay ceros
datos_test$logMP <- log(datos_test$MP + 1)
# Ajustar modelo con log
modelo_log <- lm(logMP ~ Cb + Gt + V + Es, data = datos_entrena)
summary(modelo_log)##
## Call:
## lm(formula = logMP ~ Cb + Gt + V + Es, data = datos_entrena)
##
## Residuals:
## Min 1Q Median 3Q Max
## -1.9218 -0.9689 -0.2772 0.6379 4.7248
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 7.0662642 0.1420970 49.728 < 0.0000000000000002 ***
## Cb 0.5273956 0.0435202 12.118 < 0.0000000000000002 ***
## Gt 2.2608874 0.1573646 14.367 < 0.0000000000000002 ***
## V 0.1656335 0.0605698 2.735 0.00639 **
## Es 0.0013088 0.0001566 8.356 0.000000000000000288 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 1.144 on 793 degrees of freedom
## Multiple R-squared: 0.3726, Adjusted R-squared: 0.3695
## F-statistic: 117.8 on 4 and 793 DF, p-value: < 0.00000000000000022# Predicciones (en escala log)
pred_log <- predict(modelo_log, newdata = datos_test)
# Volver a escala original
predicciones_original <- exp(pred_log) - 1
# Evaluación en escala original
comparativa_log <- data.frame(
Real = datos_test$MP,
Predicho = predicciones_original,
Error = datos_test$MP - predicciones_original,
Error_Abs = abs(datos_test$MP - predicciones_original)
)
# Métricas de desempeño
RMSE_log <- sqrt(mean(comparativa_log$Error^2))
MAE_log <- mean(comparativa_log$Error_Abs)
cat("RMSE Modelo Log:", round(RMSE_log, 2), "\n")## RMSE Modelo Log: 187383.2## MAE Modelo Log: 45729.6# R² en escala original
sst <- sum((datos_test$MP - mean(datos_test$MP))^2)
sse <- sum((datos_test$MP - predicciones_original)^2)
r2_original <- 1 - (sse/sst)
cat("R² en escala original:", round(r2_original, 4), "\n")## R² en escala original: 0.0157Interpretación: Al evaluar el modelo en la escala original, se obtuvo un \(R^2\) de solo 0.0157 y un RMSE de 187.383 COP, valores considerablemente peores que los del modelo 2 original (\(R^2\) = 0.665 en prueba).
Los resultados obtenidos indican que, aunque la transformación logarítmica es una técnica recomendada en presencia de asimetría, en este caso específico no mejoró la capacidad predictiva del modelo. Por lo tanto, se decide buscar otra alternativa de solución.
21 Tratamiento de observaciones influyentes
Los gráficos diagnósticos y la distancia de Cook identificaron observaciones que ejercen una influencia desproporcionada sobre los coeficientes del modelo. Como alternativa de solución, se opta por eliminar estas observaciones del conjunto de análisis, dado que su presencia viola los supuestos de normalidad de residuos y homocedasticidad, comprometiendo la validez de las inferencias estadísticas.
# Identificar observaciones influyentes por Distancia de Cook
cooks_dist <- cooks.distance(modelo2)
umbral_cook <- 4 / nrow(tabla_analisis)
# Observaciones influyentes
puntos_influyentes <- which(cooks_dist > umbral_cook)
cat("Número de observaciones influyentes:", length(puntos_influyentes), "\n")## Número de observaciones influyentes: 23# Crear nuevo dataset sin las observaciones influyentes
datos_limpios <- tabla_analisis[-puntos_influyentes, ]
cat("Observaciones originales:", nrow(tabla_analisis), "\n")## Observaciones originales: 1141## Observaciones después de eliminar influyentes: 1118# Volver a dividir en entrenamiento y prueba
set.seed(123)
muestra_limpia <- sample(1:nrow(datos_limpios), size = 0.7*nrow(datos_limpios))
datos_entrena_limpio <- datos_limpios[muestra_limpia, ]
datos_test_limpio <- datos_limpios[-muestra_limpia, ]
# Ajustar nuevo modelo sin observaciones influyentes
modelo_limpio <- lm(MP ~ Cb + Gt + V + Es,
data = datos_entrena_limpio)
# Predicciones y evaluación
predicciones_limpio <- predict(modelo_limpio, newdata = datos_test_limpio)
# Métricas
comparativa_limpio <- data.frame(
Real = datos_test_limpio$MP,
Predicho = predicciones_limpio,
Error = datos_test_limpio$MP - predicciones_limpio
)Interpretación: Se eliminaron 23 observaciones, quedando un total de 1.118 registros.
22 Evaluación de hipótesis del modelo limpio
22.1 Prueba F global
##
## Call:
## lm(formula = MP ~ Cb + Gt + V + Es, data = datos_entrena_limpio)
##
## Residuals:
## Min 1Q Median 3Q Max
## -527290 -28382 -23528 95 1582265
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -109778.8 14519.6 -7.561 0.000000000000113 ***
## Cb 21972.6 4666.8 4.708 0.000002958426153 ***
## Gt 50722.4 16050.0 3.160 0.00164 **
## V 80051.7 6017.2 13.304 < 0.0000000000000002 ***
## Es 991.4 38.0 26.092 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 113900 on 777 degrees of freedom
## Multiple R-squared: 0.5344, Adjusted R-squared: 0.532
## F-statistic: 222.9 on 4 and 777 DF, p-value: < 0.00000000000000022Interpretación: Dado que el valor -p de la prueba F es < 0.05 (prácticamente 0), se rechaza la hipótesis nula. Por lo tanto, el modelo en su conjunto es estadísticamente significativo.
22.2 Evaluación individual de coeficientes
22.2.1 Hipótesis
- \(H_0 : \beta_i = 0\)
- \(H_1 : \beta_i \neq 0\)
Interpretación: Todas las variables incluidas en el modelo (Cb, Gt, V y Es) presentan un valor -p < 0.05, razón por la cual se rechaza la hipótesis nula en cada caso. Todas las variables aportan significativamente a la explicación del monto a pagar.
23 Evaluación predictiva - Modelo limpio
23.1 Error Cuadrático Medio (MSE)
## [1] 3009848299223.2 Raiz del Error Cuadrático Medio (RMSE)
## RMSE: 173489.1Interpretación: El Error Cuadrático Medio (MSE) obtenido es de 30.098.482.992. Este valor representa el promedio de los errores de predicción elevados al cuadrado, en este caso no es recomendable trabajar con MSE ya que está en pesos colombianos al cuadrado, por lo que es díficil interpretarlo directamente. Por esta razón, se prefiere reportar el RMSE (Raíz del Error Cuadrático Medio) que indica, en promedio, que las predicciones del modelo se desvían en aproximadamente 173.489 COP de los valores reales.
23.3 Error Absoluto Medio (MAE)
Es el promedio de los errores absolutos (sin elevar al cuadrado) entre los valores reales y los valores predichos por el modelo.
## MAE: 53137.16Interpretación: El Error Absoluto Medio indica que, en promedio, el modelo se equivoca en $53.137 COP por cada predicción realizada sobre el conjunto de prueba.
23.4 Coeficiente R² en prueba
sst <- sum((datos_test_limpio$MP - mean(datos_test_limpio$MP))^2)
sse <- sum((datos_test_limpio$MP - predicciones_limpio)^2)
r_squared_limpio <- 1 - (sse/sst)
r_squared_limpio## [1] 0.8426113Interpretación: Según los resultados obtenidos, el modelo limpio presenta un gran poder predictivo, explicando el 84.26% de la variabilidad del monto a pagar en el conjunto de prueba.
24 Validación de supuestos - Modelo limpio
24.1 Gráficos diagnósticos
Interpretación de gráficos:
Residuals vs Fitted: Se observa una mejora notable respecto al modelo original. Los residuos se distribuyen de forma más aleatoria alrededor de cero, aunque aún se observan algunos valores extremos. Se considera que la dispersión es más uniforme en comparación con el modelo original.
Normal Q-Q: Se determina una mejora considerable en la cola izquierda, aunque todavía existen desviaciones en la cola derecha. La normalidad sigue sin cumplirse completamente, pero es mejor que en el modelo original.
Scale-Location: Al analizar el gráfico, la línea roja es más horizontal que en el modelo original, lo que indica una reducción importante en el problema de heterocedasticidad, aunque persiste una leve tendencia.
Residuals vs Leverage: Tras la eliminación de las observaciones influyentes, ya no se observan puntos con influencia extrema (Cook’s distance). Los puntos con mayor leverage (3930, 595, etc.) tienen menos impacto en el modelo.
Gracias a los resultados mencionados, se evidencia que la eliminación de las 23 observaciones influyentes generó una mejora clara en el comportamiento de los residuos y en la estabilidad del modelo. Sin embargo, los supuestos de normalidad y homocedasticidad no se cumplen a cabalidad.
24.2 Normalidad de los residuos
24.2.1 Hipótesis
- \(H_0 :\) Los residuos siguen una distribución normal.
- \(H_1 :\) Los residuos no siguen una distribución normal.
##
## Shapiro-Wilk normality test
##
## data: residuals(modelo_limpio)
## W = 0.30038, p-value < 0.00000000000000022Interpretación: Dado que el p-value < 0.05, se rechaza la hipótesis nula. Los residuos no siguen una distribución normal.
24.3 Homocedasticidad
24.3.1 Hipótesis
- \(H_0 :\) Varianza constante (homocedasticidad).
- \(H_1 :\) Existe heterocedasticidad.
## Non-constant Variance Score Test
## Variance formula: ~ fitted.values
## Chisquare = 22.52014, Df = 1, p = 0.0000020795Interpretación: Dado que el p-value < 0.05, se rechaza H0. Existe evidencia de heterocedasticidad.
24.4 Independencia de residuos
- \(H_0 :\) Los residuos son independientes.
- \(H_1 :\) Existe autocorrelación.
## lag Autocorrelation D-W Statistic p-value
## 1 -0.003457312 2.006152 0.826
## Alternative hypothesis: rho != 0Interpretación: Dado que el p-value > 0.05, no se rechaza H0. Se considera que los residuos son independientes.
25 Multicolinealidad - Modelo limpio
La multicolinealidad se evalúa mediante el Factor de Inflación de la Varianza (VIF). * VIF < 5: Baja multicolinealidad. * VIF entre 5 y 10: Multicolinealidad moderada. * VIF > 10: Multicolinealidad problemática.
## Cb Gt V Es
## 1.031631 1.038123 1.010837 1.017959Interpretación: Todos los valores de VIF están muy cercanos a 1. Esto indica que no existe problema de multicolinealidad entre las variables Cb, Gt, V y Es en el modelo limpio.
26 Observaciones influyentes - Modelo limpio
En el modelo original se identificaron múltiples observaciones con alto leverage, las cuales estaban influyendo fuertemente en el ajuste del modelo. ### Gráfico de Leverage
leverage_limpio <- hatvalues(modelo_limpio)
umbral_leverage <- 2 * (length(coef(modelo_limpio)) / nrow(datos_limpios))
puntos_alto_leverage <- which(leverage_limpio > umbral_leverage)
puntos_alto_leverage## 9 36 44 50 74 97 107 121 125 129 144 174 187 203 208 218 220 223 225 231
## 9 36 44 50 74 97 107 121 125 129 144 174 187 203 208 218 220 223 225 231
## 237 239 246 248 270 276 283 289 291 296 312 343 347 354 361 364 368 374 381 384
## 237 239 246 248 270 276 283 289 291 296 312 343 347 354 361 364 368 374 381 384
## 391 393 402 418 424 429 459 476 477 489 516 525 533 540 580 595 610 624 631 657
## 391 393 402 418 424 429 459 476 477 489 516 525 533 540 580 595 610 624 631 657
## 664 673 681 688 690 700 714 727 770 773 775 777 778 779 782
## 664 673 681 688 690 700 714 727 770 773 775 777 778 779 782plot(leverage_limpio,
main = "Valores de Leverage - Modelo Limpio",
ylab = "Leverage",
xlab = "Observaciones")
abline(h = umbral_leverage,
col = "red",
lty = 2,
lwd = 2)
Interpretación: Al analizar el gráfico, se observa que tras eliminar los 23 puntos influyentes identificados por la distancia de Cook, el gráfico de leverage del modelo limpio muestra que persisten dos observaciones con leverage extremadamente alto (~0.98 y ~0.84), ubicadas aproximadamente en las posiciones 380 y 600 del conjunto de entrenamiento limpio. Esto indica que dichas observaciones tienen combinaciones de valores en las variables independientes muy inusuales respecto al patrón general de los datos, y que su influencia no fue capturada por el criterio de distancia de Cook utilizado en la eliminación anterior, ya que este combina leverage con tamaño del residuo. A diferencia del modelo original donde ambas observaciones extremas estaban concentradas en la misma zona (~580–650), en el modelo limpio aparecen en posiciones distintas, lo que sugiere que corresponden a observaciones diferentes. La presencia persistente de estos puntos de alto leverage indica que el conjunto de datos contiene registros con características inherentemente atípicas que dificultan la obtención de un modelo perfectamente equilibrado en la escala original de MP.
26.1 Distancia de Cook - Modelo Limpio
cooks_dist <- cooks.distance(modelo_limpio)
puntos_influyentes <- which(
cooks_dist > 4 / nrow(datos_limpios)
)
puntos_influyentes## 44 144 150 214 223 289 312 342 347 368 371 393 443 448 454 552 595 631 638 743
## 44 144 150 214 223 289 312 342 347 368 371 393 443 448 454 552 595 631 638 743
## 776
## 776Interpretación: Se determina que después de la limpieza, el número de observaciones que superan el umbral de $ $ es bajo. Esto indica que el modelo limpio es mucho más robusto y que las observaciones restantes tienen una influencia razonable sobre el modelo.
27 Tabla comparativa de modelos
tabla_comparativa <- data.frame(
Métrica = c(
"R² entrenamiento",
"R² ajustado entrenamiento",
"R² prueba",
"MSE prueba",
"RMSE prueba (COP)",
"MAE prueba (COP)",
"Normalidad residuos",
"Homocedasticidad",
"Independencia residuos",
"Multicolinealidad (VIF)",
"Observaciones influyentes Cook"
),
Modelo_Original = c(
"0.7832",
"0.7821",
"0.6650",
"11.949.828.085",
"109.315",
"43.177",
"Violada (p < 2.2e-16)",
"Violada (p < 0.0001)",
"Cumple (p = 0.408)",
"< 1.04 en todas",
"23 observaciones"
),
Modelo_Limpio = c(
"0.5344",
"0.5320",
"0.8426",
"30.098.482.992",
"173.489",
"53.137",
"Violada (p < 0.05)",
"Violada (p < 0.05)",
"Cumple (p > 0.05)",
"< 1.04 en todas",
"Bajo (umbral 4/n)"
)
)
tabla_comparativaInterpretación de la comparación entre modelos:
La comparación entre el modelo original y el modelo limpio revela un hallazgo estadístico importante relacionado con el sobreajuste (overfitting).
Comportamiento del R²:
A primera vista, el descenso del R² de entrenamiento de 0.7832 a 0.5344 podría interpretarse como un empeoramiento del modelo. Sin embargo, el análisis conjunto con el R² de prueba revela lo contrario: el modelo original presentaba un sobreajuste producido por las 23 observaciones influyentes, las cuales “jalaban” la recta de regresión hacia ellas durante el entrenamiento, generando un ajuste artificialmente alto sobre los datos de entrenamiento (0.7832) que no se sostenía al generalizar a datos nuevos (0.665). Al eliminar dichas observaciones, el modelo pierde ese ajuste forzado sobre los casos extremos, pero gana una capacidad de generalización real, elevando el R² de prueba de 0.665 a 0.8426, una mejora de 17.8 puntos porcentuales.
Errores de predicción:
El RMSE aumentó de 109.315 a 173.489 COP y el MAE de 43.177 COP a 53.137 COP en el modelo limpio. Aunque estos valores son más altos en términos absolutos, deben interpretarse cuidadosamente: el conjunto de prueba del modelo limpio tiene una distribución diferente al del modelo original (sin los casos extremos que generaban predicciones negativas y errores de signo contrario), por lo que la comparación directa de los errores no es completamente equivalente entre ambos modelos.
Supuestos estadísticos:
Ambos modelos presentan violaciones en los supuestos de normalidad de residuos y homocedasticidad, lo cual es una limitación inherente a la naturaleza de la variable MP en escala original. Esta situación es consistente con lo identificado desde el análisis exploratorio: la distribución fuertemente asimétrica de MP y la presencia de permisos con características extremas hacen que el modelo lineal en escala original no pueda cumplir completamente estos supuestos. El supuesto de independencia de residuos se cumple en ambos modelos, y la multicolinealidad es prácticamente inexistente en los dos casos, con valores VIF menores a 1.04 para todas las variables.
Conclusión:
El modelo limpio es estadísticamente superior al modelo original en términos de capacidad predictiva real (R² prueba = 0.8426 vs 0.6650) y estabilidad de los coeficientes, siendo la opción preferida para la interpretación y uso del modelo en el contexto del sistema de liquidación de la tasa compensatoria por caza de fauna silvestre. El fenómeno de sobreajuste identificado en el modelo original confirma la importancia del diagnóstico y tratamiento de observaciones influyentes como etapa fundamental en el proceso de modelamiento estadístico.
28 Selección de modelos con AIC
## Start: AIC=26930.61
## MP ~ Cb + Gt + V + Es + FR
##
## Df Sum of Sq RSS AIC
## - FR 1 21464193971 20122195610998 26930
## <none> 20100731417027 26931
## - V 1 128515391147 20229246808174 26936
## - Cb 1 208619940354 20309351357380 26940
## - Gt 1 365439606901 20466171023928 26949
## - Es 1 62964793205327 83065524622353 28548
##
## Step: AIC=26929.82
## MP ~ Cb + Gt + V + Es
##
## Df Sum of Sq RSS AIC
## <none> 20122195610998 26930
## + FR 1 21464193971 20100731417027 26931
## - Gt 1 356471351483 20478666962481 26948
## - Cb 1 743276572369 20865472183367 26969
## - V 1 2248202381730 22370397992728 27049
## - Es 1 63006820240624 83129015851622 2854629 Regla de selección
modelo_cb <- lm(MP ~ Cb,
data = tabla_analisis)
modelo_cb_gt <- lm(MP ~ Cb + Gt,
data = tabla_analisis)
modelo_cb_gt_v <- lm(MP ~ Cb + Gt + V,
data = tabla_analisis)
modelo_cb_gt_v_es <- lm(MP ~ Cb + Gt + V + Es,
data = tabla_analisis)
modelo_completo <- lm(MP ~ Cb + Gt + V + Es + FR,
data = tabla_analisis)
AIC(modelo_cb,
modelo_cb_gt,
modelo_cb_gt_v,
modelo_cb_gt_v_es,
modelo_completo)Interpretación: El mejor modelo planteado por la regla de selección AIC es el modelo que utiliza las variables Cb, Gt, V, Es que coincide con el modelo 2 sugerido al inicio de los modelos de regresión. La diferencia entre los valores de AIC es muy baja pero el modelo desmuestra un buen ajuste, equilibrio estadístico y menor complejidad.