Indice de moran
dt
5/10/2020
library(clhs)
df <- data.frame(
X = runif(1000),
Y = rnorm(1000),
Hz = sample(seq(1:5), size = 1000, replace = TRUE))
dim(df)## [1] 1000 3# Returning the indices of the sampled points
muestra50 <- clhs(df, size = 50, iter = 100, progress = FALSE, simple = TRUE)
muestra50## [1] 697 341 144 14 797 133 884 997 796 905 41 889 930 919 838 466 60 957 188
## [20] 381 7 379 275 815 26 456 752 599 167 693 646 568 349 73 16 820 677 22
## [39] 68 97 813 668 625 137 3 494 790 533 319 979df50=df[muestra50,]#Dataframe con solo la muestraplot(df$X,df$Y,col=df$Hz,pch=19,cex=0.5)
points(df50$X,df50$Y,pch=13,cex=1)
#Datos de ozono
ozone <- read.table("https://stats.idre.ucla.edu/stat/r/faq/ozone.csv", sep=",", header=T)
head(ozone, n=10)## Station Av8top Lat Lon
## 1 60 7.225806 34.13583 -117.9236
## 2 69 5.899194 34.17611 -118.3153
## 3 72 4.052885 33.82361 -118.1875
## 4 74 7.181452 34.19944 -118.5347
## 5 75 6.076613 34.06694 -117.7514
## 6 84 3.157258 33.92917 -118.2097
## 7 85 5.201613 34.01500 -118.0597
## 8 87 4.717742 34.06722 -118.2264
## 9 88 6.532258 34.08333 -118.1069
## 10 89 7.540323 34.38750 -118.5347str(ozone)## 'data.frame': 32 obs. of 4 variables:
## $ Station: int 60 69 72 74 75 84 85 87 88 89 ...
## $ Av8top : num 7.23 5.9 4.05 7.18 6.08 ...
## $ Lat : num 34.1 34.2 33.8 34.2 34.1 ...
## $ Lon : num -118 -118 -118 -119 -118 ...plot(ozone$Lon,ozone$Lat,cex=0.3*ozone$Av8top,col="lightblue",pch=19)
plot(ozone$Lon,ozone$Lat,cex=1,col=0.3*ozone$Av8top,pch=19)
library(ape)
ozone.dists <- as.matrix(dist(cbind(ozone$Lon, ozone$Lat)))#Matriz de distancias
head(ozone.dists,n=5)## 1 2 3 4 5 6 7
## 1 0.0000000 0.3937326 0.4088031 0.6144127 0.1854888 0.3529452 0.1820077
## 2 0.3937326 0.0000000 0.3749446 0.2206810 0.5743590 0.2685580 0.3021020
## 3 0.4088031 0.3749446 0.0000000 0.5116772 0.4994034 0.1078698 0.2301238
## 4 0.6144127 0.2206810 0.5116772 0.0000000 0.7944601 0.4226993 0.5095533
## 5 0.1854888 0.5743590 0.4994034 0.7944601 0.0000000 0.4785934 0.3126778
## 8 9 10 11 12 13 14
## 1 0.3104545 0.1907020 0.6609031 0.5373218 0.4950538 0.5926597 0.3157506
## 2 0.1405634 0.2280588 0.3046982 0.1870421 0.2585831 0.5452384 0.5374451
## 3 0.2466955 0.2719279 0.6622189 0.3503443 0.2088277 0.8681962 0.2750202
## 4 0.3354875 0.4432559 0.1880560 0.1690395 0.3208746 0.6338945 0.7286474
## 5 0.4750001 0.3559325 0.8463845 0.7029626 0.6358365 0.7309451 0.2946199
## 15 16 17 18 19 20 21 22
## 1 0.2115536 0.5593270 0.4607347 1.409834 0.5128292 0.8036861 1.0682343 1.741442
## 2 0.4414920 0.8316243 0.6339081 1.801976 0.9042258 1.1698189 1.4606363 2.131946
## 3 0.2573998 0.5344991 0.2999415 1.644656 0.7836674 0.9514534 1.3150254 1.956648
## 4 0.6441959 1.0199047 0.8027289 2.022089 1.1243602 1.3818795 1.6808890 2.351367
## 5 0.2446098 0.4436728 0.4293460 1.227654 0.3299031 0.6220967 0.8864384 1.557873
## 23 24 25 26 27 28 29
## 1 0.7591125 0.6123126 0.2962018 0.6589204 0.4210195 0.6506294 0.7753133
## 2 1.1137395 1.0028927 0.6899308 1.0438697 0.8147263 1.0439424 1.1690234
## 3 0.8791815 0.8682796 0.6246085 1.0058960 0.7369042 0.9568847 1.0642379
## 4 1.3224455 1.2227520 0.9106118 1.2618939 1.0353954 1.2644785 1.3896906
## 5 0.5826760 0.4286014 0.1276036 0.5094913 0.2493491 0.4794728 0.6000001
## 30 31 32
## 1 0.4378266 0.7458638 0.07266991
## 2 0.8250787 1.1300249 0.46500673
## 3 0.6937060 1.0910280 0.46427212
## 4 1.0443273 1.3476445 0.68558759
## 5 0.2524623 0.5957770 0.12617659dim(ozone.dists)## [1] 32 32ozone.dists.inv <- 1/ozone.dists #inverso de la matriz de las distancias
head(ozone.dists.inv,n=5)## 1 2 3 4 5 6 7 8
## 1 Inf 2.539795 2.446165 1.627570 5.391160 2.833301 5.494272 3.221084
## 2 2.539795 Inf 2.667061 4.531428 1.741071 3.723590 3.310140 7.114228
## 3 2.446165 2.667061 Inf 1.954357 2.002389 9.270438 4.345487 4.053580
## 4 1.627570 4.531428 1.954357 Inf 1.258716 2.365748 1.962503 2.980737
## 5 5.391160 1.741071 2.002389 1.258716 Inf 2.089456 3.198180 2.105263
## 9 10 11 12 13 14 15 16
## 1 5.243784 1.513081 1.861082 2.019982 1.687309 3.167057 4.726934 1.7878628
## 2 4.384834 3.281936 5.346391 3.867229 1.834060 1.860655 2.265047 1.2024660
## 3 3.677445 1.510075 2.854335 4.788636 1.151813 3.636097 3.885007 1.8709107
## 4 2.256033 5.317565 5.915779 3.116482 1.577549 1.372406 1.552323 0.9804838
## 5 2.809521 1.181496 1.422551 1.572731 1.368092 3.394204 4.088144 2.2539131
## 17 18 19 20 21 22 23
## 1 2.170446 0.7093032 1.9499670 1.2442669 0.9361242 0.5742369 1.3173278
## 2 1.577516 0.5549462 1.1059185 0.8548332 0.6846331 0.4690551 0.8978760
## 3 3.333983 0.6080300 1.2760516 1.0510237 0.7604416 0.5110782 1.1374216
## 4 1.245751 0.4945380 0.8893947 0.7236521 0.5949233 0.4252845 0.7561748
## 5 2.329124 0.8145620 3.0311934 1.6074672 1.1281099 0.6419007 1.7162196
## 24 25 26 27 28 29 30
## 1 1.6331527 3.376076 1.5176340 2.3751868 1.5369733 1.2898011 2.2840091
## 2 0.9971157 1.449421 0.9579739 1.2274061 0.9579072 0.8554149 1.2120055
## 3 1.1517028 1.601003 0.9941386 1.3570285 1.0450580 0.9396396 1.4415329
## 4 0.8178274 1.098163 0.7924596 0.9658146 0.7908399 0.7195846 0.9575542
## 5 2.3331701 7.836768 1.9627422 4.0104418 2.0856242 1.6666665 3.9609872
## 31 32
## 1 1.3407274 13.760853
## 2 0.8849363 2.150506
## 3 0.9165668 2.153909
## 4 0.7420354 1.458603
## 5 1.6784804 7.925400diag(ozone.dists.inv) <- 0 #pasar la diagona de infinito a cero
ozone.dists.inv[1:5, 1:5]## 1 2 3 4 5
## 1 0.000000 2.539795 2.446165 1.627570 5.391160
## 2 2.539795 0.000000 2.667061 4.531428 1.741071
## 3 2.446165 2.667061 0.000000 1.954357 2.002389
## 4 1.627570 4.531428 1.954357 0.000000 1.258716
## 5 5.391160 1.741071 2.002389 1.258716 0.000000Moran.I(ozone$Av8top, ozone.dists.inv) #indice de moran## $observed
## [1] 0.2265501
##
## $expected
## [1] -0.03225806
##
## $sd
## [1] 0.03431138
##
## $p.value
## [1] 4.596323e-14set.seed(12345)
MO=rnorm(n = 150,mean = 3,sd = 0.5)
xy=expand.grid(x=seq(1,10),y=seq(1,15))
head(as.data.frame(MO),n=10)## MO
## 1 3.292764
## 2 3.354733
## 3 2.945348
## 4 2.773251
## 5 3.302944
## 6 2.091022
## 7 3.315049
## 8 2.861908
## 9 2.857920
## 10 2.540339plot(xy,col=MO,pch=19)
xymatrix<- as.matrix(dist(cbind(xy$x, xy$y)))
xymatrix.inv <- 1/xymatrix
diag(xymatrix.inv) <- 0
xymatrix.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.0000000Moran.I(MO, xymatrix.inv)## $observed
## [1] -0.009650003
##
## $expected
## [1] -0.006711409
##
## $sd
## [1] 0.007694112
##
## $p.value
## [1] 0.7025151##datos reales
library(readxl)
BD_MODELADO <- BD_MORAN_1_ <- read_excel("C:/Users/User/Downloads/BD_MORAN (1).xlsx",
n_max = 3000)lat<- BD_MODELADO $Y_WGS84
lon<- BD_MODELADO $X_WGS84
ce75cm<-BD_MODELADO $CEa_075
# por peso de lso vectores y capacidad de mi pc trabaje con 3000 datos
dfprof75<-data.frame(lat,lon,ce75cm)
matdatos75<- as.matrix(dist(cbind(dfprof75$lon, dfprof75$lat)))
dim(matdatos75)## [1] 3000 3000#Generando la matriz de pesos (inversa)
invmatdatos75<- 1/matdatos75
diag(invmatdatos75)<-0
dim(invmatdatos75)## [1] 3000 3000invmatdatos75[1:5,1:5]## 1 2 3 4 5
## 1 0.0 406473.8 183582.7 136552.7 99512.0
## 2 406473.8 0.0 334028.1 205286.9 131603.4
## 3 183582.7 334028.1 0.0 532631.2 217163.1
## 4 136552.7 205286.9 532631.2 0.0 366652.3
## 5 99512.0 131603.4 217163.1 366652.3 0.0invmatdatos75[is.infinite(invmatdatos75)] <- 0indice de Moran
library(ape)
Moran75<-Moran.I(dfprof75$ce75cm,invmatdatos75)
Moran75$observed## [1] 0.6078195el indice de Moran es de 0.6 , los datos de conductividad electrica aparente a 75 cm presentan autocorrelación espacial positiva . Se rechaza la hipotesis nula de que los datos presentan autocorrelación espacial cero
##Profundidad de 150 cm
lat2<- BD_MODELADO $Y_WGS84
lon2<- BD_MODELADO $X_WGS84
ce150cm<-BD_MODELADO$CEa_150
# por peso de lso vectores y capacidad de mi pc trabaje con 3000 datos
dfprof150<-data.frame(lat2,lon2,ce150cm)
matdatos150<- as.matrix(dist(cbind(dfprof150$lon2, dfprof150$lat2)))
dim(matdatos150)## [1] 3000 3000#Generando la matriz de pesos (inversa)
invmatdatos150<- 1/matdatos150
invmatdatos150[is.infinite(invmatdatos150)] <- 0
dim(invmatdatos75)## [1] 3000 3000invmatdatos150[1:5,1:5]## 1 2 3 4 5
## 1 0.0 406473.8 183582.7 136552.7 99512.0
## 2 406473.8 0.0 334028.1 205286.9 131603.4
## 3 183582.7 334028.1 0.0 532631.2 217163.1
## 4 136552.7 205286.9 532631.2 0.0 366652.3
## 5 99512.0 131603.4 217163.1 366652.3 0.0indice de Moran
library(ape)
Moran150<-Moran.I(dfprof150$ce150cm,invmatdatos150)
Moran150$observed## [1] 0.4883378el indice de Moran es de 0.48 , los datos de conductividad electrica aparente a 150 cm presentan autocorrelación espacial positiva . Se rechaza la hipotesis nula de que los datos presentan autocorrelación espacial cero
Mapa de Conductividad Electrica Aparente (Para un maximo de 3000 puntos)
Latitud<- BD_MODELADO$Y_WGS84
Longitud<- BD_MODELADO$X_WGS84
plot(Longitud,Latitud,col=0.3*BD_MODELADO$CEa_150, main = "Conductividad Electrica Aparente a 150 cm",cex.main=0.8,cex.lab=0.8,cex.axis=0.8)
plot(Longitud,Latitud,col=0.3*BD_MODELADO$CEa_075, main = "Conductividad Electrica Aparente a 75 cm",cex.main=0.8,cex.lab=0.8,cex.axis=0.8)
### Mapa de Conductividad Electrica Aparente
BD_MODELADO2 <- BD_MORAN_1_ <- read_excel("C:/Users/User/Downloads/BD_MORAN (1).xlsx")
par(mfrow = c(1, 2))
Latitud<- BD_MODELADO2$Y_WGS84
Longitud<- BD_MODELADO2$X_WGS84
plot(Longitud,Latitud,col=0.3*BD_MODELADO2$CEa_150, main = "Conductividad Electrica Aparente a 150 cm",cex.main=0.8,cex.lab=0.8,cex.axis=0.8)
plot(Longitud,Latitud,col=0.3*BD_MODELADO2$CEa_075, main = "Conductividad Electrica Aparente a 75 cm",cex.main=0.8,cex.lab=0.8,cex.axis=0.8)