library(spdep)
library(ape)
library(sp)
library(MVA)
library(Hmisc)
library(normtest)
library(nortest)
library(spatialreg)
library(asbio)
library(readxl)
library(plotly)
library(readxl)
library(data.table)
library(ape)
library(corrplot)
library(Hmisc)
library(psych)
library(goftest)
library(moments)
library(Metrics)
library(akima)
BD_MODELADO <- read_excel("C:/Users/User/Downloads/BD_MODELADO.xlsx")

dataBD=data.frame(BD_MODELADO)
head(dataBD,10)
##    Avg_X_MCB Avg_Y_MCE Avg_CEa_07 Avg_CEa_15     NDVI      DEM    SLOPE
## 1   843449.6  955962.0   7.237480   18.02656 0.863030 199.0000 6.385167
## 2   843454.8  955962.4   6.787250   18.02737 0.866502 197.1667 1.981082
## 3   843493.6  955951.3   6.848250   18.70444 0.874883 197.0000 0.577682
## 4   843439.6  955977.6   7.135162   18.34237 0.845838 197.0000 1.175075
## 5   843468.7  955972.8   6.826763   17.92409 0.797179 197.0000 0.210996
## 6   843500.6  955975.3   6.699966   18.39441 0.758272 197.6667 4.357386
## 7   843525.7  955982.6   6.180742   17.84332 0.763436 199.7500 6.628445
## 8   843413.4  956014.0   8.539024   18.75812 0.823320 197.1667 1.462050
## 9   843433.6  956013.4   8.869958   18.85396 0.759923 197.3333 1.663344
## 10  843468.8  956010.1   7.231308   18.34269 0.757382 197.6667 3.541936
##       Avg_z
## 1  193.0512
## 2  193.2986
## 3  193.5659
## 4  194.4116
## 5  193.9931
## 6  195.3814
## 7  196.6780
## 8  194.9936
## 9  196.1356
## 10 197.8522

##VARIABLES RESPUESTA

v_res75 <- dataBD[,3]
v_res150 <- dataBD[,4]

##VARIABLES EXPLICATVIAS (independientes)

v_exp75 <- dataBD[,4:8]
v_exp150 <- dataBD[,c(3,5:8)]
m75 <- as.matrix(v_exp75)
m150 <- as.matrix(v_exp150)

COORDENADAS

xycor <- as.matrix(dataBD[,1:2])

MATRIZ DE PESOS

MD <- as.matrix(dist(xycor, diag = T, upper = T))
inv_MD <- as.matrix(1/MD)
diag(inv_MD) <- 0
W <- as.matrix(inv_MD)

Estandarizado

sma <- apply(W, 1, sum)
WE <- W/sma  
 EST =apply(WE, 1, sum)
 head(EST,20)
##  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15 16 17 18 19 20 
##  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1

Indice de moran

IM_C75 <- Moran.I(v_res75, WE)
IM_C75$p.value
## [1] 0
IM_C150 <- Moran.I(v_res150, WE)
IM_C150$p.value
## [1] 0
IM_CNDVI <- Moran.I(dataBD[,5], WE)
IM_CNDVI$p.value
## [1] 0
IM_DEM  <- Moran.I(dataBD[,6], WE)
IM_DEM$p.value
## [1] 0
IM_SLOPE <- Moran.I(dataBD[,7], WE)
IM_SLOPE$p.value
## [1] 0
IM_Z <- Moran.I(dataBD[,8], WE)
IM_Z$p.value
## [1] 0

lo anterior lo podemos simplificar con un ciclo for

dataBD1=dataBD[,3:8]
Moran_p=list()
for(i in 1:6){
  Moran_p[i]=(Moran.I(dataBD[,i], WE))$p.value}
Moran_o=list()
for(i in 1:6){
Moran_o[i]=(Moran.I(dataBD[,i], WE))$observed}
IM_PO = data.frame(cbind("Variables"=c('CE a 70cm','CE a 150cm','NDVI','DEM','Slope',"Z (elevation)"),"Valor-o"=unlist( Moran_o), "P-valor"= unlist(Moran_p)))

head(IM_PO,6)
##       Variables            Valor.o P.valor
## 1     CE a 70cm  0.454063254612529       0
## 2    CE a 150cm  0.430058691816291       0
## 3          NDVI  0.268746828627044       0
## 4           DEM  0.160951018381243       0
## 5         Slope 0.0975040302526088       0
## 6 Z (elevation)  0.309670845207254       0

Se observa dependencia espacial en todas las variables

##Matrices de correlacion

MCP <- rcorr(as.matrix(dataBD[,3:8]), type = 'p')
MCS <- rcorr(as.matrix(dataBD[,3:8]), type = 's')
Mcorp <- MCP$r
Mcors <- MCS$r

##grafico metodo de Pearson

corrplot.mixed(Mcorp, order = 'hclust', tl.col = 'darkgreen', tl.cex = 0.5,  lower.col = "black", number.cex = .5,upper= "square")

##grafico metodo de Spearman

corrplot.mixed(Mcors, order = 'hclust', tl.col = 'darkgreen', tl.cex = 0.5, lower.col = "black", number.cex = .5,upper= "square")

Las mayores relaciones se dan entre DEM y Avg_z, con una relacion considerable de Avg_z con CEa_07. CEa_15 se relaciona conAvg_z y en menor medida con DEM. DEM y Avg_z son cordenadas de altura

describe(dataBD[,3:8])
##            vars   n   mean   sd median trimmed  mad    min    max range  skew
## Avg_CEa_07    1 313   9.77 1.53   9.65    9.76 1.40   6.17  13.60  7.44  0.10
## Avg_CEa_15    2 313  18.50 0.74  18.44   18.47 0.64  16.80  21.51  4.71  0.57
## NDVI          3 313   0.83 0.03   0.84    0.84 0.02   0.70   0.88  0.17 -1.39
## DEM           4 313 205.11 4.55 204.83  205.13 5.68 196.17 213.33 17.17  0.06
## SLOPE         5 313   4.13 2.16   3.78    3.94 2.10   0.21  12.72 12.51  0.95
## Avg_z         6 313 202.47 3.66 202.49  202.57 4.34 193.05 212.35 19.30 -0.14
##            kurtosis   se
## Avg_CEa_07    -0.28 0.09
## Avg_CEa_15     1.22 0.04
## NDVI           2.48 0.00
## DEM           -1.13 0.26
## SLOPE          1.25 0.12
## Avg_z         -0.54 0.21

##conductividad electrica 75 cm

bx_75 <- ggplot(dataBD, aes(x=0, y=Avg_CEa_07)) + 
  geom_boxplot(color="black", fill="darkgreen") +
  
  xlim(-1,1)

hist_75 <-ggplot(dataBD, aes(x=Avg_CEa_07)) + 
  geom_histogram(color="black", fill="darkgreen", binwidth=0.6)

bx_75 

hist_75

##conductividad electrica 150 cm

bx_150 <- ggplot(dataBD, aes(x=0, y=Avg_CEa_15)) + 
  geom_boxplot(color="black", fill="darkblue") +
  
  xlim(-1,1)

hist_150 <-ggplot(dataBD, aes(x=Avg_CEa_15)) + 
  geom_histogram(color="black", fill="darkblue", binwidth=0.6)

bx_150 

hist_150 

##75 cm
cv_75<-cvm.test(dataBD[,3])
cv_75
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  dataBD[, 3]
## omega2 = 104.33, p-value < 2.2e-16
##150 cm
cv_150<-cvm.test(dataBD[,4])
cv_150
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  dataBD[, 4]
## omega2 = 104.33, p-value < 2.2e-16
##75 cm
sh_75<-shapiro.test(dataBD[,3])
sh_75
## 
##  Shapiro-Wilk normality test
## 
## data:  dataBD[, 3]
## W = 0.99217, p-value = 0.09827
##150 cm
sh_150<-shapiro.test(dataBD[,4])
sh_150
## 
##  Shapiro-Wilk normality test
## 
## data:  dataBD[, 4]
## W = 0.97785, p-value = 9.283e-05
ks.test(dataBD[,3], dataBD[,4])
## 
##  Two-sample Kolmogorov-Smirnov test
## 
## data:  dataBD[, 3] and dataBD[, 4]
## D = 1, p-value < 2.2e-16
## alternative hypothesis: two-sided

###MODELACION DE LOS DATOS ##NO ESPACIAL

mod_cla75 <- lm(Avg_CEa_07~SLOPE+Avg_z+Avg_CEa_15,data= dataBD)
mod_cla150 <- lm(Avg_CEa_15~SLOPE+Avg_z+Avg_CEa_07,data= dataBD)

mod_cla75_log <- lm(Avg_CEa_07~SLOPE+Avg_z+log(Avg_CEa_15),data= dataBD)
mod_cla150_log <- lm(log(Avg_CEa_15)~SLOPE+Avg_z+log(Avg_CEa_07),data= dataBD)

summary(mod_cla75)
## 
## Call:
## lm(formula = Avg_CEa_07 ~ SLOPE + Avg_z + Avg_CEa_15, data = dataBD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.41298 -0.71115 -0.05541  0.64834  3.01718 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -73.98404    4.74765 -15.583  < 2e-16 ***
## SLOPE        -0.12525    0.02799  -4.475 1.07e-05 ***
## Avg_z         0.33680    0.01834  18.369  < 2e-16 ***
## Avg_CEa_15    0.86887    0.09270   9.373  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.05 on 309 degrees of freedom
## Multiple R-squared:  0.5315, Adjusted R-squared:  0.527 
## F-statistic: 116.9 on 3 and 309 DF,  p-value: < 2.2e-16
summary(mod_cla150)
## 
## Call:
## lm(formula = Avg_CEa_15 ~ SLOPE + Avg_z + Avg_CEa_07, data = dataBD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.50608 -0.39071  0.01281  0.36509  2.04703 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 47.55651    2.11792  22.454  < 2e-16 ***
## SLOPE        0.07595    0.01503   5.053 7.45e-07 ***
## Avg_z       -0.15733    0.01123 -14.009  < 2e-16 ***
## Avg_CEa_07   0.25480    0.02718   9.373  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.5686 on 309 degrees of freedom
## Multiple R-squared:  0.4107, Adjusted R-squared:  0.405 
## F-statistic: 71.78 on 3 and 309 DF,  p-value: < 2.2e-16
summary(mod_cla75_log)
## 
## Call:
## lm(formula = Avg_CEa_07 ~ SLOPE + Avg_z + log(Avg_CEa_15), data = dataBD)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -2.42707 -0.71701 -0.05763  0.63705  3.00432 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)     -104.75300    7.56476 -13.847  < 2e-16 ***
## SLOPE             -0.12589    0.02811  -4.479 1.06e-05 ***
## Avg_z              0.33663    0.01844  18.259  < 2e-16 ***
## log(Avg_CEa_15)   16.07153    1.74114   9.230  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.054 on 309 degrees of freedom
## Multiple R-squared:  0.5284, Adjusted R-squared:  0.5238 
## F-statistic: 115.4 on 3 and 309 DF,  p-value: < 2.2e-16
summary(mod_cla150_log)
## 
## Call:
## lm(formula = log(Avg_CEa_15) ~ SLOPE + Avg_z + log(Avg_CEa_07), 
##     data = dataBD)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.083486 -0.020916 -0.000046  0.019431  0.102855 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)      4.3365073  0.1049644  41.314  < 2e-16 ***
## SLOPE            0.0039929  0.0007972   5.009 9.24e-07 ***
## Avg_z           -0.0085806  0.0006020 -14.254  < 2e-16 ***
## log(Avg_CEa_07)  0.1330214  0.0139088   9.564  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.03024 on 309 degrees of freedom
## Multiple R-squared:  0.4194, Adjusted R-squared:  0.4137 
## F-statistic: 74.39 on 3 and 309 DF,  p-value: < 2.2e-16

##Residuos

res1 <- mod_cla75$residuals
shapiro.test(res1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res1
## W = 0.99092, p-value = 0.05044
hist(res1, col = "darkred")

res1.2 <- mod_cla150$residuals
shapiro.test(res1.2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res1.2
## W = 0.99144, p-value = 0.06673
hist(res1.2, col = "darkblue")

res2 <- mod_cla75_log$residuals
shapiro.test(res2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.99096, p-value = 0.05142
hist(res2, col = "darkgreen")

res2.1 <- mod_cla150_log$residuals
shapiro.test(res2.1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res2.1
## W = 0.99471, p-value = 0.3564
hist(res2.1, col = "yellow")

Un modelo considerado adecuado, debe cumplir con una distribucion de residuos normal es decir media 0 y varianza mínima

##simetria de los datos

skew(dataBD[,3])
## [1] 0.1002837
skewness.norm.test(dataBD[,3])
## 
##  Skewness test for normality
## 
## data:  dataBD[, 3]
## T = 0.10077, p-value = 0.4575
skew(dataBD[,4])
## [1] 0.5702493
skewness.norm.test(dataBD[,4])
## 
##  Skewness test for normality
## 
## data:  dataBD[, 4]
## T = 0.57299, p-value < 2.2e-16

Los datos en la variable CE a 150cm no son simétricos

##Simetria para los residuales

skew(res1)
## [1] 0.2441868
skewness.norm.test(res1)
## 
##  Skewness test for normality
## 
## data:  res1
## T = 0.24536, p-value = 0.079
skew(res1.2)
## [1] 0.257621
skewness.norm.test(res1.2)
## 
##  Skewness test for normality
## 
## data:  res1.2
## T = 0.25886, p-value = 0.055
skew(res2)
## [1] 0.2440019
skewness.norm.test(res2)
## 
##  Skewness test for normality
## 
## data:  res2
## T = 0.24518, p-value = 0.0785
skew(res2)
## [1] 0.2440019
skewness.norm.test(res2.1)
## 
##  Skewness test for normality
## 
## data:  res2.1
## T = 0.1436, p-value = 0.291

##Graficos de valores observados vs estimados

VS1 <- mod_cla75$fitted.values
VS1.2 <- mod_cla150$fitted.values
VS2 <- mod_cla75_log$fitted.values
VS2.2 <- mod_cla150_log$fitted.values


plot(dataBD$Avg_CEa_07,VS1,pch=16, col = 'darkgreen')

plot(dataBD$Avg_CEa_07,VS1.2,pch=16, col = 'darkgreen')

plot(dataBD$Avg_CEa_15,VS2,pch=16, col = 'darkred')

plot(dataBD$Avg_CEa_15,VS2.2,pch=16, col = 'darkred')

cor(dataBD$Avg_CEa_07,VS1)
## [1] 0.7290603
cor(dataBD$Avg_CEa_07,VS1.2)
## [1] 0.01700518
cor(dataBD$Avg_CEa_15,VS2)
## [1] 0.00948048
cor(dataBD$Avg_CEa_15,VS2.2)
## [1] 0.6480663
res_c1 <- mod_cla75$residuals
res_c2 <- mod_cla150$residuals

shapiro.test(res_c1)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_c1
## W = 0.99092, p-value = 0.05044
shapiro.test(res_c2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_c2
## W = 0.99144, p-value = 0.06673
cvm.test(res_c1)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_c1
## omega2 = 25.075, p-value < 2.2e-16
cvm.test(res_c2)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_c2
## omega2 = 33.859, p-value < 2.2e-16
moranres1 <- Moran.I(res_c1, WE)
moranres1
## $observed
## [1] 0.1580117
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004665226
## 
## $p.value
## [1] 0
moranres2 <- Moran.I(res_c2, WE)
moranres2
## $observed
## [1] 0.08362486
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004659635
## 
## $p.value
## [1] 0
AIC(mod_cla75)
## [1] 924.809
AIC(mod_cla75_log)
## [1] 926.9123
AIC(mod_cla150)
## [1] 540.8476
AIC(mod_cla150_log)
## [1] -1295.924

#Muestran dependencia espacial (p-value=0)

###Modelo autoregresivo puro \[Y = \lambda W Y + u\]

contnb <- dnearneigh(coordinates(xycor),0,854, longlat = F)
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, xycor)
dlist <- lapply(dlist, function(x) 1/x)
Wve <- nb2listw(contnb, glist = dlist, style = 'W')

##Modelo aplicado a CEa 75cm

map_75 <- spautolm(Avg_CEa_07~1, data = dataBD, listw = Wve) 
summary(map_75)
## 
## Call: spautolm(formula = Avg_CEa_07 ~ 1, data = dataBD, 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.011866 
## 
## 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
CE75E1 <- as.data.frame(map_75$fit['fitted.values'])
head(CE75E1,10)
##    fitted.values
## 1       8.630840
## 2       8.677614
## 3       8.915327
## 4       8.905593
## 5       8.823590
## 6       8.855964
## 7       8.926651
## 8       9.075587
## 9       9.005442
## 10      8.976890
CE75E <- map_75$fit$fitted.values ## Estimados
head(CE75E,10)
##        1        2        3        4        5        6        7        8 
## 8.630840 8.677614 8.915327 8.905593 8.823590 8.855964 8.926651 9.075587 
##        9       10 
## 9.005442 8.976890
df75 <- data.frame(v_res75, CE75E)
colnames(df75) <-  c('CE_obs','CE_est')
plot(df75$CE_obs, df75$CE_est, cex=0.5, pch =19, col = "darkgreen")

resmap1 <- map_75$fit$residuals
cor(df75$CE_obs, df75$CE_est)
## [1] 0.7977199
plot(xycor[,1]  ,xycor[,2], col = floor(abs(resmap1))+3, pch =20)

plot(xycor[,1]  ,xycor[,2], cex =abs(resmap1), pch =20)

plot(xycor[,1]  ,xycor[,2], cex =0.1*df75$CE_obs, pch =20)

data2 <- data.frame(xycor[,1]  ,xycor[,2], df75$CE_obs, df75$CE_est)
colnames(data2) <- c('x', 'y', 'CE_obs', 'CE_est')
plot1<-ggplot(data = data2, aes(xycor[,1]  ,xycor[,2])) +
  geom_point(cex = data2$CE_obs*0.2) +
  geom_point(color = data2$CE_est)
 
plot1

im_res_map_75 <- moran.mc(resmap1,Wve,nsim = 2000);im_res_map_75
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resmap1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.16722, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

###Modelo aplicado a CE 150cm

map_150 <- spautolm(Avg_CEa_15~1, data =dataBD, listw = Wve) 
summary(map_150)
## 
## Call: spautolm(formula = Avg_CEa_15 ~ 1, data = dataBD, 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
CE150E <- as.data.frame(map_150$fit['fitted.values'])
head(CE150E,10)
##    fitted.values
## 1       18.52785
## 2       18.52287
## 3       18.57280
## 4       18.59249
## 5       18.58692
## 6       18.58429
## 7       18.60574
## 8       18.66210
## 9       18.64986
## 10      18.64010
CE150E <- map_150$fit$fitted.values ## Estimados
head(CE150E,10)
##        1        2        3        4        5        6        7        8 
## 18.52785 18.52287 18.57280 18.59249 18.58692 18.58429 18.60574 18.66210 
##        9       10 
## 18.64986 18.64010
df150 <- data.frame(v_res150, CE150E)
colnames(df150) <-  c('CE_obs','CE_est')
plot(df150$CE_obs, df150$CE_est,col="darkblue", cex=0.8, pch =20)

resmap2 <- map_150$fit$residuals
cor(df150$CE_obs, df150$CE_est)
## [1] 0.6642671
plot(xycor[,1]  ,xycor[,2], col = floor(abs(resmap2))+3, pch =20)

plot(xycor[,1]  ,xycor[,2], cex =abs(resmap2), pch =20)

plot(xycor[,1]  ,xycor[,2], cex =0.1*df150$CE_obs, pch =20)

data2_1 <- data.frame(xycor[,1]  ,xycor[,2], df75$CE_obs, df150$CE_est)
colnames(data2_1) <- c('x', 'y', 'CE_obs', 'CE_est')
plot2<-ggplot(data = data2_1, aes(xycor[,1]  ,xycor[,2])) +
  geom_point(cex = data2_1$CE_obs*0.2) +
  geom_point(color = data2_1$CE_est)
 
plot2

im_res_map_150 <- moran.mc(resmap2,Wve,nsim = 2000);im_res_map_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resmap2 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.094365, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

###Modelo espacial del error \[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \]

##Conductividad electrica a 75cm

mod_es_er_1 <- errorsarlm(Avg_CEa_07~NDVI+SLOPE+Avg_z+Avg_CEa_15+DEM,data=dataBD,listw=Wve)

summary(mod_es_er_1)
## 
## Call:errorsarlm(formula = Avg_CEa_07 ~ NDVI + SLOPE + Avg_z + Avg_CEa_15 + 
##     DEM, data = dataBD, 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
## NDVI         -2.395368   1.907913  -1.2555  0.209301
## SLOPE        -0.073067   0.024760  -2.9510  0.003168
## Avg_z         0.257034   0.028465   9.0299 < 2.2e-16
## Avg_CEa_15    0.859898   0.083054  10.3535 < 2.2e-16
## DEM           0.036792   0.020974   1.7542  0.079402
## 
## 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)
mod_es_er_2 <- errorsarlm(Avg_CEa_07~SLOPE+Avg_z+Avg_CEa_15+DEM,data=dataBD,listw=Wve)

summary(mod_es_er_2)
## 
## Call:errorsarlm(formula = Avg_CEa_07 ~ SLOPE + Avg_z + Avg_CEa_15 + 
##     DEM, data = dataBD, 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
## SLOPE        -0.074849   0.024782  -3.0203  0.002525
## Avg_z         0.251732   0.028220   8.9203 < 2.2e-16
## Avg_CEa_15    0.871288   0.082765  10.5273 < 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)
mod_es_er_3 <- errorsarlm(Avg_CEa_07~SLOPE+Avg_z+Avg_CEa_15,data=dataBD,listw=Wve)

summary(mod_es_er_3)
## 
## Call:errorsarlm(formula = Avg_CEa_07 ~ SLOPE + Avg_z + Avg_CEa_15, 
##     data = dataBD, 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
## SLOPE        -0.079881   0.024777  -3.224  0.001264
## Avg_z         0.286926   0.021256  13.498 < 2.2e-16
## Avg_CEa_15    0.874324   0.083217  10.507 < 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)

#el mejor AIC fue en el modelo 2 ##Normalidad de residuales

res_mod_es_er <- mod_es_er_2$residuals
shapiro.test(res_mod_es_er)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod_es_er
## W = 0.99235, p-value = 0.1078
cvm.test(res_mod_es_er)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_mod_es_er
## omega2 = 26.54, p-value < 2.2e-16

Los residuales son normales

plot(df75$CE_obs, mod_es_er_2$fitted.values, cex=0.8, pch =20, col = "darkred")

cor(df75$CE_obs, mod_es_er_2$fitted.values)
## [1] 0.8246131
moran_error_75 <- moran.mc(res_mod_es_er,Wve,nsim=2000)
moran_error_75
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_mod_es_er 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.12895, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

conductividad electrica a 150cm

mod_es_er_4 <- errorsarlm(Avg_CEa_15~NDVI+SLOPE+Avg_z+Avg_CEa_07+DEM,data=dataBD,listw=Wve)
summary(mod_es_er_4)
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ NDVI + SLOPE + Avg_z + Avg_CEa_07 + 
##     DEM, data = dataBD, 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.5088606  2.7607169 16.4844 < 2.2e-16
## NDVI        -1.1627276  1.1220178 -1.0363 0.3000703
## SLOPE        0.0518339  0.0144655  3.5833 0.0003393
## Avg_z       -0.1327355  0.0171932 -7.7202 1.155e-14
## Avg_CEa_07   0.2959247  0.0286368 10.3337 < 2.2e-16
## DEM         -0.0093728  0.0123588 -0.7584 0.4482181
## 
## Lambda: 0.96877, LR test value: 49.814, p-value: 1.69e-12
## Asymptotic standard error: 0.022011
##     z-value: 44.014, 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)
mod_es_er_5 <- errorsarlm(Avg_CEa_15~SLOPE+Avg_z+Avg_CEa_07+DEM,data=dataBD,listw=Wve)
summary(mod_es_er_5)
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ SLOPE + Avg_z + Avg_CEa_07 + 
##     DEM, data = dataBD, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.477728 -0.316994 -0.014091  0.367283  2.009858 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.035550   2.729375 16.5003 < 2.2e-16
## SLOPE        0.051310   0.014481  3.5432 0.0003953
## Avg_z       -0.136386   0.016857 -8.0910 6.661e-16
## Avg_CEa_07   0.299387   0.028492 10.5079 < 2.2e-16
## DEM         -0.008240   0.012332 -0.6682 0.5040078
## 
## Lambda: 0.96901, LR test value: 50.337, p-value: 1.2949e-12
## Asymptotic standard error: 0.021843
##     z-value: 44.363, p-value: < 2.22e-16
## Wald statistic: 1968.1, p-value: < 2.22e-16
## 
## Log likelihood: -239.9531 for error model
## ML residual variance (sigma squared): 0.26594, (sigma: 0.51569)
## Number of observations: 313 
## Number of parameters estimated: 7 
## AIC: 493.91, (AIC for lm: 542.24)
mod_es_er_6 <- errorsarlm(Avg_CEa_15~SLOPE+Avg_z+Avg_CEa_07,data=dataBD,listw=Wve)
summary(mod_es_er_6)
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ SLOPE + Avg_z + Avg_CEa_07, 
##     data = dataBD, 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
## SLOPE        0.052248   0.014422   3.6227 0.0002915
## Avg_z       -0.143289   0.013322 -10.7558 < 2.2e-16
## Avg_CEa_07   0.297488   0.028368  10.4868 < 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)

#el mejor modelo por AIC 5

res_mod_es_er_2 <- mod_es_er_5$residuals
shapiro.test(res_mod_es_er_2)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mod_es_er_2
## W = 0.98823, p-value = 0.01217
cvm.test(res_mod_es_er_2)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_mod_es_er_2
## omega2 = 38.142, p-value < 2.2e-16
plot(df150$CE_obs, mod_es_er_6$fitted.values, cex=0.8, pch =20, col = 'darkblue')

cor(df150$CE_obs, mod_es_er_6$fitted.values)
## [1] 0.718156
moran_error_150 <- moran.mc(res_mod_es_er_2,Wve,nsim=2000)
moran_error_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_mod_es_er_2 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.07571, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

###Modelo lambda y Rho

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

\[u=\pWu+ \e \] ##conductividad eletrica a 75 cm

mlrho <- sacsarlm(Avg_CEa_07~SLOPE+Avg_z+Avg_CEa_15+DEM,data=dataBD,listw=Wve)
summary(mlrho)
## 
## Call:sacsarlm(formula = Avg_CEa_07 ~ SLOPE + Avg_z + Avg_CEa_15 + 
##     DEM, data = dataBD, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.090770 -0.476848 -0.031738  0.518306  2.235748 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -60.138673  15.556907 -3.8657 0.0001108
## SLOPE        -0.062334   0.021580 -2.8885 0.0038705
## Avg_z         0.187865   0.031034  6.0535 1.417e-09
## Avg_CEa_15    0.854968   0.072736 11.7543 < 2.2e-16
## DEM           0.028548   0.018830  1.5160 0.1295090
## 
## Rho: 0.97458
## Asymptotic standard error: 0.38031
##     z-value: 2.5626, p-value: 0.01039
## Lambda: 0.9722
## Asymptotic standard error: 0.41632
##     z-value: 2.3352, p-value: 0.019532
## 
## LR test value: 179.63, p-value: < 2.22e-16
## 
## Log likelihood: -366.569 for sac model
## ML residual variance (sigma squared): 0.58432, (sigma: 0.76441)
## Number of observations: 313 
## Number of parameters estimated: 8 
## AIC: 749.14, (AIC for lm: 924.77)
res_mlrho <- mlrho$residuals
summary(res_mlrho)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -2.09077 -0.47685 -0.03174  0.00000  0.51831  2.23575
shapiro.test(res_mlrho)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mlrho
## W = 0.99543, p-value = 0.4903
cvm.test(res_mlrho)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_mlrho
## omega2 = 27.131, p-value < 2.2e-16
plot(df75$CE_obs, mlrho$fitted.values, cex=0.8, pch =20,col="yellow")

cor(df75$CE_obs, mlrho$fitted.values)
## [1] 0.8703349
moran_mlrho <- moran.mc(res_mlrho,Wve,nsim=2000)
moran_mlrho
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_mlrho 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.09234, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

El modelo λ y ρ para este modelo se ajusta bien, ya que ρ es significativo (p value=0.01039) ##conductividad eletrica a 150 cm

mlrho_150 <- sacsarlm(Avg_CEa_15~SLOPE+Avg_z+Avg_CEa_07,data=dataBD,listw=Wve)
summary(mlrho_150)
## 
## Call:
## sacsarlm(formula = Avg_CEa_15 ~ SLOPE + Avg_z + Avg_CEa_07, data = dataBD, 
##     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.438935  1.3902  0.164454
## SLOPE        0.041795   0.013504  3.0949  0.001969
## Avg_z       -0.114901   0.017473 -6.5759 4.835e-11
## Avg_CEa_07   0.290649   0.032187  9.0301 < 2.2e-16
## 
## 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)
res_mlrho_150 <- mlrho_150$residuals
summary(res_mlrho_150)
##     Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
## -1.44370 -0.27891 -0.02039  0.00000  0.31097  1.95931
shapiro.test(res_mlrho_150)
## 
##  Shapiro-Wilk normality test
## 
## data:  res_mlrho_150
## W = 0.98084, p-value = 0.0003408
cvm.test(res_mlrho_150)
## 
##  Cramer-von Mises test of goodness-of-fit
##  Null hypothesis: uniform distribution
##  Parameters assumed to be fixed
## 
## data:  res_mlrho_150
## omega2 = 41.963, p-value < 2.2e-16
plot(df75$CE_obs, mlrho_150$fitted.values, cex=0.8, pch =20,col="salmon")

cor(df75$CE_obs, mlrho_150$fitted.values)
## [1] 0.06931742
moran_mlrho_150 <- moran.mc(res_mlrho_150,Wve,nsim=2000)
moran_mlrho_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_mlrho_150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057956, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater

Este modelo no puede modelar CEa 150cm ya que el p value=0.082 ###Valores observados vs. estimados ##Modelo Clásico

rmse(dataBD$Avg_CEa_07, VS1)#75cm
## [1] 1.04333
rmse(dataBD$Avg_CEa_15, VS1.2)#150cm
## [1] 0.5649946

##Modelo autorregresivo puro

rmse(df75$CE_obs, df75$CE_est)#75cm
## [1] 1.16061
rmse(df150$CE_obs, df150$CE_est)#150cm
## [1] 0.6337842

##Modelo del error

rmse(df75$CE_obs, mod_es_er_2$fitted.values)#75cm
## [1] 0.8774362
rmse(df150$CE_obs, mod_es_er_6$fitted.values)#150cm
## [1] 0.5160437

##Modelo Lambda y RHo

rmse(df75$CE_obs, mlrho$fitted.values)#75cm
## [1] 0.7644114
rmse(df150$CE_obs, mlrho_150$fitted.values)#150cm
## [1] 0.4825702

####Gráficas

fig_Z <- plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$Avg_z,data=dataBD, size = I(90),
                  marker = list(color=rgb(0.1,0.3,0.6),
                                line = list(color = rgb(0.1,0.3,0.6))))%>%
            layout(title = "Z vs XY",
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Elevation (M.A.S.L.)")
    )
  )%>%
add_markers(color = "cyan")
fig_Z
fig_DEM <- plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$DEM,data=dataBD, size = I(90),
                   marker = list(color=rgb(0.1,0.3,0.6),
                                line = list(color = rgb(0.1,0.3,0.6))))%>%
            layout(title = "DEM vs XY",
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "DEM")
    )
  )%>%
add_markers(color = "cyan")
fig_DEM

se observan diferencias en la altitud una posible causa la vegetacion de la zona

###Gráficas de conductividad electrica a 75cm con variables relacionadas

fig_Slope <- plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$SLOPE ,data=dataBD, size = I(90))%>%
            layout(title = 'CE_75 vs Slope',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Slope")
    )
  )%>%
add_markers(color = dataBD$Avg_CEa_07)
fig_Slope
fig_150 <-  plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$Avg_CEa_15 ,data=dataBD, size = I(90))%>%
            layout(
                  scene = list(title = 'CEa_150 vs Slope',
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "CEa_150")
    )
  )%>%
add_markers(color = dataBD$Avg_CEa_15)
fig_150
fig_Z_1 <- plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$Avg_z,data=dataBD, size = I(90))%>%
            layout(title = 'CEa_75 vs Altitud (z)',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Altitud")
    )
  )%>%
add_markers(color = dataBD$Avg_CEa_07)
fig_Z_1
fig_DEM_1 <- plot_ly(x=dataBD$Avg_X_MCB  ,y=dataBD$Avg_Y_MCE,z= dataBD$DEM,data=dataBD, size = I(90))%>%
            layout(title = 'CEa_150 vs Altitud (z)',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "DEM")
    )
  )%>%
add_markers(color = dataBD$Avg_CEa_15)
fig_DEM_1

###Coordenada de interés ##AUTOREGRESIVO PURO

plot(dataBD$Avg_X_MCB,dataBD$Avg_Y_MCE, main="CONTORNO CONVEXO", col="darkred",pch=20)
points(843750,956280,col="darkgreen",pch=19,cex=1)

#matriz de datos modificada
plot3 <- data.frame(843750,956280,0,0,0,0,0,0)
names(plot3) <- c("Avg_X_MCB","Avg_Y_MCE","Avg_CEa_07","Avg_CEa_15","NDVI","DEM","SLOPE","Avg_z")
nueva_W <- rbind(dataBD,plot3)
nueva_XYcor_W <- as.matrix(nueva_W[,1:2])
nueva_count <- dnearneigh(coordinates(nueva_XYcor_W),0,854,longlat = F)
nueva_dlist <- nbdists(nueva_count,nueva_XYcor_W)
nueva_dlist <- lapply(nueva_dlist,function(x)1/x)
#matroz de pesos 
nueva_Wve <- nb2listw(nueva_count,glist = nueva_dlist,style = 'W')
#matriz de indentidad
matrix_id <- diag(314)

##ARticulo en el articulo Villafañe, R., Abarca, O., Azpúrua, M., Ruiz, T., & Dugarte, J. (1999). Distribución espacial de la salinidad en los suelos de Quíbor y su relación con las limitaciones de drenaje y la calidad del agua. Bioagro, 11(2), 43-50. y Mata-Fernández, I., Rodríguez-Gamiño, M., López-Blanco, J., & Vela-Correa, G. (2014). Dinámica de la salinidad en los suelos. Revista Digital del Departamento El Hombre y su Ambiente, 1(5), 26-35. se denota claramente la relacion de la salinidad del suelo con su distribucion espacial

https://d1wqtxts1xzle7.cloudfront.net/34148267/1._Distribucion_espacial.pdf?1404817253=&response-content-disposition=inline%3B+filename%3DDISTRIBUCION_ESPACIAL_DE_LA_SALINIDAD_EN.pdf&Expires=1607408826&Signature=GVO4yTAPGYTbYaz5Ar~lh8Z-KjPSX-G1gup75KXBcY9NAMxVciwtZ0U0mhzMgnT6x42o-9j~AnY5VgP15PyFkTB9fco5mNT3eJG2NYBzqukBsF5y2bhJRl6DLi8p2AjiF8DLk9efl1oZ1wDl4fCLZFf~ljCt~PneJkAC0Zdjvt4g3HDfrMjp6P1nKAHFRWBrn9fUCkCyOQxQQFOOJyiCNxAr3hDeanKiNrCvB1dpThtr3-H0J41m0FOc8BppZODdRzQ3lXkbibr~KezAJxxFcreYPXkov94BzF1KVLGwYxAyJOsiF4nHr9ABa6DGjp60VSGmLvrxXOtN1sFgZG~eMA__&Key-Pair-Id=APKAJLOHF5GGSLRBV4ZA

http://cbs1.xoc.uam.mx/e_bios/docs/2014/05_SALINIDAD_EN_SUELOS_ESPANOL.pdf