Parte 2. Análisis en R con fitdistrplus

El presente documento ejecuta el modelamiento probabilístico de dos variables continuas (conductividad del suelo y consumo agroindustrial) mediante el uso de estimación por máxima verosimilitud en R.

Paso 1. Instalacion Y Carga De Paquetes

# --- PASO 1: INSTALACIÓN Y CARGA DE PAQUETES ---

# Creamos una lista con los paquetes necesarios
paquetes <- c("fitdistrplus", "carData", "ggplot2", "readr", "dplyr")

# Esta función instala solo los que te falten
nuevos_paquetes <- paquetes[!(paquetes %in% installed.packages()[,"Package"])]
if(length(nuevos_paquetes)) install.packages(nuevos_paquetes)

# Cargamos todas las librerías
lapply(paquetes, library, character.only = TRUE)
## Loading required package: MASS
## Loading required package: survival
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
## [[1]]
##  [1] "fitdistrplus" "survival"     "MASS"         "stats"        "graphics"    
##  [6] "grDevices"    "utils"        "datasets"     "methods"      "base"        
## 
## [[2]]
##  [1] "carData"      "fitdistrplus" "survival"     "MASS"         "stats"       
##  [6] "graphics"     "grDevices"    "utils"        "datasets"     "methods"     
## [11] "base"        
## 
## [[3]]
##  [1] "ggplot2"      "carData"      "fitdistrplus" "survival"     "MASS"        
##  [6] "stats"        "graphics"     "grDevices"    "utils"        "datasets"    
## [11] "methods"      "base"        
## 
## [[4]]
##  [1] "readr"        "ggplot2"      "carData"      "fitdistrplus" "survival"    
##  [6] "MASS"         "stats"        "graphics"     "grDevices"    "utils"       
## [11] "datasets"     "methods"      "base"        
## 
## [[5]]
##  [1] "dplyr"        "readr"        "ggplot2"      "carData"      "fitdistrplus"
##  [6] "survival"     "MASS"         "stats"        "graphics"     "grDevices"   
## [11] "utils"        "datasets"     "methods"      "base"

Paso 2 y 3. Carga de Datasets y Selección de Variables

Importamos la conductividad del suelo (Soils$Conduc) y el consumo agroindustrial (acero.csv). Para la variable operativa, aplicamos un parámetro de lectura que interpreta la coma como separador decimal, garantizando que R la reconozca como una variable numérica continua y no genere errores de lectura.

# --- Variable 1: Ingeniería Agrícola (Conductividad) ---
data("Soils", package = "carData")
variable_soils <- na.omit(Soils$Conduc)

cat("--- Resumen Estadístico: Soils$Conduc ---\n")
## --- Resumen Estadístico: Soils$Conduc ---
summary(variable_soils)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   0.670   2.790   6.635   6.589   9.852  13.320
# --- Variable 2: Ingeniería Agroindustrial (Consumo) ---
acero <- read_csv("acero.csv", locale = locale(decimal_mark = ","), show_col_types = FALSE)
variable_acero <- na.omit(acero$consumo)

cat("\n--- Resumen Estadístico: acero$consumo ---\n")
## 
## --- Resumen Estadístico: acero$consumo ---
summary(variable_acero)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   17.50   99.09  140.07  139.46  182.48  290.72

Interpretación Técnica Preliminar:

Al observar el resumen estadístico, si la media es numéricamente superior a la mediana, tenemos el primer indicio de una asimetría hacia la derecha. Esto significa que valores atípicos altos (picos de salinidad en el suelo o picos de sobreconsumo en la maquinaria) están desplazando el promedio.

Paso 4. Explorar la variable gráficamente

Generamos histogramas empíricos estandarizados a densidades de probabilidad y superponemos una estimación de densidad suavizada.

par(mfrow=c(1,2)) 

# Gráfico Soils
hist(variable_soils,
     probability = TRUE,
     main = "Histograma de Conduc - Soils",
     xlab = "Conductividad",
     col = "lightgreen")
lines(density(variable_soils), lwd = 2, col = "darkgreen")

# Gráfico Acero
hist(variable_acero,
     probability = TRUE,
     main = "Histograma de consumo - acero",
     xlab = "Consumo",
     col = "lightblue")
lines(density(variable_acero), lwd = 2, col = "darkblue")

Paso 5. Analizar la forma de los datos (Cullen y Frey)

El gráfico de Cullen y Frey es una herramienta de diagnóstico que evalúa la asimetría y el peso de las colas (curtosis) frente a distribuciones teóricas.

par(mfrow=c(1,2))

cat("Diagnóstico Soils:\n")
## Diagnóstico Soils:
descdist(variable_soils, boot = 100)
## summary statistics
## ------
## min:  0.67   max:  13.32 
## median:  6.635 
## mean:  6.588542 
## estimated sd:  3.987459 
## estimated skewness:  -0.01279401 
## estimated kurtosis:  1.552185
cat("\nDiagnóstico Acero:\n")
## 
## Diagnóstico Acero:
descdist(variable_acero, boot = 100)

## summary statistics
## ------
## min:  17.5   max:  290.72 
## median:  140.07 
## mean:  139.4565 
## estimated sd:  55.18525 
## estimated skewness:  0.003243626 
## estimated kurtosis:  2.654387

Interpretación de la forma de los datos:

La observación empírica (representada por el punto azul) demuestra que ambas variables no son simétricas, confirmando un sesgo hacia la derecha. Se ubican en zonas donde resultan ser más compatibles con distribuciones diseñadas para datos positivos y asimétricos, como lognormal, gamma o Weibull.

Paso 6. Ajustar distribuciones candidatas

Calculamos los parámetros óptimos de cuatro distribuciones teóricas mediante el algoritmo matemático de fitdist.

# Ajustes para Soils
fit_norm_soils <- fitdist(variable_soils, "norm")
fit_lnorm_soils <- fitdist(variable_soils, "lnorm")
fit_gamma_soils <- fitdist(variable_soils, "gamma")
fit_weibull_soils <- fitdist(variable_soils, "weibull")

# Ajustes para Acero
fit_norm_acero <- fitdist(variable_acero, "norm")
fit_lnorm_acero <- fitdist(variable_acero, "lnorm")
fit_gamma_acero <- fitdist(variable_acero, "gamma")
fit_weibull_acero <- fitdist(variable_acero, "weibull")

Paso 7. Comparar gráficamente los ajustes

Evaluamos el rendimiento de los modelos matemáticos utilizando cuatro perspectivas visuales para identificar cuál engloba mejor nuestra realidad física experimental.

leyenda <- c("Normal", "Lognormal", "Gamma", "Weibull")

# Comparación gráfica para Soils
par(mfrow=c(2,2))
lista_soils <- list(fit_norm_soils, fit_lnorm_soils, fit_gamma_soils, fit_weibull_soils)
denscomp(lista_soils, legendtext = leyenda, main = "Densidad (Soils)")
cdfcomp(lista_soils, legendtext = leyenda, main = "Acumulada FDA (Soils)")
qqcomp(lista_soils, legendtext = leyenda, main = "Gráfico Q-Q (Soils)")
ppcomp(lista_soils, legendtext = leyenda, main = "Gráfico P-P (Soils)")

# Comparación gráfica para Acero
par(mfrow=c(2,2))
lista_acero <- list(fit_norm_acero, fit_lnorm_acero, fit_gamma_acero, fit_weibull_acero)
denscomp(lista_acero, legendtext = leyenda, main = "Densidad (Acero)")
cdfcomp(lista_acero, legendtext = leyenda, main = "Acumulada FDA (Acero)")
qqcomp(lista_acero, legendtext = leyenda, main = "Gráfico Q-Q (Acero)")
ppcomp(lista_acero, legendtext = leyenda, main = "Gráfico P-P (Acero)")

Interpretación de los gráficos de ajuste:

En ingeniería, el gráfico Q-Q es fundamental, ya que penaliza las diferencias en los extremos (las colas). La línea teórica que mejor siga a los puntos empíricos en la parte superior derecha será la que mejor predice los valores críticos de salinidad extrema o de un alto consumo anómalo.

Paso 8. Comparar con criterios numéricos

Aplicamos la función gofstat() para obtener los estadísticos definitivos que sustentarán nuestra decisión, eliminando cualquier sesgo visual.

cat("--- Criterios de Información: Soils$Conduc ---\n")
## --- Criterios de Información: Soils$Conduc ---
gofstat(lista_soils, fitnames = leyenda)
## Goodness-of-fit statistics
##                                 Normal Lognormal     Gamma   Weibull
## Kolmogorov-Smirnov statistic 0.1311265 0.1753742 0.1621837 0.1530099
## Cramer-von Mises statistic   0.2035863 0.3274557 0.2500609 0.2358162
## Anderson-Darling statistic   1.2551940 1.9602320 1.4823249 1.4455864
## 
## Goodness-of-fit criteria
##                                  Normal Lognormal    Gamma  Weibull
## Akaike's Information Criterion 271.9903  279.1255 270.7293 267.6851
## Bayesian Information Criterion 275.7327  282.8679 274.4717 271.4275
cat("\n--- Criterios de Información: acero$consumo ---\n")
## 
## --- Criterios de Información: acero$consumo ---
gofstat(lista_acero, fitnames = leyenda)
## Goodness-of-fit statistics
##                                  Normal Lognormal     Gamma    Weibull
## Kolmogorov-Smirnov statistic 0.04222202 0.1358363 0.1043149 0.05571241
## Cramer-von Mises statistic   0.04000662 0.4946261 0.2556162 0.06139650
## Anderson-Darling statistic   0.25908130 2.9849698 1.5227200 0.39782876
## 
## Goodness-of-fit criteria
##                                  Normal Lognormal    Gamma  Weibull
## Akaike's Information Criterion 1273.530  1307.480 1286.900 1272.961
## Bayesian Information Criterion 1279.054  1313.004 1292.424 1278.486

Conclusión del Ajuste Numérico:

Para seleccionar el modelo ganador y responder a las preguntas finales de la tarea, nos basamos en los criterios de información. La regla estricta establece que el modelo con el menor valor de AIC y BIC es el definitivo, ya que representa el balance perfecto entre precisión matemática y simplicidad paramétrica.