Modelos Estadisticos

#####paquetes
library(spdep)
library(ape)
library(sp)
library(MVA)
library(Hmisc)
library(normtest)
library(nortest)
library(spatialreg)

library(ggplot2)
library(Metrics)
library(plotly)
library(DT)
library(corrplot)
library(readxl)
BD_MODELADO <- read_excel("BD_MODELADO.xlsx")

data<-data.frame('X'=BD_MODELADO$Avg_X_MCB,'Y'=BD_MODELADO$Avg_Y_MCE,'CEa75'=BD_MODELADO$Avg_CEa_07,'CEa150'=BD_MODELADO$Avg_CEa_15,'NDVI'=BD_MODELADO$NDVI,'DEM'=BD_MODELADO$DEM,'SLOPE'=BD_MODELADO$SLOPE,'Z'=BD_MODELADO$Avg_z)


XYdata=as.matrix(data[,1:2])
distancias<-as.matrix(dist(XYdata))
distancias[1:5,1:5]
##           1         2        3        4        5
## 1  0.000000  5.175854 45.29328 18.50485 21.93578
## 2  5.175854  0.000000 40.37963 21.50151 17.42547
## 3 45.293281 40.379625  0.00000 60.05426 32.89910
## 4 18.504849 21.501506 60.05426  0.00000 29.47992
## 5 21.935776 17.425472 32.89910 29.47992  0.00000
max(distancias)
## [1] 853.0099
# 1' == 111.11 km

corp =rcorr(as.matrix(data[,3:8]), type = 'p')
cors=rcorr(as.matrix(data[,3:8]), type = 's')

mcorp <- corp$r
mcors <- cors$r

par(mfrow=c(1,2))
corrplot(mcorp,order="hclust",tl.col="black",main="pearson")
corrplot(mcors,order="hclust",tl.col="black",main="spearman")

Analisando el grafico de corelacion de pearson y spearman se puede concluir que para ambos casos la conductividad electrica a 75 cm de profundidad tiene una correlacion con la altura y la pendiente. Debido a que hay dos variables que miden la altura se observa que la variable Z tiene una mejor correlacion que DEM por lo tanto se asume que es una variable que a fururo puede ser descartada de las variables explicativas.

Adicionalmente se observa como la altura del terreno tiene una correlacio negativa con la CEa150 es decir, a media que se disminuye la variable Z aumenta la CEa en a 150 cm de prfuncidad. Por el contrario, para la CEa75 entre mayor sea Z mayores seran los valores de CEa

Modelos Espaciales

MODELO 1: Modelo autoregresivo puro

contnb=dnearneigh(coordinates(XYdata),d1=0,d2=380000,longlat = F) #Distancia entre poligonos
contnb
## Neighbour list object:
## Number of regions: 313 
## Number of nonzero links: 97656 
## Percentage nonzero weights: 99.68051 
## Average number of links: 312
dlist <- nbdists(contnb, XYdata)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W") #Estandarizando matriz de pesos

map=spautolm(CEa75~1,data=data,listw=Wve)  #### modelo
summary(map) #CEa75 tiene dependencia espacial
## 
## Call: spautolm(formula = CEa75 ~ 1, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -3.258254 -0.650679 -0.071829  0.824652  3.063002 
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   5.6941     5.5177   1.032   0.3021
## 
## Lambda: 0.98811 LR test value: 162.5 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.011876 
## 
## Log likelihood: -494.8231 
## ML residual variance (sigma squared): 1.347, (sigma: 1.1606)
## Number of observations: 313 
## Number of parameters estimated: 3 
## AIC: 995.65
map1=spautolm(CEa150~1,data=data,listw=Wve)  
summary(map1) #CEa150 tiene dependencia espacial
## 
## Call: spautolm(formula = CEa150 ~ 1, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.453255 -0.397645 -0.042934  0.322283  2.953512 
## 
## Coefficients: 
##             Estimate Std. Error z value  Pr(>|z|)
## (Intercept)  19.3151     1.5515   12.45 < 2.2e-16
## 
## Lambda: 0.97691 LR test value: 86.774 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.023 
## 
## Log likelihood: -304.7918 
## ML residual variance (sigma squared): 0.40168, (sigma: 0.63378)
## Number of observations: 313 
## Number of parameters estimated: 3 
## AIC: 615.58
res_map<- map$fit$residuals #Residuales "u" del modelo
library(spdep)
im_res_map<-moran.mc(res_map,Wve, nsim = 2000) #Si hay dependencia espacial

res_map1<- map1$fit$residuals #Residuales "u" del modelo
shapiro.test(res_map1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_map1
## W = 0.95729, p-value = 6.37e-08
cvm.test(res_map1) #Residuales no son normales
## 
##  Cramer-von Mises normality test
## 
## data:  res_map1
## W = 0.33539, p-value = 0.0001306
library(spdep)
im_res_map1<-moran.mc(res_map1,Wve, nsim = 2000) #Si hay dependencia espacial

\[Y = \lambda W Y + u\] \[\hat{CEa75} = 0.98811*W*CEa75 + 5.6941\\ \hat{CEa75}-0.98811*W*CEa75 = 5.6941\\ (I-0.98811*W)CEa75 = 5.6941\\ \hat{CEa75} = 5.6941(I-0.98811*W)^{-1}\\\]

\[\hat{CEa150} = 0.97691*W*CEa150 + 19.3151\\ \hat{CEa150}-0.97691*W*CEa150 = 19.3151\\ (I-0.97691*W)CEa150 = 19.3151\\ \hat{CEa150} = 19.3151(I-0.97691*W)^{-1}\\\]

CEa75_1=map$fit$fitted.values #estimado  
DF_map=data.frame(data$CEa75,CEa75_1)
colnames(DF_map) <- c("CEa75_o","CEa75_e")

CEa150_1=map1$fit$fitted.values #estimado 
DF_map1=data.frame(data$CEa150,CEa150_1)
colnames(DF_map1) <- c("CEa150_o","CEa150_e")


RMSE1=rmse(data$CEa75,CEa75_1);RMSE1
## [1] 1.16061
RMSE1_1=rmse(data$CEa150,CEa150_1);RMSE1_1
## [1] 0.6337842

CEa75

###### CEa75


fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$CEa75_o, opacity = 0.5))
fig # CEa75 Observada
## Warning: `arrange_()` is deprecated as of dplyr 0.7.0.
## Please use `arrange()` instead.
## See vignette('programming') for more help
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_warnings()` to see where this warning was generated.
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$CEa75_e, opacity = 0.5))
fig # CEa75 Estimada con modelo 1
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map$CEa75_e - DF_map$CEa75_o)*5, opacity = 0.5))
fig # Diferencia entre CEa75 estimada y observada modelo 1
library(ggplot2)
library(akima)
# Mapa de interpolacion lineal de CEa75 esperada con el modelo 1
interpol = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map$CEa75_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol)
contour(interpol, add = T)

plot(data$CEa75, DF_map$CEa75_e, cex=0.5, pch =19) #Grafico de CEa75 E vs O

CEa150

###### CEa150

fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map1$CEa150_o*0.5, opacity = 0.5))
fig # CEa150 Observada
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map1$CEa150_e*0.5, opacity = 0.5))
fig # CEa150 Estimada con modelo 1
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map1$CEa150_e - DF_map1$CEa150_o)*10, opacity = 0.5))
fig # Diferencia entre CEa150 estimada y observada modelo 1
# Mapa de interpolacion lineal de CEa150 esperada con el modelo 1
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map1$CEa150_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa150, DF_map1$CEa150_e, cex=0.5, pch =19) #Grafico de CEa150 E vs O

MODELO 2: Modelo Autorregresivo SARAR

CEa75

##### CEa75

map_sarar=sacsarlm(CEa75~1,data=data,listw = Wve)
## Warning in sacsarlm(CEa75 ~ 1, data = data, listw = Wve): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   reciprocal condition number = 2.88892e-20 - using numerical Hessian.
summary(map_sarar) # Rho y Lambda son significativos
## 
## Call:sacsarlm(formula = CEa75 ~ 1, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -3.405867 -0.620449 -0.028054  0.640915  2.891325 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)  -1.7267     3.2488 -0.5315   0.5951
## 
## Rho: 0.97863
## Approximate (numerical Hessian) standard error: 0.021304
##     z-value: 45.936, p-value: < 2.22e-16
## Lambda: 0.97863
## Approximate (numerical Hessian) standard error: 0.02128
##     z-value: 45.988, p-value: < 2.22e-16
## 
## LR test value: 254.56, p-value: < 2.22e-16
## 
## Log likelihood: -448.7933 for sac model
## ML residual variance (sigma squared): 0.98538, (sigma: 0.99267)
## Number of observations: 313 
## Number of parameters estimated: 4 
## AIC: 905.59, (AIC for lm: 1156.2)
CEa75_e2=map_sarar$fitted.values
res2=map_sarar$residuals
moran_msarar=moran.mc(res2, Wve,nsim=2000)
moran_msarar # Residuales son dependientes
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res2 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.10493, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
shapiro.test(res2) #Los residuales son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.99417, p-value = 0.2746
RMSE2=rmse(data$CEa75,CEa75_e2);RMSE2
## [1] 0.9926651
DF_map_sarar=data.frame(map_sarar$fitted.values,data$CEa75)
colnames(DF_map_sarar) <- c("CEa75_e","CEa75_o")

fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map_sarar$CEa75_e - DF_map_sarar$CEa75_o)*5, opacity = 0.5))
fig # Diferencia entre CEa75 estimada y observada modelo 2
# Mapa interpolacion lineal datos esperados de CEa75 modelo 2
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map_sarar$CEa75_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa75, DF_map_sarar$CEa75_e, cex=0.5, pch =19) #Grafico de CEa75 E vs O

CEa150

#### CEa150

map1_sarar=sacsarlm(CEa150~1,data=data,listw = Wve)
## Warning in sacsarlm(CEa150 ~ 1, data = data, listw = Wve): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   reciprocal condition number = 5.00737e-18 - using numerical Hessian.
summary(map1_sarar) # Rho y lambda son sigificativos
## 
## Call:sacsarlm(formula = CEa150 ~ 1, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.364908 -0.357504 -0.063033  0.289858  2.878760 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   1.1253     1.1735  0.9589   0.3376
## 
## Rho: 0.95895
## Approximate (numerical Hessian) standard error: 0.040829
##     z-value: 23.487, p-value: < 2.22e-16
## Lambda: 0.95895
## Approximate (numerical Hessian) standard error: 0.040764
##     z-value: 23.524, p-value: < 2.22e-16
## 
## LR test value: 133.68, p-value: < 2.22e-16
## 
## Log likelihood: -281.3411 for sac model
## ML residual variance (sigma squared): 0.34087, (sigma: 0.58384)
## Number of observations: 313 
## Number of parameters estimated: 4 
## AIC: 570.68, (AIC for lm: 700.36)
CEa150_e2=map1_sarar$fitted.values
res2_1=map1_sarar$residuals
moran_msarar1=moran.mc(res2_1, Wve,nsim=2000)
moran_msarar1 # Residuales son dependientes
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res2_1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057087, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
shapiro.test(res2_1) #Los residuales no son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res2_1
## W = 0.94498, p-value = 2.076e-09
RMSE2_1=rmse(data$CEa150,CEa150_e2);RMSE2_1
## [1] 0.5838387
DF_map1_sarar=data.frame(map1_sarar$fitted.values,data$CEa150)
colnames(DF_map1_sarar) <- c("CEa150_e","CEa150_o")

fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map1_sarar$CEa150_e - DF_map1_sarar$CEa150_o)*10, opacity = 0.5))
fig # Diferencia entre CEa150 estimada y observada modelo 2
# Mapa interpolacion lineal datos esperados de CEa150 modelo 2
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map1_sarar$CEa150_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa150, DF_map1_sarar$CEa150_e, cex=0.5, pch =19) #Grafico de CEa150 E vs O

MODELO 3: Modelo espacial de error

CEa75

# Modelo que involucra variables explicativas

## CEa75
mser1=errorsarlm(formula=CEa75~CEa150+NDVI+DEM+SLOPE+Z, data = data,listw = Wve)
summary(mser1) # Se descarta la variable NDVI
## 
## Call:errorsarlm(formula = CEa75 ~ CEa150 + NDVI + DEM + SLOPE + Z, 
##     data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.019160 -0.540466 -0.045367  0.513314  2.592838 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -64.737579   5.752902 -11.2530 < 2.2e-16
## CEa150        0.859898   0.083054  10.3535 < 2.2e-16
## NDVI         -2.395368   1.907913  -1.2555  0.209301
## DEM           0.036792   0.020974   1.7542  0.079402
## SLOPE        -0.073067   0.024760  -2.9510  0.003168
## Z             0.257034   0.028465   9.0299 < 2.2e-16
## 
## Lambda: 0.9825, LR test value: 99.359, p-value: < 2.22e-16
## Asymptotic standard error: 0.012342
##     z-value: 79.604, p-value: < 2.22e-16
## Wald statistic: 6336.8, p-value: < 2.22e-16
## 
## Log likelihood: -406.1005 for error model
## ML residual variance (sigma squared): 0.76603, (sigma: 0.87523)
## Number of observations: 313 
## Number of parameters estimated: 8 
## AIC: 828.2, (AIC for lm: 925.56)
mser2=errorsarlm(formula=CEa75~CEa150+SLOPE+Z+DEM, data = data,listw = Wve)
summary(mser2) #Se descarta la variable DEM
## 
## Call:errorsarlm(formula = CEa75 ~ CEa150 + SLOPE + Z + DEM, data = data, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.068942 -0.573110 -0.041672  0.535538  2.620533 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -66.334322   5.621356 -11.8004 < 2.2e-16
## CEa150        0.871288   0.082765  10.5273 < 2.2e-16
## SLOPE        -0.074849   0.024782  -3.0203  0.002525
## Z             0.251732   0.028220   8.9203 < 2.2e-16
## DEM           0.039380   0.020925   1.8819  0.059845
## 
## Lambda: 0.98246, LR test value: 98.998, p-value: < 2.22e-16
## Asymptotic standard error: 0.012369
##     z-value: 79.427, p-value: < 2.22e-16
## Wald statistic: 6308.6, p-value: < 2.22e-16
## 
## Log likelihood: -406.8867 for error model
## ML residual variance (sigma squared): 0.76989, (sigma: 0.87744)
## Number of observations: 313 
## Number of parameters estimated: 7 
## AIC: 827.77, (AIC for lm: 924.77)
mser3=errorsarlm(formula=CEa75~CEa150+SLOPE+Z, data = data,listw = Wve)
summary(mser3) #Modelo definitivo
## 
## Call:errorsarlm(formula = CEa75 ~ CEa150 + SLOPE + Z, data = data, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.150527 -0.558459 -0.045187  0.540349  2.578564 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -65.325177   5.620712 -11.622 < 2.2e-16
## CEa150        0.874324   0.083217  10.507 < 2.2e-16
## SLOPE        -0.079881   0.024777  -3.224  0.001264
## Z             0.286926   0.021256  13.498 < 2.2e-16
## 
## Lambda: 0.98237, LR test value: 97.514, p-value: < 2.22e-16
## Asymptotic standard error: 0.012433
##     z-value: 79.011, p-value: < 2.22e-16
## Wald statistic: 6242.7, p-value: < 2.22e-16
## 
## Log likelihood: -408.6476 for error model
## ML residual variance (sigma squared): 0.77863, (sigma: 0.8824)
## Number of observations: 313 
## Number of parameters estimated: 6 
## AIC: 829.3, (AIC for lm: 924.81)
DF_mser=data.frame(mser3$fitted.values,data$CEa75)
colnames(DF_mser) <- c("CEa75_e","CEa75_o")
RMSE3=rmse(data$CEa75,DF_mser$CEa75_e);RMSE3
## [1] 0.8824011
res3=mser3$residuals
shapiro.test(res3) #Los residuales son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res3
## W = 0.99348, p-value = 0.1948
moran_mee=moran.mc(res3, Wve,nsim=2000)
moran_mee #Residuales tiene dependencia espacial
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res3 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.12762, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_mser$CEa75_e - DF_mser$CEa75_o)*5, opacity = 0.5))
fig # Diferencia entre CEa75 estimada y observada modelo 3
# Mapa interpolacion lineal datos esperados de CEa75 modelo 3
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_mser$CEa75_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa75, DF_mser$CEa75_e, cex=0.5, pch =19) #Grafico de CEa75 E vs O

CEa150

## CEa150
mser4=errorsarlm(formula=CEa150~CEa75+NDVI+DEM+SLOPE+Z, data = data,listw = Wve)
summary(mser4) # Se descarta la variable NDVI y DEM
## 
## Call:errorsarlm(formula = CEa150 ~ CEa75 + NDVI + DEM + SLOPE + Z, 
##     data = data, listw = Wve)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.46790 -0.32241 -0.02153  0.36060  1.99578 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.5088605  2.7607167 16.4844 < 2.2e-16
## CEa75        0.2959247  0.0286368 10.3337 < 2.2e-16
## NDVI        -1.1627276  1.1220178 -1.0363 0.3000703
## DEM         -0.0093728  0.0123588 -0.7584 0.4482181
## SLOPE        0.0518339  0.0144655  3.5833 0.0003393
## Z           -0.1327355  0.0171932 -7.7202 1.155e-14
## 
## Lambda: 0.96877, LR test value: 49.814, p-value: 1.69e-12
## Asymptotic standard error: 0.022011
##     z-value: 44.013, p-value: < 2.22e-16
## Wald statistic: 1937.2, p-value: < 2.22e-16
## 
## Log likelihood: -239.4171 for error model
## ML residual variance (sigma squared): 0.26504, (sigma: 0.51482)
## Number of observations: 313 
## Number of parameters estimated: 8 
## AIC: 494.83, (AIC for lm: 542.65)
mser5=errorsarlm(formula=CEa150~CEa75+SLOPE+Z, data = data,listw = Wve)
summary(mser5) #Modelo definitivo
## 
## Call:errorsarlm(formula = CEa150 ~ CEa75 + SLOPE + Z, data = data, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.470171 -0.318709 -0.014481  0.355224  2.014963 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) 44.749718   2.699784  16.5753 < 2.2e-16
## CEa75        0.297488   0.028368  10.4868 < 2.2e-16
## SLOPE        0.052248   0.014422   3.6227 0.0002915
## Z           -0.143289   0.013322 -10.7558 < 2.2e-16
## 
## Lambda: 0.96928, LR test value: 50.495, p-value: 1.1945e-12
## Asymptotic standard error: 0.021655
##     z-value: 44.761, p-value: < 2.22e-16
## Wald statistic: 2003.5, p-value: < 2.22e-16
## 
## Log likelihood: -240.1762 for error model
## ML residual variance (sigma squared): 0.2663, (sigma: 0.51604)
## Number of observations: 313 
## Number of parameters estimated: 6 
## AIC: 492.35, (AIC for lm: 540.85)
DF_mser1=data.frame(mser5$fitted.values,data$CEa150)
colnames(DF_mser1) <- c("CEa150_e","CEa150_o")
RMSE3_1=rmse(data$CEa150,DF_mser1$CEa150_e);RMSE3_1
## [1] 0.5160437
res3_1=mser5$residuals
shapiro.test(res3_1) #Los residuales no son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res3_1
## W = 0.98814, p-value = 0.01166
moran_mee1=moran.mc(res3_1, Wve,nsim=2000)
moran_mee1 #Residuales tiene dependencia espacial
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res3_1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.076281, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_mser1$CEa150_e - DF_mser1$CEa150_o)*10, opacity = 0.5))
fig # Diferencia entre CEa150 estimada y observada modelo 3
# Mapa interpolacion lineal datos esperados de CEa150 modelo 3
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_mser1$CEa150_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa150, DF_mser1$CEa150_e, cex=0.5, pch =19) #Grafico de CEa150 E vs O

MODELO 4: Modelo con Lambda y Rho

CEa75

# CEa75
mlrho=sacsarlm(formula=CEa75~CEa150+SLOPE+Z, data = data,listw = Wve)
summary(mlrho) #Rho y Lambda son significativos
## 
## Call:sacsarlm(formula = CEa75 ~ CEa150 + SLOPE + Z, data = data, listw = Wve)
## 
## Residuals:
##         Min          1Q      Median          3Q         Max 
## -2.11409261 -0.48405766 -0.00093027  0.51350442  2.20507301 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -59.437374  14.726420 -4.0361 5.435e-05
## CEa150        0.857117   0.073143 11.7183 < 2.2e-16
## SLOPE        -0.065991   0.021538 -3.0639  0.002185
## Z             0.213388   0.028874  7.3904 1.463e-13
## 
## Rho: 0.97498
## Asymptotic standard error: 0.37591
##     z-value: 2.5937, p-value: 0.0094956
## Lambda: 0.97212
## Asymptotic standard error: 0.41932
##     z-value: 2.3183, p-value: 0.020432
## 
## LR test value: 179.22, p-value: < 2.22e-16
## 
## Log likelihood: -367.7943 for sac model
## ML residual variance (sigma squared): 0.58887, (sigma: 0.76738)
## Number of observations: 313 
## Number of parameters estimated: 7 
## AIC: 749.59, (AIC for lm: 924.81)
res4=mlrho$residuals
shapiro.test(res4) #Los residuales son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res4
## W = 0.99548, p-value = 0.4995
moran_mlrho=moran.mc(res4, Wve,nsim=2000)
moran_mlrho #Residuales tiene dependencia espacial
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res4 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.091618, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
DF_mlrho=data.frame(mlrho$fitted.values,data$CEa75)
colnames(DF_mlrho) <- c("CEa75_e","CEa75_o")
RMSE4=rmse(data$CEa75,DF_mlrho$CEa75_e);RMSE4
## [1] 0.7673774
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_mlrho$CEa75_e - DF_mlrho$CEa75_o)*5, opacity = 0.5))
fig # Diferencia entre CEa75 estimada y observada modelo 4
# Mapa interpolacion lineal datos esperados de CEa75 modelo 4
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_mlrho$CEa75_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa75, DF_mlrho$CEa75_e, cex=0.5, pch =19) #Grafico de CEa75 E vs O

\[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \\ u=\rho W u + \epsilon\] #### CEa150

# CEa150
mlrho1=sacsarlm(formula=CEa150~CEa75+SLOPE+Z, data = data,listw = Wve)
summary(mlrho1) #Rho no es significativo y Lambda si es significativo
## 
## Call:sacsarlm(formula = CEa150 ~ CEa75 + SLOPE + Z, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.443696 -0.278907 -0.020386  0.310970  1.959313 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 21.463941  15.438933  1.3902  0.164454
## CEa75        0.290649   0.032187  9.0301 < 2.2e-16
## SLOPE        0.041795   0.013504  3.0949  0.001969
## Z           -0.114901   0.017473 -6.5759 4.835e-11
## 
## Rho: 0.94903
## Asymptotic standard error: 0.54667
##     z-value: 1.736, p-value: 0.082562
## Lambda: 0.95738
## Asymptotic standard error: 0.45839
##     z-value: 2.0886, p-value: 0.036747
## 
## LR test value: 87.934, p-value: < 2.22e-16
## 
## Log likelihood: -221.4568 for sac model
## ML residual variance (sigma squared): 0.23287, (sigma: 0.48257)
## Number of observations: 313 
## Number of parameters estimated: 7 
## AIC: 456.91, (AIC for lm: 540.85)
res4_1=mlrho1$residuals
shapiro.test(res4_1) #Los residuales no son normales
## 
##  Shapiro-Wilk normality test
## 
## data:  res4_1
## W = 0.98084, p-value = 0.0003408
moran_mlrho1=moran.mc(res4_1, Wve,nsim=2000)
moran_mlrho1 #Residuales tiene dependencia espacial
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res4_1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057956, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
DF_mlrho1=data.frame(mlrho1$fitted.values,data$CEa150)
colnames(DF_mlrho1) <- c("CEa150_e","CEa150_o")
RMSE4_1=rmse(data$CEa150,DF_mlrho1$CEa150_e);RMSE4_1
## [1] 0.4825702
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_mlrho1$CEa150_e - DF_mlrho1$CEa150_o)*10, opacity = 0.5))
fig # Diferencia entre CEa150 estimada y observada modelo 4
# Mapa interpolacion lineal datos esperados de CEa150 modelo 4
interpol1 = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_mlrho1$CEa150_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol1)
contour(interpol1, add = T)

plot(data$CEa150, DF_mlrho1$CEa150_e, cex=0.5, pch =19) #Grafico de CEa150 E vs O

\[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \\ \]

AIC

Modelos=c('Modelo1','Modelo2','Modelo3','Modelo4','Modelo1','Modelo2','Modelo3','Modelo4')
Variable=c(rep("CEa75",4),rep("CEa150",4))
AIC=c(995.65,905.59,829.3,749.59,615.58,570.68,492.35,456.91)
RMSE=c(RMSE1,RMSE2,RMSE3,RMSE4,RMSE1_1,RMSE2_1,RMSE3_1,RMSE4_1)

Tabla1=data.frame(Modelos,Variable,AIC,(round(RMSE,2)))
Tabla1_1<-datatable(Tabla1, class = 'cell-border stripe', editable = 'cell', rownames = FALSE, colnames = c("Modelos", "Variable", "AIC", "RMSE"));Tabla1_1

De acuerdo a los analisis estadisticos realizados, tanto para la CEa75 y CEa150 el modelo 4 es el que mejor se ajusta obtniendo valores mas bajos de AIC y RMSE lo que nos indica que ese modelo pude estimar de mejor forma la conductividad electrica en ambas profuncidades segun los datos obtenidos en campo. En ambos casos, se lego a la conclusion de que las variables DEM y NDVI no mostraban una dependencia espacial relacionada con la conductividad electrica por lo cual no se tomaron en cuenta al momento de analizar las variables explicativas del modelo.

En cuanto a las variables que si teneian dependencia espacila como lo es la profundidad a la que se tomó la CEa se observa que a una mayor profundida la CEa es mayor como se observo en los estudios de (Castro Franco, et al.) en Buenos Aires, donde el CEa era inversamnete proporcional con contenido de MO. DE igual forma, la altura tiene una relacion con el CEa siendo inversamente proporcional aumentanto a medida que se disminuye la altura (Cicore et al., 2016)

Estimacion de CEa en una coordenada no muestreada

cc = chull(XYdata)
cc = c(cc, cc[1])
plot(XYdata, pch = 16)
lines(XYdata[cc,], type = 'l')
### Grafico del punto en el que se quiere estimar la CEa75
points(x=843750, y=956280, col="red", pch=19)

n_data<-data.frame(843750, 956280, 0,0,0,0,0,0)
names(n_data)<-c("X","Y","CEa75", "CEa150","NDVI","DEM", "SLOPE","Z")
n<-rbind(data,n_data)
XYdata1=as.matrix(n[,1:2])

n_1<-as.matrix(dist(cbind(n$X, n$Y)))
invn_1<- 1/n_1
invn_1[is.infinite(invn_1)] <- 0
mid<-diag(314)
nsumas<-apply(invn_1,1,sum)
# Matriz de Pesos estandarizada
nnwe<-invn_1/nsumas

CEa75 estimada

p1nn<-0.98811*nnwe
p2nn<-mid-p1nn
p3nn<-solve(p2nn) # Matriz inversa
p4nn<-5.6941*p3nn #Matriz con los resultados finales
# dim(p4nn) Para comprobar las dimensiones de la matriz (314 observaciones)
p4nn[314,314]
## [1] 7.347418

CEa150 estimada

p1nn1<-0.97691*nnwe
p2nn1<-mid-p1nn1
p3nn1<-solve(p2nn1) # Matriz inversa
p4nn1<-19.3151*p3nn1 #Matriz con los resultados finales
# dim(p4nn) Para comprobar las dimensiones de la matriz (314 observaciones)
p4nn1[314,314]
## [1] 22.15947

Referencias

-Castro Franco, M., & Costa, J. L. (2012). CONDUCTIVIDAD ELÉCTRICA APARENTE Y SU RELACIÓN CON PROPIEDADES DEL SUELO, PARA LA DELIMITACIÓN DE ZONAS PARA MANEJO SITIO ESPECÍFICO.Tomado de: https://www.researchgate.net/publication/299388042_CONDUCTIVIDAD_ELECTRICA_APARENTE_Y_SU_RELACION_CON_PROPIEDADES_DEL_SUELO_PARA_LA_DELIMITACION_DE_ZONAS_PARA_MANEJO_SITIO_ESPECIFICO

-Cicore, P., Sánchez, H., Peralta, N., Aparicio, V., Castro, F., Costa, J. (2016). Utilización de la conductividad electrica aparente y la elevación para delimitar ambientes edaficos en suelos ganaderos.Tomado de: http://www.ipni.net/publication/ia-lahp.nsf/0/D658DCB73C43294E8525808E00758AB0/$FILE/Art%204.pdf