Kriging Funcional

Se realizará el análisis de los datos mediante kriging funcional usando los scores para los datos correspondientes a las 24 horas de los primeros 7 días del mes de enero del año 2018 en Seúl. Este análisis se llevará a cabo para la variable correspondiente a PM2.5.

Recordemos que las partículas en suspensión (PM2.5) son partículas en el aire que son menores a una medida de 2.5 micras. En gran medida se usan ya que están relacionadas con el aumento de emisiones de la quema de combustibles tipo Diesel. Además, al ser considerablemente pequeñas (100 veces menor en tamaño a un cabello humano) tienen una gran capacidad de alojarse e impactar el sistema respiratorio. Pueden estar conformadas por polvo, cenizas, hollín, partículas metálicas, cemento, polen, etc. Tienen efectos en la función pulmonar y recientemente se relacionan con problemas de tipo cardiovascular. Producen latidos irregulares, asma grave, infartos de miocardio no mortales. Estas partículas también contribuyen de manera negativa al medio ambiente, pues pueden hacer que los lagos se vuelvan ácidos, una disminución de nutrientes en los suelos perjudicando así la agricultura en las regiones.

El aumento de este material particulado es provocado también por fenómenos meteorológicos. En Corea del Sur, se pueden ver aumentos de dicho contaminante por el denominado “polvo amarillo” proveniente de los desiertos de China y Mongolia, lo cual ocasiona que las mediciones de estas variables se vean considerablemente más grandes.

Con una mayor cantidad de datos medidos a través del tiempo en cada una de las 25 estaciones de medición en Seúl haremos uso de estos registros para representar dichos datos a través de curvas.

Tomaremos en cuenta las variables funcionales \(\chi_1^2,...,\chi_n^2\) y nuestros datos funcionales de nuestra base de datos serán una realización \(x_1,x_2,...,x_n\) de dichas variables funcionales.

Tenemos así, 168 datos en cada estación de medida con los cuales haremos el correspondiente suavizamiento para trabajar con nuestras funciones.

Procedemos a ver los datos

library(readxl) library(rainbow) library(fda) library(sf) library(ggplot2) library(plotly) library(sp) library(geoR)

Datos<-read_excel("BaseInforme2.xlsx",sheet = "Hoja2")
head(Datos)
## # A tibble: 6 x 13
##   `Measurement date`    Dia  Hora `Station code` Address      Latitude Longitude
##   <dttm>              <dbl> <dbl>          <dbl> <chr>           <dbl>     <dbl>
## 1 2018-01-01 00:00:00     1     0            101 19, Jong-ro~     37.6      127.
## 2 2018-01-01 01:00:00     1     1            101 19, Jong-ro~     37.6      127.
## 3 2018-01-01 02:00:00     1     2            101 19, Jong-ro~     37.6      127.
## 4 2018-01-01 03:00:00     1     3            101 19, Jong-ro~     37.6      127.
## 5 2018-01-01 04:00:00     1     4            101 19, Jong-ro~     37.6      127.
## 6 2018-01-01 05:00:00     1     5            101 19, Jong-ro~     37.6      127.
## # ... with 6 more variables: SO2 <dbl>, NO2 <dbl>, O3 <dbl>, CO <dbl>,
## #   PM10 <dbl>, PM2.5 <dbl>

Representamos graficamente los diferentes datos medidos a través del tiempo en las 25 estaciones

par(mfrow=c(1,1))
Horas<-1:168

PM2.5<-matrix(Datos$PM2.5,nrow = 168,25)
colnames(PM2.5)<-c(1:25)
plot(Horas, PM2.5[,1], type="b",ylim=c(0,max(PM2.5)), xlab="Hora", ylab="PM2.5")
for(i in 2:ncol(PM2.5))
  lines(Horas, PM2.5[,i], type="b", col=i)

Procedemos con la identificación de curvas atípicas

PMTs = sfts(ts(as.numeric(PM2.5), start = c(1,1),frequency = 168), xname = "Hora",
            yname = "PM2.5 oC")
plot(PMTs)

par(mfrow=c(1,2))
fboxplot(data= PMTs, plot.type = "functional", type = "bag", projmethod = "PCAproj")
## Warning in fts(x, notchlow): Please assign column name for the data matrix.
## Warning in fts(x, notchupper): Please assign column name for the data matrix.
## Warning in fts(x, centercurve): Please assign column name for the data matrix.
fboxplot(data= PMTs, plot.type = "bivariate", type = "bag", projmethod = "PCAproj",
         ylim = c(-100,100), xlim = c(-70,60))

par(mfrow=c(1,2))
fboxplot(data= PMTs, plot.type = "functional", type = "hdr", projmethod = "PCAproj")
## Warning in fts(x, centercurve): Please assign column name for the data matrix.
fboxplot(data= PMTs, plot.type = "bivariate", type = "hdr", projmethod = "PCAproj",
         ylim = c(-100,100), xlim = c(-70,60))

Donde podemos observar que las curvas correspondientes a las estaciones 5,12,25 presentan comportamientos atípicos y en el caso de la estación 5 vemos que para todo tiempo registró un valor de cero, con lo cuál podría pensarse en que dicha estación pudo sufrir un tipo de daño o simplemente no fue utilizada. Con lo cuál, procedemos a retirar dicha información y trabajaremos con las 24 estaciones restantes.

Datas<-read_excel("BaseInforme2.xlsx",sheet = "Hoja4")

par(mfrow=c(1,1))
PM2.5b<-matrix(Datas$PM2.5,nrow = 168,24)
colnames(PM2.5b)<-c(1:24)
plot(Horas, PM2.5b[,1], type="b",ylim=c(0,max(PM2.5b)), xlab="Hora", ylab="PM2.5")
for(i in 2:ncol(PM2.5b))
  lines(Horas, PM2.5b[,i], type="b", col=i)

Haremos uso del método no paramétrico de Bsplines para el correspondiente suavizamiento de las diferentes curvas. En este método consideramos que la curva está generada por distintos vértices y cada vértice va asociado con una función. Con lo cuál, su expresión viene dada por \[ P(t)= \sum_{i=0}^n P_iN_{i,j}(t)\] donde \(P_i\) son los \(n+1\) vértices del polígono. En este caso usaremos 5 funciones en la base.

nbasis <-5
hourange <- c(1,nrow(PM2.5b))
lambda=0.1
harmaccelLfd <- vec2Lfd(c(1,10), hourange)
hourbasis_Bsplines <- create.bspline.basis(hourange,nbasis)
PM2.5_fdPar_Bspline<-fdPar(fdobj=hourbasis_Bsplines,Lfdobj=harmaccelLfd,lambda)
PM2.5_fd_Bspline <- smooth.basis(argvals=1:nrow(PM2.5b),PM2.5b,PM2.5_fdPar_Bspline)
PM2.5_fd_Bspl=PM2.5_fd_Bspline$fd
par(mfrow=c(1,1))
plot(PM2.5_fd_Bspl)
## [1] "done"
lines(PM2.5_fd_Bspl,col=rainbow(10),lwd=2,lty=1)

Realizamos una reducción de la dimensionalidad por medio del método de componentes principales y de aquí, nos quedaremos con los scores.

PCA=pca.fd(PM2.5_fd_Bspl,centerfns=T)
PCA$varprop
## [1] 0.6755073 0.1812974
head(PCA$scores)
##              [,1]        [,2]
## [1,] -34.81538582  -0.5421375
## [2,] -43.79978334  -3.8870126
## [3,] -37.73866518  -1.9606849
## [4,] -37.04741081 -16.9722118
## [5,]   0.01582169  -3.3056796
## [6,]   7.65759520  -5.7427368

Podemos observar que la columna correspondiente a la primera componenete principal capta un \(68\%\) de la variabilidad en las curvas y la segunda está captando un \(18\%\).

Continuaremos el análisis con esta primera componente.

Datos2b<-read_excel("BaseInforme2.xlsx",sheet = "Hoja5")
Datos2b<-data.frame(Datos2b,PCA$scores[,1])
names(Datos2b)
## [1] "Station.code"    "Latitude"        "Longitude"       "PCA.scores...1."
colnames(Datos2b)<-c("Station.code","Latitude","Longitude","Scores1")
head(Datos2b)
##   Station.code Latitude Longitude      Scores1
## 1          101 37.57202  127.0050 -34.81538582
## 2          102 37.56426  126.9747 -43.79978334
## 3          103 37.54003  127.0049 -37.73866518
## 4          104 37.60982  126.9348 -37.04741081
## 5          106 37.55558  126.9056   0.01582169
## 6          107 37.54186  127.0497   7.65759520

Con lo cual, consideramos el campo escalar con los scores

shp_map <- read_sf("seoul_municipalities.shp")
shp_map1<-shp_map[5]


map_variable<- function(datos, variable, map) {
  
  plot1 <- ggplot() +
    geom_sf(data = map, size = 0.3) +
    geom_point(data = datos, aes_string(x = "Longitude", y = "Latitude",
                                        color = variable)) +
    scale_color_viridis_c(direction = -1) +
    labs(
      x = "",
      y = "",
      title = paste(variable, "map Seoul")
    ) +
    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(Datos2b,
                    "Scores1",
                    shp_map1)

ggplotly(pl1)

Y transformamos las coordenadas al plano XY

###Transformamos las coordenadas al plano XY
par(mfrow=c(1,1))

deci_coordb = SpatialPoints(cbind(Datos2b$Longitude,
                                  Datos2b$Latitude),
                            proj4string = CRS("+proj=longlat"))
class(deci_coordb)
## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"
plot(deci_coordb, axes = TRUE, main = "Lat-Long Coordinates")

CRS_UTM_S = (c("+proj=utm +zone=52 +datum=WGS84 
                  +units=m +no_defs"))

utm_coordb = spTransform(deci_coordb, CRS(CRS_UTM_S))
utm_coord_dfb = as.data.frame(utm_coordb)
plot(utm_coordb, axes = TRUE, main = "UTM Coordinates", col = "red")

Datos2b<-data.frame(Datos2b,utm_coord_dfb)

Observamos los semivariagramas en diferentes direcciones

### Isotropia 

names(Datos2b)
## [1] "Station.code" "Latitude"     "Longitude"    "Scores1"      "coords.x1"   
## [6] "coords.x2"
c<-as.geodata(Datos2b[,c(5,6,4)])
str(c)
## List of 2
##  $ coords.x1                  , coords.x2                  : num [1:24, 1:2] 323823 321125 323733 317719 315002 ...
##   ..- attr(*, "dimnames")=List of 2
##   .. ..$ : NULL
##   .. ..$ : chr [1:2] "coords.x1" "coords.x2"
##  $ data                       : num [1:24] -34.8154 -43.7998 -37.7387 -37.0474 0.0158 ...
##  - attr(*, "class")= chr "geodata"
df_reg=c
df_reg$data=Datos2b$Scores

vari_0=variog(df_reg, estimator.type = "modulus", dir=0)
## variog: computing variogram for direction = 0 degrees (0 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_45=variog(df_reg, estimator.type = "modulus", dir=pi/4)
## variog: computing variogram for direction = 45 degrees (0.785 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_90=variog(df_reg, estimator.type = "modulus", dir=pi/2)
## variog: computing variogram for direction = 90 degrees (1.571 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
vari_135=variog(df_reg, estimator.type = "modulus", dir=3*pi/4)
## variog: computing variogram for direction = 135 degrees (2.356 radians)
##         tolerance angle = 22.5 degrees (0.393 radians)
par(mfrow=c(2,2))
plot(vari_0,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*0)))
plot(vari_90,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*90)))


plot(vari_45,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*45)))

plot(vari_135,main="",xlab="",ylab="Semivarianza direccional")
title(expression(hat(gamma)(h)*paste(", "*phi*"="*135)),xlab="h",ylab="Semivarianza direccional")

Ahora procederemos a realizar el adecuado ajuste del modelo de semivarianza

###Estimación Semivariograma 

coordinates(Datos2b) <-  ~ coords.x1 + coords.x2

library(gstat)
g_obj <- gstat(id = "Scores1" ,
               formula = Scores1 ~ 1,
               data = Datos2b)
var_s <- variogram(g_obj,cutoff=15000)
var_s <- var_s[-1, ]  
plot(var_s, pl = T)

show.vgms()

vgm_model1 <- vgm("Exp",
                  psill = 1200,
                  range = 7000)
plot(var_s, vgm_model1, pl = T)

par(mfrow=c(1,3))

fitted_vgm1 <- fit.variogram(object = var_s,
                             model = vgm_model1)
plot(var_s, fitted_vgm1, pl = T)

vgm_model2 <- vgm("Gau",
                  psill = 1300,
                  range = 6000)
plot(var_s, vgm_model2, pl = T)

fitted_vgm2 <- fit.variogram(object = var_s,
                             model = vgm_model2)
plot(var_s, fitted_vgm2, pl = T)

vgm_model3 <- vgm("Wav",
                  psill = 1200,
                  range = 8500)
plot(var_s, vgm_model3, pl = T)

fitted_vgm3 <- fit.variogram(object = var_s,
                             model = vgm_model3)
plot(var_s, fitted_vgm3, pl = T)

par(mfrow=c(1,3))

plot(var_s, fitted_vgm1, pl = T)

plot(var_s, fitted_vgm2, pl = T)

plot(var_s, fitted_vgm3, 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
var_s$fitted3 <- variogramLine(fitted_vgm3, dist_vector = var_s$dist)$gamma


cmr <- function(obs, fit) {
  return(mean((obs - fit)**2))
}

tabla <- data.frame(rbind(fitted_vgm1,
                          fitted_vgm2,
                          fitted_vgm3
),
CMR =
  c(cmr(var_s$gamma, var_s$fitted1),
    cmr(var_s$gamma, var_s$fitted2),
    cmr(var_s$gamma, var_s$fitted3)
  ))


pander::pander(tabla[c(1:3,10)])
model psill range CMR
Exp 1774 11507 37485
Gau 1046 4880 43838
Wav 1074 8500 49383

En donde el modelo de semivarianza Exponencial de silla 1774 y rango 11507 es el que presenta un menor error cuadrático medio y por ende trabajaremos con este.

Procedemos a realizar las predicciones basados en este modelo de semivarianza

g_obj <- gstat(id = "Scores1",
               formula = Scores1 ~ 1,
               model = fitted_vgm1, # Modelo ajustado 1
               data = Datos2b)


map <- readOGR("seoul_municipalities.shp")
## OGR data source with driver: ESRI Shapefile 
## Source: "C:\Users\AxRen\Documents\UNAL\Espacial\seoul_municipalities.shp", layer: "seoul_municipalities"
## with 25 features
## It has 6 fields
## Integer64 fields read as strings:  ESRI_PK
class(map)
## [1] "SpatialPolygonsDataFrame"
## attr(,"package")
## [1] "sp"
new <- sp::spsample(map, n = 20000, type = "regular")
## Warning in proj4string(obj): CRS object has comment, which is lost in output; in tests, see
## https://cran.r-project.org/web/packages/sp/vignettes/CRS_warnings.html
utm_coord2 = spTransform(new, CRS(CRS_UTM_S))
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]
plot1<-ggplot()+
  geom_sf(aes(col=Scores1.pred),data=predic)+
  scale_color_carto_c(palette='Prism', na.value = '#111111')

plot1  

plot12<-ggplot()+
  geom_sf(aes(col=Scores1.var),data=predic)+
  scale_fill_viridis_c(option = "magma",
                       direction = -1) +
  theme_void()


plot12

Y así, podemos observar los mapas de predicción de los scores y el respectivo mapa con su varianza asociada.

Con lo cuál, llevando estas predicciones de nuevo a los datos funcionales tenemos

coef_scores=(PCA$harmonics$coefs[,1]%*% t(Datos2b$Scores1))
vec<-PCA$meanfd$coefs
M<-NULL
MM<-NULL
for (i in 1:24) {
  M<-coef_scores[,i]+vec
  MM<-cbind(MM,M)
}

par(mfrow=c(1,1))

v1=(fd(MM, PCA$harmonics$basis))
plot(v1)

## [1] "done"
coef_scores2=PCA$harmonics$coefs[,1] %*%t(predic$Scores1.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"

En donde hacemos el producto de dichas predicciones con su respectiva función propia asociada y le sumamos la media a cada función. Permitiendo así visualizar el comportamiento de las diferentes funciones en los diferentes lugares donde se realiza la predicción.

Observemos algunas de estas predicciones

NN[,1:6]
##          mean     mean     mean     mean     mean     mean
## [1,] 22.56890 22.56669 22.56899 22.57128 22.57359 22.57595
## [2,] 15.94401 15.94143 15.94410 15.94678 15.94948 15.95223
## [3,] 21.08966 21.09016 21.08964 21.08912 21.08859 21.08806
## [4,] 36.27531 36.27029 36.27549 36.28070 36.28595 36.29130
## [5,] 34.62425 34.62204 34.62434 34.62663 34.62895 34.63131
NN[,1000:1005]
##          mean     mean     mean     mean     mean     mean
## [1,] 22.40058 22.41141 22.42435 22.43918 22.45565 22.47349
## [2,] 15.74764 15.76028 15.77537 15.79267 15.81189 15.83270
## [3,] 21.12790 21.12544 21.12250 21.11913 21.11539 21.11133
## [4,] 35.89340 35.91798 35.94733 35.98098 36.01835 36.05882
## [5,] 34.45580 34.46665 34.47959 34.49443 34.51092 34.52877
NN[,10000:10005]
##          mean     mean     mean     mean     mean     mean
## [1,] 23.17947 23.19983 23.22352 23.25353 23.30848 23.41818
## [2,] 16.65629 16.68005 16.70769 16.74270 16.80680 16.93478
## [3,] 20.95095 20.94632 20.94094 20.93412 20.92164 20.89672
## [4,] 37.66061 37.70681 37.76057 37.82866 37.95332 38.20224
## [5,] 35.23526 35.25564 35.27935 35.30938 35.36437 35.47416

MAPA DE PREDICCIONES

Aquí podemos visualizar el mapa de predicciones para el tiempo 100

hourange1<-seq(1:168)
Todas = eval.fd(hourange1,v2)
Todas1= Todas[100,]
head(Todas1)
##     mean     mean     mean     mean     mean     mean 
## 26.61450 26.61234 26.61458 26.61682 26.61908 26.62139
S<-data.frame(Todas1,predic$geometry)
S1<-st_as_sf(S)

plot14<-ggplot()+
  geom_sf(aes(col=Todas1),data=S1)+
  scale_color_carto_c(palette='Prism', na.value = '#111111')

plot14

Referencias