Kriging funcional con datos de la velocidad del viento a diferentes horas del día del mes de febrero de 2023 en el departamento de Caldas-Colombia
Jessica López Mejía y Nicolás Felipe Romero Solano
Departamento de Caldas
En el departamento de Caldas, Colombia, la velocidad del viento varía dependiendo de la zona geográfica y las condiciones climáticas. En general, Caldas no se caracteriza por ser una región con vientos extremadamente fuertes o constantes.
La velocidad promedio del viento en Caldas suele oscilar entre 3 y 6 metros por segundo (m/s), con algunas variaciones según la altitud y la ubicación geográfica específica. La velocidad media en los días que se tomaron los datos es de 0.99 m/s. Adicionalmente, en las zonas montañosas y elevadas, como la Sierra Nevada del Ruiz, es posible que se registren velocidades ligeramente más altas debido a los efectos del relieve.
Es importante tener en cuenta que los vientos en Caldas pueden experimentar cambios estacionales y diarios. Durante la temporada de lluvias, es posible que se registren vientos más fuertes asociados a tormentas y sistemas meteorológicos. Sin embargo, en comparación con otras regiones del país, Caldas no es conocido por ser propenso a eventos climáticos extremos relacionados con vientos intensos.
Kriging espacio tiempo con datos funcionales
Los datos funcionales son un tipo de datos en el que cada observación se representa como una función o curva en lugar de un valor único. El análisis de datos funcionales se basa en el estudio de la función que describe la variabilidad de un conjunto de datos en un espacio de muestras1. Esto permite aprovechar la estructura y dependencias presentes en los datos para realizar análisis y predicciones más precisas.
Análisis gráfico de los datos funcionales espaciales
Veamos un resumen de cada variable que compone la base de datos:
## Warning: package 'rgdal' was built under R version 4.2.3## Warning: package 'geoR' was built under R version 4.2.3## Warning: package 'raster' was built under R version 4.2.3## Warning: package 'interp' was built under R version 4.2.3## Warning: package 'akima' was built under R version 4.2.3## Warning: package 'sf' was built under R version 4.2.3## Warning: package 'gstat' was built under R version 4.2.3## Warning: package 'MASS' was built under R version 4.2.3## Warning: package 'lattice' was built under R version 4.2.3## Warning: package 'maptools' was built under R version 4.2.3## Warning: package 'plotly' was built under R version 4.2.3## Warning: package 'sqldf' was built under R version 4.2.3## Warning: package 'gsubfn' was built under R version 4.2.3## Warning: package 'proto' was built under R version 4.2.3## Warning: package 'RSQLite' was built under R version 4.2.3## Warning: package 'fda' was built under R version 4.2.3## Warning: package 'fds' was built under R version 4.2.3## Warning: package 'rainbow' was built under R version 4.2.3## Warning: package 'pcaPP' was built under R version 4.2.3## Warning: package 'RCurl' was built under R version 4.2.3## Warning: package 'deSolve' was built under R version 4.2.3## Warning: package 'fda.usc' was built under R version 4.2.3## Latitud Longitud Region Departamento
## Min. :4.986 Min. :-75.87 Length:2241 Length:2241
## 1st Qu.:5.082 1st Qu.:-75.65 Class :character Class :character
## Median :5.297 Median :-75.51 Mode :character Mode :character
## Mean :5.283 Mean :-75.44
## 3rd Qu.:5.424 3rd Qu.:-75.16
## Max. :5.610 Max. :-74.67
## Municipio Fecha Hora Velocidad_Viento
## Length:2241 Length:2241 Length:2241 Min. :0.0000
## Class :character Class :character Class :character 1st Qu.:0.6000
## Mode :character Mode :character Mode :character Median :1.0000
## Mean :0.9963
## 3rd Qu.:1.3000
## Max. :2.5000## [,1]
## Media 0.9962517
## Mediana 1.0000000
## Varianza 0.2212627
## Desviacion_Estandar 0.4703857
## Cv 47.2155517
## MIN 0.0000000
## MAX 2.5000000Veamos los municipios a los que se les hizo la medición de la velocidad del viento desde las 0:00 hasta las 23:00
## [1] "AGUADAS" "ANSERMA" "ARANZAZU" "BELALCAZAR" "CHINCHINA"
## [6] "FILADELFIA" "LA DORADA" "LA MERCED" "MANIZALES" "MANZANARES"
## [11] "MARMATO" "MARQUETALIA" "MARULANDA" "NEIRA" "NORCASIA"
## [16] "PACORA" "PALESTINA" "PENSILVANIA" "RIOSUCIO" "RISARALDA"
## [21] "SALAMINA" "SAMANA" "SAN JOSE" "SUPIA" "VICTORIA"
## [26] "VILLAMARIA" "VITERBO"Gráfico descriptivo de la variable velocidad del viento en el tiempo por cada ubicación espacial
## named integer(0)Método Bag
En este caso, se utilizará un gráfico de caja “bag”, que muestra la dispersión de las funciones en la dimensión funcional. Con el método de proyección “PCAproj”, que es una técnica de reducción de dimensionalidad basada en el análisis de componentes principales.
par(mfrow=c(1,1))
fboxplot(data= TempTs, plot.type = "bivariate", type = "bag", projmethod = "PCAproj",
ylim = c(-3.5,2), xlim = c(-3,3))
Hay curva atipica en el municipio 7
Metodo HDR
En este caso, se utiliza el tipo “hdr” (heterogeneidad y dispersión relativa), que muestra la heterogeneidad y la dispersión relativa de las funciones en la dimensión funcional. Este tipo de gráfico resalta las diferencias en la forma, amplitud y posición de las funciones.
par(mfrow=c(1,2))
fboxplot(data= TempTs, plot.type = "functional", type = "hdr", projmethod = "PCAproj")## Warning in fts(x, centercurve): Please assign column name for the data matrix.fboxplot(data= TempTs, plot.type = "bivariate", type = "hdr", projmethod = "PCAproj")
Vemos que se detecta una curva atípica en el municipio 7 y 22.
Método BD2
El método BD2 se basa en la idea de que los valores atípicos pueden afectar las derivadas de las funciones en lugar de las funciones originales. En lugar de trabajar directamente con las funciones originales, se calculan las derivadas de las funciones y se aplica el enfoque de bagging para construir múltiples modelos predictivos de las derivadas.
par(mfrow=c(1,1))
fbplot(fit= C, method = "BD2", xlab = "Hora", ylab = "Velocidad Viento", outliercol = 2)
## $depth
## [1] 0.08760684 0.07763533 0.11680912 0.07407407 0.08689459 0.18803419
## [7] 0.07407407 0.13817664 0.07407407 0.11396011 0.07407407 0.10541311
## [13] 0.07407407 0.09401709 0.11111111 0.14886040 0.13888889 0.08903134
## [19] 0.07763533 0.07407407 0.09259259 0.07407407 0.07407407 0.11111111
## [25] 0.07407407 0.07549858 0.08689459
##
## $outpoint
## integer(0)
##
## $medcurve
## [1] 6Este método detecta outliers.
Al comparar la detección de datos atípicos utilizando las tres herramientas mencionadas (bagplot, hdr plot y fbplot) en la muestra de las 27 curvas de velocidad del viento, se evidencia que el hdr plot resulta ser la herramienta mas estricta o exigente en cuanto a que detecta mejor los outliers.
Otras medidas resumen
- El parámetro “trim” especifica la fracción de datos que se eliminaran.
- Alpha representa el nivel de confianza utilizado en ciertos cálculos de profundidad funcional.
- Trimmed mean, es una medida de tendencia central que se calcula eliminando un cierto porcentaje de valores extremos o atípicos en los extremos de un conjunto de datos antes de calcular la media.
par(mfrow=c(1,1))
MedFM = fdepth(data = TempTs, type = "FM", trim = 0.1)
plot(MedFM)
Este tipo de profundidad funcional utiliza la distancia de Mahalanobis para calcular la profundidad de un punto de datos en relación con la distribución multivariante general de los datos. Los puntos que tienen una mayor profundidad funcional están más cerca del centro de la distribución y se consideran más “típicos”.
MedRP = fdepth(data = TempTs, type = "RP", trim = 0.1)
plot(MedRP)
Esta medida de profundidad funcional se basa en proyectar los datos en una dirección aleatoria y evaluar la posición relativa de los puntos de datos en función de la dispersión de las proyecciones.
MedRPD = fdepth(data = TempTs, type = "RPD", trim = 0.1)## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.
## Warning in fts(1:dim(vproject)[2], t(vproject)): Please assign column name for
## the data matrix.plot(MedRPD)
Similar al RP depth, pero específico para datos funcionales.
Parámetros estimados - Moda
MODA = fdepth(data = TempTs, type = "mode", trim = 0.1)
plot(MODA)
Esta medida de profundidad se basa en el concepto de encontrar los modos de las curvas de contorno generadas a partir de los datos. Los puntos que se encuentran en las regiones más densas de los datos tienen una mayor profundidad funcional.
RADIUS = fdepth(data = TempTs, type = "radius", trim = 0.1, alpha=0.05, weight = "hard")
plot(RADIUS)
RADIUS = fdepth(data = TempTs, type = "radius", trim = 0.1, alpha=0.05, weight = "soft")
plot(RADIUS)
Esta medida calcula la profundidad de un punto en relación con el radio de un contenedor convexo que contiene la mayoría de los datos. Los puntos que se encuentran mas cerca del centro del contenedor convexo tienen una mayor profundidad funcional.
Finalmente, a través de procedimientos bootstrap es posible estimar intervalos de confianza para las medidas de resumen.
Boostrap Funcional
Este método se basa en la idea de generar multiples muestras (con reemplazo) a partir de la muestra original y calcular las estimaciones de interes en cada una de estas muestras.
BSpl <- create.bspline.basis(rangeval=c(0,24),nbasis=8, norder=4)
Ftemp <- Data2fd(C2,basisobj=BSpl,argvals=seq( 0 , 24, length= 24 ))
plot(Ftemp, main="Datos suavizados de la Velocidad del viento B-splines(4,k=8)", xlab = "Horas", ylab = "Velocidad del Viento")## [1] "done"legend("topright", legend = colnames(C2), col = 2:ncol(C), lty = 1, cex = 0.3)
out.boot1=fdata.bootstrap(Ftemp,statistic=func.mean,nb=200,draw=TRUE)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.FM,nb=100,draw=TRUE)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.mode,nb=100,draw=TRUE, trim = 0.1)
out.boot2=fdata.bootstrap(Ftemp,statistic=func.trim.RP,nb=100,draw=TRUE, trim = 0.1)
ACPF
PCA para datos funcionales
El Análisis de Componentes Principales (PCA) es una técnica estadística utilizada para reducir la dimensionalidad de un conjunto de datos y encontrar la variabilidad en el mismo1. En el caso de los datos funcionales, el objetivo del análisis en componentes principales funcionales es reducir la dimensión del espacio, ya que las variables aleatorias funcionales se encuentran en espacios de dimensión infinita. Esto permite representar los datos de manera más compacta y eficiente, facilitando su análisis y modelado.
par(mfrow=c(2,3))
PCA = pca.fd(Ftemp, 6)
plot(PCA)
Tomaremos una varianza del 90%. Y así, con dos funciones propias se explica este porcentaje de la varianza.
Gráficos de los scores
par(mfrow=c(1,1))
plot(x = PCA$scores[,1], y = PCA$scores[,2])
text(x = PCA$scores[,1], y = PCA$scores[,2] , seq(14,16,1), cex = 1.5)
En la literatura existen varias funciones base que permiten desarrollar la expansión. La selección del tipo de funciones de la base depende de las características que cumpla que fenómeno de estudio. Por ejemplo, las series de Fourier se utilizan para funciones con comportamientos cíclicos o periódicos y las bases monomiales se utilizan cuando existen tendencias simples que pueden ser ajustadas mediante líneas rectas, polinomios cuadráticos, polinomios de orden superior, etc. Las bases más usuales son aquellas basadas en splines debido a que son más flexibles y se ajustan de mejor manera a diferentes comportamientos. Dentro de estas se encuentran B-splines, M-splines, I-splines, y funciones de potencia truncadas.
Elegir BASES
En general, el número de bases utilizado puede variar según varios factores, como la complejidad del problema, la cantidad de datos disponibles y las limitaciones computacionales. Aquí hay algunas pautas generales que puedes considerar al seleccionar el número de bases.
Complejidad del problema: Si el problema es complejo y se espera que requiera una representación más rica de las funciones, puede ser necesario aumentar el número de bases. Esto permitira capturar mejor las variaciones y las formas más complejas de las funciones que estás modelando.
Cantidad de datos: Si tienes una gran cantidad de datos disponibles, puedes permitirte utilizar un número mayor de bases. Sin embargo, ten en cuenta que un exceso de bases puede llevar a un sobre ajuste (overfitting) si no se tiene cuidado.
Limitaciones computacionales: Si tienes restricciones de recursos computacionales, como memoria o tiempo de cálculo limitados, debes considerar limitar el número de bases para evitar problemas de rendimiento.
Bases monomiales
Las bases monomiales requieren el dominio y el número de base. Por ejemplo, la base monomial con K=7 funciones de base definidas en el intervalo [0,3] se puede construir con:
bbasis_obj = create.monomial.basis(rangeval=c(0,24), nbasis = 8)
Ftemp <- Data2fd(argvals=seq( 0 , 24 , length= 24 ),C2, basisobj=bbasis_obj)
plot(Ftemp, main="Datos suavizados Velocidad_viento Monomiales (k=8)", xlab = "Horas", ylab = "Velocidad Viento")## [1] "done"legend("topright", legend = colnames(C2), col = 2:ncol(C2), lty = 1, cex = 0.3)
Bases de fourier
Para utilizarlas bases de Fourier se requiere que se defina el dominio, el periodo de oscilación y el número de funciones de base.
fbasis_obj <- create.fourier.basis(rangeval=c(0,24), nbasis=8, period = 24)
Ftemp <- Data2fd(seq(0,24,length= 24),C2, basisobj=fbasis_obj)
plot(Ftemp, main="Datos suavizados de la Velocidad del viento Fourier(k=8), periodo 24", xlab = "Horas", ylab = "Velocidad Viento")## [1] "done"legend("topright", legend = colnames(C2), col = 2:ncol(C2), lty = 1, cex = 0.3)
B-splines Las bases B-spline requieren el dominio, el número de funciones de la base y el orden.
bsbasis_obj <- create.bspline.basis(rangeval=c(0,24),nbasis=8, norder=4)
plot(bsbasis_obj)
Ftemp <- Data2fd(seq(0,24,length= 24),C2, basisobj=bsbasis_obj)
plot(Ftemp, main="Datos suavizados de la Velocidad del viento B-splines(4,k=8)", xlab = "Horas", ylab = "Velocidad Viento")## [1] "done"legend("topright", legend = colnames(C2), col = 2:ncol(C2), lty = 1, cex = 0.3)
Estimación parámetros
loglam = seq(0,2,0.001)
nlam = length(loglam)
dfsave = rep(NA,nlam)
gcvsave = rep(NA,nlam)
for (ilam in 1:nlam) {
lambda = loglam[ilam]
fdParobj = fdPar(bsbasis_obj, Lfdobj=NULL, lambda= lambda)
smoothlist = smooth.basis(seq(0,24,length= 24), C, fdParobj)
dfsave[ilam] = smoothlist$df
gcvsave[ilam] = sum(smoothlist$gcv)
}
par(mfrow=c(1,1))
plot(loglam,gcvsave,xlab=expression(lambda),ylab=expression(GCV(lambda)),
main="Parametros de suavizamiento versus GCV",type="b",cex=0.7)
Tomando como referencia los 27 municipios en la base de datos de interés, se evidencia que el punto que minimiza la estadística GCV se encuentra alrededor de v=0.001 valores superiores antes de 0.6 y después de 0.7. Se destaca que dicho parámetro minimiza la suma del criterio GCV de todas las curvas que se encuentran dentro de la muestra, es decir, el valor de dicho parámetro resulta ser el óptimo en relación con el proceso de suavizamiento de todas las curvas dentro de la muestra.
Estadística descriptiva
Ftemp es usado con b-splines
meanfdh <- mean.fd(Ftemp)
varfdh <-var.fd(Ftemp)
stdvfdh <- stddev.fd(Ftemp)
plot(Ftemp,col=8, lty=1, xlab = "Hora", ylab = "Velocidad del viento")## [1] "done"lines(meanfdh,col=2,lwd=2)
par(mfrow=c(1,2))
plot(varfdh, main=" Superficie de varianza", xlab = "v", ylab = "s")
plot(stdvfdh, main="Desviacion estándar", xlab = "Hora", ylab = "Velocidad del viento")
## [1] "done"ACP sin las curvas atípicas 7 y 22
CC <- C
CC <- CC[,-7]
CC <- CC[,-21]
CC2 <-CC
Listado_Municpios2 <- Listado_Municpios
Listado_Municpios2 <- Listado_Municpios[-7]
Listado_Municpios2 <- Listado_Municpios2[-21]
colnames(CC2) <-c(Listado_Municpios2)
par(mfrow=c(1,1))
horas <-0:23
plot(horas, CC[,1], type="b", ylim=c(0, 2.5), xlab="Horas", ylab="Velocidad Viento",col=1)
for(i in 2:ncol(CC)){
lines(horas, CC[,i], type="b", col=i)}
legend("topright", legend = colnames(CC2), col = 2:ncol(CC2), lty = 1, cex = 0.3)
BSpl <- create.bspline.basis(rangeval=c(0,24),nbasis=8, norder=4)
FtempCC <- Data2fd(CC2,argvals=seq( 0 , 24 , length= 24 ), basisobj=BSpl)
plot(FtempCC, main="Datos suavizados de la Velocidad del viento B-splines(4,k=8)", xlab = "Horas", ylab = "Velocidad Viento")## [1] "done"legend("topright", legend = colnames(CC2), col = 2:ncol(CC2), lty = 1, cex = 0.3)
PCA = pca.fd(FtempCC, 6)
par(mfrow=c(2,3))
plot(PCA)
(PCA$varprop[1]+PCA$varprop[2])*100## [1] 90.16705Tomaremos una varianza del 90%. Con dos funciones propias se explica el 90,1% de la varianza.
Gráfico de los scores sin las curvas atípicas
par(mfrow=c(1,1))
plot(x = PCA$scores[,1], y = PCA$scores[,2], cex=0.001,pch=16)
text(x = PCA$scores[,1], y = PCA$scores[,2] , seq(14,16,1), cex = 1)
Organizando la tabla de datos
CALDAS3 <- sqldf("select Municipio, Latitud, Longitud
from CALDAS
group by Municipio ")
CALDAS3 <- CALDAS3[-7,]
CALDAS3 <- CALDAS3[-21,]
CALDAS3## Municipio Latitud Longitud
## 1 AGUADAS 5.610248 -75.45487
## 2 ANSERMA 5.236476 -75.78434
## 3 ARANZAZU 5.271196 -75.49129
## 4 BELALCAZAR 4.993821 -75.81190
## 5 CHINCHINA 4.985813 -75.60748
## 6 FILADELFIA 5.297098 -75.56247
## 8 LA MERCED 5.396372 -75.54657
## 9 MANIZALES 5.057687 -75.49105
## 10 MANZANARES 5.255708 -75.15284
## 11 MARMATO 5.474002 -75.60025
## 12 MARQUETALIA 5.297525 -75.05309
## 13 MARULANDA 5.284335 -75.25940
## 14 NEIRA 5.166706 -75.51984
## 15 NORCASIA 5.574794 -74.88954
## 16 PACORA 5.527175 -75.45962
## 17 PALESTINA 5.017839 -75.62446
## 18 PENSILVANIA 5.383279 -75.16030
## 19 RIOSUCIO 5.423672 -75.70207
## 20 RISARALDA 5.164192 -75.76727
## 21 SALAMINA 5.403019 -75.48722
## 23 SAN JOSE 5.082311 -75.79206
## 24 SUPIA 5.446834 -75.64966
## 25 VICTORIA 5.317418 -74.91116
## 26 VILLAMARIA 5.043807 -75.51376
## 27 VITERBO 5.062663 -75.87061CALDAS2 <- sqldf("select Municipio, Latitud, Longitud
from CALDAS
group by Municipio ")
CALDAS2 <- CALDAS2[-7,]
CALDAS2 <- CALDAS2[-21,]
CALDAS2## Municipio Latitud Longitud
## 1 AGUADAS 5.610248 -75.45487
## 2 ANSERMA 5.236476 -75.78434
## 3 ARANZAZU 5.271196 -75.49129
## 4 BELALCAZAR 4.993821 -75.81190
## 5 CHINCHINA 4.985813 -75.60748
## 6 FILADELFIA 5.297098 -75.56247
## 8 LA MERCED 5.396372 -75.54657
## 9 MANIZALES 5.057687 -75.49105
## 10 MANZANARES 5.255708 -75.15284
## 11 MARMATO 5.474002 -75.60025
## 12 MARQUETALIA 5.297525 -75.05309
## 13 MARULANDA 5.284335 -75.25940
## 14 NEIRA 5.166706 -75.51984
## 15 NORCASIA 5.574794 -74.88954
## 16 PACORA 5.527175 -75.45962
## 17 PALESTINA 5.017839 -75.62446
## 18 PENSILVANIA 5.383279 -75.16030
## 19 RIOSUCIO 5.423672 -75.70207
## 20 RISARALDA 5.164192 -75.76727
## 21 SALAMINA 5.403019 -75.48722
## 23 SAN JOSE 5.082311 -75.79206
## 24 SUPIA 5.446834 -75.64966
## 25 VICTORIA 5.317418 -74.91116
## 26 VILLAMARIA 5.043807 -75.51376
## 27 VITERBO 5.062663 -75.87061CALDAS3<-data.frame(CALDAS3,PCA$scores[,1])
CALDAS2<-data.frame(CALDAS2,PCA$scores[,2])
shp_map <- read_sf("C:/Users/nicol/OneDrive/Documentos/U. Nacional/MATERIAS/Estadistica Espacial/Datos funcionales/CaldasM.shp")
shp_map1<-shp_map[5]
Cal <- shp_map1
ggplot() +
geom_sf(data=Cal)
names(CALDAS2)## [1] "Municipio" "Latitud" "Longitud" "PCA.scores...2."names(CALDAS3)## [1] "Municipio" "Latitud" "Longitud" "PCA.scores...1."colnames(CALDAS2)<-c("Municipio","Latitud","Longitud","Score2")
colnames(CALDAS3)<-c("Municipio","Latitud","Longitud","Score1")Mapa del departamento de Caldas con sus respectivos scores
library(rcartocolor)## Warning: package 'rcartocolor' was built under R version 4.2.3library(plotly)
map_variable<- function(datos, variable, map) {
plot1 <- ggplot() +
geom_sf(data = map, size = 0.3) +
geom_point(data = datos, aes_string(x = "Longitud", y = "Latitud",
color = variable)) +
scale_color_viridis_c(direction = -1) +
labs(
x = "",
y = "",
title = paste(variable, "Mapa Dept Caldas")
) +
theme_void() +
theme(
plot.title = element_text(size = 14,
face = "bold",
colour = "black"),
plot.margin = unit(c(1, 1, 1.5, 1.2), "cm"))
return(plot1)
}
pl1 <- map_variable(CALDAS3,
"Score1",
shp_map1)## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation ideoms with `aes()`ggplotly(pl1)pl2 <- map_variable(CALDAS2,
"Score2",
shp_map1)
ggplotly(pl2)Kriging Score 2
Organizando los datos…
## Longitud Latitud Score2
## 1 -75.45487 5.610248 -0.4540961
## 2 -75.78434 5.236476 0.6943075
## 3 -75.49129 5.271196 -0.3766098
## 4 -75.81190 4.993821 0.0131652
## 5 -75.60748 4.985813 -0.5679960
## 6 -75.56247 5.297098 -0.3079728
## [1] "Municipio" "Latitud" "Longitud" "Score2" "Longitud.1"
## [6] "Latitud.1"## List of 2
## $ Longitud.1 , Latitud.1 : num [1:25, 1:2] 449624 413082 445561 409993 432655 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "1" "2" "3" "4" ...
## .. ..$ : chr [1:2] "Longitud.1" "Latitud.1"
## $ data : num [1:25] -0.4541 0.6943 -0.3766 0.0132 -0.568 ...
## - attr(*, "class")= chr "geodata"Explorando anisotropía
## 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)
Por lo que de acuerdo a los resultados obtenidos mostrados en los diferentes plots, se puede decir que la velocidad del viento en el departamento de Caldas es Anisotrópica pues la estructura de los semivariogramas varía conforme cambia la dirección.
Explorando dos posibles modelos de semivariograma
par(mfrow=c(2,3))
coordinates(CALDAS2) <- ~ Longitud.1 + Latitud.1
g_obj <- gstat(id = "Score2" ,
formula = Score2 ~ 1,
data = CALDAS2)
var_s <- variogram(g_obj)
var_s <- var_s[-1, ]
plot(var_s, pl = T)
show.vgms()
vgm_model1 <- vgm("Wav",
psill = 0.35,
range = 31000)
plot(var_s, vgm_model1, pl = T)
vgm_model2 <- vgm("Gau",
psill = 0.35,
range = 22000)
plot(var_s, vgm_model2, pl = T)
par(mfrow=c(1,2))
fitted_vgm1 <- fit.variogram(object = var_s,
model = vgm_model1)## Warning in fit.variogram(object = var_s, model = vgm_model1): No convergence
## after 200 iterations: try different initial values?plot(var_s, fitted_vgm1, pl = T)
fitted_vgm2 <- fit.variogram(object = var_s,
model = vgm_model2)
plot(var_s, fitted_vgm2, pl = T)
var_s$fitted1 <- variogramLine(fitted_vgm1, dist_vector = var_s$dist)$gamma
var_s$fitted2 <- variogramLine(fitted_vgm2, dist_vector = var_s$dist)$gamma
cmr <- function(obs, fit) {
return(mean((obs - fit)**2))
}
tabla <- data.frame(rbind(fitted_vgm1,
fitted_vgm2
),
CMR =
c(cmr(var_s$gamma, var_s$fitted1),
cmr(var_s$gamma, var_s$fitted2)
))
pander::pander(tabla[c(1:3,10)])| model | psill | range | CMR |
|---|---|---|---|
| Wav | 0.2739 | 31000 | 0.00654 |
| Gau | 0.4496 | 32292 | 0.007263 |
De acuerdo a lo presentado en la tabla resumen escogimos el modelo Wave con silla 0.4184 y rango 31000 para modelar la semivarianza, ya que tiene un CMR mas bajo.
Predicción de los datos funcionales
g_obj <- gstat(id = "Score2",
formula = Score2 ~ 1,
model = fitted_vgm1, # Modelo ajustado 1
data = CALDAS2)
st_write(Cal, "Cal.shp",append=FALSE)## Deleting layer `Cal' using driver `ESRI Shapefile'
## Writing layer `Cal' to data source `Cal.shp' using driver `ESRI Shapefile'
## Writing 27 features with 12 fields and geometry type Multi Polygon.map <- readOGR("Cal.shp")## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\nicol\OneDrive\Documentos\U. Nacional\MATERIAS\Estadistica Espacial\Datos funcionales\Cal.shp", layer: "Cal"
## with 27 features
## It has 12 fieldsclass(map)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"new <- sp::spsample(map, n = 20000, type = "regular")
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
plot(utm_coord, axes = TRUE, main = "UTM Coordinates", col = "green",cex = 0.7)
utm_coord_df = as.data.frame(utm_coord)
utm_coord2 = spTransform(new, CRS(CRS_UTM_NY))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))
predic <- predict(g_obj, newdata = utm_coord2)## [using ordinary kriging]Mapa de predicciones
par(mfrow=c(1,1))
plot1 <- ggplot() +
geom_sf(aes(col = Score2.pred), data = predic) +
scale_color_gradientn(colors = c("red", "orange", "blue"), na.value = "#111111")
plot1 
Mapa de varianza de las predicciones
plot12<- ggplot() +
geom_sf(aes(col = Score2.var), data = predic) +
scale_color_gradientn(colors = c("red", "orange", "blue"), na.value = "#111111")
plot12 
Y así, podemos observar los mapas de predicción de los scores y el respectivo mapas con su varianza asociada, con lo cual podemos decir que las predicciones son bastante buenas dado que e general las respectivas varianzas están entre 0.04 y 0.16.
Ahora,llevando estas predicciones de nuevo a los datos funcionales tenemos:
coef_scores=(PCA$harmonics$coefs[,2]%*% t(CALDAS2$Score2))
vec<-PCA$meanfd$coefs
M<-NULL
MM<-NULL
for (i in 1:25) {
M<-coef_scores[,i]+vec
MM<-cbind(MM,M)
}
par(mfrow=c(1,1))
v1=(fd(MM, PCA$harmonics$basis))
plot(v1)## [1] "done"legend("topright", legend = colnames(CC2), col = 2:ncol(CC2), lty = 1, cex = 0.3)
coef_scores2=PCA$harmonics$coefs[,2] %*%t(predic$Score2.pred)
N<-NULL
NN<-NULL
for (j in 1:ncol(coef_scores2)) {
N<-coef_scores2[,j]+vec
NN<-cbind(NN,N)
}
v2=(fd(NN, PCA$harmonics$basis))
plot(v2)
## [1] "done"MAPA DE PREDICCIONES Aquí podemos visualizar el mapa de predicciones para la hora 7
Eval.fd — evalúa un objeto de datos funcionales en valores de argumento especificados, o evalúe una derivada o el resultado de aplicar un operador diferencial lineal al objeto funcional.
hourange1<-seq(1:24)
Todas = eval.fd(hourange1,v2)
Velocidad= Todas[7,]
head(Velocidad)## mean mean mean mean mean mean
## 0.9131519 0.8984195 0.8836731 0.9557453 0.9417752 0.9277081S<-data.frame(Velocidad,predic$geometry)
S1<-st_as_sf(S)
plot14<-ggplot()+
geom_sf(aes(col=Velocidad),data=S1)+
scale_color_carto_c(palette='Prism', na.value = '#111111')
plot14 
Validación cruzada
La validación cruzada en gstat está implementada en la función krige.cv() por lo que se emplea esa función para realizarla. Veamos:
CALDAS22 <- sqldf("select Municipio, Latitud, Longitud
from CALDAS
group by Municipio ")
CALDAS22 <- CALDAS22[-7,]
CALDAS22 <- CALDAS22[-21,]
CALDAS22<-data.frame(CALDAS22,PCA$scores[,2])
deci_coord = SpatialPoints(cbind(CALDAS22$Longitud, CALDAS22$Latitud),proj4string = CRS("+proj=longlat"))
Randina_geod <- data.frame(deci_coord, Score2)
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
utm_coord_df = as.data.frame(utm_coord)
CALDAS22<-data.frame(CALDAS22,utm_coord_df)
CALDAS_sf <- st_as_sf(CALDAS22, coords = c("coords.x1" , "coords.x2"), remove = FALSE, agr = "constant")
valcru <-krige.cv(formula=PCA.scores...2. ~1,locations=CALDAS_sf, model=vgm("Wav", psill = 0.27, range = 31000))Realizada la validación es posible calcular estadísticos resumen como el error medio (ME), error cuadrático medio (MSE), media del cociente de la desviación cuadrática (mean squared deviation ratio, MSDR), raíz del error cuadrático medio (RMSE), la RMSE relativa a la media de los observados (RMSE_rel)
ME <- mean(valcru$residual)
MSE <- mean(valcru$residual ^ 2)
MSDR <- mean(valcru$zscore ^ 2)
RMSE <- sqrt(mean(valcru$residual ^ 2))
RMSE_rel <-
sqrt(mean(valcru$residual ^ 2)) /
mean(valcru$observed) * 100
r <- cor(valcru$observed,
valcru$observed - valcru$residual)
tabla <- data.frame(ME, MSE, RMSE,
RMSE_rel, MSDR, r)
tabla## ME MSE RMSE RMSE_rel MSDR r
## 1 -0.2885333 4.881488 2.209409 -2.93085e+19 2746.976 0.54027En donde podemos decir que:
ME (Error medio): El valor ME representa el promedio de los residuos (diferencia entre los valores observados y predichos). Un ME de -0.2885333 indica que, en promedio, los residuos tienden a ser ligeramente negativos, lo que sugiere que el modelo tiende a subestimar ligeramente los valores observados.
MSE (Error cuadrático medio): El MSE es una medida del promedio de los errores al cuadrado. Un MSE de 4.881488 indica que el modelo tiene un error cuadrático medio de aproximadamente 4.88. Cuanto menor sea el valor del MSE, mejor será el ajuste del modelo.
RMSE (Raíz del error cuadrático medio): El RMSE es simplemente la raíz cuadrada del MSE. Un RMSE de 2.209409 indica que, en promedio, el modelo tiene un error absoluto de alrededor de 2.21 unidades en la misma escala que la variable observada. Cuanto menor sea el valor del RMSE, mejor será el ajuste del modelo.
RMSE_rel (RMSE relativo): El RMSE relativo se calcula dividiendo el RMSE por la media de los valores observados y multiplicándolo por 100 para expresarlo como un porcentaje. Un RMSE_rel de -2.93085e+19 sugiere que hay un problema con el cálculo de esta métrica, ya que el valor resultante es extremadamente grande en valor absoluto y negativo. Esto podría ser causado por una división por cero o por problemas en los datos utilizados.
MSDR (Dispersión media de los residuos estandarizados): El MSDR representa la media de los residuos estandarizados al cuadrado. Un MSDR de 2746.976 indica que hay una dispersión media relativamente alta en los residuos estandarizados. Cuanto menor sea el valor del MSDR, mejor será el ajuste del modelo.
r (Coeficiente de correlación): El valor r es el coeficiente de correlación entre los valores observados y los residuos (diferencia entre los valores observados y predichos). Un valor de r de 0.54027 indica una correlación positiva moderada entre los valores observados y los residuos. Esto sugiere que el modelo puede capturar parte de la variabilidad en los datos, pero aún existe una parte no explicada por el modelo.
Kriging Score 1
Organizando los datos…
## Longitud Latitud Score1
## 1 -75.45487 5.610248 -0.2627150
## 2 -75.78434 5.236476 -1.9686645
## 3 -75.49129 5.271196 1.5560901
## 4 -75.81190 4.993821 -2.7556619
## 5 -75.60748 4.985813 0.2347540
## 6 -75.56247 5.297098 0.4083734
## [1] "Municipio" "Latitud" "Longitud" "Score1" "Longitud.1"
## [6] "Latitud.1"## List of 2
## $ Longitud.1 , Latitud.1 : num [1:25, 1:2] 449624 413082 445561 409993 432655 ...
## ..- attr(*, "dimnames")=List of 2
## .. ..$ : chr [1:25] "1" "2" "3" "4" ...
## .. ..$ : chr [1:2] "Longitud.1" "Latitud.1"
## $ data : num [1:25] -0.263 -1.969 1.556 -2.756 0.235 ...
## - attr(*, "class")= chr "geodata"Explorando anisotropía
## 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)
Por lo que de acuerdo a los resultados obtenidos mostrados en los diferentes plots, se puede decir que la velocidad del viento en el departamento de Caldas es Anisotrópica pues la estructura de los semivariogramas varía conforme cambia la dirección.
Explorando dos posibles modelos de semivariograma
par(mfrow=c(2,3))
coordinates(CALDAS3) <- ~ Longitud.1 + Latitud.1
g_obj <- gstat(id = "Score1" ,
formula = Score1 ~ 1,
data = CALDAS3)
var_s <- variogram(g_obj)
var_s <- var_s[-1, ]
plot(var_s, pl = T)
show.vgms()
vgm_model1 <- vgm("Wav",
psill = 2.8,
range = 31000)
plot(var_s, vgm_model1, pl = T)
vgm_model2 <- vgm("Gau",
psill = 2.8,
range = 22000)
plot(var_s, vgm_model2, pl = T)
par(mfrow=c(1,2))
fitted_vgm1 <- fit.variogram(object = var_s,
model = vgm_model1)
plot(var_s, fitted_vgm1, pl = T)
fitted_vgm2 <- fit.variogram(object = var_s,
model = vgm_model2)
plot(var_s, fitted_vgm2, pl = T)
var_s$fitted1 <- variogramLine(fitted_vgm1, dist_vector = var_s$dist)$gamma
var_s$fitted2 <- variogramLine(fitted_vgm2, dist_vector = var_s$dist)$gamma
cmr <- function(obs, fit) {
return(mean((obs - fit)**2))
}
tabla <- data.frame(rbind(fitted_vgm1,
fitted_vgm2
),
CMR =
c(cmr(var_s$gamma, var_s$fitted1),
cmr(var_s$gamma, var_s$fitted2)
))
pander::pander(tabla[c(1:3,10)])| model | psill | range | CMR |
|---|---|---|---|
| Wav | 3.244 | 31000 | 0.6931 |
| Gau | 6.898 | 39496 | 0.8045 |
De acuerdo a lo presentado en la tabla resumen escogimos el modelo Wave con silla 3.379 y rango 31000 para modelar la semivarianza, ya que tiene un CMR mas bajo.
Predicción de los datos funcionales
g_obj <- gstat(id = "Score1",
formula = Score1 ~ 1,
model = fitted_vgm1, # Modelo ajustado 1
data = CALDAS3)
st_write(Cal, "Cal.shp",append=FALSE)## Deleting layer `Cal' using driver `ESRI Shapefile'
## Writing layer `Cal' to data source `Cal.shp' using driver `ESRI Shapefile'
## Writing 27 features with 12 fields and geometry type Multi Polygon.map <- readOGR("Cal.shp")## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## Warning: OGR support is provided by the sf and terra packages among others## OGR data source with driver: ESRI Shapefile
## Source: "C:\Users\nicol\OneDrive\Documentos\U. Nacional\MATERIAS\Estadistica Espacial\Datos funcionales\Cal.shp", layer: "Cal"
## with 27 features
## It has 12 fieldsclass(map)## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"new <- sp::spsample(map, n = 20000, type = "regular")
CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
par(mfrow=c(1,1))
plot(utm_coord, axes = TRUE, main = "UTM Coordinates", col = "green",cex = 0.7)
utm_coord_df = as.data.frame(utm_coord)
utm_coord2 = spTransform(new, CRS(CRS_UTM_NY))
utm_coord2<-as.data.frame(utm_coord2)
utm_coord2<-st_as_sf(utm_coord2,coords = c("x1","x2"))
predic <- predict(g_obj, newdata = utm_coord2)## [using ordinary kriging]Mapa de predicciones
Mapa de varianza de las predicciones
plot12<-ggplot() +
geom_sf(aes(col = Score1.var), data = predic) +
scale_color_gradientn(colors = c("red", "orange", "blue"), na.value = "#111111")
plot12 
Y así, podemos observar los mapas de predicción de los scores y el respectivo mapas con su varianza asociada, con lo cual podemos decir que las predicciones son bastante buenas dado que e general las respectivas varianzas están entre 0.05 y 1.5.
Ahora,llevando estas predicciones de nuevo a los datos funcionales tenemos:
coef_scores=(PCA$harmonics$coefs[,1]%*% t(CALDAS3$Score1))
vec<-PCA$meanfd$coefs
M<-NULL
MM<-NULL
for (i in 1:25) {
M<-coef_scores[,i]+vec
MM<-cbind(MM,M)
}
par(mfrow=c(1,1))
v1=(fd(MM, PCA$harmonics$basis))
plot(v1)## [1] "done"legend("topright", legend = colnames(CC2), col = 2:ncol(CC2), lty = 1, cex = 0.3)
coef_scores1=PCA$harmonics$coefs[,1] %*%t(predic$Score1.pred)
N<-NULL
NN<-NULL
for (j in 1:ncol(coef_scores1)) {
N<-coef_scores1[,j]+vec
NN<-cbind(NN,N)
}
v2=(fd(NN, PCA$harmonics$basis))
plot(v2)
## [1] "done"MAPA DE PREDICCIONES Aquí podemos visualizar el mapa de predicciones para la hora 7
Eval.fd -> evalúa un objeto de datos funcionales en valores de argumento especificados, o evalúe una derivada o el resultado de aplicar un operador diferencial lineal al objeto funcional.
hourange1<-seq(1:24)
Todas = eval.fd(hourange1,v2)
Velocidad= Todas[7,]
head(Velocidad)## mean mean mean mean mean mean
## 1.762511 1.808723 1.817420 1.822845 1.739968 1.759502S<-data.frame(Velocidad,predic$geometry)
S1<-st_as_sf(S)
plot14<-ggplot()+
geom_sf(aes(col=Velocidad),data=S1)+
scale_color_carto_c(palette='Prism', na.value = '#111111')
plot14 
Validación cruzada
La validación cruzada en gstat está implementada en la función krige.cv() por lo que se emplea esa función para realizarla. Veamos:
CALDAS33 <- sqldf("select Municipio, Latitud, Longitud
from CALDAS
group by Municipio ")
CALDAS33 <- CALDAS33[-7,]
CALDAS33 <- CALDAS33[-21,]
CALDAS33<-data.frame(CALDAS33,PCA$scores[,1])
deci_coord = SpatialPoints(cbind(CALDAS33$Longitud, CALDAS33$Latitud),proj4string = CRS("+proj=longlat"))
Randina_geod <- data.frame(deci_coord, Score2)
head(Randina_geod)## coords.x1 coords.x2 Score2
## 1 -75.45487 5.610248 -0.4540961
## 2 -75.78434 5.236476 0.6943075
## 3 -75.49129 5.271196 -0.3766098
## 4 -75.81190 4.993821 0.0131652
## 5 -75.60748 4.985813 -0.5679960
## 6 -75.56247 5.297098 -0.3079728CRS_UTM_NY = (c("+proj=utm +zone=18 +datum=NAD83
+units=m +no_defs"))
utm_coord = spTransform(deci_coord, CRS(CRS_UTM_NY))
utm_coord_df = as.data.frame(utm_coord)
CALDAS33<-data.frame(CALDAS33,utm_coord_df)
CALDAS_sf <- st_as_sf(CALDAS33, coords = c("coords.x1" , "coords.x2"), remove = FALSE, agr = "constant")
valcru <-krige.cv(formula=PCA.scores...1. ~1,locations=CALDAS_sf, model=vgm("Wav", psill = 2.8, range = 31000))Realizada la validación es posible calcular estadísticos resumen como el error medio (ME), error cuadrático medio (MSE), media del cociente de la desviación cuadrática (mean squared deviation ratio, MSDR), raíz del error cuadrático medio (RMSE), la RMSE relativa a la media de los observados (RMSE_rel)
ME <- mean(valcru$residual)
MSE <- mean(valcru$residual ^ 2)
MSDR <- mean(valcru$zscore ^ 2)
RMSE <- sqrt(mean(valcru$residual ^ 2))
RMSE_rel <-
sqrt(mean(valcru$residual ^ 2)) /
mean(valcru$observed) * 100
r <- cor(valcru$observed,
valcru$observed - valcru$residual)
tabla <- data.frame(ME, MSE, RMSE,
RMSE_rel, MSDR, r)
tabla## ME MSE RMSE RMSE_rel MSDR r
## 1 -0.4692697 18.02398 4.245466 2.823099e+18 666.3627 0.4542892En donde podemos decir que:
ME (Error medio): El valor ME representa el promedio de los residuos (diferencia entre los valores observados y predichos). Un ME de -0.4692697 indica que, en promedio, los residuos tienden a ser ligeramente negativos, lo que sugiere que el modelo tiende a subestimar ligeramente los valores observados.
MSE (Error cuadrático medio): El MSE es una medida del promedio de los errores al cuadrado. Un MSE de 18.02398 indica que el modelo tiene un error cuadrático medio de aproximadamente 18.02. Cuanto menor sea el valor del MSE, mejor será el ajuste del modelo.
RMSE (Raíz del error cuadrático medio): El RMSE es simplemente la raíz cuadrada del MSE. Un RMSE de 4.245466 indica que, en promedio, el modelo tiene un error absoluto de alrededor de 4.245 unidades en la misma escala que la variable observada. Cuanto menor sea el valor del RMSE, mejor será el ajuste del modelo.
RMSE_rel (RMSE relativo): El RMSE relativo se calcula dividiendo el RMSE por la media de los valores observados y multiplicándolo por 100 para expresarlo como un porcentaje. Un RMSE_rel de 2.823099e+18 sugiere que hay un problema con el cálculo de esta métrica, ya que el valor resultante es extremadamente grande en valor absoluto y negativo. Esto podría ser causado por una división por cero o por problemas en los datos utilizados.
MSDR (Dispersión media de los residuos estandarizados): El MSDR representa la media de los residuos estandarizados al cuadrado. Un MSDR de 666.3627 indica que hay una dispersión media relativamente alta en los residuos estandarizados. Cuanto menor sea el valor del MSDR, mejor será el ajuste del modelo.
r (Coeficiente de correlación): El valor r es el coeficiente de correlación entre los valores observados y los residuos (diferencia entre los valores observados y predichos). Un valor de r de 0.4542892 indica una correlación positiva moderada entre los valores observados y los residuos. Esto sugiere que el modelo puede capturar parte de la variabilidad en los datos, pero aún existe una parte no explicada por el modelo.
Conclusiones
Se puede observar que usando menos varianza explicado, el caso del score1, no tenemos una predicción con mas varianza. Mientras que al aumentar la varianza explicada, el caso del score2, se obtiene una predicción mas acertada.
En adición, se puede observar que los las funciones predichas son mucho mejor cuando tenemos mayor varianza explicada.
El kriging funcional es una herramienta de predicción de la Geoestadística funcional que aprovecha las dependencias espacio-temporales existentes en los registros medioambientales para realizar predicciones en localizaciones e instantes temporales no observados. En el caso del modelamiento de la velocidad del viento en el departamento de Caldas, el uso de datos funcionales y kriging funcional permitió obtener predicciones más precisas en todas las ubicaciones y en cada hora del día de una manera computacionalmente más eficiente.
Sobre las validaciones cruzadas se puede pensar en ajustar otro modelo que ayude a mejorar la predicciones. Este valor tan “bajo” de la validación cruzada se puede pensar que es debido a no haber hecho modelacion de la tendencia.