Bitácora
Marlon Ruiz
13/10/2021
Análisis Geoestadístico para la variable CS137
Apartir de los datos entregados, para la zona que cubre parcialmente al país de Ucrania, se carga el conjunto de datos el cual ya tiene una limpieza previa de los mismos, y una zona de exclusión donde se concentran la mayor cantidad de puntos de la muestra.
CS137.table <- read.csv("~/Geoestadistica/Datos/Chernobyl_ExclusionZone_SInRepetidos.txt")
## Crear DB de Rgeostats
library(RGeostats)## Warning: package 'RGeostats' was built under R version 4.0.0## Loading required package: Rcpp## Loading RGeostats - Version:12.0.0cs137.db <- db.create(CS137.table, flag.grid = F, ndim = 2, autoname = F)
z<-CS137.table$MAX_Cs137
#Cambio de rol de las variables
cs137.db <- db.locate(cs137.db, "MAX_Cs137", "z",1)
cs137.db <- db.locate(cs137.db, "Este", "x",1)Se puede observar al sur oeste la falta de datos, por lo tanto al realizar las interpolaciones será el área de estudio de menor interés, además, se observa que en la parte norte, aproximadamente en la coordenada E=720000, y N=5700000, se encuentran los valores más altos o de mayor concentración para la variable Cs137, esto corresponde a la cercania que tiene la muestra con el reactor nuclear.
# Mapeo de los puntos
plot(cs137.db, pch=21, bg="green", title="Muestra de Cs137")
Posterior a la carga del dataset y selección de la variable Cs137, la cual va a ser objeto de estudio del analisis geoestadístico univariado, se procede con el análisis exploratorio de datos espaciales
AEDE
Se procede con el análisis univariante para la variable , mediante estadísticas descriptivas y métodos gráficos como el histograma y diagrama de cajas y bigotes para observar el comportamiento de la variable objeto de estudio.
Medidas de tendencia central
Moda <- function(x) {
uni <- unique(x)
uni[which.max(tabulate(match(x, uni)))]
}
tabla.mtc= data.frame(stringsAsFactors=FALSE,
Media = mean(z),
Mediana = median(z),
Moda=Moda(z),
Minimo=min(z),
Maximo=max(z)
)Se puede observar en la siguiente tabla que los datos de la media y mediana se encuentran muy alejados, la media es superior a la mediana (correspondiente al segundo cuartil o Q50%), Media\(>>\)Mediana, por lo que podemos inferir dos cosas inicialmente, la primera, se espera un histograma asimetrico positivo, con el 50% de datos entre 6 y 248, en donde el mayor recuento de observaciones se esperaria para el valor de Cs137 igual a 60.
La segunda, es que la variable no se distribuye normalmente, esto puede darse por los valores tan altos que pueden llegar a ser atipicos y que influencian a la media de los datos, como se observa el valor máximo de la muestra corresponde a 110000.
tabla.mtc## Media Mediana Moda Minimo Maximo
## 1 1202.541 248 60 6 110000Dentro de los métodos gráficos podemos observar las medidas de tendencia dentro del histograma de frecuencias de la variable en estudio, como se muestra acontinuación:
hist(z, freq = T,nclass = 20)
abline(v=mean(z),col='red',lwd=2)
abline(v=median(z),col='blue',lwd=2)
abline(v=Moda(z),col='green',lwd=2)
Medidas de dispersión
Dentro de las medidas de dispersión, podemos ver que la desviación estandar es bastante alta, por lo que se esperaria un coeficiente de variación también de la misma caracteristica
tabla.md= data.frame(stringsAsFactors=FALSE,
Desviación = sd(z),
Varianza = var(z),
Coeficiente_Variación=(sd(z)/mean(z))*100
)
tabla.md## Desviación Varianza Coeficiente_Variación
## 1 5030.296 25303874 418.3054Como se menciono previamente, el coeficiente de variación es muy alto por lo que el conjunto de datos tiende a ser heterogeneo y por lo tanto no se puede inferir sobre la población que la media de los datos de Cs137 sea igual a 1202.541, dada su gran variabilidad.
Medidas de distribución
Dentro de las medidas de distribución podemos tener el skewness, el cual verifica el comportamiento de la función de distribución de probabilidad, en donde esta, corrobora lo previamente inferido por medio de las medidas de tendencia central. Aquí determinamos que la distribución es asímetrica positiva si el estadístico del skewness es mayor a cero, si este es inferior a cero, se tiene asimetría negativa, por el contrario si el valor del estadístico es cero o tiende a cero, la distribución es simetrica.
También se puede determinar la forma de la curva de la función de distribución de los datos, en donde se tiene en cuenta la curtosis de los datos, aquí, si el estadistico es igual a 3, la función es mesocúrtica, si es menor a 3 platicúrtica y si por el contrario es superior a 3 tiende a ser leptocúrtica.
library(EnvStats)##
## Attaching package: 'EnvStats'## The following object is masked from 'package:RGeostats':
##
## print## The following objects are masked from 'package:stats':
##
## predict, predict.lm## The following object is masked from 'package:base':
##
## print.defaultskewness(z)## [1] 14.88686kurtosis(z, excess = FALSE)## [1] 272.3134Como se habia inferido de las medidas de tendencia central, la distribución es asimetrica positiva dado que el skewness es mayor a cero, por otro lado podemos observar que la curtosis es muy superior a 3, por lo tanto es leptocúrtica. Para poder observar mejor esto, observaremos las pruebas de normalidad de la variable y las representacones gráficas expuestas a continuación.
Normalidad
En búsqueda de establecer si el comportamiento de la variable de estudio es Gaussiana, se utilizaron algunas ayudas gráficas como lo son el histograma y la función de densidad de probabilidad empírica contra una distribución Normal teorica así mismo se establecieron algunas pruebas estadísticas.
hist(z, freq = F,col = "blue")
lines(density(z),col="red",lwd=2)
qqnorm(z)
qqline(z,col="blue")
La variable de estudio tiene un histograma asimétrico positivo tal como se habia mencionado anteriormente, donde su media es mayor que su mediana, esto implica la posible No normalidad de la variable como se puede identificar con el gráfico anterior, una posible causa podría ser la presencia de presuntos datos atípicos. Para corroborar la normalidad de la variable se procede a utilizar los test estadísticos de y bajo la siguiente prueba de hipótesis:
\[H_{0}:\:la\:variable\:sigue\:una\:distribución\: normal\]
library(nortest)
#Prueba de normalidad Kolmogorov-Smirnov-Lillefors
#P_value>5%, se rechaza la H0:Variable se distribuye normal.
lillie.test(z)##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: z
## D = 0.40599, p-value < 2.2e-16#Prueba de shapiro wilks
shapiro.test(z)##
## Shapiro-Wilk normality test
##
## data: z
## W = 0.17729, p-value < 2.2e-16Como el P_value) es muy inferior al \(5\%\), se rechaza la hipótesis nula, por lo tanto corroboramos que la variable de estudio no se comporta bajo una distribución normal, por lo tanto procedemos a transformar la variable.
Transformación de la variable respuesta Cs137
Se realiza una transformación a la variable respuesta Cs137, dentro de las opciones podemos efectuar una transformación de la familia boxcox, la cual pertenece a una tranformación de tipo logaritmica, también podemos realizar transformaciones como la anamorfosis gaussiana, sin embargo, observamos que se tiene una distribución de punto de masa de acuerdo con los histogramas iniciales, por lo tanto dentro de la primera aproximación, se realizará una tranformación de la siguiente manera: \[z_{trans}=\frac{Cs137^L-1}{L-1}\:\:\:L \ne 0\] \[z_{trans}=log(Cs137)\:\:\:L = 0\] Se realiza la transformación más sencilla y se grafica su histograma junto con su respectiva media y mediana, para determinar su asimetria.
library(EnvStats)
boxcox.list<-boxcox(cs137.db[,5],objective.name = "Shapiro-Wilk",optimize=T)
### Transformacion Ln(Cs137)
cs137.db = db.add(cs137.db,
ln.cs137=ifelse(cs137.db[,5]!=0,
(((cs137.db[,5]^boxcox.list$lambda)-1)/(boxcox.list$lambda-1)),
log(cs137.db[,5])))
hist(cs137.db[,6], freq = F,nclass = 20)
abline(v=mean(cs137.db[,6]),col='red',lwd=2)
abline(v=median(cs137.db[,6]),col='blue',lwd=2)
lines(density(cs137.db[,6]))
Como se puede apreciar en la figura anterior, la transformación a la variable, tuvo un efecto significativo sobre esta, en donde el histograma se convierte un poco más simetrico y se vislumbra que la diferencia entre la media y la mediana no es muy grande, sin embargo, al dibujar la función de distribución de probabilidad, esta no es completamente normal, pese a esto se aplican los test estadísticos para validar la normalidad de la variable transformada. Posterior a esto validamos las medidas de distribución para saber si la variable en cuestión ya cumple con el críterio para poder modelar.
\[H_{0}:\:la\:variable\:sigue\:una\:distribución\: normal\]
library(nortest)
lillie.test(cs137.db[,6])##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: cs137.db[, 6]
## D = 0.028949, p-value = 0.02224shapiro.test(cs137.db[,6])##
## Shapiro-Wilk normality test
##
## data: cs137.db[, 6]
## W = 0.99557, p-value = 0.001751# Medidas de distribución
skewness(cs137.db[,6])## [1] 0.01052493kurtosis(cs137.db[,6], excess = FALSE)## [1] 2.57139Como se puede observar, los test de normalidad, aún siguen sin ser significativos al \(5\%\), sin embargo, las medidas de distribución mejoraron, dando así una distribución asímetrica positiva y que tiende a ser mesocúrtica sin aún serlo completamente, pese a esto se debe realizar otro tratamiento a los datos, que sea más definitivo para poder continuar con la parte del modelado del semivariograma.
Anamorfosis Gaussiana
### Polinomios de Hermite ###
cs137.herm<-anam.fit(cs137.db, name = "ln.cs137", type = "gaus")
### Transformacion Raw to Gaussian - Z2Y
cs137.anam<-anam.z2y(anam = cs137.herm, cs137.db, names = "ln.cs137", flag.bound = T, radix = "Gaus" )
hist(cs137.anam[,7], freq = F,nclass = 20)
abline(v=mean(cs137.anam[,7]),col='red',lwd=2)
abline(v=median(cs137.anam[,7]),col='blue',lwd=2)
lines(density(cs137.anam[,7]))
Como se puede apreciar en la figura anterior, la transformación gaussiana que se realizo a la variable ya transformada bajo box cox, tiende a ser una distribución normal estandarizada con media cero y desviación estandar uno, la media y la mediana tienden al mismo valor y por esto se tiene una distribución simétrica, se realizan los test de normalidad y forma para validar que la variable sea normal.
\[H_{0}:\:la\:variable\:sigue\:una\:distribución\: normal\]
library(nortest)
lillie.test(cs137.anam[,7])##
## Lilliefors (Kolmogorov-Smirnov) normality test
##
## data: cs137.anam[, 7]
## D = 0.0093471, p-value = 0.9989shapiro.test(cs137.anam[,7])##
## Shapiro-Wilk normality test
##
## data: cs137.anam[, 7]
## W = 0.99964, p-value = 0.9996# Medidas de distribución
skewness(cs137.anam[,7])## [1] -0.0004222879kurtosis(cs137.anam[,7], excess = FALSE)## [1] 2.964812Como se observa, los test de normalidad fueron aprobados dado que el \(P\_Value>5\%\), además las medidas de distribución, como el skewness \(\sim0\), y la kurtosis \(\sim3\). por lo que podemos afirmar que esta variable cumple con el supuesto de normalidad bajo la tranformación de anamorfosis Gaussiana a la transformación box cox de la variable de estudio Cs137.
Análisis estructural
Dentro del análisis geoestadístico, se debe realizar una parte que corresponde al análisis estructural donde se evidencia el mayor rigor para poder determinar la forma funcional de los parámetros que intervienen en la interpolación, aquí se realiza una identificación de la tendencia que tienen los datos respecto a su ubicación espacial. Dentro de este se debe también verificar el grado de tendencia que existe, dado que es un parámetro que tendrá peso en cuanto a la interpolación se refiere.
También se debe realizar un análisis al variograma en diferentes direcciones y así saber si los datos tienen anisotropía. La isotropía implica que la dependencia espacial es una función de solo la distancia y no la dirección de la separación espacial entre las ubicaciones de muestreo.
Tendencia
Para el análisis de tendencia, utilizamos la variable ya transformada que cumple con el supuesto de normalidad, en este caso se realiza una distinción para la tendencia para cada drección (Estes, Nortes), además de la combinación de las mismas.
Tendencia Este
Se observa la nube de puntos entre la variable respuesta normalizada y la componente este de su ubicación, \(Z\_trans\sim Este\), allí mostramos la media de la variable transformada para saber si los puntos se distribuyen de forma aleatoria sobre esta, o si se puede observar un comportamiento definido.
#Tendencia grado 1 orden 1
fit1orden<-lm(cs137.anam[,7]~cs137.anam[,4])
#Tendencia orden 2
fit2orden<-lm((formula=cs137.anam[,7]~poly(cs137.anam[,4],2,raw = T)))
eq = function(x){fit2orden$coefficients[3]*x*x+fit2orden$coefficients[2]*x+fit2orden$coefficients[1]}
#Nuevamente se verifica la tendencia 1 en el gráfico
plot(cs137.anam[,7]~cs137.anam[,4])
abline(h=mean(cs137.anam[,7]),col='red',lwd=2)
lines(smooth.spline(cs137.anam[,7]~cs137.anam[,4]),lty = 2, col = "cyan")#Tendencia Orden 3
abline(fit1orden,col='green',lwd=2)
curve(eq,col='purple', add=T)
# Add a legend
legend(673600, 3.4, legend=c("Media", "Tendencia 1", "Tendencia 2", "Tendencia 3"),
col=c("red", "green","purple","cyan"), lty=1:3, cex=0.8)
En la gráfica anterior se observan las tendencias de orden 1, 2 y 3, además de la media de la variable Gaussiana. Allí no se encuentra una tendencia clara en donde si se aumenta o disminuye la dirección Oeste-Este, afecte directamente o inversamente la variable de estudio, todas las tendencias siguen en cierta medida la media de la variable transformada, por lo tanto se puede inferir que no existe una tendencia en la dirección de los Estes.
summary(fit1orden)##
## Call:
## lm(formula = cs137.anam[, 7] ~ cs137.anam[, 4])
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.2726 -0.6064 0.0296 0.6593 3.2472
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.524e+00 1.508e+00 2.999 0.00276 **
## cs137.anam[, 4] -6.379e-06 2.126e-06 -3.000 0.00276 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9962 on 1169 degrees of freedom
## Multiple R-squared: 0.00764, Adjusted R-squared: 0.006791
## F-statistic: 9 on 1 and 1169 DF, p-value: 0.002758summary(fit2orden)##
## Call:
## lm(formula = (formula = cs137.anam[, 7] ~ poly(cs137.anam[, 4],
## 2, raw = T)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.3031 -0.5773 0.0372 0.6354 3.1322
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.474e+02 7.388e+01 -4.702 2.88e-06 ***
## poly(cs137.anam[, 4], 2, raw = T)1 9.856e-04 2.082e-04 4.733 2.48e-06 ***
## poly(cs137.anam[, 4], 2, raw = T)2 -6.988e-10 1.467e-10 -4.764 2.13e-06 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9871 on 1168 degrees of freedom
## Multiple R-squared: 0.02656, Adjusted R-squared: 0.02489
## F-statistic: 15.93 on 2 and 1168 DF, p-value: 1.49e-07fit3orden<-lm((formula=cs137.anam[,7]~poly(cs137.anam[,4],3,raw = T)))
summary(fit3orden)##
## Call:
## lm(formula = (formula = cs137.anam[, 7] ~ poly(cs137.anam[, 4],
## 3, raw = T)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.0918 -0.6616 0.0662 0.6580 2.9810
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.618e+04 3.377e+03 7.751 1.97e-14 ***
## poly(cs137.anam[, 4], 3, raw = T)1 -1.112e-01 1.428e-02 -7.786 1.51e-14 ***
## poly(cs137.anam[, 4], 3, raw = T)2 1.575e-07 2.013e-08 7.821 1.16e-14 ***
## poly(cs137.anam[, 4], 3, raw = T)3 -7.430e-14 9.458e-15 -7.856 8.94e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9624 on 1167 degrees of freedom
## Multiple R-squared: 0.07545, Adjusted R-squared: 0.07308
## F-statistic: 31.75 on 3 and 1167 DF, p-value: < 2.2e-16Previamente se habia inferido a partir de la grafica que no se tenia una tendencia de ningún orden, sin embargo, con el ajuste de un modelo para cada orden de la tendencia se evalua cual de estos explica mejor la variable de estudio transformada. Aqui, en los resultados previos, ninguna de las tres tendencias tien un \(R^2_{Adjus}\) significativo, ya que todos estos se encuentran por debajo del \(10\%\).
Tendencia Norte
Se observa la nube de puntos entre la variable respuesta normalizada y la componente norte de su ubicación, \(Z\_trans\sim Norte\), allí mostramos la media de la variable transformada para saber si los puntos se distribuyen de forma aleatoria sobre esta, o si se puede observar un comportamiento definido.
#Tendencia grado 1 orden 1
fit1orden<-lm(cs137.anam[,7]~cs137.anam[,3])
#Tendencia orden 2
fit2orden<-lm((formula=cs137.anam[,7]~poly(cs137.anam[,3],2,raw = T)))
eq = function(x){fit2orden$coefficients[3]*x*x+fit2orden$coefficients[2]*x+fit2orden$coefficients[1]}
plot(cs137.anam[,7]~cs137.anam[,3])
abline(h=mean(cs137.anam[,7]),col='red',lwd=2)
lines(smooth.spline(cs137.anam[,7]~cs137.anam[,3]),lty = 2, col = "cyan")#Tendencia Orden 3
abline(fit1orden,col='green',lwd=2)
curve(eq,col='purple', add=T)
# Add a legend
legend(5660000, 3.4, legend=c("Media", "Tendencia 1", "Tendencia 2", "Tendencia 3"),
col=c("red", "green","purple","cyan"), lty=1:3, cex=0.8)
A diferencia de la tendencia en dirección de las Estes, mediante los métodos graficos se puede observar que entre más aumenta la dirección de la Norte, aumenta la variable transformada, allí se observa una notoria tendencia de orden 1 y 2, sin embargo, para poder identificar cual de estas es la más significativa, se determinará por medio del \(R^2_{Adjus}\).
summary(fit1orden)##
## Call:
## lm(formula = cs137.anam[, 7] ~ cs137.anam[, 3])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.66972 -0.46847 0.05837 0.46193 2.64218
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.199e+02 1.105e+01 -28.96 <2e-16 ***
## cs137.anam[, 3] 5.624e-05 1.942e-06 28.96 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7631 on 1169 degrees of freedom
## Multiple R-squared: 0.4178, Adjusted R-squared: 0.4173
## F-statistic: 839 on 1 and 1169 DF, p-value: < 2.2e-16summary(fit2orden)##
## Call:
## lm(formula = (formula = cs137.anam[, 7] ~ poly(cs137.anam[, 3],
## 2, raw = T)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.65064 -0.46409 0.07183 0.45679 2.64411
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.330e+03 4.944e+03 1.887 0.0594 .
## poly(cs137.anam[, 3], 2, raw = T)1 -3.337e-03 1.738e-03 -1.919 0.0552 .
## poly(cs137.anam[, 3], 2, raw = T)2 2.983e-10 1.528e-10 1.952 0.0512 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7621 on 1168 degrees of freedom
## Multiple R-squared: 0.4197, Adjusted R-squared: 0.4187
## F-statistic: 422.4 on 2 and 1168 DF, p-value: < 2.2e-16fit3orden<-lm((formula=cs137.anam[,7]~poly(cs137.anam[,3],3,raw = T)))
summary(fit3orden)##
## Call:
## lm(formula = (formula = cs137.anam[, 7] ~ poly(cs137.anam[, 3],
## 3, raw = T)))
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.65064 -0.46409 0.07183 0.45679 2.64411
##
## Coefficients: (1 not defined because of singularities)
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 9.330e+03 4.944e+03 1.887 0.0594 .
## poly(cs137.anam[, 3], 3, raw = T)1 -3.337e-03 1.738e-03 -1.919 0.0552 .
## poly(cs137.anam[, 3], 3, raw = T)2 2.983e-10 1.528e-10 1.952 0.0512 .
## poly(cs137.anam[, 3], 3, raw = T)3 NA NA NA NA
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7621 on 1168 degrees of freedom
## Multiple R-squared: 0.4197, Adjusted R-squared: 0.4187
## F-statistic: 422.4 on 2 and 1168 DF, p-value: < 2.2e-16Dentro de los \(R^2_{Adjus}\), se observa que para las 3 tendencias es del \(40\%\), sin embargo, la tendencia de orden 1, es la única que tiene significativo al \(5\%\) sus covariables (intercepto y orden de la variable). Por lo tanto se puede inferir que existe una tendencia de orden 1 en cuanto a la variación de la coordenada Norte.
Tendencia Este - Norte
Para establecer si la tendencia se da únicamente en una dirección, se debe realizar el anaílis conjunto entre las variables, dado que puede existir una correlación entre las componentes Este y Norte de la muestra. Para esto se ejecuta el siguiente anális y se evalúa el \(R^2_{Adjus}\) del modelo.
##Tendencia De orden 1 en X,Y
fit.xy<-lm(cs137.anam[,7]~cs137.anam[,4]+cs137.anam[,3])
#VISUALIZACIÓN DE LA SUPERFICIE EN 3D
library(plot3D)
#1. Grilla
grid.lines = 25
x.pred <- seq(min(cs137.anam[,4]),max(cs137.anam[,4]), length.out= grid.lines)
y.pred <- seq(min(cs137.anam[,3]),max(cs137.anam[,3]), length.out= grid.lines)
grid.xy <- expand.grid(X=x.pred,Y=y.pred)
z.pred.xy <- matrix(predict(fit.xy,newdata = grid.xy), nrow = grid.lines,ncol = grid.lines)
scatter3D(cs137.anam[,4], cs137.anam[,3], cs137.anam[,7],
pch = 8, cex = 0.2, theta = 50, phi = 20,
ticktype = "detailed",xlab = "Este", ylab = "Norte",
zlab = "CS137", surf = list(x = x.pred, y = y.pred,
z = z.pred.xy, facets = NA), main = "Tendencia xy orden 1")
En el ScatterPlot en 3D, se observa la tendencia notoria de que el aumento en la componente Norte, se da mayor variación de la variable transformada, además de que se mantiene el comportamiento de que la dirección Este permanece invariante en cuanto a su relación con la variable de estudio.
#2. Imagen
image2D(z.pred.xy,x=x.pred,y=y.pred)
summary(fit.xy)##
## Call:
## lm(formula = cs137.anam[, 7] ~ cs137.anam[, 4] + cs137.anam[,
## 3])
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.78121 -0.47280 0.05778 0.44559 2.64190
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -3.383e+02 1.179e+01 -28.694 < 2e-16 ***
## cs137.anam[, 4] 7.137e-06 1.682e-06 4.244 2.37e-05 ***
## cs137.anam[, 3] 5.858e-05 2.005e-06 29.217 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.7576 on 1168 degrees of freedom
## Multiple R-squared: 0.4267, Adjusted R-squared: 0.4257
## F-statistic: 434.6 on 2 and 1168 DF, p-value: < 2.2e-16Finalmente podemos apreciar que la variación se da únicamente en el aumento de la dirección Norte, sin embargo, al evaluar el modelo respecto a la tendecia de orden 1 en ambas direcciones, observamos que el \(R^2_{adjus}\) es superior a los analizados previamente, además de que ambas componentes resultan ser significativas y directamente proporcionales a la variable transformada, por lo tanto se puede apreciar una tendencia de orden 1 en ambas direcciones. Esto se puede dar debido a la relación que tienen ambas componentes.
Análisis Variografico
Se calcula la diferencia al cuadrado entre los valores de cada par de valor de la muestra separado por una distancia h. Las diferencias resultantes se grafican frente a la separación de pares de muestras en el espacio geográfico y forman la nube de variogramas. La nube de variogramas en sí misma es una herramienta poderosa para explorar características de datos.
#Nube variografica [Z(x)-Z(x+h)]^2
cloud_cs137<-cloud.calc(cs137.anam, lagnb = 2000)
plot(cloud_cs137, xlab="h", ylab="[Z(x)-Z(x+h)]^2",pos.legend=1)
De acuerdo con la nube variografica, podemos observar que la distancia mayor a la que se encuentra el variograma es cercana a los 80000, además observamos como la tendencia cerca a la distancia \(h=\) 30000 tiende a bajar, a partir de esto se va a utilizar la mitad de la distancia maxima para poder interpretar y ajustar los variogramas.
dmax<-(((max(CS137.table$Este)-min(CS137.table$Este))^2+
(max(CS137.table$Norte)-min(CS137.table$Norte))^2)^(1/2))/2
dmax## [1] 42635.41#variograma omnidireccional
cs137.anam <- db.locate(cs137.anam, "Gaus.ln.cs137", "z",1)
cs137.vario<-vario.calc(cs137.anam,lag =2000 ,nlag = 20)
plot(cs137.vario,npairdw=T)
Como se puede observar la mitad de la distancia máxima ronda los 40000, por lo tanto se establecieron los parámetros de distancia y cantidad de separaciones que como resultado nos dieran esos 40000, en este caso fueron 20 intervalos, cada uno a una separación de 2000.
Como se puede apreciar en el semivariograma empírico, la meseta está muy cercano a los 30000, el valor de esta es aproximadamente 1.2. Al inicio podemos observar que el punto inicial es muy cercano a cero, por lo tanto no se observaría un efecto pepita.
Sin embargo, se debe realizar los variogramas direccionales para poder establecer que la variable de estudio, en este caso la variable con doble transformación sea isotrópica, por tal motivo, se calculan los variogramas direccionales, a 0, 45, 90, 135 grados respectivamente.
#variograma direccional
trans.cs137.variodir<-vario.calc(cs137.anam,lag =2000 ,nlag = 20,dir=c(0,45,90,135))
plot(trans.cs137.variodir,npairdw=T,pos.legend=5)
Dentro de los variogramas direccionales, se puede apreciar, que existen dos diferentes comportamientos para el posible modelado de la variable transformada, esto puede corresponder a una Anisotopía zonal o sencillamente a un modelo no estacionario para dos direcciones, se observa que al rededor de una distancia igual a 10000 metros, se separan en dos grupos, los variogramas con dirección 0 y 135 grados, se estabilizan con un valor aproximado de 0.6 de meseta, por otro lado, los variogramas direccionales para 45 y 90 grados no se alcanza a previsualizar la meseta, puede que sea un modelo no estacionario el que explique este fenomeno, dado que hasta aquí corresponde la parte estructurada de la variable.
Como se habia mencionado, se establece un modelo sin efecto pepita, puesto que todos los variogramas direccionales salen muy cercanos al origen, además, se tiene un modelo anidado puesto que se observa una parte no estructurada. Estos modelos que se escogieron fue, el modelo estable, dado que su comportamiento tiende a ser parabólico hasta que encuentra la meseta y se estabiliza, por otra parte se anida con un modelo esférico dado que su comportamiento tiende a ser lineal para modelos no estacionarios, cabe aclarar que esto sucede cuando el rango tiende a ser muy grande.
# melem.name()
trans.cs137.vardirteo<- model.auto(trans.cs137.variodir,
struct = c("Stable","Spherical"),
title="Modelo Direccional Transformacion(Cs137)")
#trans.cs137.vardirteoComo resultado del modelo ajustado se encuentran ciertas particularidades, donde el modelo estable, su el rango mínimo es sobre los 2600 metros allí alcanza una meseta muy pequeña de 0.3 aproximadamente, además, tiene un angulo de Anisotropía cercano a 45 grados. Por otro lado, se tiene un modelo esférico, el cual alcanza la meseta en 2.1 a una distancia aproximada de 54000 metros, superior a la distancia donde se encuentra la parte estructurada de la variable (40000 metros), además, un ángulo de Anisotropía a 345 grados aproximadamente.
Para validar otro tipo de modelo, se realiza una simuación mediante la libreria Rgeostats, la cual, establece un número de iteraciones bajos los parámetros ingresados y determina cual sería el mejor modelo a establecer.
trans.cs137.vardirteo.auto<-model.auto(trans.cs137.variodir,
maxiter=10000,
title="Modelo Direccional automatico Transformacion(Cs137)")
#trans.cs137.vardirteo.autoDe acuerdo con las simulaciones, se genera un modelo anidado, la primera parte es un modelo Gaussiano, con una meseta aproximada a 0.34, muy similar al del modelo estable, además con un rango en 443230.588 metros, y un ángulo de 397.279 grados, esto varia drasticamente con el modelo inicialmente ajustado. La segunda parte del modelo corresponde a un modelo esférico, al igual que la alternativa seleccionada previamente, sin embargo, la meseta de este, llega hasta 2.5, en un rango de 53000 metros con un ángulo de anisotropía de 342.132 grados, este puede ser más certero puesto que modela un comportamiento no estacionario.
Validación cruzada
Para determinar cual de los modelos se debe seleccionar, se realiza una validación cruzada a estos, donde a partir de sus errores, se determinará cual modelo es más apto para realizar la interpolación de la variable transformada.
Se realiza el calculo de los errores de estimación y los errores estandarizados para cada observación, sin embargo, se calcula la media para ambos valores y la varianza de los errores estandarizados.
####Vecindad movil
trans.cs137.neigh.mv<-neigh.create(ndim = 2,type = 2,nmaxi = 10,radius = 35000,
flag.sector = T,nsect = 8,nsmax = 2)
####Validación cruzada modelo 1, estable esferico
cs137.anam<-xvalid(cs137.anam,trans.cs137.vardirteo,trans.cs137.neigh.mv,
radix ="model1")
#Cambio de rol de las variables
cs137.anam <- db.locate(cs137.anam, "Gaus.ln.cs137", "z",1)
cs137.anam <- db.locate(cs137.anam, 9)
####Validación cruzada modelo 2, Gaussiano esferico
cs137.anam<-xvalid(cs137.anam,model=trans.cs137.vardirteo.auto,trans.cs137.neigh.mv,
radix ="model2")
tabla.xvalid= data.frame(stringsAsFactors=FALSE,
row.names = c("Estable-esférico","Gaussiano-esférico"),
Media_Error_Estimacion = c(mean(cs137.anam@items$model1.Gaus.ln.cs137.esterr),
mean(cs137.anam@items$model2.Gaus.ln.cs137.esterr)),
Media_Error_Estandar=c(mean(cs137.anam@items$model1.Gaus.ln.cs137.stderr),
mean(cs137.anam@items$model2.Gaus.ln.cs137.stderr)),
Desviacion_Error_sd=c(sd(cs137.anam@items$model1.Gaus.ln.cs137.stderr),
sd(cs137.anam@items$model2.Gaus.ln.cs137.stderr))
)
tabla.xvalid## Media_Error_Estimacion Media_Error_Estandar
## Estable-esférico 0.004213305 0.01102078
## Gaussiano-esférico 0.002548425 0.03489675
## Desviacion_Error_sd
## Estable-esférico 1.440527
## Gaussiano-esférico 5.448022Si se analiza la tabla anterior, el modelo, Estable-Esférico, es quién debería ser seleccionado para la interpolación dado que este mantiene una desviación de los errores estandarizados, muy cercano a uno. Sin embargo, se analizarán los tres parámetros.
La media del error de estimación para ambos modelos es muy cercana a cero, esto corresponde a que ambos modelos llegan a predecir el comportamiento de la variabilidad de la variable transformada, lo cual es apto para interpolar, pero, cuando se observa los errores esntadarizados, ya se marca una diferencia a favor del primer modelo de la tabla, en donde por una potencia de 10, el modelo estable esférico tiene una media de error mucho más cercana a cero, por lo cual se observa que su desvicación tiende a 1, tal cual como se puede apreciar el histograma de la variable transformada.
Estimación Kriging
Para la interpolación por Kriging, se utiliza el Krigeado de tipo Universal, dado que se encontró tendencia dentro de los datos.
Se realiza la interpolación a la variable con doble transformación (Boxcox y anamorfosis), además, con tendencia de orden 1 en las coordenadas x, y.
### !!!! Cambio de rol de variable de estudio Gaus.ln.cs137
cs137.anam<-db.locate(cs137.anam, "Gaus.ln.cs137", "z",1)
cs137.anam<-db.locate(cs137.anam, c(10,11))
### Crear grilla
cs137.grid <- db.grid.init(cs137.anam, nodes = c(100, 100))
### Kriging Universal
cs137.grid<-kriging(dbin=cs137.anam, dbout=cs137.grid,model = trans.cs137.vardirteo,
neigh =trans.cs137.neigh.mv, radix = "KU",uc = c("1","xy"))
summary(cs137.grid@items$KU.Gaus.ln.cs137.stdev)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 0.1062 0.3228 0.4122 0.5965 0.7285 2.3404par(mar = c(4, 4, .1, .1))
plot(cs137.grid, name.image="KU.Gaus.ln.cs137.stdev",pos.legend=1, title="Desviación KU para Cs137")
arrows( 718000,5700000, 735000,5703000)
text(735000,5709000,"Reactor",pos = 1,cex = 1.2)
plot(cs137.db, bg="green",pch=21, title="Muestra de Cs137")
A partir de los errores de estimación del kriging universal, se puede apreciar que estos tienden a ser mínimos en la misma superficie donde se espaciliza la muestra inicial, allí se puede detectar una pequeña alteración dentro de está área de influencia que corresponde al cuerpo de agua cercano al reactor nuclear, en donde, no se tiene una muestra para el valor de Cesio 137. Por otro lado los errores tienden a ser homogeneos dentro de está área y se observa que aumenta en la periferia.
Posterior al Krigeado se debe realizar la doble anti transformación, primero se debe pasar los valores de las predicciones a la inversa de la anamorfosis, posterior a esto se debe realizar la inversa de boxcox, tal como se presenta a continuación.
##Anti transformar para volver a valores reales
cs137.grid<-anam.y2z(anam=cs137.herm, cs137.grid, names = "KU.Gaus.ln.cs137.estim",
radix = "bxcx.raw")
### Transfomacion inversa - (boxcox)^-1
cs137.grid = db.add(cs137.grid,
raw.cs137=(cs137.grid@items$bxcx.raw.KU.Gaus.ln.cs137.estim*
(boxcox.list$lambda-1)+1)^(1/boxcox.list$lambda))Si observamos el mapa de predicciones resultado del Kriging Universal, tiene un valor mínimo de 6 y un valor superior cercano a 80000, si recordamos las medidas de tendencia central de la variable inicial, este tenia un valor mínimo muy similar, sin embargo, el valor máximo oscilaba los 110000.
Si comparamos el mapa de predicciones con la muestra inicial, podemos observar una tendencia de los datos predecidos a tener un comportamiento muy similar al muestreo, concentrando los mayores valores en las cercanias al reactor nuclear.
par(mar = c(4, 4, .1, .1))
plot(cs137.grid, name.image="raw.cs137",pos.legend=1, title="Interpolación por KU para Cs137")
plot(cs137.db, bg="green",pch=21, title="Muestra de Cs137")
Simulaciones
Si se quiere identificar con mayor precisión el valor de la superficie más probable, se realizan las simulaciones para poder establecerla. Para esto se debe utilizar la variable con doble transformación, aquella que sigue un comportamiento normal, además del variograma teorico que mejor ajusto y el vecindacio movil, además, se debe seleccionar un valor ara realizar las iteraciones/simulaciones en este caso se seleccionaron 1000 simulaciones.
cs137.simu <- db.grid.init(cs137.anam, nodes = c(250, 250))
##Simulacion
cs137.simu<-simtub(dbin=cs137.anam, dbout = cs137.simu, trans.cs137.vardirteo, trans.cs137.neigh.mv, nbsimu = 1000, nbtuba = 100)
##varianza de superfice + probable
cs137.simu <- db.compare(cs137.simu, names = seq(4,1003), fun = "var")
#cs137.simu<-db.delete(cs137.simu, names=c("desviacion","bxcx.raw.desviacion","sd.raw.cs137"))
cs137.simu = db.add(cs137.simu, desviacion=(cs137.simu@items$var)^(1/2))Lo primero que se hace es analizar la variabilidad de la superficie creada a partir de las simulaciones y/o realizaciones. Seleccionamos la función de varianza para las 1000 realizaciones y posterior sacamos su raíz para tener la desviación de las realizaciones.
par(mar = c(4, 4, .1, .1))
plot(cs137.simu, name.image="desviacion",pos.legend=1, title="Desviación Estandar para la Simulación de Cs137")
arrows( 718000,5700000, 735000,5703000)
text(735000,5709000,"Reactor",pos = 1,cex = 1.2)
plot(cs137.db, bg="green",pch=21, title="Muestra de Cs137")
Como se observa, al igual que el Kriging universal, los valores de mayor variabilidad o incertidumbre en este caso, son los valores extremos, o en donde la variable no cuenta con un muestreo definido, también se observa que el comportamiento de la incertidumbre asemeja un mapa de calor para la variable de estudio Cs137 en donde los valores mínimos se establecen en donde se encuentran muestras de el mismo. Si realizamos una comparativa general entre los errores del Kriging Universal y de las 1000 realizaciones, que sus valores no son muy disimilares, en cuanto a comportamiento se refiere, se tienen pequeñas variaciones o mayor precisión en este caso dado que las simulaciones no suavizan la superficie.
## Superficie más probable
cs137.simu <- db.compare(cs137.simu, names = seq(4,1003), fun = "mean")
cs137.simu<-anam.y2z(anam=cs137.herm, cs137.simu, names = "mean", radix = "bxcx.raw")
cs137.simu = db.add(cs137.simu, raw.cs137=(cs137.simu@items$bxcx.raw.mean*(boxcox.list$lambda-1)+1)^(1/boxcox.list$lambda))Si se analiza el mapa de la superficie más probable realizado por las simulaciones, se tiene un valor mínimo de 7.5, la muestra inicialmente tenia un valor mínimo de 6. Además su valor máximo es de 90000 y el de la muestra inicial era de 110000. Sin embargo, al igual que en el Kriging Universal el mapa de la superficie más probable sigue el mismo comportamiento que el de la muestra, dando como resultado un mapa con un comportamiento similar al del suavizado por Kriging.
par(mar = c(4, 4, .1, .1))
plot(cs137.simu, name.image="raw.cs137",pos.legend=1, title="Simulación para Cs137")
plot(cs137.db, bg="green",pch=21, title="Muestra de Cs137")