Pulimiento de medianas

Simulación de un campo aleatorio.

Cargar librerias

Lista de librerías con link a la documentación.

gstat
sp

Lista de librerías con link a la documentación.

Grilla de las ubicaciones espaciales.

n_x <- 4
n_y <- 6
x <- seq(0, 1, len = n_x)
y <- seq(0, 1, len = n_y)
coordenadas <- as.data.frame(expand.grid(x, y))
names(coordenadas) <- c("X", "Y")

Encabezado coordenadas

X	Y
0.0000000	0.0
0.3333333	0.0
0.6666667	0.0
1.0000000	0.0
0.0000000	0.2
0.3333333	0.2

Definición de objeto VGM

Esto define un objeto vgm que es el tipo de objeto que usa el paquete gstat para los modelos teóricos de variograma. Con este objeto se pueden definir modelos anidados.

vario <- vgm(10, # Punto de silla
             "Exp", # Modelo, ver documentación
             0.5)  # Rango
print(vario)

##   model psill range
## 1   Exp    10   0.5

Matriz de varianza dadas coordenadas.

coordinates(coordenadas) <- ~X + Y
class(coordenadas) # Cambio de objedto dataframe a sp

## [1] "SpatialPoints"
## attr(,"package")
## [1] "sp"

cov_mat <- vgmArea(coordenadas, # Matriz de ubiaciones SP
        vgm = vario) # VGM object

print(dim(cov_mat))

## [1] 24 24

Simulación.

Simulación dada la media y la matriz de varianza

mu  <- rep(0, n_x * n_y) # Media del proceso
simu <- rmvnorm(1,
                mean = mu,
                sigma = cov_mat)
print(simu[1:5])

## [1]  0.3202784  1.8696771 -0.3138076  1.3413568  1.6079527

Pulimiento de medianas

Unir las coordenadas con la columna de simulación

data <- as.data.frame(cbind(coordenadas@coords,
                            Simula = t(simu)))
names(data) <- c("X", "Y", "Var")
print(head(data))

X	Y	Var
0.0000000	0.0	0.3202784
0.3333333	0.0	1.8696771
0.6666667	0.0	-0.3138076
1.0000000	0.0	1.3413568
0.0000000	0.2	1.6079527
0.3333333	0.2	0.7935526

Reshape para matriz, esto transforma la tabla de datos en matriz

tabla <- reshape2::dcast(data,
                         X ~ Y,
                         value.var = "Var")
rownames(tabla) <- tabla[, 1]
tabla <- tabla[, c(-1)]
print(tabla)

	0	0.2	0.4	0.6	0.8	1
0	0.3202784	1.6079527	0.9830152	7.519702	1.7294144	-1.1107318
0.333333333333333	1.8696771	0.7935526	4.8470537	3.284085	-0.6781032	3.4614707
0.666666666666667	-0.3138076	0.4327380	5.2294823	4.273177	3.0492405	-0.1208445
1	1.3413568	3.9792647	3.9830081	7.672783	4.4182895	2.1628881

Pulimiento de medianas de la tabla

med <- medpolish(tabla)

## 1: 27.06093
## Final: 26.97376

geo_data <- reshape2::melt(med$residuals)
print(med)

## 
## Median Polish Results (Dataset: "tabla")
## 
## Overall: 1.880851
## 
## Row Effects:
##                 0 0.333333333333333 0.666666666666667                 1 
##        -0.4260350         0.4260350        -0.5304161         2.0022273 
## 
## Column Effects:
##          0        0.2        0.4        0.6        0.8          1 
## -1.3993903 -0.4107555  1.3200484  3.3562228  0.4049045 -1.5957351 
## 
## Residuals:
##                          0      0.2     0.4      0.6      0.8        1
## 0                  0.26485  0.56389 -1.7918  2.70866 -0.13031 -0.96981
## 0.333333333333333  0.96218 -1.10258  1.2201 -2.37902 -3.38989  2.75032
## 0.666666666666667 -0.26485 -0.50694  2.5590 -0.43348  1.29390  0.12446
## 1                 -1.14233  0.50694 -1.2201  0.43348  0.13031 -0.12446

Reshape de los datos, con efecto de la fila y la columna

tabla_residuales <- as.data.frame(med$residuals)
names(tabla_residuales) <- med$col
rownames(tabla_residuales) <- med$row
geo_data <- reshape2::melt(as.matrix(tabla_residuales))

geo_data <- cbind(data,
                  geo_data,
                  med$overall)
names(geo_data) <- c("X",
                     "Y",
                     "Var",
                     "Efecto fila",
                     "Efecto columa",
                     "Residual",
                     "Efecto Global")
print(geo_data)

X	Y	Var	Efecto fila	Efecto columa	Residual	Efecto Global
0.0000000	0.0	0.3202784	-0.4260350	-1.3993903	0.2648525	1.880851
0.3333333	0.0	1.8696771	0.4260350	-1.3993903	0.9621811	1.880851
0.6666667	0.0	-0.3138076	-0.5304161	-1.3993903	-0.2648525	1.880851
1.0000000	0.0	1.3413568	2.0022273	-1.3993903	-1.1423315	1.880851
0.0000000	0.2	1.6079527	-0.4260350	-0.4107555	0.5638920	1.880851
0.3333333	0.2	0.7935526	0.4260350	-0.4107555	-1.1025782	1.880851
0.6666667	0.2	0.4327380	-0.5304161	-0.4107555	-0.5069416	1.880851
1.0000000	0.2	3.9792647	2.0022273	-0.4107555	0.5069416	1.880851
0.0000000	0.4	0.9830152	-0.4260350	1.3200484	-1.7918495	1.880851
0.3333333	0.4	4.8470537	0.4260350	1.3200484	1.2201189	1.880851
0.6666667	0.4	5.2294823	-0.5304161	1.3200484	2.5589987	1.880851
1.0000000	0.4	3.9830081	2.0022273	1.3200484	-1.2201189	1.880851
0.0000000	0.6	7.5197021	-0.4260350	3.3562228	2.7086630	1.880851
0.3333333	0.6	3.2840853	0.4260350	3.3562228	-2.3790238	1.880851
0.6666667	0.6	4.2731768	-0.5304161	3.3562228	-0.4334812	1.880851
1.0000000	0.6	7.6727825	2.0022273	3.3562228	0.4334812	1.880851
0.0000000	0.8	1.7294144	-0.4260350	0.4049045	-0.1303064	1.880851
0.3333333	0.8	-0.6781032	0.4260350	0.4049045	-3.3898941	1.880851
0.6666667	0.8	3.0492405	-0.5304161	0.4049045	1.2939008	1.880851
1.0000000	0.8	4.4182895	2.0022273	0.4049045	0.1303064	1.880851
0.0000000	1.0	-1.1107318	-0.4260350	-1.5957351	-0.9698130	1.880851
0.3333333	1.0	3.4614707	0.4260350	-1.5957351	2.7503195	1.880851
0.6666667	1.0	-0.1208445	-0.5304161	-1.5957351	0.1244554	1.880851
1.0000000	1.0	2.1628881	2.0022273	-1.5957351	-0.1244554	1.880851

Validación de la descomposición

valida <- cbind(geo_data$Var,
                geo_data[["Efecto fila"]] +
                geo_data[["Efecto columa"]] +
                geo_data[["Residual"]] +
                geo_data[["Efecto Global"]])
valida <- as.data.frame(valida)
names(valida) <- c("datos", "suma")
print(valida)

datos	suma
0.3202784	0.3202784
1.8696771	1.8696771
-0.3138076	-0.3138076
1.3413568	1.3413568
1.6079527	1.6079527
0.7935526	0.7935526
0.4327380	0.4327380
3.9792647	3.9792647
0.9830152	0.9830152
4.8470537	4.8470537
5.2294823	5.2294823
3.9830081	3.9830081
7.5197021	7.5197021
3.2840853	3.2840853
4.2731768	4.2731768
7.6727825	7.6727825
1.7294144	1.7294144
-0.6781032	-0.6781032
3.0492405	3.0492405
4.4182895	4.4182895
-1.1107318	-1.1107318
3.4614707	3.4614707
-0.1208445	-0.1208445
2.1628881	2.1628881