#preparando librerias
library(spdep)
## Loading required package: sp
## Loading required package: spData
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
## Loading required package: sf
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(ape)
## Registered S3 method overwritten by 'ape':
##   method   from 
##   plot.mst spdep
library(sp)
library(MVA)
## Loading required package: HSAUR2
## Loading required package: tools
library(Hmisc)
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
## Loading required package: ggplot2
## 
## Attaching package: 'Hmisc'
## The following object is masked from 'package:ape':
## 
##     zoom
## The following objects are masked from 'package:base':
## 
##     format.pval, units
library(normtest)
library(nortest)
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:Hmisc':
## 
##     src, summarize
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(spatialreg)
## Loading required package: Matrix
## Registered S3 methods overwritten by 'spatialreg':
##   method                   from 
##   residuals.stsls          spdep
##   deviance.stsls           spdep
##   coef.stsls               spdep
##   print.stsls              spdep
##   summary.stsls            spdep
##   print.summary.stsls      spdep
##   residuals.gmsar          spdep
##   deviance.gmsar           spdep
##   coef.gmsar               spdep
##   fitted.gmsar             spdep
##   print.gmsar              spdep
##   summary.gmsar            spdep
##   print.summary.gmsar      spdep
##   print.lagmess            spdep
##   summary.lagmess          spdep
##   print.summary.lagmess    spdep
##   residuals.lagmess        spdep
##   deviance.lagmess         spdep
##   coef.lagmess             spdep
##   fitted.lagmess           spdep
##   logLik.lagmess           spdep
##   fitted.SFResult          spdep
##   print.SFResult           spdep
##   fitted.ME_res            spdep
##   print.ME_res             spdep
##   print.lagImpact          spdep
##   plot.lagImpact           spdep
##   summary.lagImpact        spdep
##   HPDinterval.lagImpact    spdep
##   print.summary.lagImpact  spdep
##   print.sarlm              spdep
##   summary.sarlm            spdep
##   residuals.sarlm          spdep
##   deviance.sarlm           spdep
##   coef.sarlm               spdep
##   vcov.sarlm               spdep
##   fitted.sarlm             spdep
##   logLik.sarlm             spdep
##   anova.sarlm              spdep
##   predict.sarlm            spdep
##   print.summary.sarlm      spdep
##   print.sarlm.pred         spdep
##   as.data.frame.sarlm.pred spdep
##   residuals.spautolm       spdep
##   deviance.spautolm        spdep
##   coef.spautolm            spdep
##   fitted.spautolm          spdep
##   print.spautolm           spdep
##   summary.spautolm         spdep
##   logLik.spautolm          spdep
##   print.summary.spautolm   spdep
##   print.WXImpact           spdep
##   summary.WXImpact         spdep
##   print.summary.WXImpact   spdep
##   predict.SLX              spdep
## 
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
## 
##     anova.sarlm, as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
##     as_dsCMatrix_IrW, as_dsTMatrix_listw, bptest.sarlm, can.be.simmed,
##     cheb_setup, coef.gmsar, coef.sarlm, coef.spautolm, coef.stsls,
##     create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm,
##     deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
##     errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm,
##     fitted.SFResult, fitted.spautolm, get.ClusterOption,
##     get.coresOption, get.mcOption, get.VerboseOption,
##     get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone,
##     gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts,
##     Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
##     lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm,
##     LR1.sarlm, LR1.spautolm, LU_prepermutate_setup, LU_setup,
##     Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
##     mom_calc_int2, moments_setup, powerWeights, predict.sarlm,
##     predict.SLX, print.gmsar, print.ME_res, print.sarlm,
##     print.sarlm.pred, print.SFResult, print.spautolm, print.stsls,
##     print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
##     print.summary.stsls, residuals.gmsar, residuals.sarlm,
##     residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup,
##     SE_interp_setup, SE_whichMin_setup, set.ClusterOption,
##     set.coresOption, set.mcOption, set.VerboseOption,
##     set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup,
##     SpatialFiltering, spautolm, spBreg_err, spBreg_lag, spBreg_sac,
##     stsls, subgraph_eigenw, summary.gmsar, summary.sarlm,
##     summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlm
library(readxl)

df= read_excel("C:/Users/Valentina Osorio/Documents/U.N/Computacion Estadistica/BD_MODELADO.xlsx")

Visualizacion de datos

View(df)

Mapeo de conductividad electrica aparente a 150 cm

library(ggplot2)
df %>% ggplot(aes(x = Avg_X_MCB, y=Avg_Y_MCE, colour=Avg_CEa_15))+
   geom_point(size = 5,shape=15)+
   scale_color_continuous(type = 'viridis')

Se observa dependencia espacial de la conductividad electrica a una profundidad de 150 cm debido a que existen regiones de un mismo color

Matriz de pesos

library(ape)
# Matriz de distancias 
df.dists <- as.matrix(dist(cbind(df$Avg_X_MCB, df$Avg_Y_MCE)))
# Inversa de las matriz 
df.dists.inv <- 1/df.dists
# Asignar ceros a la diagonal 
diag(df.dists.inv) <- 0
df.dists.inv <- round(df.dists.inv,3)
#Matriz de pesos estandarizada 
We<-df.dists.inv/rowSums(df.dists.inv)
View(We)

Indice de Moran

# Indice de Moran para CE a 150 cm
IM=Moran.I(df$Avg_CEa_15,df.dists.inv)
IM$p.value
## [1] 0

Regresión lineal no espacial

mod1=lm(df$Avg_CEa_15~df$NDVI)
summary(mod1)
## 
## Call:
## lm(formula = df$Avg_CEa_15 ~ df$NDVI)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.71991 -0.45614 -0.02841  0.38067  2.87345 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)   22.227      1.179   18.85  < 2e-16 ***
## df$NDVI       -4.461      1.412   -3.16  0.00173 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7268 on 311 degrees of freedom
## Multiple R-squared:  0.0311, Adjusted R-squared:  0.02799 
## F-statistic: 9.984 on 1 and 311 DF,  p-value: 0.001735
#Extrayendo los residuales 
df$residuales = mod1$residuals
# Mapeo de  los residuales del modelo
df %>% ggplot(aes(x = Avg_X_MCB, y=Avg_Y_MCE, colour=residuales))+
   geom_point(size = 5,shape=15)+
   scale_color_continuous(type = 'viridis')

En los modelos de regresion clasicos se requiere que que los residuales sean independientes pero en este caso como los residuales presentan dependencia espacial ajustar un modelo como estos no es conveniente

Analisis exploratorio

Separacion de la informacion georeferenciada de la informacion de las variables explicativas 

df_xy=df[,c(1,2)] #Coordenadas
X= df[,-c(1,2)] #Variables explicativas

Indice de Moran para todas las variables en X

p.values=c()
Moran_observado =c()
for(i in 1:ncol(X)){
p.values[i]=Moran.I(unlist(df[,i]),We)$p.value
Moran_observado[i]=Moran.I(unlist(df[,i]),We)$observed
}
colnames(X)
## [1] "Avg_CEa_07" "Avg_CEa_15" "NDVI"       "DEM"        "SLOPE"     
## [6] "Avg_z"      "residuales"
p.values
## [1] 0 0 0 0 0 0 0
Moran_observado
## [1] 0.45160403 0.42704518 0.26648613 0.15986601 0.09699894 0.30877314 0.06945855

Los p-valores < 0,05 indican que todas las variables explicativas parecen tener dependencia espacial, sin embargo, al evidenciar los valores del indice de Moran se puede afirmar que la mayor dependencia espacial parece estar en la CE a 75 cm. Todo lo anterior indica que se debe realizar un modelado espacial.

library(corrplot)
## corrplot 0.89 loaded
mcp=rcorr(as.matrix(X),type="pearson")
mcs=rcorr(as.matrix(X),type="spearman")
mcorp=mcp$r
mcors=mcs$r
par(mfrow=c(1,2))
corrplot(mcorp, order="hclust", tl.col="black", tl.srt=45,main="Pearson")
corrplot(mcors, order="hclust", tl.col="black", tl.srt=45,main="Spearman")

El analisis exploratorio de correlaciones permite intuir que variables se pueden usar en el modelo para explicar la dependencia espacial. En esta caso se evidencia que las variables con alta correlacion positiva son el DEM con la altura (z), Z con la CE a 75 cm mientras que hay una alta correlacion negativa de Z y DEM con la CE a 150 cm.

Comportamiento de las variables

for(i in 1:ncol(X)){
boxplot(X[,i],main=colnames(X)[i])  
}

En la CE a 150 cm, el NDVI y la pendiente se evidencian algunos datos atipicos

Prueba de normalidad para las variables

pruebas_normales=c()
for(i in 1:ncol(X)){
pruebas_normales[i]=shapiro.test(unlist(X[,i]))$p.value  
}
pruebas_normales
## [1] 9.826780e-02 9.282551e-05 1.257288e-13 1.311242e-07 5.362136e-09
## [6] 2.194832e-04 4.556127e-04

Se evidencia que la CE a 75 cm es la unica variable que cumple con los supuestos de normalidad, por tanto esta es la mas facil de modelar

Matriz de pesos para modelado espacial

contnb=dnearneigh(coordinates(df_xy),0,380000,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
class(contnb)
## [1] "nb"
#Matriz de distancias
df_xy=as.matrix(df_xy)
dlist <- nbdists(contnb, df_xy)
#Inversa de la matriz 
dlist <- lapply(dlist, function(x) 1/x)
#Matriz de pesos estandarizada
Wve=nb2listw(contnb,glist=dlist,style = "W")

Modelos de regresion espacial

Modelo autoregresivo puro (SAR)

\[Y=\lambda W Y + \alpha 1_n +\epsilon\]

SAR= spautolm(Avg_CEa_15~1, data= X, listw= Wve, family="SAR")
summary(SAR)
## 
## Call: spautolm(formula = Avg_CEa_15 ~ 1, data = X, listw = Wve, family = "SAR")
## 
## 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

Residuales del modelo

residuales_mSAR =SAR$fit$residuals
#Prueba de normalidad 
shapiro.test(residuales_mSAR)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_mSAR
## W = 0.95729, p-value = 6.37e-08

Los residuales no cumplen los supuestos de normalidad de acuerdo al p-valor

Indice de Moran para el modelo SAR

Moran.I(residuales_mSAR,We)
## $observed
## [1] 0.09349941
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004635041
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no sirve. En este caso no es suficiente explicar la CE con ella misma, se requiere colocar informacion de otras variables.

Modelo autoregresivo dos veces (SARAR) \[Y=\lambda W Y + \alpha 1_n + u\] \[u=\rho W u + \epsilon\]

SARAR= sacsarlm(Avg_CEa_15~1, data= X, listw= Wve)
## Warning in sacsarlm(Avg_CEa_15 ~ 1, data = X, listw = Wve): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16 
##   número de condición recíproco = 1.77227e-18 - using numerical Hessian.
summary(SARAR)
## 
## Call:sacsarlm(formula = Avg_CEa_15 ~ 1, data = X, 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.1713  0.9607   0.3367
## 
## Rho: 0.95895
## Approximate (numerical Hessian) standard error: 0.040792
##     z-value: 23.508, p-value: < 2.22e-16
## Lambda: 0.95895
## Approximate (numerical Hessian) standard error: 0.040755
##     z-value: 23.53, 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)

El p valor en ambos casos indica que \(\lambda W\) y \(\rho W\) son significativos, es decir, que al parecer el modelo sirve sin embargo para poder afirmar esto, se requiere calcular el indice de Moran

Residuales del modelo

residuales_SARAR =SARAR$residuals
shapiro.test(residuales_SARAR)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SARAR
## W = 0.94498, p-value = 2.076e-09

Los residuales no cumplen los supuestos de normalidad de acuerdo al p-valor

Indice de Moran para el modelo SARAR

Moran.I(residuales_SARAR,We)
## $observed
## [1] 0.0563745
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004629181
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo general espacial anidado (GNS) \[Y= \lambda WY + \alpha1_n +X\beta_{(1)}+ WX \beta_{(2)}+u \\ |\lambda| <1 \\ u=\rho W u + \epsilon \\ |\rho| <1\]

GNS= sacsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(GNS)
## 
## Call:sacsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + DEM + SLOPE + 
##     Avg_z, data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.259920 -0.257064 -0.024628  0.256195  1.934415 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     26.084548  36.251116  0.7196   0.47180
## Avg_CEa_07       0.362195   0.034530 10.4893 < 2.2e-16
## NDVI            -0.433013   1.115198 -0.3883   0.69781
## DEM              0.025591   0.016294  1.5706   0.11628
## SLOPE            0.010971   0.013341  0.8224   0.41087
## Avg_z           -0.112125   0.022705 -4.9384 7.875e-07
## lag.Avg_CEa_07  -0.486187   0.219517 -2.2148   0.02677
## lag.NDVI       -29.360420  11.602286 -2.5306   0.01139
## lag.DEM         -0.125137   0.091592 -1.3663   0.17186
## lag.SLOPE        0.883938   0.212106  4.1674 3.080e-05
## lag.Avg_z        0.205119   0.155215  1.3215   0.18633
## 
## Rho: 0.9037
## Asymptotic standard error: 0.58886
##     z-value: 1.5347, p-value: 0.12487
## Lambda: 0.94686
## Asymptotic standard error: 0.32846
##     z-value: 2.8827, p-value: 0.0039429
## 
## LR test value: 149.65, p-value: < 2.22e-16
## 
## Log likelihood: -189.5009 for sacmixed model
## ML residual variance (sigma squared): 0.19093, (sigma: 0.43695)
## Number of observations: 313 
## Number of parameters estimated: 14 
## AIC: 407, (AIC for lm: 542.65)

Residuales del modelo

residuales_GNS =GNS$residuals
shapiro.test(residuales_GNS)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_GNS
## W = 0.97364, p-value = 1.68e-05

Los residuales no cumplen el supuesto de normalidad

Indice de Moran para el modelo GNS

Moran.I(residuales_GNS,We)
## $observed
## [1] 0.04788866
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004639083
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo GNS 2 En este modelo no se tiene en cuenta el DEM debido a que como anteriormente se obtuvo un p-valor mayor a 0.05 tanto de manera natural como rezagada es necesario eliminar esta variable, la cual no esta relacionada con la CE a 150 cm

GNS2= sacsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(GNS2)
## 
## Call:sacsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + SLOPE + Avg_z, 
##     data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.242977 -0.273753 -0.027046  0.262976  1.942185 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     17.956728  33.982429  0.5284  0.597213
## Avg_CEa_07       0.373281   0.034028 10.9697 < 2.2e-16
## NDVI            -0.543373   1.116766 -0.4866  0.626571
## SLOPE            0.011352   0.013389  0.8479  0.396517
## Avg_z           -0.105276   0.022422 -4.6952 2.664e-06
## lag.Avg_CEa_07  -0.559944   0.212632 -2.6334  0.008454
## lag.NDVI       -28.742233  11.599813 -2.4778  0.013219
## lag.SLOPE        0.884282   0.203207  4.3516 1.351e-05
## lag.Avg_z        0.138273   0.128293  1.0778  0.281128
## 
## Rho: 0.90567
## Asymptotic standard error: 0.5636
##     z-value: 1.6069, p-value: 0.10807
## Lambda: 0.94735
## Asymptotic standard error: 0.31799
##     z-value: 2.9792, p-value: 0.0028904
## 
## LR test value: 147.85, p-value: < 2.22e-16
## 
## Log likelihood: -190.7997 for sacmixed model
## ML residual variance (sigma squared): 0.19248, (sigma: 0.43873)
## Number of observations: 313 
## Number of parameters estimated: 12 
## AIC: 405.6, (AIC for lm: 541.44)

En este modelo \(\rho\) no es significativo

Residuales del modelo

residuales_GNS2 =GNS2$residuals
shapiro.test(residuales_GNS2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_GNS2
## W = 0.97534, p-value = 3.293e-05

Los residuales no cumplen los supuestos de normalidad

Indice de Moran para el modelo GNS 2

Moran.I(residuales_GNS2,We)
## $observed
## [1] 0.04841028
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004640051
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo espacial del error (SEM) \[Y= \alpha 1_n + X\beta + u\] \[u=\rho W u + \epsilon\]

SEM=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(SEM)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + DEM + SLOPE + 
##     Avg_z, data = X, listw = Wve)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.4682926 -0.3245699  0.0049751  0.3215294  1.9926400 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 22.7597381  2.0758704 10.9639 < 2.2e-16
## Avg_CEa_07   0.2479677  0.0247839 10.0052 < 2.2e-16
## NDVI        -1.2527443  1.0391615 -1.2055    0.2280
## DEM         -0.0035254  0.0106002 -0.3326    0.7394
## SLOPE        0.0603502  0.0137383  4.3928 1.119e-05
## Avg_z       -0.1133123  0.0145926 -7.7650 8.216e-15
## 
## Rho: 0.96209, LR test value: 51.297, p-value: 7.9403e-13
## Asymptotic standard error: 0.02667
##     z-value: 36.074, p-value: < 2.22e-16
## Wald statistic: 1301.3, p-value: < 2.22e-16
## 
## Log likelihood: -238.6759 for lag model
## ML residual variance (sigma squared): 0.26412, (sigma: 0.51393)
## Number of observations: 313 
## Number of parameters estimated: 8 
## AIC: 493.35, (AIC for lm: 542.65)
## LM test for residual autocorrelation
## test value: 197.83, p-value: < 2.22e-16

Residuales del modelo

residuales_SEM =SEM$residuals
shapiro.test(residuales_SEM)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SEM
## W = 0.98401, p-value = 0.001482

Los residuales no cumplen el supuesto de normalidad

Indice de Moran para el modelo SEM

Moran.I(residuales_SEM,We)
## $observed
## [1] 0.06272672
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645596
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SEM 2 En este modelo no se tiene en cuenta el NDVI debido a que como anteriormente se obtuvo un p-valor mayor a 0.05, quiere decir que no esta relacionada con la variable de CE a 150 cm

SEM2=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(SEM2)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + DEM + SLOPE + Avg_z, 
##     data = X, listw = Wve)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.4816823 -0.3104710 -0.0029262  0.3123820  2.0086327 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 22.2840690  2.0415406 10.9153 < 2.2e-16
## Avg_CEa_07   0.2511915  0.0246965 10.1712 < 2.2e-16
## DEM         -0.0023501  0.0105797 -0.2221    0.8242
## SLOPE        0.0590076  0.0137250  4.2993 1.713e-05
## Avg_z       -0.1174587  0.0142130 -8.2642 2.220e-16
## 
## Rho: 0.96224, LR test value: 51.442, p-value: 7.3763e-13
## Asymptotic standard error: 0.026571
##     z-value: 36.214, p-value: < 2.22e-16
## Wald statistic: 1311.5, p-value: < 2.22e-16
## 
## Log likelihood: -239.4008 for lag model
## ML residual variance (sigma squared): 0.26534, (sigma: 0.51511)
## Number of observations: 313 
## Number of parameters estimated: 7 
## AIC: 492.8, (AIC for lm: 542.24)
## LM test for residual autocorrelation
## test value: 204.24, p-value: < 2.22e-16

Residuales del modelo

residuales_SEM2 =SEM2$residuals
shapiro.test(residuales_SEM2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SEM2
## W = 0.98402, p-value = 0.001489

Los residuales no cumplen con los supuestos de normalidad Indice de Moran para el modelo SEM 2

Moran.I(residuales_SEM2,We)
## $observed
## [1] 0.06381437
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645202
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SEM 3 En este modelo no se tiene en cuenta el NDVI y DEM debido a que como anteriormente se obtuvo un p-valor mayor a 0.05, quiere decir que no estan relacionadas con la variable de CE a 150 cm

SEM3=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+SLOPE+Avg_z, data= X, listw= Wve)
summary(SEM3)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + SLOPE + Avg_z, data = X, 
##     listw = Wve)
## 
## Residuals:
##        Min         1Q     Median         3Q        Max 
## -1.4784823 -0.3104058 -0.0036333  0.3134775  2.0108734 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) 22.244004   2.033354  10.9396 < 2.2e-16
## Avg_CEa_07   0.250818   0.024640  10.1794 < 2.2e-16
## SLOPE        0.059387   0.013617   4.3611 1.294e-05
## Avg_z       -0.119649   0.010216 -11.7118 < 2.2e-16
## 
## Rho: 0.96243, LR test value: 51.997, p-value: 5.56e-13
## Asymptotic standard error: 0.026429
##     z-value: 36.416, p-value: < 2.22e-16
## Wald statistic: 1326.1, p-value: < 2.22e-16
## 
## Log likelihood: -239.4255 for lag model
## ML residual variance (sigma squared): 0.26537, (sigma: 0.51514)
## Number of observations: 313 
## Number of parameters estimated: 6 
## AIC: 490.85, (AIC for lm: 540.85)
## LM test for residual autocorrelation
## test value: 204.57, p-value: < 2.22e-16

Residuales del modelo

residuales_SEM3 =SEM3$residuals
shapiro.test(residuales_SEM3)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SEM3
## W = 0.98381, p-value = 0.001347

Los residuales no cumplen con los supuestos de normalidad

Indice de Moran para el modelo SEM 3

Moran.I(residuales_SEM3,We)
## $observed
## [1] 0.06388428
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645126
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo espacial en rezago (SLM) \[Y=\lambda W Y + \alpha 1_n + X\beta + \epsilon\]

SML=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(SML)
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + DEM + SLOPE + 
##     Avg_z, data = X, 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.5088603  2.7607166 16.4844 < 2.2e-16
## Avg_CEa_07   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
## Avg_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)

Residuales del modelo

residuales_SML =SML$residuals
shapiro.test(residuales_SML)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SML
## W = 0.98846, p-value = 0.01377

Los residuales no cumplen los supuestos de normalidad Indice de Moran para el modelo SML

Moran.I(residuales_SML,We)
## $observed
## [1] 0.07550844
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004647281
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SML 2

En este modelo no se tiene en cuenta el NDVI debido a que como anteriormente se obtuvo un p-valor mayor a 0.05, quiere decir que no esta relacionada con la variable de CE a 150 cm

SML1=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(SML1)
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + DEM + SLOPE + 
##     Avg_z, data = X, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.477728 -0.316994 -0.014091  0.367283  2.009859 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.035550   2.729375 16.5003 < 2.2e-16
## Avg_CEa_07   0.299387   0.028492 10.5079 < 2.2e-16
## DEM         -0.008240   0.012332 -0.6682 0.5040078
## SLOPE        0.051310   0.014481  3.5432 0.0003953
## Avg_z       -0.136386   0.016857 -8.0910 6.661e-16
## 
## 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)

Residuales del modelo

residuales_SML1 =SML1$residuals
shapiro.test(residuales_SML1)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SML1
## W = 0.98823, p-value = 0.01217

Los residuales no cumplen con los supuestos de normalidad

Indice de Moran para el modelo SML 2

Moran.I(residuales_SML1,We)
## $observed
## [1] 0.07620439
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004646962
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SML 3 En este modelo no se tiene en cuenta el NDVI y DEM debido a que como anteriormente se obtuvo un p-valor mayor a 0.05, quiere decir que no estan relacionadas con la variable de CE a 150 cm

SML2=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+SLOPE+Avg_z, data= X, listw= Wve)
summary(SML2) 
## 
## Call:errorsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + SLOPE + Avg_z, 
##     data = X, 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
## Avg_CEa_07   0.297488   0.028368  10.4868 < 2.2e-16
## SLOPE        0.052248   0.014422   3.6227 0.0002915
## Avg_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)

Residuales del modelo

residuales_SML2 =SML2$residuals
shapiro.test(residuales_SML2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SML2
## W = 0.98814, p-value = 0.01166

Los residuales no cumplen los supuestos de normalidad

Indice de Moran para el modelo SML 3

Moran.I(residuales_SML2,We)
## $observed
## [1] 0.0767392
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004646816
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SDE \[Y= \alpha 1_n + X\beta_1 + WX\beta_2 + u\] \[u=\rho W u + \epsilon\]

SDE=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(SDE)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + DEM + SLOPE + 
##     Avg_z, data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.356577 -0.269279 -0.016311  0.269347  1.954889 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     28.265444  13.304190  2.1246  0.033624
## Avg_CEa_07       0.318816   0.033502  9.5163 < 2.2e-16
## NDVI            -0.463244   1.170570 -0.3957  0.692295
## DEM              0.026816   0.016575  1.6179  0.105686
## SLOPE            0.010237   0.014082  0.7269  0.467260
## Avg_z           -0.127758   0.021328 -5.9901 2.098e-09
## lag.Avg_CEa_07  -0.148588   0.197032 -0.7541  0.450771
## lag.NDVI       -33.518708  10.304158 -3.2529  0.001142
## lag.DEM         -0.135244   0.080945 -1.6708  0.094759
## lag.SLOPE        1.004251   0.157352  6.3822 1.746e-10
## lag.Avg_z        0.216661   0.099651  2.1742  0.029690
## 
## Rho: 0.91953, LR test value: 20.527, p-value: 5.8806e-06
## Asymptotic standard error: 0.056295
##     z-value: 16.334, p-value: < 2.22e-16
## Wald statistic: 266.81, p-value: < 2.22e-16
## 
## Log likelihood: -202.5621 for mixed model
## ML residual variance (sigma squared): 0.21072, (sigma: 0.45905)
## Number of observations: 313 
## Number of parameters estimated: 13 
## AIC: 431.12, (AIC for lm: 449.65)
## LM test for residual autocorrelation
## test value: 128.64, p-value: < 2.22e-16

Residuales del modelo

residuales_SDE=SDE$residuals
shapiro.test(residuales_SDE)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SDE
## W = 0.97531, p-value = 3.25e-05

Los residuales no cumplen el supuesto de normalidad

Indice de Moran para el modelo SDE

Moran.I(residuales_SDE,We)
## $observed
## [1] 0.04750286
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004640154
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Modelo SDE 2 En este modelo no se tiene en cuenta el DEM debido a que como anteriormente se obtuvo un p-valor mayor a 0.05 tanto de manera natural como rezagada es necesario eliminar esta variable, la cual no esta relacionada con la CE a 150 cm

SDE2=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(SDE2)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + NDVI + SLOPE + Avg_z, 
##     data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.333622 -0.275363 -0.029168  0.282362  1.964752 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                  Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     19.647922  12.226997  1.6069  0.108070
## Avg_CEa_07       0.330155   0.032674 10.1046 < 2.2e-16
## NDVI            -0.589440   1.173067 -0.5025  0.615332
## SLOPE            0.010556   0.014147  0.7461  0.455591
## Avg_z           -0.121066   0.020727 -5.8409 5.193e-09
## lag.Avg_CEa_07  -0.223421   0.188506 -1.1852  0.235932
## lag.NDVI       -32.761470  10.341965 -3.1678  0.001536
## lag.SLOPE        1.008515   0.155431  6.4885 8.669e-11
## lag.Avg_z        0.143016   0.082839  1.7264  0.084270
## 
## Rho: 0.92005, LR test value: 20.705, p-value: 5.3568e-06
## Asymptotic standard error: 0.055914
##     z-value: 16.455, p-value: < 2.22e-16
## Wald statistic: 270.76, p-value: < 2.22e-16
## 
## Log likelihood: -204.0637 for mixed model
## ML residual variance (sigma squared): 0.21275, (sigma: 0.46124)
## Number of observations: 313 
## Number of parameters estimated: 11 
## AIC: 430.13, (AIC for lm: 448.83)
## LM test for residual autocorrelation
## test value: 131.98, p-value: < 2.22e-16

Residuales del modelo

residuales_SDE2=SDE2$residuals
shapiro.test(residuales_SDE2)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_SDE2
## W = 0.97702, p-value = 6.545e-05

Los residuales no cumplen los supuestos de normalidad

Indice de Moran para el modelo SDE 2

Moran.I(residuales_SDE2,We)
## $observed
## [1] 0.04829465
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004641113
## 
## $p.value
## [1] 0

El p-valor indica que los residuales todavia presentan dependencia espacial, lo que quiere decir que este modelo no es util

Conclusion

En este caso ocurrio una mala especificacion del modelo, es decir, ninguno de los modelos pudo explicar la dependencia espacial presentada debido a que las variables usadas no son las adecuadas, es necesario ir a campo a tomar datos de variables que esten relacionadas con la CE como las propiedades fisicas o quimicas del suelo