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.5.1, GDAL 2.2.2, PROJ 4.9.2
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(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_dgRMatrix_listw, as_dsCMatrix_I, as_dsCMatrix_IrW,
##     as_dsTMatrix_listw, as.spam.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)

Luego de realizar la instalación de los diferentes paquetes para trabajar, realizamos la debida importaión de la base de datos, para nuestro caso utilizaremos datos de Conductividad eléctrica aparente a dos profundidades distantes, las coordenadas de cada punto muestreado, un índice de vegetación de diferencia normalizada, un modelo de elevación digital, pendiente y altura.

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

Podemos generar una matriz de correlación y graficar para tener una idea de la relación de las variables trabajadas:

MCP <- rcorr(as.matrix(data[,3:8]), type = 'p')
MCS <- rcorr(as.matrix(data[,3:8]), type = 's')
Mcorp <- MCP$r
Mcors <- MCS$r
#################3     GRÁFICOS       ####################
library(corrplot)
## corrplot 0.84 loaded
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')

Generamos una matriz de pesos para realizar los 4 modelos propuestos:

XYdata=as.matrix(cbind(data$X,data$Y))
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
dim(distancias)
## [1] 313 313

PRIMER MODELO ESPACIAL: Autoregresivo puro

\[Y = \lambda W Y + u\] ### Conductividad Eléctrica a 75 cm de profundidad

contnb=dnearneigh(coordinates(XYdata),0,380000)
dlist <- nbdists(contnb, XYdata)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W")

map=spautolm(CEA75~1,data=data,listw=Wve)  #### map=Modelo autoregresivo puro
summary(map)
## 
## 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
#ESTIMADOS DEL MODELO
CEA75_e=map$fit$fitted.values 

#####################    RESIDUALES  ######################################
res_map = map$fit$residuals 
#Normalidad
shapiro.test(res_map) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  res_map
## W = 0.99351, p-value = 0.1971
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:
# Residuales "u" del modelo

im_res_map = moran.mc(res_map, Wve, nsim = 2000) # Moran - MonteCarlo
im_res_map
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_map 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.16722, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e
CEA75=data[,3]
DF_map=data.frame(CEA75,CEA75_e)
colnames(DF_map) <- c("Cea75_o","Cea75_e")
cor(CEA75, CEA75_e)
## [1] 0.7977199
plot(CEA75,CEA75_e,main="Observados vs predichos de modelo MAP a 75 cm")

###############               GRÁFICOS                 #################################
### Mapeo de Cea a 75 cm de profundidad con estimados del primer modelo espacial

library(ggplot2)
library(akima)
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)

#Valores de CEa observados a 75 cm de profundidad
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:Hmisc':
## 
##     subplot
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea75_o*0.6, opacity = 2))
fig
## 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.
#### Valores de CEa estimados a 75 cm de profundidad


fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea75_e*0.5, opacity = 2))
fig
### Comparación estimados y observados
graE_O<-ggplot(data = data, aes(x = data$X, y = data$Y)) +
  geom_point(cex = data$CEA75*0.2) +
  geom_point(color = CEA75_e)
graE_O
## Warning: Use of `data$X` is discouraged. Use `X` instead.
## Warning: Use of `data$Y` is discouraged. Use `Y` instead.
## Warning: Use of `data$X` is discouraged. Use `X` instead.
## Warning: Use of `data$Y` is discouraged. Use `Y` instead.

#### Residuales del Modelo MAP para CEa a 75 cm de profundidad


fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map$Cea75_e - DF_map$Cea75_o)*5, opacity = 1))
fig

De acuerdo a lo anterior tenemos:

Se obtiene el siguiente modelo: \[Y = \lambda W Y + u\\ Y = 0,98811* W Y + u\\\]

GRÁFICOS

Gráficamos la CEa observada contra los valores estimados:

Valores de CEa observados a 75 cm de profundidad

library(plotly)
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea75_o*0.6, opacity = 2))
fig

Valores de CEa estimados a 75 cm de profundidad

fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea75_e*0.5, opacity = 2))
fig

Residuales del Modelo MAP para CEa a 75 cm de profundidad

fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map$Cea75_e - DF_map$Cea75_o)*5, opacity = 1))
fig

Conductividad Eléctrica a 150 cm de profundidad

map150=spautolm(CEA150~1,data=data,listw=Wve)  #### map=Modelo autoregresivo puro
summary(map150)
## 
## 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
#ESTIMADOS DEL MODELO
CEA150_e=map$fit$fitted.values 

#####################    RESIDUALES  ######################################
res_map150 = map150$fit$residuals 
#Normalidad
shapiro.test(res_map150) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  res_map150
## W = 0.95729, p-value = 6.37e-08
cvm.test(res_map150)
## 
##  Cramer-von Mises normality test
## 
## data:  res_map150
## W = 0.33539, p-value = 0.0001306
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:
# Residuales "u" del modelo

im_res_map150 = moran.mc(res_map150, Wve, nsim = 2000) # Moran - MonteCarlo
im_res_map150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res_map150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.094365, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e
CEA150=data[,3]
DF_map150=data.frame(CEA150,CEA150_e)
colnames(DF_map150) <- c("Cea150_o","Cea150_e")
cor(CEA150, CEA150_e)
## [1] 0.7977199
plot(CEA150,CEA150_e,main="Observados vs predichos de modelo MAP a 150 cm")

###############               GRÁFICOS                 #################################
### Mapeo de Cea a 150 cm de profundidad con estimados del primer modelo espacial

library(ggplot2)
library(akima)
interpol = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map150$Cea150_e,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol)
contour(interpol, add = T)

#Valores de CEa observados a 150 cm de profundidad
library(plotly)
fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea150_o*0.6, opacity = 2))
fig
#### Valores de CEa estimados a 150 cm de profundidad


fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF_map$Cea150_e*0.5, opacity = 2))
fig
### Comparación estimados y observados
graE_O<-ggplot(data = data, aes(x = data$X, y = data$Y)) +
  geom_point(cex = data$CEA150*0.2) +
  geom_point(color = CEA150_e)
graE_O
## Warning: Use of `data$X` is discouraged. Use `X` instead.
## Warning: Use of `data$Y` is discouraged. Use `Y` instead.
## Warning: Use of `data$X` is discouraged. Use `X` instead.
## Warning: Use of `data$Y` is discouraged. Use `Y` instead.

#### Residuales del Modelo MAP para CEa a 150 cm de profundidad


fig <- plot_ly(as.data.frame(XYdata), 
               x = ~data$X,
               y = ~data$Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF_map$Cea150_e - DF_map$Cea150_o)*5, opacity = 1))
fig

De acuerdo a lo anterior tenemos:

  • Con el resultado del Test de Normalidad de Shapiro-Wilk, se rechaza la hipótesis nula concluyendo así que no hay normalidad en los residuales del modelo. Y según Cramer-von Mises normality test con un p-value = 0.0001306 no hay normalidad en los residuales.
  • Lambda=0,97691
  • AIC=615,58
  • Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0,094365, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

Se obtiene el siguiente modelo: \[Y = \lambda W Y + u\\ Y = 0,97691* W Y + u\\\]

SEGUNDO MODELO ESPACIAL: SARAR / Modelo de autocorrelación espacial

Sin variables explicativas, solo Rho y Lambda.

\[Y = \lambda W Y + \alpha1_{n}+ u\\ u=\rho W u+ \epsilon\]

Conductividad Eléctrica a 75 cm de profundidad

# map = Modelo Autorregresivo- SARAR
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)
## 
## 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)
#ESTIMADOS DEL MODELO
CEA75_e2=map_sarar$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
res2=map_sarar$residuals
#Normalidad
shapiro.test(res2) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  res2
## W = 0.99417, p-value = 0.2746
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_msarar = moran.mc(res2, Wve, nsim = 2000) # Mora-MonteCarlo
moran_msarar
## 
##  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
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e2
DF_map=data.frame(CEA75,CEA75_e2)
colnames(DF_map) <- c("Cea75_o","Cea75_e2")

plot(CEA75,CEA75_e2,main="Observados vs predichos de modelo MAP a 75 cm")

De lo anterior modelo SAC se concluye: * Intercepto= -1,7267 * Error=3,2488 * Rho=0,9786 y p-value <2,22e-16 * Lambda= 0,97863 y p-value <2,22e-16 * AIC: 905.59, (AIC para lm: 1156,2) * Con el resultado del Test de Normalidad de Shapiro-Wilk, no se rechaza la hipótesis nula concluyendo así que si hay normalidad en los residuales del modelo. p-value: 0,2746. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0,10493, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

Conductividad Eléctrica a 150 cm de profundidad

# map = Modelo Autorregresivo- SARAR
map_sarar150=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(map_sarar150)
## 
## 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)
#ESTIMADOS DEL MODELO
CEA150_e2=map_sarar150$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
res2_150=map_sarar150$residuals
#Normalidad
shapiro.test(res2_150) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  res2_150
## W = 0.94498, p-value = 2.076e-09
cvm.test(res2_150)
## 
##  Cramer-von Mises normality test
## 
## data:  res2_150
## W = 0.41988, p-value = 1.637e-05
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_msarar150 = moran.mc(res2_150, Wve, nsim = 2000) # Mora-MonteCarlo
moran_msarar150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  res2_150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057087, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e2
DF_map150=data.frame(CEA150,CEA150_e2)
colnames(DF_map) <- c("Cea150_o","Cea150_e2")

plot(CEA150,CEA150_e2,main="Observados vs predichos de modelo MAP a 150 cm")

De lo anterior modelo SAC se concluye: * Intercepto= 1.1253 * Error=3,2488 * Rho=0.95895 y p-value <2,22e-16 * Lambda= 0,95895 y p-value <2,22e-16 * AIC: 570.68, (AIC para lm: 700.36) * Con el resultado del Test de Normalidad de Shapiro-Wilk, se rechaza la hipótesis nula concluyendo así que no hay normalidad en los residuales del modelo. p-value: 2.076e-09. De igual forma, al correr el test de Normalidad de Cramer-von Mises, con un p-value = 1.637e-05, se rechaza la hipótesis nula por tanto los residuales no son normales. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0.057087, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

TERCER MODELO ESPACIAL: Autoregresivo para el error

Conductividad Eléctrica a 75 cm de profundidad

QUITAR DEM Y NDVI

mee=errorsarlm(formula=CEA75~Z+NDVI+DEM+SLOPE+CEA150,data=data,listw=Wve)
summary(mee)
## 
## Call:errorsarlm(formula = CEA75 ~ Z + NDVI + DEM + SLOPE + CEA150, 
##     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
## Z             0.257034   0.028465   9.0299 < 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
## CEA150        0.859898   0.083054  10.3535 < 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)

Según lo anterior, no hay diferencias signicativas entre NVDI y DEM, por tanto probaremos el mismo modelo pero solo tomando como variables explicativas Z, SLOPE y la conductividad eléctrica aparente a 150 cm de profundidad.

mser2=errorsarlm(formula=CEA75~Z+SLOPE+CEA150,data=data,listw=Wve)
summary(mser2)
## 
## Call:errorsarlm(formula = CEA75 ~ Z + SLOPE + CEA150, 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
## Z             0.286926   0.021256  13.498 < 2.2e-16
## SLOPE        -0.079881   0.024777  -3.224  0.001264
## CEA150        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)

Al correr el modelo 2, disminuye el Criterio de información de Akaike, el cual es una medida de la calidad relativa de un modelo estadístico y al disminuir nos indica que el modelo 2 es más acertado. El estadístico AIC es una medida deldesajuste de la distribución predictiva al especificar el modelo real, y que la minimización de este estadístico (AIC) supone minimizar el desajuste y obtener el mejor modelo predictivo. Probando solo con las variables explicativas Z y CEA150, se tiene:

mser3=errorsarlm(formula=CEA75~Z+CEA150,data=data,listw=Wve)
summary(mser3)
## 
## Call:errorsarlm(formula = CEA75 ~ Z + CEA150, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.960809 -0.651968  0.005767  0.560672  2.621649 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -64.240228   5.804254 -11.068 < 2.2e-16
## Z             0.282741   0.021573  13.107 < 2.2e-16
## CEA150        0.840157   0.083883  10.016 < 2.2e-16
## 
## Lambda: 0.98359, LR test value: 106.95, p-value: < 2.22e-16
## Asymptotic standard error: 0.011573
##     z-value: 84.988, p-value: < 2.22e-16
## Wald statistic: 7223, p-value: < 2.22e-16
## 
## Log likelihood: -413.7575 for error model
## ML residual variance (sigma squared): 0.8041, (sigma: 0.89672)
## Number of observations: 313 
## Number of parameters estimated: 5 
## AIC: 837.52, (AIC for lm: 942.47)

Aumenta el AIC con respecto al primero y segundo, así que nos quedamos con el segundo modelo que tiene la pendiente, la altura y la CEA a 150 cm de profundidad como variables explicativas.

#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA75_e3=mser2$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm1=mser2$residuals
#Normalidad
shapiro.test(resm1) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm1
## W = 0.99348, p-value = 0.1948
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mser2 = moran.mc(resm1, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mser2
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.12762, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map=data.frame(CEA75,CEA75_e3)
colnames(DF_map) <- c("Cea75_o","Cea75_e3")

plot(CEA75,CEA75_e3,main="Observados vs predichos del modelo Autoregresivo para el error a 75 cm")

Así pues, se tiene lo siguiente del modelo Autoregresivo para el error a 75 cm de profundidad: * AIC: 827,77 * Lambda: 0.98246 * Con el resultado del Test de Normalidad de Shapiro-Wilk para los residuales del modelo, no se rechaza la hipótesis nula concluyendo así que si hay normalidad en los residuales del modelo. p-value: 0,07491. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0,12913, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

Conductividad Eléctrica a 150 cm de profundidad

mee150=errorsarlm(formula=CEA150~Z+NDVI+DEM+SLOPE+CEA75,data=data,listw=Wve)
summary(mee150)
## 
## Call:errorsarlm(formula = CEA150 ~ Z + NDVI + DEM + SLOPE + CEA75, 
##     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
## Z           -0.1327355  0.0171932 -7.7202 1.155e-14
## 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
## CEA75        0.2959247  0.0286368 10.3337 < 2.2e-16
## 
## 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)

Según lo anterior, hay diferencias signicativas entre Cea a 75 cm de profundidad y Z, por tanto probaremos el mismo modelo pero solo tomando como variables explicativas Z y la conductividad eléctrica aparente a 75 cm de profundidad.

mser2_150=errorsarlm(formula=CEA150~Z+CEA75,data=data,listw=Wve)
summary(mser2_150)
## 
## Call:errorsarlm(formula = CEA150 ~ Z + CEA75, data = data, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.606421 -0.342975 -0.009752  0.333663  1.921523 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) 44.479245   2.810630  15.8254 < 2.2e-16
## Z           -0.140173   0.013585 -10.3183 < 2.2e-16
## CEA75        0.288265   0.028836   9.9968 < 2.2e-16
## 
## Lambda: 0.97373, LR test value: 62.516, p-value: 2.6645e-15
## Asymptotic standard error: 0.01852
##     z-value: 52.576, p-value: < 2.22e-16
## Wald statistic: 2764.3, p-value: < 2.22e-16
## 
## Log likelihood: -246.5914 for error model
## ML residual variance (sigma squared): 0.27716, (sigma: 0.52646)
## Number of observations: 313 
## Number of parameters estimated: 5 
## AIC: 503.18, (AIC for lm: 563.7)

Sin embargo, al correr el modelo 2, aumenta el Criterio de información de Akaike, el cual es una medida de la calidad relativa de un modelo estadístico y al aumentar nos indica que el modelo 1 es más acertado. El estadístico AIC es una medida deldesajuste de la distribución predictiva al especificar el modelo real, y que la minimización de este estadístico (AIC) supone minimizar el desajuste y obtener el mejor modelo predictivo.

mser2_150=errorsarlm(formula=CEA150~Z+SLOPE+CEA75,data=data,listw=Wve)
summary(mser2_150)
## 
## Call:errorsarlm(formula = CEA150 ~ Z + SLOPE + CEA75, 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
## Z           -0.143289   0.013322 -10.7558 < 2.2e-16
## SLOPE        0.052248   0.014422   3.6227 0.0002915
## CEA75        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)

Al correr el 3 modelo y ver que disminuye el AIC con respecto al primero, nos quedamos con el último modelo que tiene como variables explicativas la altura, la pendiente y la CEA a 75 cm de profundidad.

#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA150_e3=mser2_150$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm1_150=mser2_150$residuals
#Normalidad
shapiro.test(resm1_150) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm1_150
## W = 0.98814, p-value = 0.01166
cvm.test(resm1_150)
## 
##  Cramer-von Mises normality test
## 
## data:  resm1_150
## W = 0.065621, p-value = 0.3182
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mser2_150 = moran.mc(resm1_150, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mser2_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm1_150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.076281, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map150=data.frame(CEA150,CEA150_e3)
colnames(DF_map150) <- c("Cea150_o","Cea150_e3")

plot(CEA150,CEA150_e3,main="Observados vs predichos del modelo Autoregresivo para el error a 150 cm")

De esta forma, se tiene lo siguiente del modelo Autoregresivo para el error a 150 cm de profunfdidad: * AIC: 492,35 * Lambda: 0.96928, p-value: 1.1945e-12 * Con el resultado del Test de Normalidad de Shapiro-Wilk para los residuales del modelo, se rechaza la hipótesis nula concluyendo así que no hay normalidad en los residuales del modelo. p-value: 0.01166. Sin embargo, al correr el test de Normalidad de Cramer-von Mises, con un p-value = 0.3182, no se rechaza la hipótesis nula por tanto los residuales son normales. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0.076281, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

CUARTO MODELO ESPACIAL: Lambda y Rho

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

Conductividad Eléctrica a 75 cm de profundidad

mlrho1=sacsarlm(formula=CEA75~Z+NDVI+DEM+SLOPE+CEA150,data=data,listw=Wve)
summary(mlrho1)
## 
## Call:sacsarlm(formula = CEA75 ~ Z + NDVI + DEM + SLOPE + CEA150, data = data, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.095552 -0.452042 -0.016313  0.479376  2.211327 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -58.705272  16.043512 -3.6591 0.0002531
## Z             0.192558   0.031300  6.1520 7.651e-10
## NDVI         -2.117899   1.662599 -1.2738 0.2027170
## DEM           0.026261   0.018843  1.3936 0.1634245
## SLOPE        -0.060744   0.021561 -2.8173 0.0048425
## CEA150        0.844893   0.072944 11.5828 < 2.2e-16
## 
## Rho: 0.97457
## Asymptotic standard error: 0.38121
##     z-value: 2.5565, p-value: 0.010573
## Lambda: 0.97219
## Asymptotic standard error: 0.41724
##     z-value: 2.33, p-value: 0.019805
## 
## LR test value: 180.04, p-value: < 2.22e-16
## 
## Log likelihood: -365.758 for sac model
## ML residual variance (sigma squared): 0.58131, (sigma: 0.76243)
## Number of observations: 313 
## Number of parameters estimated: 9 
## AIC: 749.52, (AIC for lm: 925.56)

\[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \] Este modelo tiene un AIC mucho menor de 749.52 correspondiente al modelo anterior.

#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA75_e4=mlrho1$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm4=mlrho1$residuals
#Normalidad
shapiro.test(resm4) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm4
## W = 0.99436, p-value = 0.3016
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mlrho1 = moran.mc(resm4, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mlrho1
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm4 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.092291, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map=data.frame(CEA75,CEA75_e4)
colnames(DF_map) <- c("Cea75_o","Cea75_e4")

plot(CEA75,CEA75_e4,main="Observados vs predichos del modelo Autoregresivo para las variables y el error a 75cm")

Modelo de menor AIC - Eliminando a NDVI y DEM como variable explicativa

#Modelo de mejor AIC - Borrando a NDVI como variable explicativa

mlrho2=sacsarlm(formula=CEA75~Z+SLOPE+CEA150,data=data,listw=Wve)
summary(mlrho2)
## 
## Call:sacsarlm(formula = CEA75 ~ Z + SLOPE + CEA150, 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.726424 -4.0361 5.435e-05
## Z             0.213388   0.028874  7.3904 1.463e-13
## SLOPE        -0.065991   0.021538 -3.0639  0.002185
## CEA150        0.857117   0.073143 11.7183 < 2.2e-16
## 
## 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)
#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA75_e4.1=mlrho2$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm4.1=mlrho2$residuals
#Normalidad
shapiro.test(resm4.1) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm4.1
## W = 0.99548, p-value = 0.4995
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mlrho2 = moran.mc(resm4.1, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mlrho2
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm4.1 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.091618, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map=data.frame(CEA75,CEA75_e4.1)
colnames(DF_map) <- c("Cea75_o","Cea75_e4.1")

plot(CEA75,CEA75_e4.1,main="Observados vs predichos del modelo Autoregresivo para las variables y el error a 75cm")

Así pues, se tiene lo siguiente del modelo Autoregresivo para las variables explicativas y el error: * AIC: 749,14 * Lambda: 0.97219 * Rho: 0.97457 * Con el resultado del Test de Normalidad de Shapiro-Wilk para los residuales del modelo, no se rechaza la hipótesis nula concluyendo así que si hay normalidad en los residuales del modelo. p-value: 0,3016. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0,092291, con un p-valor=0,0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

Conductividad Eléctrica a 150 cm de profundidad

mlrho1_150=sacsarlm(formula=CEA150~Z+NDVI+DEM+SLOPE+CEA75,data=data,listw=Wve)
summary(mlrho1_150)
## 
## Call:sacsarlm(formula = CEA150 ~ Z + NDVI + DEM + SLOPE + CEA75, data = data, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.435255 -0.287126 -0.017044  0.306288  1.947100 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 21.886264  16.212844  1.3499    0.1770
## Z           -0.111676   0.018912 -5.9051 3.524e-09
## NDVI        -1.006362   1.049559 -0.9588    0.3376
## DEM         -0.001063   0.012093 -0.0879    0.9300
## SLOPE        0.042252   0.013536  3.1215    0.0018
## CEA75        0.287666   0.033326  8.6318 < 2.2e-16
## 
## Rho: 0.94878
## Asymptotic standard error: 0.57536
##     z-value: 1.649, p-value: 0.099142
## Lambda: 0.95692
## Asymptotic standard error: 0.4852
##     z-value: 1.9722, p-value: 0.048585
## 
## LR test value: 86.653, p-value: < 2.22e-16
## 
## Log likelihood: -220.9976 for sac model
## ML residual variance (sigma squared): 0.23222, (sigma: 0.48189)
## Number of observations: 313 
## Number of parameters estimated: 9 
## AIC: 460, (AIC for lm: 542.65)

\[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \] Este modelo tiene un AIC mucho menor de 495.

#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA150_e4=mlrho1_150$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm4_150=mlrho1_150$residuals
#Normalidad
shapiro.test(resm4_150) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm4_150
## W = 0.98154, p-value = 0.0004664
cvm.test(resm4_150)
## 
##  Cramer-von Mises normality test
## 
## data:  resm4_150
## W = 0.11491, p-value = 0.06995
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mlrho1_150 = moran.mc(resm4_150, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mlrho1_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm4_150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057283, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map_150=data.frame(CEA150,CEA150_e4)
colnames(DF_map_150) <- c("Cea150_o","Cea150_e4")

plot(CEA150,CEA150_e4,main="Observados vs predichos del modelo Autoregresivo para las variables y el error a 150cm")

Modelo de mejor AIC - Borrando a NDVI y DEM como variable explicativa

#Modelo de mejor AIC - Borrando a NDVI como variable explicativa

mlrho2_150=sacsarlm(formula=CEA150~Z+SLOPE+CEA75,data=data,listw=Wve)
summary(mlrho2_150)
## 
## Call:sacsarlm(formula = CEA150 ~ Z + SLOPE + CEA75, 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
## Z           -0.114901   0.017473 -6.5759 4.835e-11
## SLOPE        0.041795   0.013504  3.0949  0.001969
## CEA75        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)
#ESTIMADOS DEL MODELO AUTOREGRESIVO PARA EL ERROR
CEA150_e4.1=mlrho2_150$fitted.values

#####################    RESIDUALES  ######################################
# Residuales "u" del modelo
resm4.1_150=mlrho2_150$residuals
#Normalidad
shapiro.test(resm4.1_150) #/No se debería rechazar
## 
##  Shapiro-Wilk normality test
## 
## data:  resm4.1_150
## W = 0.98084, p-value = 0.0003408
cvm.test(resm4.1_150)
## 
##  Cramer-von Mises normality test
## 
## data:  resm4.1_150
## W = 0.13217, p-value = 0.04088
#Ahora, realizamos una simulación MonteCarlo para ver si queda u o epsilon, es decir, residuos con o sin dependencia espacial:

moran_mlrho2_150 = moran.mc(resm4.1_150, Wve, nsim = 2000) # Mora-MonteCarlo
moran_mlrho2_150
## 
##  Monte-Carlo simulation of Moran I
## 
## data:  resm4.1_150 
## weights: Wve  
## number of simulations + 1: 2001 
## 
## statistic = 0.057956, observed rank = 2001, p-value = 0.0004998
## alternative hypothesis: greater
###############     ANÁLISIS DE PREDICCIÓN DEL MODELO        ###########################
# Data CEa_o y CEa_e3
DF_map_150=data.frame(CEA150,CEA150_e4.1)
colnames(DF_map_150) <- c("Cea150_o","Cea150_e4.1")

plot(CEA150,CEA150_e4.1,main="Observados vs predichos del modelo Autoregresivo para las variables y el error a 150cm")

Así pues, se tiene lo siguiente del modelo Autoregresivo para las variables explicativas y el error a 150 cm de profundidad: * AIC: 456,91 * Lambda: 0.95737 * Rho: 0.94901 pero con un p-valor:0.096902 * Con el resultado del Test de Normalidad de Shapiro-Wilk para los residuales del modelo, se rechaza la hipótesis nula concluyendo así que no hay normalidad en los residuales del modelo. p-value: 0.0003407. * Con la simulación Monte-Carlo del Índice de Moran, obtenemos un I de Morán de 0.05795, con un p-valor=0.0004998 lo que nos indica que hay dependencia espacial positiva de los residuales.

ANÁLISIS DE MODELOS ESPACIALES

Para saber cual es el mejor modelo espacial tabulamos todos los AIC y se obtiene lo siguiente:

Tabla Resumen

Modelo Profundidad CEA Coeficiente de Información Kaike Lambda Rho
Autoregresivo puro 75 cm 995,65 0,98811
150 cm 615,58 0,97691
SARAR 75 cm 905,59 0,97863 0,97863
150 cm 570,68 0,95895 0,0,95895
MEE 75 cm 827,77 0,98246
150 cm 492.35 0.96928
MLRHO 75 cm 749,52 0,97219 0,97457
150 cm 460 0,95692 0,94878*
MLRHO(-NDVI) 75 cm 749,14 0,9722 0,97458
150 cm 458,91 0,95737 0,94901*

También podemos hacer uso de la librería metrics para comparar los valores observados y estimados RMSE (más chico mejor).

library(Metrics)
rmse(CEA75,CEA75_e)
## [1] 1.16061
rmse(CEA75,CEA75_e2)
## [1] 0.9926651
rmse(CEA75,CEA75_e3)
## [1] 0.8824011
rmse(CEA75,CEA75_e4)
## [1] 0.7624345
rmse(CEA75,CEA75_e4.1)
## [1] 0.7673774
rmse(CEA150,CEA150_e)
## [1] 1.16061
rmse(CEA150,CEA150_e2)
## [1] 8.892897
rmse(CEA150,CEA150_e3)
## [1] 8.871811
rmse(CEA150,CEA150_e4)
## [1] 8.874602
rmse(CEA150,CEA150_e4.1)
## [1] 8.8745

LUEGO DE AJUSTADO EL MEJOR MODELO Vamos a mapear las variables que evidencien relación estadística con CEA75:

Mapeo de Cea a 75 cm de profundidad

library(ggplot2)
library(akima)
interpol = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = DF_map$Cea75_e4.1,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol)
contour(interpol, add = T)

puntos_muestra = sample(1:nrow(XYdata), nrow(XYdata)*0.8, F)
muestra = XYdata[puntos_muestra,]
plot(XYdata, col = (1:nrow(XYdata)%in%puntos_muestra)+2, pch = 16)

cc = chull(XYdata)
cc = c(cc, cc[1])
plot(XYdata, pch = 16)
lines(XYdata[cc,], type = 'l')

Predicción de un valor de CEa en el espacio

Si se quiere encontrar un valor estimado conociendo las coordenadas del punto se tiene:

\[Y = \lambda W Y + \epsilon\] \[Y - \lambda W Y = \epsilon\\ ( I - \lambda W) Y = \epsilon\\ Y = ( I - \lambda W)^{-1}\epsilon\] \[E(Y) = E((I-\lambda W)^{-1}\epsilon)\\ E(Y) = E((I-\lambda W)^{-1})E(\epsilon)\\ \require{cancel} E(Y) = E((I-\lambda W)^{-1})\cancelto{0}{E(\epsilon)}\\\]

Modelo Espacial Autorregresivo Puro final

\[\hat{Cea_{75}} = \alpha + \lambda W Cea_{75}\]

Realizamos el gráfico con el punto seleccionado:

mapa = chull(XYdata)
mapeo = c(cc, cc[1])
plot(XYdata, pch = 16)
lines(XYdata[cc,], type = 'l')
### Grafico del punto seleccionado
points(x=843560, y=956126, col="red", pch=19)

Para la digresión es necesario crear una nueva matriz de pesos incluyendo el valor de las nuevas coordenadas de manera que procedemos a realizar la matriz de pesos y a reemplazar con la nueva matriz y el despeje del modelo autoregresivo puro inicial, obteniendo el valor predicho o estimado de Cea para una profundidad de 75 cm, para las coordenadas x=843560, y=956126.

point<-data.frame(843560,956126,0,0,0,0,0,0)
names(point)<-c("X","Y","CEA75","CEA150","NDVI","DEM","SLOPE","Z")
w2<-rbind(data,point)
matriz<-as.matrix(dist(cbind(w2$X, w2$Y)))
inversa<- 1/matriz
inversa[is.infinite(inversa)] <- 0
matriz2<-diag(314)
# Estandarizando la  nueva matriz
matriz3<-apply(inversa,1,sum)
# Matriz de Pesos estandarizada "we"
Wfin<-inversa/matriz3

A<-0.98811*Wfin
B<-matriz2-A
C<-solve(B) # Matriz inversa
D<-5.6941*C #Matriz con los resultados finales
D[314,314] # Valor estimado de CEA a 75cm en la coordenada del punto seleccionado
## [1] 7.296308
######TRASNFORMACIÓN
#dat_subset_sp <- BD_PARA_REGRESION
#coordinates(dat_subset_sp) <- ~ X + Y
#str(dat_subset_sp)
#dat_subset_sp
#proy <- CRS("+proj=longlat +datum=WGS84")
#dat_subset_sp@proj4string <- proy

FUENTES BIBLIOGRÁFICAS