índice de Moran

Autocorrelación espacial de materia organica

Para el estudio de variables, como lo es para la materia organica en este caso, se dan autocorrelaciones según el espacio en la misma, esto puede ser determinado a través del índice de Moran.

el índice de Moran se peude definir de la siguiente forma

El índice de moran se peude definir de la sighuiente forma: \[I=\frac{N}{S_0}\frac{\sum_{(2)}W_{ij}(Y_i-\bar Y)(Y_j-\bar Y)}{\sum_{i=1}^N(Y_i-\bar Y)^2}\\Donde: W_{ij}=elementos~de~la~matriz~de~pesos\\espaciales~correspondientes~al~par~(i,j)\\Y=Valor~medio~o~esperado~de~la~variable~y\\N=número~de~observaciones~o~tamaño~muestral\\S_0=\sum~_i\sum~_jW_{ij}=\sum~_{(2)}W_{ij}\]

Con el siguiente grafico se puede mostrar la forma distribucional de los contenidos de materia organica en un parcela idealmente rectangular para facilitar la ejemplificación, donde sus datos se distribuyen de forma normal.

Materia.O=rnorm(n=150, mean = 3, sd = 0.5)
xy=expand.grid(x=seq(1,10), y=seq(1,15))
plot(xy,col=Materia.O, pch=16)

Para determinar las autocorrelaciones se es necesario dar las medidas de distancia de cada dato en la matirz de datos utilizando una medida especifica. posteriormente estos son invertidos, lo que quiere decir que cada dato va a aser divisor de 1 y se regresan los datos diagonales, como se muestra a continuación.

mo.dist<-as.matrix(dist(cbind(xy)))
mo.dist.inv=1/mo.dist
diag(mo.dist.inv)<-0
mo.dist.inv[1:10, 1:10]
##            1         2         3         4         5         6         7
## 1  0.0000000 1.0000000 0.5000000 0.3333333 0.2500000 0.2000000 0.1666667
## 2  1.0000000 0.0000000 1.0000000 0.5000000 0.3333333 0.2500000 0.2000000
## 3  0.5000000 1.0000000 0.0000000 1.0000000 0.5000000 0.3333333 0.2500000
## 4  0.3333333 0.5000000 1.0000000 0.0000000 1.0000000 0.5000000 0.3333333
## 5  0.2500000 0.3333333 0.5000000 1.0000000 0.0000000 1.0000000 0.5000000
## 6  0.2000000 0.2500000 0.3333333 0.5000000 1.0000000 0.0000000 1.0000000
## 7  0.1666667 0.2000000 0.2500000 0.3333333 0.5000000 1.0000000 0.0000000
## 8  0.1428571 0.1666667 0.2000000 0.2500000 0.3333333 0.5000000 1.0000000
## 9  0.1250000 0.1428571 0.1666667 0.2000000 0.2500000 0.3333333 0.5000000
## 10 0.1111111 0.1250000 0.1428571 0.1666667 0.2000000 0.2500000 0.3333333
##            8         9        10
## 1  0.1428571 0.1250000 0.1111111
## 2  0.1666667 0.1428571 0.1250000
## 3  0.2000000 0.1666667 0.1428571
## 4  0.2500000 0.2000000 0.1666667
## 5  0.3333333 0.2500000 0.2000000
## 6  0.5000000 0.3333333 0.2500000
## 7  1.0000000 0.5000000 0.3333333
## 8  0.0000000 1.0000000 0.5000000
## 9  1.0000000 0.0000000 1.0000000
## 10 0.5000000 1.0000000 0.0000000

con estos datos ya es posible hacer el calculo del índice de moran de la siguiente manera

library(ape)
Moran.I((Materia.O), mo.dist.inv)
## $observed
## [1] -0.01161115
## 
## $expected
## [1] -0.006711409
## 
## $sd
## [1] 0.007678122
## 
## $p.value
## [1] 0.5233806

Para este caso, se puede observar que los datos presentan cierta dispersión en su distribución pero que se pueden clasificar como datos sin dependencia espacial o que la relacion es mínima.

Ejercicio de ejemplificación

Otra manera de demostar su fincionalidad es con un gran numero de datos, que para ete caso serán de conductividad eléctrica.

library(readxl)
BD_MORAN_1<- read_excel("C:/Users/Maiker/Downloads/BD_MORAN (1).xlsx", 
    col_types = c("numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric", 
        "numeric"))
View(BD_MORAN_1)
head(BD_MORAN_1, n=10)
## # A tibble: 10 x 7
##        z X_WGS84 Y_WGS84 CEa_075 CEa_150   X_MCB   Y_MCE
##    <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>   <dbl>
##  1  194.   -72.5    4.20    6.17    18.1 843499. 955943.
##  2  194.   -72.5    4.20    6.21    17.7 843499. 955943.
##  3  194.   -72.5    4.20    6.37    18.6 843499. 955944.
##  4  194.   -72.5    4.20    6.33    18.2 843498. 955944.
##  5  194.   -72.5    4.20    6.41    18.2 843498. 955944.
##  6  194.   -72.5    4.20    6.37    18.8 843498. 955944.
##  7  194.   -72.5    4.20    6.33    18.0 843498. 955944.
##  8  193.   -72.5    4.20    6.72    18.6 843497. 955944.
##  9  193.   -72.5    4.20    6.48    18.7 843497. 955944.
## 10  193.   -72.5    4.20    6.60    18.4 843497. 955944.
dim(BD_MORAN_1)
## [1] 18526     7

En esto tenemos que se trata de cerca de 18526 datos de conductividad electrica a los cuales va a ser necesario determinar si hay algún tipo de dependencia o relación entre ellos mismos, que para este caso son datos a 75cm y 150 cm.

sin embargo, para este caso se requiere dar una froma real de la parcela por lo que se hará lo siguiente:

plot(BD_MORAN_1$Y_WGS84, BD_MORAN_1$X_WGS84, cex=1, col=0.3*BD_MORAN_1$CEa_075, pch=16, ylab = "", xlab = "")

plot(BD_MORAN_1$Y_WGS84, BD_MORAN_1$X_WGS84, cex=1, col=0.3*BD_MORAN_1$CEa_150, pch=19, ylab="", xlab="")

donde las coordenadas de los datos nos permiten determinar los sitios donde fueron tomados y asi lograr hacer una presunción de la forma del lote.

Sin embargo el alto numero de datos no permite que se pueda hacer un calculo directo, por lo que es necesario dividirlos en 3 partes, de la siguiente manera.

library(readxl)
parte1 <- read_excel("C:/Users/Maiker/Downloads/BD_MORAN (1).xlsx", 
    sheet = "Hoja2", col_types = c("numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))

parte2 <- read_excel("C:/Users/Maiker/Downloads/BD_MORAN (1).xlsx", 
    sheet = "Hoja3", col_types = c("numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))

parte3 <- read_excel("C:/Users/Maiker/Downloads/BD_MORAN (1).xlsx", 
    sheet = "Hoja4", col_types = c("numeric", 
        "numeric", "numeric", "numeric", 
        "numeric", "numeric", "numeric"))

grafica1<-plot(parte1$Y_WGS84, parte1$X_WGS84, cex=1, col=0.3*parte1$CEa_075, pch=16, ylab = "", xlab = "")

grafica2<-plot(parte2$Y_WGS84, parte2$X_WGS84, cex=1, col=0.3*parte2$CEa_075, pch=16, ylab = "", xlab = "")

grafica3<-plot(parte3$Y_WGS84, parte3$X_WGS84, cex=1, col=0.3*parte3$CEa_075, pch=16, ylab = "", xlab = "")

para 75 cm

grafica1<-plot(parte1$Y_WGS84, parte1$X_WGS84, cex=1, col=0.3*parte1$CEa_150, pch=16, ylab = "", xlab = "")

grafica2<-plot(parte2$Y_WGS84, parte2$X_WGS84, cex=1, col=0.3*parte2$CEa_150, pch=16, ylab = "", xlab = "")

grafica3<-plot(parte3$Y_WGS84, parte3$X_WGS84, cex=1, col=0.3*parte3$CEa_150, pch=16, ylab = "", xlab = "")

de esta manera, podrtemos obtener 3 resultadods diferentes del índice de Moran como ya se habia estudiado anteriormente.

Primero se dan los valores invertidos espaciales, que en este caso podemos usar las coordenadas de los datos, para cada cojunto de datos.

ce.dists1<-as.matrix(dist(cbind(parte1$X_WGS84, parte1$Y_WGS84)))
ce.dists.inv1 <-1/ce.dists1
ce.dists.inv1[is.infinite(ce.dists.inv1)] <- 0
diag(ce.dists.inv1)<-0
ce.dists.inv1[1:6, 1:6]
##           1        2        3        4        5         6
## 1      0.00 406473.8 183582.7 136552.7  99512.0  80767.51
## 2 406473.84      0.0 334028.1 205286.9 131603.4 100690.59
## 3 183582.68 334028.1      0.0 532631.2 217163.1 144140.71
## 4 136552.66 205286.9 532631.2      0.0 366652.3 197620.35
## 5  99512.00 131603.4 217163.1 366652.3      0.0 428663.73
## 6  80767.51 100690.6 144140.7 197620.4 428663.7      0.00
ce.dists2<-as.matrix(dist(cbind(parte2$X_WGS84, parte2$Y_WGS84)))
ce.dists.inv2 <-1/ce.dists2
ce.dists.inv2[is.infinite(ce.dists.inv2)] <- 0
diag(ce.dists.inv2)<-0
ce.dists.inv2[1:6, 1:6]
##           1        2         3         4         5         6
## 1      0.00 222810.0 130001.76  91892.63  66665.19  54501.54
## 2 222810.01      0.0 312102.51 156393.25  95127.50  72149.30
## 3 130001.76 312102.5      0.00 313473.46 136833.86  93841.09
## 4  91892.63 156393.2 313473.46      0.00 242832.20 133927.38
## 5  66665.19  95127.5 136833.86 242832.20      0.00 298442.32
## 6  54501.54  72149.3  93841.09 133927.38 298442.32      0.00
ce.dists3<-as.matrix(dist(cbind(parte3$X_WGS84, parte3$Y_WGS84)))
ce.dists.inv3<-1/ce.dists3
ce.dists.inv3[is.infinite(ce.dists.inv3)] <- 0
diag(ce.dists.inv3)<-0
ce.dists.inv3[1:6, 1:6]
##           1         2         3         4         5         6
## 1      0.00 241412.12  88965.96  66895.51  47194.97  38773.39
## 2 241412.12      0.00 140885.18  92537.63  58663.31  46192.33
## 3  88965.96 140885.18      0.00 269655.43 100518.15  68725.52
## 4  66895.51  92537.63 269655.43      0.00 160256.00  92232.20
## 5  47194.97  58663.31 100518.15 160256.00      0.00 217288.13
## 6  38773.39  46192.33  68725.52  92232.20 217288.13      0.00

con esto calculoamos nuestro indice de Moran para cada profundidad para cada conjunto de datos

I.Moran1_75<-Moran.I(parte1$CEa_075, ce.dists.inv1); I.Moran1_75
## $observed
## [1] 0.5137397
## 
## $expected
## [1] -0.0001666944
## 
## $sd
## [1] 0.001569696
## 
## $p.value
## [1] 0
I.Moran2_75<-Moran.I(parte2$CEa_075, ce.dists.inv2); I.Moran2_75
## $observed
## [1] 0.6147387
## 
## $expected
## [1] -0.0001666667
## 
## $sd
## [1] 0.001341977
## 
## $p.value
## [1] 0
I.Moran3_75<-Moran.I(parte3$CEa_075, ce.dists.inv3); I.Moran3_75
## $observed
## [1] 0.4491044
## 
## $expected
## [1] -0.0001532802
## 
## $sd
## [1] 0.001303272
## 
## $p.value
## [1] 0

Y para los datos de 150 cm.

I.Moran1_150<-Moran.I(parte1$CEa_150, ce.dists.inv1); I.Moran1_150
## $observed
## [1] 0.3310082
## 
## $expected
## [1] -0.0001666944
## 
## $sd
## [1] 0.001569439
## 
## $p.value
## [1] 0
I.Moran2_150<-Moran.I(parte2$CEa_150, ce.dists.inv2); I.Moran2_150
## $observed
## [1] 0.3012206
## 
## $expected
## [1] -0.0001666667
## 
## $sd
## [1] 0.001341892
## 
## $p.value
## [1] 0
I.Moran3_150<-Moran.I(parte3$CEa_150, ce.dists.inv3); I.Moran3_150
## $observed
## [1] 0.4800928
## 
## $expected
## [1] -0.0001532802
## 
## $sd
## [1] 0.001303416
## 
## $p.value
## [1] 0

finalmente se hace un computo de los datos obtenidos en el p-valor a 75 cm y 10 cm para obetener conclusiones.

library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
bind_rows(I.Moran1_75, I.Moran2_75, I.Moran3_75, id=NULL)
## # A tibble: 3 x 4
##   observed  expected      sd p.value
##      <dbl>     <dbl>   <dbl>   <dbl>
## 1    0.514 -0.000167 0.00157       0
## 2    0.615 -0.000167 0.00134       0
## 3    0.449 -0.000153 0.00130       0
bind_rows(I.Moran1_150, I.Moran2_150, I.Moran3_150, id=NULL)
## # A tibble: 3 x 4
##   observed  expected      sd p.value
##      <dbl>     <dbl>   <dbl>   <dbl>
## 1    0.331 -0.000167 0.00157       0
## 2    0.301 -0.000167 0.00134       0
## 3    0.480 -0.000153 0.00130       0

Finalmente se puede asegurar que en todos los datose da un p-valor de 0