Planteamiento del problema

La península de Yucatán se extiende entre el Golfo de México al oeste y norte, y el Mar Caribe al este. Los sitios arqueológicos que en esta se encuentran, como Chichén Itza y Uxmal (designados como Patrimonio de la Humanidad por la UNESCO) sumado a las playas de Cancún hacen de esta península uno de los puntos más turísticos del país. Esto resalta la importancia de analizar sus condiciones climáticas, pues es una variable a tener en cuenta al momento de planear una visita al territorio.

El clima de Yucatán es uno de los más calurosos de México; esto se debe a la ubicación geográfica. Las altas temperaturas azotan el estado durante todo el año, dando a la zona un clima semihúmedo y caluroso. Aproximadamente el 85% del territorio recibe entre 24° C y 30° C durante todo el año. El estado de Yucatán puede llegar a los 42° C en verano.

Descripción de los datos

Los datos usados provienen de la información histórica de las estaciones climatológicas convencionales que conforman la Red Nacional de la Comisión Nacional de Agua (CONAGUA). Se tomó una muestra de 33 estaciones ubicadas al sureste de México en la península de Yucatán, concretamente en los estados de Yucatán, Quintana Roo y Campeche. Esta zona es una de las más famosas del país por sus playas en el golfo de México y las ruinas mayas, resaltando la importancia de analizar las condiciones climáticas de esta.

Se seleccionaron las variables Temperatura máxima diaria y Precipitación acumulada diaria desde 2002-01-01 hasta 2024-12-31, resultando en 8401 datos. Aunque se cuenta con datos históricos anteriores, se toma este periodo de tiempo para no tener tantos datos faltantes en la base de datos.

Los datos se pueden consultar en: https://smn.conagua.gob.mx/es/climatologia/informacion-climatologica/informacion-estadistica-climatologica

Las unidades de las variables son:

Para este primer análisis solamente espacial, se toman los datos recolectados del día 2018-06-16.

Carga de librerías

library(readxl)
## Warning: package 'readxl' was built under R version 4.5.2
library(geoR)
## Warning: package 'geoR' was built under R version 4.5.2
## --------------------------------------------------------------
##  Analysis of Geostatistical Data
##  For an Introduction to geoR go to http://www.leg.ufpr.br/geoR
##  geoR version 1.9-6 (built on 2025-08-29) is now loaded
## --------------------------------------------------------------
library(gstat)
## Warning: package 'gstat' was built under R version 4.5.2
library(spacetime)
## Warning: package 'spacetime' was built under R version 4.5.2
library(sp)
## Warning: package 'sp' was built under R version 4.5.2
library(sf)
## Warning: package 'sf' was built under R version 4.5.2
## Linking to GEOS 3.13.1, GDAL 3.11.4, PROJ 9.7.0; sf_use_s2() is TRUE
library(tidyr)
library(dplyr)
## 
## Adjuntando el paquete: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(tmap)
## Warning: package 'tmap' was built under R version 4.5.2
library(rnaturalearth)
## Warning: package 'rnaturalearth' was built under R version 4.5.2
library(rnaturalearthdata)
## Warning: package 'rnaturalearthdata' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'rnaturalearthdata'
## The following object is masked from 'package:rnaturalearth':
## 
##     countries110
library(rnaturalearthhires)
library(ggplot2)

setwd("C:/Users/CAMILA/OneDrive - Universidad Nacional de Colombia/Escritorio/Universidad/Maestría Estadística/2do semestre/Espacial")

Lecura de datos

base_PRECIP <- read_excel("base_PRECIP_2.xlsx", sheet = "Sheet1")
base_TMAX <- read_excel("base_TMAX_2.xlsx", sheet = "Sheet1")
ubicaciones <- read_excel("base_TMAX_2.xlsx", sheet = "Ubicaciones")
## New names:
## • `` -> `...1`
fecha <- as.POSIXct("2018-06-16", tz="UTC")
dato_TMAX <- base_TMAX[base_TMAX$FECHA==fecha,]
dato_PRECIP <- base_PRECIP[base_PRECIP$FECHA==fecha,]

df_clima <- cbind(Temp = as.numeric(dato_TMAX[-1]), Precip = as.numeric(dato_PRECIP[-1]))
df_clima <- as.data.frame(df_clima)
head(df_clima)
##   Temp Precip
## 1   30    9.0
## 2   32   28.0
## 3   34   29.0
## 4   32   36.2
## 5   30   40.0
## 6   30   40.0
par(mfrow=c(1, 2))
boxplot(df_clima$Temp, main = "Boxplot Temperatura máxima")
boxplot(df_clima$Precip, main = "Boxplot Precipitación acumulada")

Convertir a coordenadas planas

# Crear un objeto SpatialPoints
ori_coords <- SpatialPoints(cbind(ubicaciones$LONGITUD, ubicaciones$LATITUD),
                        proj4string = CRS("+proj=longlat"))

# Convertirlo a UTM
utm_coords <- spTransform(ori_coords, CRS("+proj=utm +zone=16 +datum=WGS84 +units=m +no_defs"))
plot(utm_coords, axes = TRUE, main = "Coordenadas UTM")

utm_coords_df <- as.data.frame(utm_coords)
colnames(utm_coords_df) <- c("Easting", "Northing")

Mapa de la zona de interés

# Convertir tus puntos de sp → sf
utm_coords_sf <- st_as_sf(utm_coords)

# 2. Descargar el polígono de México
mexico <- ne_countries(scale = "medium", country = "Mexico", returnclass = "sf")

# 3. Transformar México a UTM zona 16N
mexico_utm <- st_transform(mexico, 32616)

# 4. Graficar
ggplot() +
  geom_sf(data = mexico_utm, fill = "gray90", color = "black") +
  geom_sf(data = utm_coords_sf, color = "red", size = 2) +
  theme_minimal()

ggplot() +
  geom_sf(data = mexico_utm, fill = "gray90", color = "black") +
  geom_sf(data = utm_coords_sf, color = "red", size = 2) +
  coord_sf(xlim = range(st_coordinates(utm_coords_sf)[,1]) + c(-60000, 15000),
           ylim = range(st_coordinates(utm_coords_sf)[,2]) + c(-30000, 30000)) +
  theme_minimal()

# GEOESTADÍSTICA ## Análisis de tendencia

df_clima <- cbind(utm_coords_df, df_clima)
head(df_clima)
##    Easting Northing Temp Precip
## 1 212518.7  2214240   30    9.0
## 2 173715.3  2213196   32   28.0
## 3 162137.3  2210460   34   29.0
## 4 168065.1  2229804   32   36.2
## 5 172492.9  2234861   30   40.0
## 6 173684.1  2236378   30   40.0
pairs(df_clima[,-4])

pairs(df_clima[,-3])

Los graficos de dispersión dejan ver que no se cuenta con ningún tipo de tendencia en los datos ni de Temperatura ni de Precipitación. Es por esta razón que se puede asumir que las medias de los procesos no dependen de la ubicación espacial, luego son de media constante.

# Convertimos los datos a un objeto de la clase geodata
temp.g <- as.geodata(df_clima, data.col = 3)
precip.g <- as.geodata(df_clima, data.col = 4)

# Resumen para la variable Temperatura
summary(temp.g)
## Number of data points: 33 
## 
## Coordinates summary
##      Easting Northing
## min 162137.3  2210460
## max 296321.7  2285992
## 
## Distance summary
##        min        max 
##   1067.368 145794.046 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## 26.00000 28.50000 30.00000 30.93636 33.00000 38.00000
plot(temp.g)

# Resumen para la variable Precipitación 
summary(precip.g)
## Number of data points: 33 
## 
## Coordinates summary
##      Easting Northing
## min 162137.3  2210460
## max 296321.7  2285992
## 
## Distance summary
##        min        max 
##   1067.368 145794.046 
## 
## Data summary
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
##  0.00000 10.00000 25.00000 24.74545 36.20000 55.00000
plot(precip.g)

Una vez finalizado este análisis de tendencia, concluimos que dado que la media del proceso espacial es constante y que podemos asumir que la función de covarianza \(C(\mathbf{h})\) existe y es función única del vector de separación \(\mathbf{h}\) tenemos un proceso estacionario débil.

Explorar anisotropía

Un supuesto importante a revisar es el de isotropía. Este es, corroborar si la función de covarianza o semivarianza dependen únicamete de la magnitud de la separación \(\mathbf{h}\). Para esto se llevan a cabo estimaciones empíricas de cuatro semivariogramas direccionales, como se observa en las siguientes gráficas.

var_direc_temp <- variog4(as.geodata(df_clima, data.col = 3))
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
plot(var_direc_temp, main = "Semivariogramas direccionales para la Temperatura máxima diaria")

var_direc_precip <- variog4(as.geodata(df_clima, data.col = 4))
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
## variog: computing omnidirectional variogram
plot(var_direc_precip, main = "Semivariogramas direccionales para la Precipitación acumulada diaria")

Note que los semivariogramas se ven diferentes depeniendo de la dirección en la que se consideren, en este sentido los semivariogramas serían anisotrópicos.

Estimación del semivariograma

El semivariograma se define como la varianza de la variable incrementos. Este se estima empíricamente como \[\hat{\gamma}(\mathbf{h})=\frac{1}{2\mid N(\mathbf{h})\mid}\sum_{N(\mathbf{h})}[Z(\mathbf{s+h})-Z(\mathbf{s})]^2\]

vari_temp <- variog(as.geodata(df_clima, data.col = 3), max.dist = 100000)
## variog: computing omnidirectional variogram
vari_precip <- variog(as.geodata(df_clima, data.col = 4), max.dist = 100000)
## variog: computing omnidirectional variogram
plot(vari_temp)

plot(vari_precip, main = "Semivariograma empírico para la Precipitación acumulada diaria")

rob_vari_temp <- variog(as.geodata(df_clima, data.col = 3), 
                        estimator.type = "modulus", max.dist = 100000)
## variog: computing omnidirectional variogram
rob_vari_precip <- variog(as.geodata(df_clima, data.col = 4), 
                          estimator.type = "modulus", max.dist = 100000)
## variog: computing omnidirectional variogram
plot(rob_vari_temp)

plot(rob_vari_precip, main = "Semivariograma empírico robusto para la Precipitación acumulada diaria")

La función de semivarianza empírica robusta muestra un comportamiento más sinusoidal que la estimada por el método clásico, además de un menor efecto pepita. En ambos casos no se observa que la función de semivarianza tienda concretamente hacia un punto de silla \(C(\mathbf{0})=\sigma^2\), aunque si se evidencia su comportamiento creciente.

Modelos para los semivariogramas

Modelos para la variable: Temperatura

Para la variable Temperatura se prueban los modelos Matern, exponencial y powered exponencial dado que fueron los que mostraron un mejor acople en el análisis inicial a sentimiento. Para cada uno de estos se presentan las estimaciones mediante MCo y MCP

MSE <- function(x) mean(x**2)

#eyefit(vari_temp)
#    cov.model           sigmasq      phi tausq kappa kappa2   practicalRange
#1         exponential   11.14 30997.59  1.39  <NA>   <NA> 92860.4807632263
#2 powered.exponential   11.14 31997.59     0  2.57   <NA> 49037.3503979688
#3              matern   11.14 27538.39     0  0.21   <NA> 56375.4783222614

# - - - - - Modelos exponenciales - - - - - - - - 
ini1 <- c(11.14, 30997.59)
fitvar_temp_exp1 <- variofit(vari_temp, 
                            cov.model = "exponential", 
                            ini.cov.pars = ini1,
                            wei= "equal",
                            fix.nugget = TRUE,
                            nugget = 0)
## variofit: covariance model used is exponential 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
fitvar_temp_exp2 <- variofit(vari_temp, 
                            cov.model = "exponential", 
                            ini.cov.pars = ini1,
                            wei= "npairs",
                            fix.nugget = TRUE,
                            nugget = 0)
## variofit: covariance model used is exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
fitvar_temp_exp3 <- variofit(vari_temp, 
                            cov.model = "exponential", 
                            ini.cov.pars = ini1,
                            wei= "cressie",
                            fix.nugget = TRUE,
                            nugget = 0)
## variofit: covariance model used is exponential 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
fitvar_temp_exp4 <- likfit(temp.g, 
                           cov.model = "exponential", 
                           ini.cov.pars = ini1,
                           fix.nugget = TRUE, 
                           nugget =  0)
## kappa not used for the exponential correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optimize.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
## 
## WARNING: estimated range is less than 1 tenth of the minimum distance between two points. Consider re-examine the model excluding spatial dependence
mse1 <- MSE(vari_temp$v-fitvar_temp_exp1$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_exp1$cov.pars))
mse2 <- MSE(vari_temp$v-fitvar_temp_exp2$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_exp2$cov.pars))
mse3 <- MSE(vari_temp$v-fitvar_temp_exp3$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_exp3$cov.pars))
mse4 <- MSE(vari_temp$v-fitvar_temp_exp4$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_exp4$cov.pars))
mse1;mse2;mse3;mse4
## [1] 5.318534
## [1] 5.324614
## [1] 7.46398
## [1] 6.909518
plot(vari_temp, main ="Modelos exponenciales")
lines(fitvar_temp_exp1, col = "red")
lines(fitvar_temp_exp2, col = "green")
lines(fitvar_temp_exp3, col = "blue")
lines(fitvar_temp_exp4, col = "black")
legend("bottomright", legend = c("MCO", "MCP con npairs", "MCP con cressie", "MLE"), col =c("red", "green", "blue", "black"), lty = rep(1, 4))

# El modelo exponencial con menor MSE es el de MCO : 5.3185
# Modelos powered exponenciales
ini1 <- c(11.14, 31997.59)
fitvar_temp_pow1 <- variofit(vari_temp, 
                            cov.model = "powered.exponential", 
                            ini1,
                            wei= "equal",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 2.57)
## variofit: covariance model used is powered.exponential 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
fitvar_temp_pow2 <- variofit(vari_temp, 
                            cov.model = "powered.exponential", 
                            ini1,
                            wei= "npairs",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 2.57)
## variofit: covariance model used is powered.exponential 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
fitvar_temp_pow3 <- variofit(vari_temp, 
                            cov.model = "powered.exponential", 
                            ini1,
                            wei= "cressie",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 2.57)
## variofit: covariance model used is powered.exponential 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
fitvar_temp_pow4 <- likfit(temp.g, 
                           cov.model = "powered.exponential", 
                           ini.cov.pars = ini1,
                           fix.nugget = TRUE, 
                           nugget =  0)
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optimize.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
mse1 <- MSE(vari_temp$v-fitvar_temp_pow1$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_pow1$cov.pars))
mse2 <- MSE(vari_temp$v-fitvar_temp_pow2$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_pow2$cov.pars))
mse3 <- MSE(vari_temp$v-fitvar_temp_pow3$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_pow3$cov.pars))
mse4 <- MSE(vari_temp$v-fitvar_temp_pow4$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_pow4$cov.pars))
mse1;mse2;mse3;mse4
## [1] 6.723288
## [1] 6.850782
## [1] 6.647058
## [1] 6.46428
# EL modelo powered exponencial de menor MSE es el de Máxima verosimilitud: 6.4642

plot(vari_temp, main ="Modelos powered exponenciales")
lines(fitvar_temp_pow1, col = "red")
lines(fitvar_temp_pow2, col = "green")
lines(fitvar_temp_pow3, col = "blue")
lines(fitvar_temp_pow4, col = "black")
legend("bottomright", legend = c("MCO", "MCP con npairs", "MCP con cressie", "MLE"), col =c("red", "green", "blue", "black"), lty = rep(1, 4))

# Modeos matern
ini1 <- c(11.14, 30997.59)
fitvar_temp_mat1 <- variofit(vari_temp, 
                            cov.model = "matern", 
                            ini1,
                            wei= "equal",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 0.21)
## variofit: covariance model used is matern 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
fitvar_temp_mat2 <- variofit(vari_temp, 
                            cov.model = "matern", 
                            ini1,
                            wei= "npairs",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 0.21)
## variofit: covariance model used is matern 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
fitvar_temp_mat3 <- variofit(vari_temp, 
                            cov.model = "matern", 
                            ini1,
                            wei= "cressie",
                            fix.nugget = TRUE,
                            nugget = 0,
                            kappa = 0.21)
## variofit: covariance model used is matern 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
fitvar_temp_mat4 <- likfit(temp.g, 
                           cov.model = "matern", 
                           ini.cov.pars = ini1,
                           fix.nugget = TRUE, 
                           nugget =  0)
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optimize.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optimize.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
## 
## WARNING: estimated range is less than 1 tenth of the minimum distance between two points. Consider re-examine the model excluding spatial dependence
mse1 <- MSE(vari_temp$v-fitvar_temp_mat1$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_mat1$cov.pars))
mse2 <- MSE(vari_temp$v-fitvar_temp_mat2$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_mat2$cov.pars))
mse3 <- MSE(vari_temp$v-fitvar_temp_mat3$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_mat3$cov.pars))
mse4 <- MSE(vari_temp$v-fitvar_temp_mat4$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_temp_mat4$cov.pars))
mse1;mse2;mse3;mse4
## [1] 9.040737
## [1] 8.128857
## [1] 7.094529
## [1] 6.909518
# El modelo con meno MSE es el de Máxima verosimilitud: 6.9095

plot(vari_temp, main ="Modelos matern")
lines(fitvar_temp_mat1, col = "red")
lines(fitvar_temp_mat2, col = "green")
lines(fitvar_temp_mat3, col = "blue")
lines(fitvar_temp_mat4, col = "black")
legend("bottomright", legend = c("MCO", "MCP con npairs", "MCP con cressie", "MLE"), col =c("red", "green", "blue", "black"), lty = rep(1, 4))

plot(vari_temp, main ="Modelos para la Temperatura")
lines(fitvar_temp_exp1, col = "red")
lines(fitvar_temp_pow4, col = "green")
lines(fitvar_temp_mat4, col = "blue")
legend("bottomright", legend = c("Exponencial - MCO", "Powered - MLE", "Matern - MLE"), col =c("red", "green", "blue"), lty = rep(1, 3))

Modelos para la variable: Precipitación

#eyefit(vari_precip)
#  cov.model sigmasq    phi tausq kappa kappa2  practicalRange
#1      wave     210 4500.8    50  <NA>   <NA>  13463.94711203
#2      wave     250 4000.8     0  <NA>   <NA> 11968.218895711

ini1 <- c(250, 4000.8)

fitvar_precip1 <- variofit(vari_precip, 
                        cov.model = "wave", 
                        ini1,
                        fix.nugget = TRUE, 
                        nugget =  88, 
                        wei= "equal")
## variofit: covariance model used is wave 
## variofit: weights used: equal 
## variofit: minimisation function used: optim
fitvar_precip2 <- variofit(vari_precip, 
                        cov.model = "wave", 
                        ini1,
                        fix.nugget = TRUE, 
                        nugget =  88, 
                        wei= "npairs")
## variofit: covariance model used is wave 
## variofit: weights used: npairs 
## variofit: minimisation function used: optim
fitvar_precip3 <- variofit(vari_precip, 
                        cov.model = "wave", 
                        ini1,
                        fix.nugget = TRUE, 
                        nugget =  88, 
                        wei= "cressie")
## variofit: covariance model used is wave 
## variofit: weights used: cressie 
## variofit: minimisation function used: optim
fitvar_precip4 <- likfit(precip.g, 
                         cov.model = "wave", 
                         ini.cov.pars = ini1,
                         fix.nugget = TRUE, 
                         nugget =  88)
## kappa not used for the wave correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
ini1 <- c(210, 4500.8)

fitvar_precip5 <- likfit(precip.g, 
                         cov.model = "wave", 
                         ini.cov.pars = ini1,
                         fix.nugget = TRUE, 
                         nugget =  88)
## kappa not used for the wave correlation function
## ---------------------------------------------------------------
## likfit: likelihood maximisation using the function optim.
## likfit: Use control() to pass additional
##          arguments for the maximisation function.
##         For further details see documentation for optim.
## likfit: It is highly advisable to run this function several
##         times with different initial values for the parameters.
## likfit: WARNING: This step can be time demanding!
## ---------------------------------------------------------------
## likfit: end of numerical maximisation.
mse1 <- MSE(vari_temp$v-fitvar_precip1$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_precip1$cov.pars))
mse2 <- MSE(vari_temp$v-fitvar_precip2$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_precip2$cov.pars))
mse3 <- MSE(vari_temp$v-fitvar_precip3$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_precip3$cov.pars))
mse4 <- MSE(vari_temp$v-fitvar_precip4$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_precip4$cov.pars))
mse5 <- MSE(vari_temp$v-fitvar_precip5$cov.pars[1] +cov.spatial(vari_temp$u, cov.model = "exponential",
                                                                  cov.pars = fitvar_precip5$cov.pars))
mse1;mse2;mse3;mse4;mse5
## [1] 25267.29
## [1] 28224.4
## [1] 30634.36
## [1] 25129.32
## [1] 30297.76
# El menor MSE se obtiene con el primer modelo de Máxima verosimiltud: 25129.32
# El segundo mejor es el de MCO: 25267.29

plot(vari_precip)
lines(fitvar_precip1, col = "red")
lines(fitvar_precip2, col = "green")
lines(fitvar_precip3, col = "blue")
lines(fitvar_precip4, col = "black")
lines(fitvar_precip5, col = "brown")
legend("bottomright", legend = c("MCO", "MCP con npairs", "MCP con cressie", "MLE", "MLE con otro ini"), col =c("red", "green", "blue", "black", "brown"), lty = rep(1, 4))

plot(vari_precip, main = "Mejor modelo Precipitación")
lines(fitvar_precip1, col = "red")
lines(fitvar_precip4, col = "black")

legend("bottomright", legend = c("MCO", "MLE"), col =c("red", "black"), lty = rep(1, 2))

Modelos de semivariograma comparativos

A continuación se explora la programación “a pedal” del modelo teórico exponencial con ayuda de la función optim. Se exploran el MCO, MCP y MLE.

Posteriormente estos resultados se compararan en términos de MSE con los obtenidos con el paquete geoR

ini_exp <- c(11.14, 30997.59)

datos <- data.frame(estim = vari_temp$v,
              h = vari_temp$u,
              n = vari_temp$n)

# Secuencia suave de h para dibujar las curvas
h_seq <- seq(min(datos$h), max(datos$h), length.out = 1000)

## - - - - Estimación MCO - - - - - 
mco <- nls(estim ~ sigma2*(1 - exp(-3*h / a)), data = datos, 
           start = c(sigma2 = ini_exp[1], a = ini_exp[2]))

# Predicciones del modelo
mco_teorico <- mco$m$predict()

## - - - - Estimación MCP npairs - - - - - 
W1 <- diag(x = 1/datos$n)

mcp_npairs <- nls(estim ~ sigma2*(1 - exp(-3*h / a)), data = datos, 
                  start = c(sigma2 = ini_exp[1], a = ini_exp[2]), weights = diag(W1))

# Predicciones del modelo
mcp_teorico <- mcp_npairs$m$predict()

# - - -  Estimación MCP Cressie - - - - 
cressie <- function(par, vari_temp){
  sigma2 <- par[1]
  a <- par[2]
  
  h = vari_temp$u
  # Modelo exponencial
  gamma_h = sigma2*(1 - exp(-3*h/a))
  W2 = diag((8*gamma_h^2)/vari_temp$n)
  
  # Lo que quiero minimizar
  f_gamma <- t(vari_temp$v-gamma_h)%*%solve(W2)%*%(vari_temp$v-gamma_h)
  return(as.numeric(f_gamma))
}

mcp_cressie_par <- optim(par = ini_exp, fn = cressie, vari_temp = vari_temp, method = "BFGS")
mcp_cressie_par$par # sigma y a estimados
## [1]    11.06908 30997.58955
mcp_cressie <- mcp_cressie_par$par[1]*(1 - exp(-3*datos$h/mcp_cressie_par$par[2]))

## - - -  Estimación MLE - - - - ##
dist_s <- as.matrix(dist(temp.g$coords))
temp_obs <- scale(df_clima$Temp, center = TRUE, scale = FALSE)

negloglik <- function(par, h, z) {
  sigma2 <- par[1]
  a      <- par[2]
  # Modelo exponencial
  gamma <- sigma2*(1 - exp(-3*h/a))  # matriz de semivariogramas
  Sigma <- sigma2 *1 - gamma            # C(h) = sigma2 - gamma(h)
  
  # Función a minimizar
  # Cholesky estable para log-det y solución
  cholS <- tryCatch(chol(Sigma), error = function(e) NULL)
  if (is.null(cholS)) return(1e12)  # penaliza si no PD
  logdet <- 2 * sum(log(diag(cholS)))
  y <- backsolve(cholS, z, transpose = TRUE)
  quad <- sum(y^2)
  negloglik <- 0.5 * (length(z) * log(2*pi) + logdet + quad)
  return(as.numeric(negloglik))
}

mle <- optim(par = ini_exp,
             fn = negloglik,
             h = dist_s, 
             z = temp_obs,
             method = "L-BFGS-B",
             control = list(maxit = 2000))
mle$par
## [1]    14.17377 30997.58219
mle_fit <- mle$par[1]*(1 - exp(-3*datos$h/mle$par[2]))

# Ya tengo las estimaciones de los modelos con MCO, MCP y MLE. Solo me falta graficarlos y contrastarlos mediante rmse con los obtenidos mediante geoR
par(mfrow = c(2, 2))

plot(vari_temp, main ="Modelos MCO")
lines(fitvar_temp_exp1, col = "red")
lines(vari_temp$u, mco_teorico, col = "blue", lwd = 2)
legend("bottomright", legend = c("Modelo geoR", "Modelo 'a pedal'"), col =c("red", "blue"), lty = rep(1, 2))

plot(vari_temp, main ="Modelos MCP npairs")
lines(fitvar_temp_exp2, col = "red")
lines(vari_temp$u, mcp_teorico, col = "blue", lwd = 2)
legend("bottomright", legend = c("Modelo geoR", "Modelo 'a pedal'"), col =c("red", "blue"), lty = rep(1, 2))

plot(vari_temp, main ="Modelos MCP cressie")
lines(fitvar_temp_exp3, col = "red")
lines(vari_temp$u, mcp_cressie, col = "blue", lwd = 2)
legend("bottomright", legend = c("Modelo geoR", "Modelo 'a pedal'"), col =c("red", "blue"), lty = rep(1, 2))

plot(vari_temp, main ="Modelos MLE")
lines(fitvar_temp_exp4, col = "red")
lines(vari_temp$u, mle_fit, col = "blue", lwd = 2)
legend("bottomright", legend = c("Modelo geoR", "Modelo 'a pedal'"), col =c("red", "blue"), lty = rep(1, 2))

De los gráficos anteriores se evidencia que para la familia de modelos exponenciales, el modelo de Mínimos cuadrados ordinarios es el único cuya estimación mediante nls coincide con la empleada por el paquete geoR. Los modelos MCP con matriz de pesos Cressie y MLE programados resultan ser más suaves, razón por la que se ajustan mejor a la curvatura del semivariograma empírico.

Kriging espacial

# Descargar los límites de los estados de México
mexico <- ne_states(country = "Mexico", returnclass = "sf")
# Ver qué estados tiene
unique(mexico$name)
##  [1] "Sonora"              "Baja California"     "Chihuahua"          
##  [4] "Coahuila"            "Tamaulipas"          "Nuevo León"         
##  [7] "Quintana Roo"        "Campeche"            "Tabasco"            
## [10] "Chiapas"             "Colima"              "Nayarit"            
## [13] "Baja California Sur" "Sinaloa"             "Yucatán"            
## [16] "Veracruz"            "Jalisco"             "Michoacán"          
## [19] "Guerrero"            "Oaxaca"              NA                   
## [22] "México"              "Puebla"              "Morelos"            
## [25] "Querétaro"           "Hidalgo"             "Guanajuato"         
## [28] "San Luis Potosí"     "Zacatecas"           "Aguascalientes"     
## [31] "Durango"             "Tlaxcala"            "Distrito Federal"
yucatan_region <- mexico %>%
  dplyr::filter(name %in% c("Yucatán", "Campeche"))

# Transformar a UTM zona 16N (EPSG:32616)
yucatan_region_utm <- st_transform(yucatan_region, crs = 32616)

Kriging Temperatura

x.seq <- seq(min(temp.g$coords[,1]), max(temp.g$coords[,1]), length = 100)
y.seq <- seq(min(temp.g$coords[,2]), max(temp.g$coords[,2]), length = 100)
grilla_full <- expand.grid(Var1 = x.seq, Var2 = y.seq)

kriging_temp <- krige.conv(temp.g, locations = grilla_full, 
                           krige = krige.control(type.krige = "SK", beta =30.93636, 
                                                 cov.pars = c(10.5146,13341.4454),
                                                 cov.model="exponential"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
krig_temp_df <- data.frame(
  Easting = grilla_full$Var1,
  Northing = grilla_full$Var2,
  pred = kriging_temp$predict,
  var = kriging_temp$krige.var
)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 

ggplot() +
  geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
  geom_tile(data = krig_temp_df, aes(x = Easting, y = Northing, fill = pred)) +
  geom_point(data = as.data.frame(temp.g$coords),
             aes(x = Easting, y = Northing),
             color = "black", size = 1) +
  scale_fill_viridis_c(name = "Temperatura (°C)") +
  coord_sf(xlim=c(120000, 350000), ylim = c(2180000, 2400000)) +
  labs(title = "Kriging de Temperatura - Norte de la Península de Yucatán",
       x = "Easting (m)", y = "Northing (m)") +
  theme_minimal()

ggplot() +
  geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
  geom_tile(data = krig_temp_df, aes(x = Easting, y = Northing, fill = var)) +
  scale_fill_viridis_c(name = "Varianza de predicción") +
  coord_sf(xlim=c(120000, 350000), ylim = c(2180000, 2400000)) +
  labs(
    title = "Varianza del error de predicción - Kriging",
    x = "Easting (m)", y = "Northing (m)"
  ) +
  theme_minimal()

# Leer shapefile nivel municipal
mex_estados <- st_read("gadm41_MEX_shp/gadm41_MEX_1.shp")
## Reading layer `gadm41_MEX_1' from data source 
##   `C:\Users\CAMILA\OneDrive - Universidad Nacional de Colombia\Escritorio\Universidad\Maestría Estadística\2do semestre\Espacial\gadm41_MEX_shp\gadm41_MEX_1.shp' 
##   using driver `ESRI Shapefile'
## Simple feature collection with 32 features and 11 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -118.3665 ymin: 14.53507 xmax: -86.71074 ymax: 32.71863
## Geodetic CRS:  WGS 84
# Filtrar los estados de la península
yucatan_region <- mex_estados[mex_estados$NAME_1 %in% c("Yucatán", "Campeche"), ]

yucatan_region_utm <- st_transform(yucatan_region, 32616)

# Crear una grilla regular dentro del polígono
grilla_sf <- st_make_grid(
  yucatan_region_utm,
  cellsize = 1000,       # tamaño de celda en metros (ajústalo según la resolución deseada)
  what = "centers",
  square = TRUE
)

# Convertir a objeto sf (para poder manejarlo fácilmente)
grilla <- st_as_sf(data.frame(geometry = grilla_sf))

grilla_sf <- st_as_sf(data.frame(geometry = grilla))
grilla_sf <- grilla_sf[yucatan_region_utm, ]

grilla_coords <- st_coordinates(grilla_sf)
kriging_temp <- krige.conv(temp.g, locations = grilla_coords, 
                           krige = krige.control(type.krige = "SK", beta =30.93636, 
                                                 cov.pars = c(10.5146,13341.4454),
                                                 cov.model="exponential"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
krig_temp_df <- data.frame(
  X = grilla_coords[, 1],
  Y = grilla_coords[, 2],
  pred = kriging_temp$predict,
  var = kriging_temp$krige.var
)

ggplot() +
  # superficie continua del kriging
  geom_tile(data = krig_temp_df,aes(x = X, y = Y, fill = pred)) +
  scale_fill_viridis_c(name = "Temperatura (°C)") +
  # límites de los estados
  geom_sf(data = yucatan_region_utm,fill = NA,color = "black",linewidth = 0.6) +
  # puntos observados
  geom_point(data = as.data.frame(temp.g$coords),aes(x = Easting, y = Northing),
             color = "black",size = 1) +
  # recorte manual para centrarse en la región
  coord_sf(xlim = c(140000, 350000),ylim = c(2170000, 2350000)) +
  labs(
    title = "Predicciones Kriging de Temperatura",
    x = "Easting (m)",
    y = "Northing (m)"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    plot.background = element_rect(fill = "white"),
    legend.position = "right"
  )

ggplot() +
  # superficie continua del kriging
  geom_tile(data = krig_temp_df,aes(x = X, y = Y, fill = var)) +
  scale_fill_viridis_c(name = "VaR E.P.") +
  # límites de los estados
  geom_sf(data = yucatan_region_utm,fill = NA,color = "black",linewidth = 0.6) +
  # puntos observados
  #geom_point(data = as.data.frame(temp.g$coords),aes(x = Easting, y = Northing),
  #           color = "black",size = 1) +
  # recorte manual para centrarse en la región
  coord_sf(xlim = c(140000, 350000),ylim = c(2170000, 2350000)) +
  labs(
    title = "Varianza del error de predicción - Temperatura",
    x = "Easting (m)",
    y = "Northing (m)"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    plot.background = element_rect(fill = "white"),
    legend.position = "right"
  )

Kriging Precipitación

x.seq <- seq(min(precip.g$coords[,1]), max(precip.g$coords[,1]), length = 100)
y.seq <- seq(min(precip.g$coords[,2]), max(precip.g$coords[,2]), length = 100)
grilla_full <- expand.grid(Var1 = x.seq, Var2 = y.seq)

kriging_precip <- krige.conv(precip.g, locations = grilla_full, 
                             krige = krige.control(type.krige = "SK", beta =24.42, 
                                                   cov.pars = c(172.24 ,3566.44),
                                                   cov.model="wave"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
krig_precip_df <- data.frame(
  Easting = grilla_full$Var1,
  Northing = grilla_full$Var2,
  pred = kriging_precip$predict,
  var = kriging_precip$krige.var
)

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

ggplot() +
  geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
  geom_tile(data = krig_precip_df, aes(x = Easting, y = Northing, fill = pred)) +
  geom_point(data = as.data.frame(temp.g$coords),
             aes(x = Easting, y = Northing),
             color = "black", size = 1) +
  scale_fill_viridis_c(name = "Temperatura (°C)") +
  coord_sf(xlim=c(120000, 350000), ylim = c(2180000, 2400000)) +
  labs(title = "Kriging de Precipitación - Norte de la Península de Yucatán",
       x = "Easting (m)", y = "Northing (m)") +
  theme_minimal()

ggplot() +
  geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
  geom_tile(data = krig_precip_df, aes(x = Easting, y = Northing, fill = var)) +
  scale_fill_viridis_c(name = "Varianza de predicción") +
  coord_sf(xlim=c(120000, 350000), ylim = c(2180000, 2400000)) +
  labs(
    title = "Varianza del error de predicción - Precipitación",
    x = "Easting (m)", y = "Northing (m)"
  ) +
  theme_minimal()

kriging_precip <- krige.conv(precip.g, locations = grilla_coords, 
                             krige = krige.control(type.krige = "SK", beta =24.42, 
                                                   cov.pars = c(172.24 ,3566.44),
                                                   cov.model="wave"))
## krige.conv: model with constant mean
## krige.conv: Kriging performed using global neighbourhood
krig_precip_df <- data.frame(
  X = grilla_coords[, 1],
  Y = grilla_coords[, 2],
  pred = kriging_precip$predict,
  var = kriging_precip$krige.var
)

ggplot() +
  # superficie continua del kriging
  geom_tile(data = krig_precip_df,aes(x = X, y = Y, fill = pred)) +
  scale_fill_viridis_c(name = "Precipitación") +
  # límites de los estados
  geom_sf(data = yucatan_region_utm,fill = NA,color = "black",linewidth = 0.6) +
  # puntos observados
  geom_point(data = as.data.frame(temp.g$coords),aes(x = Easting, y = Northing),
             color = "black",size = 1) +
  # recorte manual para centrarse en la región
  coord_sf(xlim = c(140000, 350000),ylim = c(2170000, 2350000)) +
  labs(
    title = "Predicciones Kriging de Precipitación",
    x = "Easting (m)",
    y = "Northing (m)"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    plot.background = element_rect(fill = "white"),
    legend.position = "right"
  )

ggplot() +
  # superficie continua del kriging
  geom_tile(data = krig_precip_df,aes(x = X, y = Y, fill = var)) +
  scale_fill_viridis_c(name = "Var E.P.") +
  # límites de los estados
  geom_sf(data = yucatan_region_utm,fill = NA,color = "black",linewidth = 0.6) +
  # puntos observados
  geom_point(data = as.data.frame(temp.g$coords),aes(x = Easting, y = Northing),
             color = "black",size = 1) +
  # recorte manual para centrarse en la región
  coord_sf(xlim = c(140000, 350000),ylim = c(2170000, 2350000)) +
  labs(
    title = "Varianza del error de predicción - Precipitación",
    x = "Easting (m)",
    y = "Northing (m)"
  ) +
  theme_bw() +
  theme(
    panel.grid = element_blank(),
    panel.background = element_rect(fill = "white"),
    plot.background = element_rect(fill = "white"),
    legend.position = "right"
  )

# GEOESTADÍSTICA ESPACIO_TIEMPO

Para el análsis espacio-temporal se elige la variable de Temperatura máxima diaria (ºC) en las 33 estaciones ubicadas al sureste de México en la península de Yucatán para el periodo de tiempo comprendido entre 01-11-2021 y 2l 20-12-2021, resultando en 50 momentos consecutivos en el tiempo.

base_TMAX <- base_TMAX %>% 
  filter(FECHA >= as.POSIXct("2021-11-01", tz="UTC") & FECHA <= as.POSIXct("2021-12-20", tz="UTC"))

t <- 1:50

Análisis de medias

Pimero, se analiza el comportamiento de la media temporal con ayuda del siguiente gráfico

media_temporal <- base_TMAX %>%
  select(-FECHA) %>%
  mutate(media = rowMeans(., na.rm = TRUE))

plot(base_TMAX$FECHA, media_temporal$media, type = "l",
     xlab = "Fecha", ylab = "Temperatura promedio",
     main = "Media temporal de TMAX (promedio sobre estaciones)")

El gráfico nos deja ver que no se evidencia una tendencia en las medias en los tiempos estudiados. Estas parecen más tener un comportamiento de ruido blanco: media constante y homogeneidad de varianza.

Ahora veremos el comportamiento de la varable Temperatura para las diferenctes estaciones de monitoreo, con el fin de buscar por relaciones espaciales de la media.

media_espacial <- base_TMAX %>%
  select(-FECHA) %>%
  summarise(across(everything(), \(x) mean(x, na.rm = TRUE))) %>%
  t()

# Gráfico
boxplot(base_TMAX %>% select(-FECHA),
        main = "Distribución de TMAX por estación",
        las = 2)

df_media_tmax <- cbind(utm_coords_df, Media_TMax = media_espacial)
pairs(df_media_tmax)

Como se vio en el caso espacial, la variable promedio de la Temperatura máxima diaria no presenta una tendencia en relación con las coordenadas espaciales.

matplot(base_TMAX$FECHA,
        base_TMAX %>% select(-FECHA),
        type = "l", lty = 1,
        xlab = "Fecha", ylab = "TMAX",
        main = "Series temporales por estación")

Las series temporales para todas las ubicaciones espaciales tienen comportamientos bastante similares, hay unas cuantas cuya media parece variar levemente de la media global pero nada que de indicios de que la media tiene estructura de tipo espacial.

Todo este análisis conclute que estamos bajo un análsisis espacio-temporal de media constante. En este sentido, si posteriormente quisieramos centrar los datos lo que haremos será restar la media global.

Estimación del semivariograma espacio-temporal

times <- base_TMAX$FECHA                     # 50 fechas
temp_matrix <- as.data.frame(base_TMAX[,-1]) # quitar columna FECHA

temp_matrix_long <- as.vector(t(temp_matrix))  
df_temp <- data.frame(TMAX = temp_matrix_long)

temp_stf <- STFDF(
  sp = utm_coords,
  time = times,
  data = df_temp
)
sp_cuts  <- 15  # número de bins espaciales
tm_cuts  <- 10  # número de bins temporales

variST_temp <- variogramST(
  TMAX ~ 1,
  data = temp_stf,
  assumeRegular = FALSE
)

df_vario <- variST_temp %>%
  select(spacelag, timelag, gamma, np) %>%
  rename(
    h = spacelag,
    u = timelag,
    gamma_emp = gamma,
    N = np
  ) %>%
  filter(!is.na(gamma_emp), N > 0)
plot(variST_temp, wireframe = TRUE)

plot(variST_temp, map = TRUE)

Modelos para los semivariogramas

Para la estimación de los modelos teóricos se decidió explorar las funciones de covarianza Cressie y Exponencial. La priemra de estas entra dentro de la categoría de funciones no separables, mientras que la Exponencial si lo es.

# - - - -   Organizar los datos en formato largo    - - - - #

# Creamos índice de tiempo (1 a 50)
base_TMAX <- base_TMAX %>% mutate(t = 1:n())

# Pasamos a formato largo
long_TMAX <- base_TMAX %>%
  pivot_longer(cols = starts_with("s"),
               names_to = "station",
               values_to = "TMAX") %>%
  mutate(station = as.integer(gsub("s", "", station))) %>%
  left_join(utm_coords_df %>% mutate(station = 1:nrow(.)), by = "station")

head(long_TMAX)
## # A tibble: 6 × 6
##   FECHA                   t station  TMAX Easting Northing
##   <dttm>              <int>   <int> <dbl>   <dbl>    <dbl>
## 1 2021-11-01 00:00:00     1       1    31 212519. 2214240.
## 2 2021-11-01 00:00:00     1       2    35 173715. 2213196.
## 3 2021-11-01 00:00:00     1       3    35 162137. 2210460.
## 4 2021-11-01 00:00:00     1       4    30 168065. 2229804.
## 5 2021-11-01 00:00:00     1       5    30 172493. 2234861.
## 6 2021-11-01 00:00:00     1       6    30 173684. 2236378.
# - - - - - - -  Matrices de distancias - - - - - - - - - - #

coords <- long_TMAX %>% select(Easting, Northing)
time_vec <- long_TMAX$t

dist_s <- as.matrix(dist(coords))
dist_t <- as.matrix(dist(time_vec))
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - -  Funciones de covarianza - - - - - - - - - - - - - -  #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

# Función de covarianza Cressie
cressie1 <- function(h, u, param){
  u <- as.numeric(u)  # evitar errores con difftime

  sigma <- param[1] # sigma = la silla
  at    <- param[2] # parametro de escala temporal ------> mean(apply(base_TMAX, 2, var), na.rm = TRUE)
  as    <- param[3] # parametro de escala espacial

  (sigma^2 / (1 + (at^2) * u^2)) * exp(-(as^2) * h^2 / (1 + (at^2) * u^2))
}

gamma_cressie <- function(h, u, param){
  sigma <- param[1]
  sigma^2 - cressie1(h, u, param)
}


# Función de covarianza Exp
exp_esp_temp <- function(h, u, param){
  u <- as.numeric(u)

  sigma <- param[1] # (silla o varianza total)
  a_s   <- param[2] # parámetro de rango espacial
  a_t   <- param[3] # parámetro de rango temporal

  sigma^2 * exp(-h/a_s - u/a_t)
}

gamma_exp <- function(h, u, param){
  sigma <- param[1]
  sigma^2 - exp_esp_temp(h, u, param)
}
z_obs <- scale(long_TMAX$TMAX, center = TRUE, scale = FALSE)  # CENTRAR LOS DATOS

Cressie

# Log verosimilitud negativa a minimizar con optim
negloglik <- function(p, z, h, u) {
  # p: vector de parámetros positivos
  if (any(p <= 0)) return(Inf)
  Sigma <- cressie1(h, u, p)
  # Para estabilidad numérica
  Sigma <- Sigma + diag(1e-8, nrow(Sigma))
  # Log-verosimilitud multivariada normal
  detS <- determinant(Sigma, logarithm = TRUE)$modulus
  quad <- t(z) %*% solve(Sigma, z)
  0.5 * (detS + quad + length(z) * log(2 * pi))
}

# Estimación con la función optim
est_mle_cressie <- optim(
    par = c(11.14,3.060346 ,30997.59), 
    fn = negloglik,
    z = z_obs,
    h = dist_s,
    u = dist_t,
    method = "L-BFGS-B",
    lower = c(0.001, 0.001, 0.001),
    upper = c(20, 10, 50000)
)

param_mle_cressie <- est_mle_cressie$par
param_mle_cressie
## [1] 2.134305e+00 9.007268e-01 3.099759e+04
# Función objetivo a minimizar
obj_ols_cressie <- function(param) {
  with(df_vario,
       sum((gamma_emp - gamma_cressie(h, u, param))^2)
  )
}

est_ols_cressie <- optim(
  par = c(11.14, 3.06, 30997.6),
  fn = obj_ols_cressie,
  method = "L-BFGS-B",
  lower = c(0.001, 0.001, 0.001),
  upper = c(50, 20, 100000)
)

param_ols_cressie <- est_ols_cressie$par
vg <- df_vario %>%
  rename(h = h, u = u, N = N, v = gamma_emp
  )

# Función objetivo
obj_mcp <- function(param, vg, modelo_gamma){

  gamma_theo <- modelo_gamma(vg$h, vg$u, param)
  gamma_emp  <- vg$v

  # pesos diagonales
  w <- 8*(gamma_theo^2)/vg$N

  # evitar pesos nulos
  w[w < 1e-12] <- 1e-12

  r <- 2*gamma_emp - 2*gamma_theo

  sum( (r^2) / w )
}

fit_mcp_cressie <- optim(
  par = c(11.14, 3.06, 30997.6),
  fn = obj_mcp,
  vg = vg,
  modelo_gamma = gamma_cressie,
  method = "L-BFGS-B",
  lower = c(1e-6, 1e-6, 1e-6)
)

param_mcp_cressie <- fit_mcp_cressie$par

Exponencial

# Función negativa de log verosimilitud 
negloglik_exp <- function(p, z, h, u) {
  if (any(p <= 0)) return(Inf)
  
  Sigma <- exp_esp_temp(h, u, p)  # ← aquí cambia
  
  # Para estabilidad numérica
  Sigma <- Sigma + diag(1e-8, nrow(Sigma))
  
  detS <- determinant(Sigma, logarithm = TRUE)$modulus
  quad <- t(z) %*% solve(Sigma, z)
  
  0.5 * (detS + quad + length(z) * log(2 * pi))
}

# Estimación de sigma con optim
est_exp <- optim(
  par = c(11.14, 30997.59, 3.060346),     # σ, rango espacial, rango temporal
  fn = negloglik_exp,
  z = z_obs,
  h = dist_s,
  u = dist_t,
  method = "L-BFGS-B",
  lower = c(0.001, 1000, 0.1),
  upper = c(50, 100000, 20)
)

param_mle_exp <- est_exp$par
param_mle_exp
## [1]     3.674489 30997.589218     2.430244
# Función objetivo a minimizar
obj_ols_exp <- function(param) {
  with(df_vario,
       sum((gamma_emp - gamma_exp(h, u, param))^2)
  )
}

est_ols_exp <- optim(
  par = c(11.14, 30997.59, 3.060346),
  fn = obj_ols_exp,
  method = "L-BFGS-B",
  lower = c(0.001, 1000, 0.1),
  upper = c(50, 100000, 20)
)

param_ols_exp <- est_ols_exp$par
param_ols_exp
## [1]     2.707631 31007.532050    18.983064
fit_mcp_exp <- optim(
  par = c(11.14, 30997.59, 3.060346),
  fn = obj_mcp,
  vg = vg,
  modelo_gamma = gamma_exp,
  method = "L-BFGS-B",
  lower = c(0.001, 1000, 0.1),
  upper = c(50, 100000, 20)
)

param_mcp_exp <- fit_mcp_exp$par
param_mcp_exp
## [1]     2.585668 30997.694681     4.663507

Habiendo generado todos los modelos anteriores se procede a compararlos mediante RMSE para ver cual se ajusta mejor al semivariograma espacio-tiempo empírico, dada la imposibilidad de contrastarlos gráficamente.

rmse_model <- function(param, df, gamma_fun){
  gamma_theo <- gamma_fun(df$h, df$u, param)
  sqrt(mean((df$gamma_emp - gamma_theo)^2))
}

rmse_cressie_ols <- rmse_model(param_ols_cressie, df_vario, gamma_cressie)
rmse_cressie_mcp <- rmse_model(param_mcp_cressie, df_vario, gamma_cressie)
rmse_exp_ols     <- rmse_model(param_ols_exp, df_vario, gamma_exp)
rmse_exp_mcp     <- rmse_model(param_mcp_exp, df_vario, gamma_exp)

rmse_cressie_ols; rmse_cressie_mcp; rmse_exp_ols; rmse_exp_mcp
## [1] 1.466216
## [1] 1.661467
## [1] 1.284865
## [1] 1.718555
k <- 3
AIC_cressie <- 2*k - 2*(-est_mle_cressie$value)
AIC_exp <- 2*k - 2*(-est_exp$value)
AIC_cressie; AIC_exp # Gana el modelo cressie en el metodo MLE
## [1] 6582.995
## [1] 6804.825
BIC_cressie <- k*log(length(z_obs)) - 2*(-est_mle_cressie$value)
BIC_exp <- k*log(length(z_obs)) - 2*(-est_exp$value)
BIC_cressie; BIC_exp
## [1] 6599.22
## [1] 6821.051

En términos de AIC y BIC el modelo de máxima verosimilitud cressie es el mejor se desempeña.

Kriging espacio-temporal

# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #
# - - - - - - - - - Kriging espacio-temporal - - - - - - - - - - - - -  - #
# - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - #

## Seleccionar aleatoriamente ciertas filas de mis ubicaciones
set.seed(19)
coords_muestra <- utm_coords_df[sample(nrow(utm_coords_df),10),]

pred_grid <- expand.grid(Easting = coords_muestra$Easting, 
                         Northing = coords_muestra$Northing, 
                         t = 1:50)

# Unir los puntos a predecir con los observados
coords_all <- rbind(coords, pred_grid[, c("Easting", "Northing")])
time_all <- c(time_vec, pred_grid$t)

# Nuevas distancias
dist_s_all <- as.matrix(dist(coords_all))
dist_t_all <- as.matrix(dist(time_all))

# Estimar la covarianza
Sigma_all <- cressie1(dist_s_all, dist_t_all, param_mle_cressie)

n_obs <- nrow(coords)
n_pred <- nrow(pred_grid)

Sigma_oo <- Sigma_all[1:n_obs, 1:n_obs]
Sigma_po <- Sigma_all[(n_obs + 1):(n_obs + n_pred), 1:n_obs]
Sigma_pp <- Sigma_all[(n_obs + 1):(n_obs + n_pred), (n_obs + 1):(n_obs + n_pred)]

# Predicciones
lambda_cressie <- Sigma_po %*% solve(Sigma_oo)
z_pred_cressie <- as.numeric(lambda_cressie %*% z_obs)
Var_pred <- diag(Sigma_pp - Sigma_po %*% solve(Sigma_oo) %*% t(Sigma_po))

media_TMAX <- attr(z_obs, "scaled:center")

# Volver a la escala original
z_pred_cressie_original <- z_pred_cressie + media_TMAX
Sigma_all_exp <- exp_esp_temp(dist_s_all, dist_t_all, param_mle_exp)

n_obs <- nrow(coords)
n_pred <- nrow(pred_grid)

Sigma_oo_exp <- Sigma_all_exp[1:n_obs, 1:n_obs]
Sigma_po_exp <- Sigma_all_exp[(n_obs + 1):(n_obs + n_pred), 1:n_obs]
Sigma_pp_exp <- Sigma_all_exp[(n_obs + 1):(n_obs + n_pred), (n_obs + 1):(n_obs + n_pred)]

lambda_exp <- Sigma_po_exp %*% solve(Sigma_oo_exp)
z_pred_exp <- as.numeric(lambda_exp %*% z_obs)
Var_pred_exp <- diag(Sigma_pp_exp - Sigma_po_exp %*% solve(Sigma_oo_exp) %*% t(Sigma_po_exp))

# Volver a la escala original
media_TMAX <- attr(z_obs, "scaled:center")
z_pred_original_exp <- z_pred_exp + media_TMAX
hist(z_pred_cressie_original, main = "Histograma de predicciones Cressie")

hist(z_pred_original_exp, main = "Histograma de predicciones Exponencial")

hist(Var_pred, main = "Histograma de varianzas del error de predición Cressie")

hist(Var_pred_exp, main = "Histograma de varianzas del error de predición Exponencial")

Los histogramas dejan ver que a pesar de que le modelo Cressie se desempeó mejor en la parte explicativa del modelo de semivarianza, es muy malo prediciendo ya que casi todas sus predicciones se van a un valor constante, teniendo frecuencias de varia de error de predicción altísimas. Caso contrario, el, modelo exponencial que si parece predecir adecuadamente todo el espectro de datos.

Ahora, se seleccionan los 5 horizontes de tiempo más cercanos para llevar a cabo mapas de predicción en estos.

horizontes <- c(1, 2, 3, 4, 5)   # días adelante
t0 <- max(time_vec)

tiempos_futuros <- t0 + horizontes
tiempos_futuros
## [1] 51 52 53 54 55
grid_pred <- expand.grid(Easting = utm_coords_df$Easting, 
                         Northing = utm_coords_df$Northing)
pred_kriging_horizonte <- function(t_pred, tipo = "mle") {

  # grid espacial
  grid_pred <- expand.grid(
    Easting  = seq(min(coords$Easting),  max(coords$Easting),  length.out = 60),
    Northing = seq(min(coords$Northing), max(coords$Northing), length.out = 60)
  )
  grid_pred$t <- t_pred
  
  # unir coords
  coords_all <- rbind(coords, grid_pred[, c("Easting", "Northing")])
  time_all   <- c(time_vec, rep(t_pred, nrow(grid_pred)))
  
  # distancias
  dist_s_all <- as.matrix(dist(coords_all))
  dist_t_all <- as.matrix(dist(time_all))
  
  # elegir tipo de modelo
  if (tipo == "mle") {
    Sigma_all <- exp_esp_temp(dist_s_all, dist_t_all, param_mle_exp)
  } else if (tipo == "mco") {
    Sigma_all <- exp_esp_temp(dist_s_all, dist_t_all, param_ols_exp)
  } else if (tipo == "mcp"){
    Sigma_all <- exp_esp_temp(dist_s_all, dist_t_all, param_mcp_exp)
  }
  
  n_obs  <- nrow(coords)
  n_pred <- nrow(grid_pred)
  
  Sigma_oo <- Sigma_all[1:n_obs, 1:n_obs]
  Sigma_po <- Sigma_all[(n_obs+1):(n_obs+n_pred), 1:n_obs]
  Sigma_pp <- Sigma_all[(n_obs+1):(n_obs+n_pred), (n_obs+1):(n_obs+n_pred)]
  
  # kriging
  lambda <- Sigma_po %*% solve(Sigma_oo)
  z_pred <- as.numeric(lambda %*% z_obs)
  
  variance_pred <- diag(Sigma_pp - Sigma_po %*% solve(Sigma_oo) %*% t(Sigma_po))
  
  # re-escalar
  z_pred_original <- z_pred + media_TMAX
  
  # data.frame final
  krig_df <- data.frame(
    Easting  = grid_pred$Easting,
    Northing = grid_pred$Northing,
    prediction = z_pred_original,
    variance_prediction = variance_pred
  )
  
  return(krig_df)
}

Se usa el modelo Exponencial de tipos MCO, MCP y MLE

lista_pred_mle <- lapply(tiempos_futuros, pred_kriging_horizonte, tipo = "mle")
names(lista_pred_mle) <- paste0("t+", horizontes)

lista_pred_mco <- lapply(tiempos_futuros, pred_kriging_horizonte, tipo = "mco")
names(lista_pred_mco) <- paste0("t+", horizontes)

lista_pred_mcp <- lapply(tiempos_futuros, pred_kriging_horizonte, tipo = "mcp")
names(lista_pred_mcp) <- paste0("t+", horizontes)
library(viridis)
## Cargando paquete requerido: viridisLite
plot_mapa_pred <- function(df_pred, titulo = "Mapa de Predicción") {
  
  ggplot() +
    geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
    geom_tile(
      data = df_pred,
      aes(x = Easting, y = Northing, fill = prediction)
    ) +
    geom_point(
      data = as.data.frame(temp.g$coords),
      aes(x = Easting, y = Northing),
      color = "black",
      size = 1
    ) +
    scale_fill_viridis_c(
      name = "Temperatura (°C)",
      limits = c(min_global, max_global)  # ← AQUÍ
    ) +
    coord_sf(xlim = c(120000, 350000),
             ylim = c(2180000, 2300000)) +
    labs(title = titulo, x = "Easting (m)", y = "Northing (m)") +
    theme_minimal()
}

min_global <- min(sapply(lista_pred_mle, function(df) min(df$prediction)))
max_global <- max(sapply(lista_pred_mle, function(df) max(df$prediction)))

for (i in seq_along(lista_pred_mle)) {
  horizonte <- names(lista_pred_mle)[i]
  df <- lista_pred_mle[[i]]
  
  p <- plot_mapa_pred(df, titulo = paste("Predicción Kriging Modelo Exp MLE –", horizonte))
  print(p)
}

for (i in seq_along(lista_pred_mco)) {
  horizonte <- names(lista_pred_mco)[i]
  df <- lista_pred_mco[[i]]
  
  p <- plot_mapa_pred(df, titulo = paste("Predicción Kriging Modelo Exp MCO –", horizonte))
  print(p)
}

for (i in seq_along(lista_pred_mcp)) {
  horizonte <- names(lista_pred_mcp)[i]
  df <- lista_pred_mcp[[i]]
  
  p <- plot_mapa_pred(df, titulo = paste("Predicción Kriging Modelo Exp MCP –", horizonte))
  print(p)
}

Como esperado, con el aumento del horizonte de pronóstico las predicciones tienden hacia el valor promedio de Temperatura diaria máxima, y este comportamiento se presenta sin diferencias con el tipo de modelo usado. Así mismo se espera que la varianza del error de predicción aumente a medida que nos alejamos en el horizont de pronóstico.

GEOESTADÍSTICA FUNCIONAL

Otra forma de abordar los datos espacio-tiempo es ajustando una curva a la serie temporal, para cada una de las ubicaciones espaciales. Es así como entramos en el mundo de los Datos funcionales.

Lo primero que se hará será completar los datos de las series faltantes (NA en momentos del tiempo en los que no hubo medición de Temperatura máxima) como sigue:

Completar los datos faltantes

library(fda)
## Warning: package 'fda' was built under R version 4.5.2
## Cargando paquete requerido: splines
## Cargando paquete requerido: fds
## Warning: package 'fds' was built under R version 4.5.2
## Cargando paquete requerido: rainbow
## Warning: package 'rainbow' was built under R version 4.5.2
## Cargando paquete requerido: MASS
## 
## Adjuntando el paquete: 'MASS'
## The following object is masked from 'package:dplyr':
## 
##     select
## Cargando paquete requerido: pcaPP
## Cargando paquete requerido: RCurl
## Warning: package 'RCurl' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'RCurl'
## The following object is masked from 'package:tidyr':
## 
##     complete
## Cargando paquete requerido: deSolve
## Warning: package 'deSolve' was built under R version 4.5.2
## 
## Adjuntando el paquete: 'fda'
## The following object is masked from 'package:graphics':
## 
##     matplot
## The following object is masked from 'package:datasets':
## 
##     gait
library(fda.usc)
## Warning: package 'fda.usc' was built under R version 4.5.2
## Cargando paquete requerido: mgcv
## Cargando paquete requerido: nlme
## 
## Adjuntando el paquete: 'nlme'
## The following object is masked from 'package:dplyr':
## 
##     collapse
## This is mgcv 1.9-3. For overview type 'help("mgcv-package")'.
## Cargando paquete requerido: knitr
##  fda.usc is running sequentially usign foreach package
##  Please, execute ops.fda.usc() once to run in local parallel mode
##  Deprecated functions: min.basis, min.np, anova.hetero, anova.onefactor, anova.RPm
##  New functions: optim.basis, optim.np, fanova.hetero, fanova.onefactor, fanova.RPm
## ----------------------------------------------------------------------------------
library(splines)
library(SpatFD)
## Warning: package 'SpatFD' was built under R version 4.5.2
library(RColorBrewer)
base_TMAX <- read_excel("base_TMAX_2.xlsx", sheet = "Sheet1")

base_TMAX_filled <- base_TMAX %>%
  rowwise() %>%
  mutate(
    mean_day = round(mean(c_across(starts_with("s")), na.rm = TRUE),1),
    across(starts_with("s"), ~ ifelse(is.na(.x), mean_day, .x))
  ) %>%
  ungroup() 
base_TMAX_filled <- base_TMAX_filled[,-35]

base_TMAX_filled <- base_TMAX_filled %>% 
  filter(FECHA >= as.POSIXct("2018-11-01", tz="UTC") & FECHA <= as.POSIXct("2021-12-20", tz="UTC"))

Visualizar los datos discretizados

Una vez completadas las series temporales se visualizan los datos

n_puntos = 1146
n_curvas = 33

colors <- rainbow(n_curvas) 

colors <- brewer.pal(min(n_curvas, 12), "Paired") # Usar Set3 para tonos suaves
if (n_curvas > 12) { # Si hay más filas que colores en Set3, repetir los colores
  colors <- colorRampPalette(colors)(n_curvas)
}


plot(1:n_puntos, t(base_TMAX_filled[,-1])[1,], type = "n", # Crear un marco vacío
     xlab = "Emisiones de onda (nm)", ylab = "Intensidad",
     main = "Curvas de la Temperatura máxima",
     ylim = c(15,50))
grid()
for (i in 1:n_curvas) {
  lines(1:n_puntos, t(base_TMAX_filled[,-1])[i,], col = colors[i], lwd = 2) # Añadir cada curva
}

Suavizamiento

Bajo el hecho de que las funciones infinito-dimensionales son observadas solo en un numero finito de periodos, se ajusta un modelo para reconstruir la variable funcional. Este modelo no es más que una combinación lineal de funciones base ortonormales.

Dado que los datos parecen presentar un leve comportamiento estacional (sin ser completamente periódicos) y además dan cuenta de cambios locales, se opta por un suavizamiento de tipo BSpline.

Con el fin de establecer el mejor parámetro de suavizado \(\lambda\) se emplea validación cruzada, como se muestra en el siguiente código

# Extraer matriz de datos funcionales (tiempo x estaciones)
Y <- as.matrix(base_TMAX_filled[,-1])
argvals <- 1:nrow(base_TMAX_filled)

# Número de bases (suficiente para capturar variación diaria)
nbasis <- 40
rangeval <- range(argvals)

# Crear base B-spline
basis_bs <- create.bspline.basis(rangeval = rangeval, nbasis = nbasis)

# Penalizar segunda derivada
Lfd <- int2Lfd(2)

# Vector de lambdas para búsqueda GCV
lambdas <- 10^seq(-6, 6, length.out = 40)
gcv_values <- numeric(length(lambdas))

for(i in seq_along(lambdas)){
  fdPar_obj <- fdPar(basis_bs, Lfd, lambdas[i])
  res <- smooth.basis(argvals = argvals, y = Y, fdParobj = fdPar_obj)
  gcv_values[i] <- mean(res$gcv)
}

best_lambda <- lambdas[which.min(gcv_values)]
best_lambda
## [1] 3455.107
fdPar_best <- fdPar(basis_bs, Lfd, best_lambda)
tmax_fd <- smooth.basis(argvals, Y, fdPar_best)$fd

plot(tmax_fd, lwd=2)
## [1] "done"
title("Suavizamiento B-spline penalizado de TMAX")

matriz_curva <- eval.fd(1:1146, tmax_fd)

Función media

mean_temp <- mean.fd(tmax_fd)

plot(tmax_fd, main="Curvas de la Temperatura diaria máxima",
     col = colors, xlab = "Días", ylab = "Temperatura")
## [1] "done"
lines(mean_temp, col = 1, lwd = 2)
legend("topright", legend = c("Función Media"), fill = c("black"))

Note que la función media parece ser una buena función resumen del comportamiento de los datos, ya que parece que las demás curvas están construidas alrededor de esta.

Boxplot Funcional

Una parte importante en el análisis descriptivo de los datos funcionales es la identificación de posibles datos atípicos (en este caso de magnitud). Una forma intuitiva y sencilla de llevar esto a cabo es mediante la extensión del boxplot a FDA: El boxplot funcional, el cual usa las profundidades como “medidas de orden”. Se dice que una curva es atípica si se encuentra por fuera de la región central de curvas más profundas.

La diferencia entre MBD y BD es la función de profundidad que usan.

## - - - - - - - - - -  BOXPLOT FUNCIONAL - - - - - - - - - - - 
# Metodo MBD
tmax_fbplot1 <- fbplot(matriz_curva, method = "MBD", ylim = c(23, 45), fullout = TRUE, 
                         main= "Boxplot funcional Temperatura maxima \n Método MBD", 
                         ylab = "Intensidad", xlab = "Emisiones de onda (nm)") 

tmax_fbplot1$outpoint
## [1] 11 17 24 26
# Metodo BD
tmax_fbplot2 <- fbplot(matriz_curva, method = "BD2", ylim = c(23, 45), fullout = TRUE, 
                         main= "Boxplot funcional Temperatura maxima \n Método BD", 
                         ylab = "Intensidad", xlab = "Emisiones de onda (nm)") 

tmax_fbplot2$outpoint
## [1] 10 11 24 32

Ambos boxplot nos dejan ver que contamos con 4 curvas outliers, coincidiendo en la identificación de las estaciones de monitoreo 11 y 24 como atípicas.

FPCA

Antes de continuar con la predicción es importane llevar a cabo un análisis de componentes principales funcional, ya que este fundamenta varios procedimientos posteriores.

Con FPCA se busca proyectar la estructura de covarianza \(C\) sobre unas curvas \(\xi\) que garanticen que dicha proyección tenga varianza máxima. Mediante la maximización de la covarianza se encuentra que estas curvas no son más que las funciones propias del operador de covarianza.

Así, los puntajes/scores se definen como \[f_{\mathbf{s}_i}^k=\langle \mathcal{X}_{\mathbf{s}_i}, \xi^k\rangle \qquad k = 1, ..., K\qquad i =1, ...,n\]

estos son proyecciones de los datos funcionales sobre la respectiva función propia. Estos puntajes son escalares, lo que hace sencillo su análisis.

fpca_temp <- pca.fd(tmax_fd, nharm = 4, centerfns = TRUE)

plot(fpca_temp$harmonics, main = "Funciones propias para la Temperatura máxima",
     ylab = "Temperatura máxima", xlab = "Días", xaxt = "n") # Grafica de las funciones propias

## [1] "done"
# Proporcion de varianza explicada por cada lambda
prop_fpca <- fpca_temp$varprop; prop_fpca
## [1] 0.50493900 0.13344385 0.09120628 0.06219653
# porcentaje de variabilidad explicada
cumsum(prop_fpca) # Con dos componentes ya se obtiene el 93.98% de variabilidad explicada
## [1] 0.5049390 0.6383828 0.7295891 0.7917857

Note que el 1er componente ya explica más del 50% de la variabilidad total. La proporción de varianza explicada muestra que con 4 componentes principales (funciones propias) se alcanza un 79.1% de variabilidad explicada.

Scores

Dado que los scores son escalares, estudiamos su estructura de covarianza estimando los semivariogramas empíricos.

score1 <- as.matrix(fpca_temp$scores[,1])
score2 <- as.matrix(fpca_temp$scores[,2])
score3 <- as.matrix(fpca_temp$scores[,3])
score4 <- as.matrix(fpca_temp$scores[,4])
vari_score1 <- variog(as.geodata(cbind(score1, utm_coords_df), data.col = 1), max.dist = 100000)
## variog: computing omnidirectional variogram
vari_score2 <- variog(as.geodata(cbind(score2, utm_coords_df), data.col = 1), max.dist = 100000)
## variog: computing omnidirectional variogram
vari_score3 <- variog(as.geodata(cbind(score3, utm_coords_df), data.col = 1), max.dist = 100000)
## variog: computing omnidirectional variogram
vari_score4 <- variog(as.geodata(cbind(score4, utm_coords_df), data.col = 1), max.dist = 100000)
## variog: computing omnidirectional variogram
plot(vari_score1,main = "Semivariograma empírico para el score1")

plot(vari_score2,main = "Semivariograma empírico para el score2")

plot(vari_score3,main = "Semivariograma empírico para el score3")

plot(vari_score4,main = "Semivariograma empírico para el score4")

Kriging funcional

colnames(utm_coords_df) <- c("X", "Y")

sfd_temp <- SpatFD(base_TMAX_filled[,-1], 
                   coords = utm_coords_df, 
                   basis = "Fourier", nbasis = 40,
                   lambda = 0.0005, nharm = 4)

plot(sfd_temp$`base_TMAX_filled[, -1]`$data_fd, lwd = 2)

## [1] "done"

Kriging scores y lambda

Mediante la observación de los semivariogramas empiricos es que se establecen los valores iniciales de la silla y el rango espacial para los modelos teóricos de semivarianza que se muestran a continuación

colnames(utm_coords_df) <- c("Easting", "Northing")
pred_grid <- expand.grid(
  Easting  = seq(min(utm_coords_df$Easting),  max(utm_coords_df$Easting),  length.out = 70),
  Northing = seq(min(utm_coords_df$Northing), max(utm_coords_df$Northing), length.out = 70))

modelos <- list(vgm(psill = 1237.49, "Exp", range = 38735.74, nugget =  0),
                vgm(psill = 395.48, "Exp", range = 85218.6, nugget = 0),
                vgm(psill = 192.49, "Exp", range = 15494.28, nugget =  0),
                vgm(psill = 200.79, "Exp", range = 100712.9, nugget =  0))



fkrig_lambdas <- KS_scores_lambdas(SFD = sfd_temp, newcoords = pred_grid, model = modelos, method = "both")
## Using first variable by default
## Using fill.all = TRUE by default
## [using simple kriging]
## [using simple kriging]
## [using simple kriging]
## [using simple kriging]

Recordemos que el kriging lambda lo que hace es estimar los coeficientes de la expresión \[\mathcal{\tilde{X}}_{\mathbf{s}_0}(t)=\sum_{i=1}^n \lambda_i \mathcal{X}_{\mathbf{s}_i}(t)\]

Por otro lado, el kriging scores dice que la predicción de la curva \(\mathcal{X}_{\mathbf{s}_0}(t)\) es \[\mathcal{X}^*_{\mathbf{s}_0}(t)=\boldsymbol{\xi}^T(t)\mathbf{f}_{\mathbf{s}_0}^*\]

con \(\boldsymbol{\xi}^T(t)\) como el vector de las primeras k funciones propias. En este caso solo se debe encontrar la predicción \(\mathbf{f}_{\mathbf{s}_0}^*\) mediante el método de cokriging.

ggplot_KS(fkrig_lambdas, show.varpred = FALSE,
              main = "Plot 1 - Using Scores",
              main2 = "Plot 2 - Using Lambda",
              ylab = "Temperatura maxima diaria")
## Warning: `aes_()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`
## ℹ The deprecated feature was likely used in the SpatFD package.
##   Please report the issue to the authors.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## [[1]]

## 
## [[2]]

Karhunen-Loeve

Con la descomposición Karhunen-Loeve se pueden reconstuir las curvas usando las funciones propias y los scores. Ya que el kriging funcional nos devuelve las predicciones de los scores, es posible encontrar las predicciones de las curvas en las ubicaciones espaciales deseadas. Estas se guardan en la matriz A.

scores <- fkrig_lambdas$KS_scores$scores_pred
head(scores)
##          V1         V2        V3          V4
## 1  1.359761 -4.2746340 -31.67806 -1.28995947
## 2  5.886724 -0.9280889 -26.10608 -0.87882482
## 3 10.535773  2.6309573 -21.02176 -0.42704176
## 4 15.260706  6.3777163 -16.36813  0.06246341
## 5 19.915840 10.2257201 -12.12432  0.57519941
## 6 24.062046 13.9146091  -8.37538  1.06524624
names(scores) <- c("score1", "score2", "score3", "score4")

mu_temp <- eval.fd(1:n_puntos, mean_temp)


eigen_fun1 <- eval.fd(1:n_puntos, fkrig_lambdas$SFD$`base_TMAX_filled[, -1]`$fpca$harmonics)[,1] 
eigen_fun2 <- eval.fd(1:n_puntos, fkrig_lambdas$SFD$`base_TMAX_filled[, -1]`$fpca$harmonics)[,2] 
eigen_fun3 <- eval.fd(1:n_puntos, fkrig_lambdas$SFD$`base_TMAX_filled[, -1]`$fpca$harmonics)[,3] 
eigen_fun4 <- eval.fd(1:n_puntos, fkrig_lambdas$SFD$`base_TMAX_filled[, -1]`$fpca$harmonics)[,4] 


A <- NULL
for (i in 1:nrow(scores)) {
  KL_decom <- mu_temp + scores[i,1]*eigen_fun1 + scores[i,2]*eigen_fun2 + 
  scores[i,3]*eigen_fun3 + scores[i,4]*eigen_fun4
  
  A <- cbind(A, KL_decom)
}
dim(A)
## [1] 1146 4900

Mapas finales

Con estos se construyen los mapas de predicción para tiempos específicos. En este caso se obta por los tiempos 100, 300, 600 y 900.

tiempos_mapa <- c(100, 300, 600, 900)

colnames(A) <- paste0("s", 1:ncol(A))
A_sel <- A[tiempos_mapa, ] 

pred_plot_df <- A_sel %>%
  as.data.frame() %>%
  mutate(time = tiempos_mapa) %>%
  pivot_longer(
    cols = starts_with("s"), 
    names_to = "id",
    values_to = "value"
  ) %>%
  mutate(id = as.integer(gsub("s", "", id)))

pred_plot_df <- pred_plot_df %>%
  left_join(
    pred_grid %>% mutate(id = row_number()),
    by = "id"
  )

colnames(utm_coords_df) <- c("Easting", "Northing")


ggplot() +
  geom_sf(data = yucatan_region_utm, fill = NA, color = "black") +
  geom_tile(data = pred_plot_df,
            aes(x = Easting, y = Northing, fill = value)) +
  geom_point(data = utm_coords_df,
             aes(x = Easting, y = Northing),
             color = "black", size = 1) +
  scale_fill_viridis_c(name = "Predicción") +
  coord_sf(xlim = c(120000, 350000),
             ylim = c(2180000, 2300000)) +
  facet_wrap(~ time) +
  labs(
    title = "Kriging Funcional – Predicciones en 4 tiempos",
    x = "Easting (m)",
    y = "Northing (m)"
  ) +
  theme_minimal()