1 Introducción

La geoestadística surge como respuesta a un tipo particular de problema: cuando se tiene una variable que se mide en un conjunto conocido de sitios dentro de una región, pero el interés real está en conocer su comportamiento en todo el dominio espacial, incluyendo los puntos donde no hubo medición directa, para el caso del presente estudio, en un cultivo de aguacate se instalaron sensores microclimáticos sobre 535 árboles georreferenciados, y aunque se dispone de la temperatura puntual en cada árbol, lo que realmente interesa para fines agronómicos es entender cómo varía la temperatura de forma continua a lo largo de toda la parcela.

Tiene sentido aplicar geoestadística aquí porque la variable de interés (temperatura) no es producto del azar puro en el espacio, depende de fenómenos físicos continuos (radiación solar, sombra proyectada por los árboles vecinos, ventilación, topografía local) que generan un patrón de variación gradual. Si dos sensores están a poca distancia entre sí, es razonable esperar que sus lecturas sean parecidas; pero si están en extremos opuestos de la parcela, esa similitud se pierde. Esta dependencia entre la cercanía espacial y la similitud de los valores (autocorrelación espacial) es precisamente lo que la geoestadística modela de forma explícita a través del variograma, y es lo que permite hacer algo que un promedio simple o una interpolación lineal no puede hacer con la misma claridad, predecir el valor de la variable en cualquier punto no muestreado junto con una medida honesta de la incertidumbre asociada a esa predicción.

Lo particular de esta variable y en general de cualquier variable regionalizada apta para geoestadística es que cumple tres condiciones:

  • Está definida sobre un dominio espacial fijo y continuo (la parcela del cultivo).
  • Su valor en un punto puede tratarse como la realización de una variable aleatoria correlacionada con sus vecinos.
  • Esa correlación decae con la distancia de una manera que se puede estimar empíricamente a partir de los datos disponibles.

Con los resultados de este análisis se espera, en primer lugar, confirmar estadísticamente que existe esa autocorrelación espacial, en segundo lugar, caracterizarla mediante un modelo teórico de semivariograma (identificando pepita, meseta y rango); y en tercer lugar, construir un mapa de predicción de temperatura para toda la parcela junto con su mapa de error, validando la calidad de esa predicción mediante validación cruzada.

2 Marco teórico

2.1 Variable regionalizada

Una variable regionalizada es un proceso estocástico \(\{Z(x): x \in D \subset \mathbb{R}^m\}\) definido sobre un dominio fijo y continuo \(D\). En nuestro caso \(m=2\) (longitud y latitud) y \(Z(x)\) representa la temperatura medida en la ubicación \(x\). A diferencia de una variable aleatoria clásica, aquí importa tanto el valor como la posición en la que se observa.

2.2 Estacionariedad

Para que el variograma sea estimable a partir de una sola realización de los datos, se asume estacionariedad de segundo orden, la dependencia entre \(Z(x)\) y \(Z(x+h)\) depende únicamente de la distancia (y eventualmente la dirección) \(h\) entre los sitios, no de su ubicación absoluta.

2.3 El semivariograma

El semivariograma cuantifica cómo cambia la disimilitud entre observaciones a medida que aumenta la distancia entre ellas:

\[2\gamma(h) = V[Z(x) - Z(x+h)] = E[Z(x)-Z(x+h)]^2\]

Su estimador empírico (estimador de momentos o de Matheron) a partir de una muestra de \(n(h)\) pares de sitios separados aproximadamente por la distancia \(h\) es:

\[2\hat{\gamma}(h) = \frac{\sum (z(x+h) - z(x))^2}{n(h)}\]

De la forma del semivariograma se leen tres elementos clave: la pepita (\(\tau^2\) o nugget), que es la discontinuidad en el origen asociada a variabilidad de microescala o error de medición; la meseta (sill), el valor al cual se estabiliza la semivarianza una vez se pierde la correlación espacial; y el rango, la distancia a partir de la cual los datos dejan de estar correlacionados.

2.4 Modelos teóricos de semivariograma

El semivariograma empírico es una estimación ruidosa basada en pares de datos; para usarlo en kriging es necesario ajustarle un modelo teórico válido (positivo-definido). Los tres modelos comparados en este trabajo son:

  • Exponencial: \(\gamma(h) = c_0 + c_1\left(1 - \exp\left(-\dfrac{h}{a}\right)\right)\)
  • Gaussiano: \(\gamma(h) = c_0 + c_1\left(1 - \exp\left(-\dfrac{h^2}{a^2}\right)\right)\)
  • Esférico: \(\gamma(h) = c_0 + c_1\left(\dfrac{3h}{2a} - \dfrac{h^3}{2a^3}\right)\) para \(h \le a\), y \(c_0+c_1\) en caso contrario

donde \(c_0\) es la pepita, \(c_1\) la meseta parcial y \(a\) el parámetro de rango.

2.5 Kriging ordinario

El kriging ordinario (KO) es un predictor lineal insesgado de varianza mínima (BLUP):

\[Z^*(x_0) = \sum_{i=1}^n \lambda_i Z(x_i), \qquad \sum_{i=1}^n \lambda_i = 1\]

Los ponderadores \(\lambda_i\) se obtienen minimizando la varianza del error de predicción \(V[Z^*(x_0)-Z(x_0)]\) sujeta a la restricción de insesgamiento, lo cual se resuelve mediante un sistema de ecuaciones lineales con un multiplicador de Lagrange \(\mu\):

\[\sum_{j=1}^n \lambda_j C_{ij} - \mu = C_{i0}, \quad i=1,\dots,n \qquad \text{sujeto a} \quad \sum_{i=1}^n \lambda_i = 1\]

A diferencia de un interpolador determinístico como el inverso de la distancia (IDW), donde los pesos dependen únicamente de la distancia geométrica, en KO los pesos \(\lambda_i\) dependen de la estructura de autocorrelación estimada en el semivariograma, lo que permite además obtener la varianza del error de predicción en cada punto del dominio.

3 Descripción de los datos

library(readxl)
library(dplyr)
library(geoR)
library(leaflet)
library(ggplot2)

El conjunto de datos completo contiene mediciones repetidas en el tiempo (entre 10 y 49 visitas) para 535 árboles de aguacate georreferenciados, con variables microclimáticas (temperatura, humedad relativa, viento, presión, etc.) y variables agronómicas (estado fenológico predominante y número de frutos afectados). Para este análisis se trabaja únicamente con el corte transversal correspondiente al 01 de octubre de 2020, donde cada árbol cuenta con exactamente una medición, lo cual produce un conjunto de datos puramente espacial, ideal para el flujo de trabajo geoestadístico clásico. La variable de interés es “Temperature”, interpolada en función de “Longitude”” y “Latitude”.

# Ajustar la ruta según corresponda en su equipo
datos_completos <- read_excel("Datos_Completos_Aguacate.xlsx")

datos_dia <- datos_completos %>%
  filter(startsWith(FORMATTED_DATE_TIME, "01/10/2020")) %>%
  select(id_arbol, Longitude, Latitude, Temperature)

glimpse(datos_dia)
## Rows: 534
## Columns: 4
## $ id_arbol    <chr> "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "…
## $ Longitude   <dbl> -76.71124, -76.71120, -76.71113, -76.71119, -76.71121, -76…
## $ Latitude    <dbl> 2.393549, 2.393573, 2.393541, 2.393503, 2.393486, 2.393441…
## $ Temperature <dbl> 23.9, 23.5, 24.5, 25.9, 26.0, 24.5, 25.5, 25.7, 26.0, 25.9…
summary(datos_dia$Temperature)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   22.20   24.50   25.80   25.83   27.18   29.70

4 Análisis exploratorio espacial

Antes de construir cualquier modelo conviene visualizar la distribución espacial de los puntos y de la variable, tal como se hizo en el ejemplo de referencia revisado previamente.

leaflet() %>%
  addTiles() %>%
  addCircleMarkers(lng = datos_dia$Longitude, lat = datos_dia$Latitude,
                    radius = 3, color = "black", fillOpacity = 0.6)
ggplot(datos_dia, aes(x = Longitude, y = Latitude, color = Temperature)) +
  geom_point(size = 2) +
  scale_color_gradient2(low = "blue", mid = "yellow", high = "red",
                         midpoint = mean(datos_dia$Temperature)) +
  coord_fixed() +
  labs(title = "Distribución espacial de la temperatura - 01/10/2020",
       color = "Temp. (°C)") +
  theme_minimal()

hist(datos_dia$Temperature, breaks = 20, col = "steelblue",
     main = "Distribución de la temperatura", xlab = "Temperatura (°C)")

La parcela tiene una extensión muy pequeña (del orden de 150-170 metros de lado), por lo que se trabaja con las coordenadas en grados decimales sin proyectar a un sistema métrico (UTM); para un área de este tamaño la distorsión introducida es prácticamente despreciable.

5 Construcción de la variable regionalizada

geo_temp <- as.geodata(datos_dia,
                        coords.col = match(c("Longitude","Latitude"), names(datos_dia)),
                        data.col   = match("Temperature", names(datos_dia)))
plot(geo_temp)

distancias <- dist(datos_dia[, c("Longitude","Latitude")])
summary(distancias)
##      Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
## 1.712e-05 4.051e-04 6.408e-04 6.827e-04 9.178e-04 1.959e-03

6 Variograma empírico y prueba de significancia

Se calculó el variograma omnidireccional con bins hasta el tercer cuartil de la distribución de distancias (criterio habitual para evitar estimaciones poco confiables en las colas, donde hay pocos pares de puntos), y se acompaña de una prueba de Monte Carlo por permutación: si se reordena aleatoriamente la temperatura entre las ubicaciones (rompiendo cualquier estructura espacial real) y se repite el cálculo del variograma cientos de veces, se obtiene una envolvente de lo que cabría esperar si no hubiera autocorrelación espacial. Si el variograma empírico observado se sale de esa envolvente, es evidencia de autocorrelación espacial genuina.

cutoff <- as.numeric(quantile(distancias, 0.75))

variograma <- variog(geo_temp, option = "bin", uvec = seq(0, cutoff, length.out = 10))
## variog: computing omnidirectional variogram
plot(variograma, main = "Semivariograma empírico")

set.seed(123)
variograma_mc <- variog.mc.env(geo_temp, obj = variograma)
## variog.env: generating 99 simulations by permutating data values
## variog.env: computing the empirical variogram for the 99 simulations
## variog.env: computing the envelops
plot(variograma, envelope.obj = variograma_mc,
     main = "Semivariograma empírico con envolvente Monte Carlo (95%)")

En el ejercicio realizado para verificar esta metodología se encontró que, a distancias cortas, la semivarianza empírica cae claramente por debajo de la envolvente de permutación, lo cual es evidencia estadística de autocorrelación espacial fuerte: los árboles cercanos presentan temperaturas mucho más similares entre sí de lo que cabría esperar por azar. A distancias mayores el variograma empírico se acerca o incluso supera la envolvente, lo cual es consistente con la pérdida progresiva de correlación espacial a medida que aumenta la distancia.

7 Ajuste de modelos teóricos de semivariograma

Para una comparación justa entre los tres modelos, se estima la pepita libremente en los tres casos (a diferencia del ejemplo de referencia, donde cada modelo se ajustó con una restricción distinta sobre la pepita, lo que dificulta comparar directamente la suma de cuadrados entre ellos). Se usa el mismo criterio de ponderación por número de pares (wei = "npairs") en los tres ajustes.

ini.vals <- expand.grid(seq(0.5, 3, length.out = 8), seq(0.0001, 0.001, length.out = 8))

modelo_exp <- variofit(variograma, ini.cov.pars = ini.vals, cov.model = "exponential",
                        nugget = 1, fix.nugget = FALSE, wei = "npairs", min = "optim")
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2.29"  "0"   "1"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 2728.87463825682
modelo_gauss <- variofit(variograma, ini.cov.pars = ini.vals, cov.model = "gaussian",
                          nugget = 1, fix.nugget = FALSE, wei = "npairs", min = "optim")
## variofit: covariance model used is gaussian 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2.29"  "0"   "1"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 7644.15875110439
modelo_sph <- variofit(variograma, ini.cov.pars = ini.vals, cov.model = "spherical",
                        nugget = 1, fix.nugget = FALSE, wei = "npairs", min = "optim")
## variofit: covariance model used is spherical 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim 
## variofit: searching for best initial value ... selected values:
##               sigmasq phi   tausq kappa
## initial.value "2.29"  "0"   "1"   "0.5"
## status        "est"   "est" "est" "fix"
## loss value: 5700.33848459119
comparacion <- data.frame(
  modelo = c("Exponencial", "Gaussiano", "Esférico"),
  pepita = c(modelo_exp$nugget, modelo_gauss$nugget, modelo_sph$nugget),
  meseta_parcial = c(modelo_exp$cov.pars[1], modelo_gauss$cov.pars[1], modelo_sph$cov.pars[1]),
  rango = c(modelo_exp$cov.pars[2], modelo_gauss$cov.pars[2], modelo_sph$cov.pars[2]),
  SCE = c(modelo_exp$value, modelo_gauss$value, modelo_sph$value)
)
comparacion
plot(variograma, main = "Comparación de modelos teóricos ajustados")
lines(modelo_exp, col = "blue", lwd = 2)
lines(modelo_gauss, col = "red", lwd = 2)
lines(modelo_sph, col = "purple", lwd = 2)
legend("bottomright", legend = c("Exponencial","Gaussiano","Esférico"),
       col = c("blue","red","purple"), lwd = 2)

modelos <- list(exponential = modelo_exp, gaussian = modelo_gauss, spherical = modelo_sph)
sce_valores <- sapply(modelos, function(m) m$value)
mejor_modelo_nombre <- names(which.min(sce_valores))
mejor_modelo <- modelos[[mejor_modelo_nombre]]
cat("Modelo seleccionado:", mejor_modelo_nombre, "\n")
## Modelo seleccionado: exponential
summary(mejor_modelo)
## $pmethod
## [1] "WLS (weighted least squares)"
## 
## $cov.model
## [1] "exponential"
## 
## $spatial.component
##      sigmasq          phi 
## 2.2206792133 0.0002065834 
## 
## $spatial.component.extra
## kappa 
##   0.5 
## 
## $nugget.component
##    tausq 
## 1.088749 
## 
## $fix.nugget
## [1] FALSE
## 
## $fix.kappa
## [1] TRUE
## 
## $practicalRange
## [1] 0.0006037078
## 
## $sum.of.squares
##    value 
## 1986.628 
## 
## $estimated.pars
##        tausq      sigmasq          phi 
## 1.0887487077 2.2206792133 0.0002065834 
## 
## $weights
## [1] "npairs"
## 
## $call
## variofit(vario = variograma, ini.cov.pars = ini.vals, cov.model = "exponential", 
##     fix.nugget = FALSE, nugget = 1, weights = "npairs", minimisation.function = "optim")
## 
## attr(,"class")
## [1] "summary.variomodel"

Al verificar este procedimiento se encontró que el modelo exponencial obtiene la menor suma de cuadrados ponderada, lo cual es consistente con la forma del variograma empírico: un crecimiento marcado a distancias muy cortas seguido de una aproximación más gradual a la meseta, un comportamiento que el modelo exponencial con su curvatura cóncava característica cerca del origen captura mejor que el gaussiano (más suave en el origen) o el esférico (lineal cerca del origen).

8 Predicción espacial: Kriging Ordinario

rango_lon <- range(datos_dia$Longitude)
rango_lat <- range(datos_dia$Latitude)

grilla <- expand.grid(
  Longitude = seq(rango_lon[1], rango_lon[2], length.out = 100),
  Latitude  = seq(rango_lat[1], rango_lat[2], length.out = 100)
)
ko_pred <- krige.conv(geo_temp, loc = grilla,
                       krige = krige.control(obj.model = mejor_modelo,
                                              trend.d = "cte", trend.l = "cte"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
par(mfrow = c(1,2))
image(ko_pred, loc = grilla, main = "Predicción Kriging - Temperatura",
      xlab = "Longitud", ylab = "Latitud")
contour(ko_pred, loc = grilla, add = TRUE, drawlabels = TRUE)

image(ko_pred, loc = grilla, val = sqrt(ko_pred$krige.var),
      main = "Desviación estándar del Kriging", xlab = "Longitud", ylab = "Latitud")
contour(ko_pred, loc = grilla, val = sqrt(ko_pred$krige.var), add = TRUE, drawlabels = TRUE)

El mapa de predicción revela una estructura espacial nada trivial: una franja diagonal de temperaturas más bajas (entre 23.5 y 24.5 °C) atraviesa la parcela desde la zona centro-superior hacia el sector inferior derecho, separando dos focos claramente más cálidos (por encima de 27.5-28.5 °C), uno en la esquina inferior izquierda y otro hacia el centro-derecha de la parcela. Que esta franja fría tenga una orientación diagonal y no se distribuya de forma aleatoria es información agronómicamente valiosa: es poco probable que un patrón tan organizado se deba al azar de la instrumentación, y resulta más plausible que esté asociado a algún factor estructural del terreno o del cultivo que sigue esa misma orientación, por ejemplo un surco o hilera de siembra con mayor densidad de copa que genera sombra continua, un canal de drenaje o escorrentía que favorece la acumulación de aire frío (cold air pooling), o un corredor de viento dominante alineado en esa dirección que favorece la pérdida de calor por convección. Los dos focos cálidos, en cambio, podrían corresponder a bordes de la parcela con menor cobertura de copa (mayor exposición a radiación solar directa) o a zonas con suelo más expuesto. Esta lectura es, por supuesto, una hipótesis a contrastar con información adicional del cultivo (mapa de alturas de árbol, orientación de las hileras de siembra, topografía), pero es precisamente el tipo de pregunta que un mapa de kriging permite formular y que un promedio simple de temperatura jamás revelaría.

En cuanto a la incertidumbre, el mapa de desviación estándar del kriging se comporta exactamente como predice la teoría: es mínima en el núcleo densamente muestreado de la parcela (donde la interpolación se apoya en muchos vecinos cercanos) y crece de forma sostenida hacia los bordes, llegando a casi 1.9 °C en las esquinas más alejadas de cualquier árbol con sensor, frente a valores cercanos a 1.1-1.2 °C en el centro. Esto no es un detalle menor: significa que las predicciones en el borde de la parcela deben tomarse con considerablemente menos confianza que las del centro, algo que un mapa de predicción sin su correspondiente mapa de error simplemente ocultaría.

9 Validación cruzada

valida <- xvalid(geodata = geo_temp, model = mejor_modelo)
## xvalid: number of data locations       = 534
## xvalid: number of validation locations = 534
## xvalid: performing cross-validation at location ... 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 60, 61, 62, 63, 64, 65, 66, 67, 68, 69, 70, 71, 72, 73, 74, 75, 76, 77, 78, 79, 80, 81, 82, 83, 84, 85, 86, 87, 88, 89, 90, 91, 92, 93, 94, 95, 96, 97, 98, 99, 100, 101, 102, 103, 104, 105, 106, 107, 108, 109, 110, 111, 112, 113, 114, 115, 116, 117, 118, 119, 120, 121, 122, 123, 124, 125, 126, 127, 128, 129, 130, 131, 132, 133, 134, 135, 136, 137, 138, 139, 140, 141, 142, 143, 144, 145, 146, 147, 148, 149, 150, 151, 152, 153, 154, 155, 156, 157, 158, 159, 160, 161, 162, 163, 164, 165, 166, 167, 168, 169, 170, 171, 172, 173, 174, 175, 176, 177, 178, 179, 180, 181, 182, 183, 184, 185, 186, 187, 188, 189, 190, 191, 192, 193, 194, 195, 196, 197, 198, 199, 200, 201, 202, 203, 204, 205, 206, 207, 208, 209, 210, 211, 212, 213, 214, 215, 216, 217, 218, 219, 220, 221, 222, 223, 224, 225, 226, 227, 228, 229, 230, 231, 232, 233, 234, 235, 236, 237, 238, 239, 240, 241, 242, 243, 244, 245, 246, 247, 248, 249, 250, 251, 252, 253, 254, 255, 256, 257, 258, 259, 260, 261, 262, 263, 264, 265, 266, 267, 268, 269, 270, 271, 272, 273, 274, 275, 276, 277, 278, 279, 280, 281, 282, 283, 284, 285, 286, 287, 288, 289, 290, 291, 292, 293, 294, 295, 296, 297, 298, 299, 300, 301, 302, 303, 304, 305, 306, 307, 308, 309, 310, 311, 312, 313, 314, 315, 316, 317, 318, 319, 320, 321, 322, 323, 324, 325, 326, 327, 328, 329, 330, 331, 332, 333, 334, 335, 336, 337, 338, 339, 340, 341, 342, 343, 344, 345, 346, 347, 348, 349, 350, 351, 352, 353, 354, 355, 356, 357, 358, 359, 360, 361, 362, 363, 364, 365, 366, 367, 368, 369, 370, 371, 372, 373, 374, 375, 376, 377, 378, 379, 380, 381, 382, 383, 384, 385, 386, 387, 388, 389, 390, 391, 392, 393, 394, 395, 396, 397, 398, 399, 400, 401, 402, 403, 404, 405, 406, 407, 408, 409, 410, 411, 412, 413, 414, 415, 416, 417, 418, 419, 420, 421, 422, 423, 424, 425, 426, 427, 428, 429, 430, 431, 432, 433, 434, 435, 436, 437, 438, 439, 440, 441, 442, 443, 444, 445, 446, 447, 448, 449, 450, 451, 452, 453, 454, 455, 456, 457, 458, 459, 460, 461, 462, 463, 464, 465, 466, 467, 468, 469, 470, 471, 472, 473, 474, 475, 476, 477, 478, 479, 480, 481, 482, 483, 484, 485, 486, 487, 488, 489, 490, 491, 492, 493, 494, 495, 496, 497, 498, 499, 500, 501, 502, 503, 504, 505, 506, 507, 508, 509, 510, 511, 512, 513, 514, 515, 516, 517, 518, 519, 520, 521, 522, 523, 524, 525, 526, 527, 528, 529, 530, 531, 532, 533, 534, 
## xvalid: end of cross-validation
MAE  <- mean(abs(valida$error))
RMSE <- sqrt(mean(valida$error^2))
sesgo <- mean(valida$error)

data.frame(MAE = MAE, RMSE = RMSE, Sesgo = sesgo)
plot(valida$data, valida$predicted, xlab = "Valor observado", ylab = "Valor predicho (LOO)",
     main = "Validación cruzada: observado vs. predicho")
abline(0, 1, col = "red", lwd = 2)

razon_pepita_meseta <- mejor_modelo$nugget / (mejor_modelo$nugget + mejor_modelo$cov.pars[1])
varianza_global <- var(datos_dia$Temperature)
pseudo_R2 <- 1 - (RMSE^2 / varianza_global)

data.frame(
  Razon_pepita_meseta = round(razon_pepita_meseta, 3),
  Pseudo_R2 = round(pseudo_R2, 3)
)

El error absoluto medio (MAE ≈ 0.88 °C) y el sesgo prácticamente nulo (≈ -0.004 °C) son dos lecturas complementarias: la magnitud del error indica qué tan lejos está, en promedio, la predicción del valor real, mientras que el sesgo cercano a cero confirma que el kriging ordinario está cumpliendo su propiedad teórica de insesgamiento, no hay una tendencia sistemática a sobreestimar ni subestimar la temperatura. Puesto en contexto, la desviación estándar global de la temperatura en este corte es de aproximadamente 1.77 °C, de modo que un MAE de 0.88 °C equivale a cerca de la mitad de esa variabilidad total, y el pseudo-R² calculado a partir del RMSE de validación cruzada (1 − RMSE²/Var(Z)) ronda el 61%. Es decir, el modelo espacial ajustado logra explicar algo más de la mitad de la variabilidad observada en la temperatura usando exclusivamente la posición de cada árbol, sin ninguna covariable adicional, lo cual es un desempeño razonablemente bueno y respalda que la estructura espacial detectada en el semivariograma es real y aprovechable para predicción, no solo una curiosidad estadística.

La razón pepita/meseta (nugget-to-sill ratio) calculada arriba, de aproximadamente 33%, es otro indicador útil de qué tan fuerte es la dependencia espacial: siguiendo la clasificación propuesta por Cambardella et al. (1994) —ampliamente usada en estudios de variabilidad espacial de suelos y microclima—, razones por debajo de 25% se consideran dependencia espacial fuerte, entre 25% y 75% dependencia moderada, y por encima de 75% dependencia débil (la variable se comporta casi como ruido blanco espacial). El valor obtenido aquí (~33%) ubica a la temperatura de este cultivo en el rango de dependencia moderada: hay una componente de variación de microescala (pepita) no despreciable —que puede deberse a precisión del sensor, microvariaciones puntuales de sombra o ventilación que no se capturan a la resolución de muestreo usada—, pero la mayor parte de la variabilidad total (~67%, la meseta parcial) sigue estando estructurada espacialmente y es la que el kriging logra explotar.

10 Discusión y conclusiones

10.1 Sobre la autocorrelación espacial y el modelo ajustado

La prueba de Monte Carlo confirmó algo que es la base de todo el análisis posterior: la temperatura de los árboles no se distribuye al azar en la parcela. En las distancias más cortas evaluadas, la semivarianza empírica cayó claramente por debajo de la envolvente de permutación, lo que significa que pares de árboles cercanos presentan diferencias de temperatura sistemáticamente menores que las que produciría una asignación aleatoria de valores a las ubicaciones. A medida que la distancia entre árboles aumenta, esa diferencia se acerca y eventualmente se mezcla con la envolvente, justo el comportamiento que se espera de un proceso espacialmente autocorrelacionado: la similitud entre vecinos se diluye con la distancia hasta volverse estadísticamente indistinguible del azar.

El modelo exponencial resultó ser el de mejor ajuste entre los tres evaluados, con una suma de cuadrados ponderada (≈1987) notablemente menor que la del esférico (≈4887) y la del gaussiano (≈5424). Esto no es un detalle puramente técnico: la forma del modelo ajustado tiene una interpretación física. Un modelo exponencial, con su ascenso pronunciado cerca del origen, describe procesos en los que la correlación se pierde relativamente rápido al alejarse de un punto, consistente con un fenómeno de microclima dominado por factores muy locales (sombra puntual de un árbol, un parche de suelo más húmedo) más que por gradientes regionales suaves. El rango práctico estimado (~0.0006 grados, del orden de 60-70 metros si se aproxima 1 grado a 111 km) indica que, más allá de esa distancia, dos árboles cualesquiera de la parcela aportan información prácticamente independiente sobre la temperatura del otro. Ese número tiene una lectura práctica directa: si en el futuro se quisiera diseñar una red de sensores más económica que la actual (535 árboles instrumentados), espaciarlos a más de 60-70 metros implicaría perder la posibilidad de interpolar con la misma confianza entre ellos.

10.2 Sobre el mapa de predicción y su lectura agronómica

El valor del kriging frente a alternativas más simples se aprecia mejor en el propio mapa de predicción: la franja diagonal de menor temperatura y los dos focos más cálidos identificados no serían visibles en un simple promedio o en una tabla de valores por árbol. Esa estructura organizada -y no un mosaico de manchas aleatorias- es evidencia indirecta adicional de que el proceso subyacente tiene una causa física identificable (orientación de hileras de siembra, topografía o un patrón de sombra), y abre una pregunta concreta y útil para el equipo agronómico del cultivo: ¿qué hay físicamente en esa diagonal que la distingue del resto de la parcela? Esa es, en buena medida, la utilidad real de la geoestadística en un contexto productivo, no solo interpolar, sino señalar patrones espaciales que orientan decisiones de manejo diferenciado por zonas (por ejemplo, riego o manejo fitosanitario ajustado a las condiciones de cada sector del lote).

10.3 Conclusiones generales

Retomando las preguntas que enmarcaron este trabajo y respondiéndolas ya con los resultados obtenidos: la geoestadística tuvo sentido aplicarla aquí porque el objetivo nunca fue resumir la temperatura de la parcela con un único número, sino predecirla en cualquier punto del cultivo, y la prueba de Monte Carlo demostró que existe una autocorrelación espacial real y estadísticamente significativa que un análisis no espacial simplemente pasaría por alto. La temperatura, como variable regionalizada, mostró exactamente las propiedades que la hacen apta para este enfoque: continuidad sobre un dominio fijo, una dependencia espacial que decae con la distancia de forma consistente con un modelo exponencial, y una razón pepita/meseta (~33%) que confirma una dependencia espacial de fuerza moderada, ni tan débil como para que el muestreo sea inútil, ni tan fuerte como para que un solo sensor central bastara para toda la parcela.

En términos de resultados concretos, el análisis entregó tres productos con valor práctico directo: un semivariograma caracterizado (pepita ≈1.09, meseta parcial ≈2.22, rango práctico ≈60-70 m) que cuantifica la escala espacial del fenómeno; un mapa de predicción que revela una estructura térmica organizada dentro de la parcela con valor agronómico potencial; y una validación cruzada (MAE≈0.88 °C, pseudo-R²≈61%) que respalda que el modelo ajustado tiene poder predictivo real y no es un sobreajuste a los datos observados. En conjunto, estos resultados confirman que, incluso en una parcela de extensión reducida, la variabilidad microclimática no es despreciable ni aleatoria, y que la geoestadística ofrece tanto las herramientas de predicción como las de caracterización necesarias para describirla rigurosamente.

11 Referencias

  • Cambardella, C. A., Moorman, T. B., Novak, J. M., Parkin, T. B., Karlen, D. L., Turco, R. F., & Konopka, A. E. (1994). Field-scale variability of soil properties in central Iowa soils. Soil Science Society of America Journal, 58(5), 1501-1511.
  • Cressie, N. (1993). Statistics for Spatial Data. Wiley.
  • Diggle, P. J., & Ribeiro, P. J. (2007). Model-based Geostatistics. Springer.
  • Isaaks, E. H., & Srivastava, R. M. (1989). An Introduction to Applied Geostatistics. Oxford University Press.
  • Material de clase: Análisis de Datos Espaciales — Tema: Geoestadística.