Asignación:

Realice las regresiones espaciales vistas en clase para la variable CEA a 150 cm y publiquelo en Rpubs

Introducir datos
library(readxl)
df <- read_excel("~/university/Computacion/Clase (22 April)/BD_MODELADO.xlsx")
View(df)

Matriz de distancias

library(ape)
df.dists <- as.matrix(dist(cbind(df$Avg_X_MCB, df$Avg_Y_MCE))) # me calcula las distancias

Inverso de la distancias

df.dists.inv <- 1/df.dists
diag(df.dists.inv) <- 0  # Asignar ceros a la diagonal donde sale inf
df.dists.inv <- round(df.dists.inv,3) # redondeo 
df.dists.inv[1:5, 1:5] # visualizar
##       1     2     3     4     5
## 1 0.000 0.193 0.022 0.054 0.046
## 2 0.193 0.000 0.025 0.047 0.057
## 3 0.022 0.025 0.000 0.017 0.030
## 4 0.054 0.047 0.017 0.000 0.034
## 5 0.046 0.057 0.030 0.034 0.000
We<-df.dists.inv/rowSums(df.dists.inv) # Matriz de pesos estandarizada
df_xy=df[,c(1,2)] # Coords
X= df[,-c(1,2)]  # Explicativas

Modelo lineal no espacial

plot(df$Avg_CEa_15,df$NDVI,pch=16,cex=0.6)

plot(df$Avg_CEa_15,df$DEM,pch=16,cex=0.6)

plot(df$Avg_CEa_15,df$SLOPE,pch=16,cex=0.6)

De manera exploratoria no se aprecia relacion entre la CEa a 150 con el NDVI ni con altura ni pendiente

mod1.1=lm(df$Avg_CEa_15~df$NDVI)
summary(mod1.1)
## 
## 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
hist(mod1.1$residuals)

mod1.2=lm(df$Avg_CEa_15~df$DEM)
summary(mod1.2)
## 
## Call:
## lm(formula = df$Avg_CEa_15 ~ df$DEM)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.58433 -0.44587 -0.05758  0.35063  2.84724 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 31.270887   1.740187  17.970  < 2e-16 ***
## df$DEM      -0.062246   0.008482  -7.339 1.89e-12 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.6817 on 311 degrees of freedom
## Multiple R-squared:  0.1476, Adjusted R-squared:  0.1449 
## F-statistic: 53.86 on 1 and 311 DF,  p-value: 1.894e-12
hist(mod1.2$residuals)

mod1.3=lm(df$Avg_CEa_15~df$SLOPE)
summary(mod1.3)
## 
## Call:
## lm(formula = df$Avg_CEa_15 ~ df$SLOPE)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.58422 -0.46796 -0.03566  0.42734  2.91769 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 18.24028    0.08838 206.375  < 2e-16 ***
## df$SLOPE     0.06375    0.01897   3.361 0.000874 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.7253 on 311 degrees of freedom
## Multiple R-squared:  0.03504,    Adjusted R-squared:  0.03194 
## F-statistic: 11.29 on 1 and 311 DF,  p-value: 0.0008744
hist(mod1.3$residuals)

Se aprecia que para estos modelos no espaciales, el coeficiente de determinación es muy bajo, posiblemente explicado por la dependencia espacial que puede tener el CEa

Indice de Moran para CEa 150
Moran.I(df$Avg_CEa_15,df.dists.inv)
## $observed
## [1] 0.159866
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004644691
## 
## $p.value
## [1] 0

La CEa a 150, parece tener dependencia espacial, por lo cual es necesario utilizar modelos que me expliquen dicha dependencia

Debo estandarizar los valores de la matriz de pesos para que lamda este entre 0 y 1

Matriz de pesos para modelo espacial

#preparando librerias
library(spdep)
library(ape)
library(sp)
library(MVA)
library(Hmisc)
library(normtest)
library(nortest)
library(dplyr)
library(spatialreg)
library(corrplot)
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"
df_xy=as.matrix(df_xy)
dlist <- nbdists(contnb, df_xy)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W") # matriz de peso estandarizada

Modelo autoregresivo puro SAR

map= spautolm(Avg_CEa_15~1, data= X, listw= Wve, family="SAR")
summary(map)
## 
## 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
Normalidad de map
residuales_map =map$fit$residuals
shapiro.test(residuales_map) 
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_map
## W = 0.95729, p-value = 6.37e-08
Dependencia
Moran.I(residuales_map,We)
## $observed
## [1] 0.09349941
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004635041
## 
## $p.value
## [1] 0

Este modelo, sin otras variables no explica la dependencia espacial

Normalidad datos CEa
shapiro.test(df$Avg_CEa_15)
## 
##  Shapiro-Wilk normality test
## 
## data:  df$Avg_CEa_15
## W = 0.97785, p-value = 9.283e-05

Modelo SARAR

msarar= 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 
##   reciprocal condition number = 1.77227e-18 - using numerical Hessian.
summary(msarar)
## 
## 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 modelo SARAR es un mejor modelo que el SAR pues su ACI es menor

Residuos SARAR
residuos_msarar <- msarar$residuals
Pruebas SARAR
shapiro.test(residuos_msarar)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_msarar
## W = 0.94498, p-value = 2.076e-09
Moran.I(residuos_msarar, We)
## $observed
## [1] 0.0563745
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004629181
## 
## $p.value
## [1] 0

El modelo SARAR no explica la dependenca espacial de los datos

Modelo SLM

mslm=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(mslm)
## 
## 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)

El modelo SLM que tienen variables, es mejor que los modelos de autoregresión

Residuos SLM
residuos_mslm <- mslm$residuals
Pruebas SLM
shapiro.test(residuos_mslm)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_mslm
## W = 0.98846, p-value = 0.01377
Moran.I(residuos_mslm, We)
## $observed
## [1] 0.07550844
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004647281
## 
## $p.value
## [1] 0

El modelo SLM con todas las variables no explica la dependencia espacial

mslm.b=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(mslm.b)
## 
## 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_mslm.b <- mslm.b$residuals
shapiro.test(residuales_mslm.b)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_mslm.b
## W = 0.98823, p-value = 0.01217
Moran.I(residuales_mslm.b, We)
## $observed
## [1] 0.07620439
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004646962
## 
## $p.value
## [1] 0

Sin NDVI el modelo SLM mejoro un poco pero no explica la dependencia espacial

mslm.c=errorsarlm(formula=Avg_CEa_15~Avg_CEa_07+SLOPE+Avg_z, data= X, listw= Wve)
summary(mslm.c)
## 
## 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)
residuos_mslm.c <- mslm.c$residuals
shapiro.test(residuos_mslm.c)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_mslm.c
## W = 0.98814, p-value = 0.01166
Moran.I(residuos_mslm.c, We)
## $observed
## [1] 0.0767392
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004646816
## 
## $p.value
## [1] 0

A pesar de remover el DEM el modelo no explica la dependencia espacial

Modelo SEM

msem=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(msem)
## 
## 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
residuos_msem <- msem$residuals
shapiro.test(residuos_msem)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_msem
## W = 0.98401, p-value = 0.001482
Moran.I(residuos_msem, We)
## $observed
## [1] 0.06272672
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645596
## 
## $p.value
## [1] 0
msem.b=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve)
summary(msem.b)
## 
## 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_msem.b <- msem.b$residuals
shapiro.test(residuales_msem.b)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_msem.b
## W = 0.98402, p-value = 0.001489
Moran.I(residuales_msem.b,We)
## $observed
## [1] 0.06381437
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645202
## 
## $p.value
## [1] 0
msem.c=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+SLOPE+Avg_z, data= X, listw= Wve)
summary(msem.c)
## 
## 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
residuos_msem.c <-msem.c$residuals
shapiro.test(residuos_msem.c)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuos_msem.c
## W = 0.98381, p-value = 0.001347
Moran.I(residuos_msem.c, We)
## $observed
## [1] 0.06388428
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004645126
## 
## $p.value
## [1] 0

El modelo SEM no explica la dependencia espacial, en cambio el modelo SLM sin NDVI y sin DEM tiene menor AIC

Modelo SDE

msde=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(msde)
## 
## 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_msde <- msde$residuals
shapiro.test(residuales_msde)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_msde
## W = 0.97531, p-value = 3.25e-05
Moran.I(residuales_msde, We)
## $observed
## [1] 0.04750286
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004640154
## 
## $p.value
## [1] 0
msde.b=lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(msde.b)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + DEM + SLOPE + Avg_z, 
##     data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.379820 -0.276005 -0.029616  0.282895  2.145603 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)    15.893764  13.013420  1.2213   0.22196
## Avg_CEa_07      0.338295   0.034000  9.9498 < 2.2e-16
## DEM             0.026634   0.016976  1.5689   0.11666
## SLOPE           0.014220   0.014423  0.9859   0.32416
## Avg_z          -0.118428   0.021396 -5.5350 3.112e-08
## lag.Avg_CEa_07 -0.331487   0.197441 -1.6789   0.09317
## lag.DEM        -0.133819   0.083016 -1.6120   0.10697
## lag.SLOPE       0.627814   0.131129  4.7878 1.686e-06
## lag.Avg_z       0.141067   0.098780  1.4281   0.15327
## 
## Rho: 0.93587, LR test value: 26.617, p-value: 2.4806e-07
## Asymptotic standard error: 0.044993
##     z-value: 20.801, p-value: < 2.22e-16
## Wald statistic: 432.66, p-value: < 2.22e-16
## 
## Log likelihood: -211.0934 for mixed model
## ML residual variance (sigma squared): 0.2222, (sigma: 0.47138)
## Number of observations: 313 
## Number of parameters estimated: 11 
## AIC: 444.19, (AIC for lm: 468.8)
## LM test for residual autocorrelation
## test value: 204.42, p-value: < 2.22e-16
residuales_msde.b <- msde.b$residuals
shapiro.test(residuales_msde.b)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_msde.b
## W = 0.96557, p-value = 8.953e-07
Moran.I(residuales_msde.b, We)
## $observed
## [1] 0.06140014
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004635677
## 
## $p.value
## [1] 0

El AIC del modelo SED disminuye si se remueve el NDVI

msde.c <- lagsarlm(formula=Avg_CEa_15~Avg_CEa_07+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(msde.c)
## 
## Call:lagsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + SLOPE + Avg_z, data = X, 
##     listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.356185 -0.291875 -0.027334  0.266255  2.154409 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)     7.568151  11.915822  0.6351   0.52534
## Avg_CEa_07      0.349842   0.033080 10.5757 < 2.2e-16
## SLOPE           0.014569   0.014485  1.0058   0.31450
## Avg_z          -0.112156   0.020824 -5.3860 7.205e-08
## lag.Avg_CEa_07 -0.405714   0.188414 -2.1533   0.03129
## lag.SLOPE       0.634660   0.128416  4.9422 7.723e-07
## lag.Avg_z       0.070200   0.081047  0.8662   0.38640
## 
## Rho: 0.93603, LR test value: 26.717, p-value: 2.3559e-07
## Asymptotic standard error: 0.044869
##     z-value: 20.861, p-value: < 2.22e-16
## Wald statistic: 435.19, p-value: < 2.22e-16
## 
## Log likelihood: -212.4981 for mixed model
## ML residual variance (sigma squared): 0.2242, (sigma: 0.4735)
## Number of observations: 313 
## Number of parameters estimated: 9 
## AIC: 443, (AIC for lm: 467.71)
## LM test for residual autocorrelation
## test value: 206.76, p-value: < 2.22e-16
residuales_msde.c <- msde.c$residuals
shapiro.test(residuales_msde.c)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_msde.c
## W = 0.96734, p-value = 1.649e-06
Moran.I(residuales_msde.c, We)
## $observed
## [1] 0.06194672
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.00463645
## 
## $p.value
## [1] 0

El SDE no mejora al remover el DEM

Modelo GNS

mgns= sacsarlm(formula=Avg_CEa_15~Avg_CEa_07+NDVI+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(mgns)
## 
## 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_mgns <- mgns$residuals
shapiro.test(residuales_mgns)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_mgns
## W = 0.97364, p-value = 1.68e-05
Moran.I(residuales_mgns, We)
## $observed
## [1] 0.04788866
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004639083
## 
## $p.value
## [1] 0

Este modelo GNS es el que tiene un AIC menor, el mejor de los modelos ensayados

mgns.b= sacsarlm(formula=Avg_CEa_15~Avg_CEa_07+DEM+SLOPE+Avg_z, data= X, listw= Wve,type="mixed")
summary(mgns.b)
## 
## Call:sacsarlm(formula = Avg_CEa_15 ~ Avg_CEa_07 + DEM + SLOPE + Avg_z, 
##     data = X, listw = Wve, type = "mixed")
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.270871 -0.268895 -0.043227  0.251359  2.068286 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                 Estimate Std. Error z value  Pr(>|z|)
## (Intercept)    14.256309  35.671498  0.3997  0.689410
## Avg_CEa_07      0.376570   0.034707 10.8500 < 2.2e-16
## DEM             0.025822   0.016516  1.5635  0.117944
## SLOPE           0.012529   0.013519  0.9268  0.354028
## Avg_z          -0.107267   0.022501 -4.7673 1.868e-06
## lag.Avg_CEa_07 -0.611325   0.220653 -2.7705  0.005597
## lag.DEM        -0.123820   0.093118 -1.3297  0.183615
## lag.SLOPE       0.604869   0.189089  3.1989  0.001380
## lag.Avg_z       0.144362   0.158025  0.9135  0.360957
## 
## Rho: 0.91887
## Asymptotic standard error: 0.59159
##     z-value: 1.5532, p-value: 0.12037
## Lambda: 0.95352
## Asymptotic standard error: 0.34156
##     z-value: 2.7917, p-value: 0.0052437
## 
## LR test value: 140.96, p-value: < 2.22e-16
## 
## Log likelihood: -194.644 for sacmixed model
## ML residual variance (sigma squared): 0.19691, (sigma: 0.44375)
## Number of observations: 313 
## Number of parameters estimated: 12 
## AIC: 413.29, (AIC for lm: 542.24)
residuales_mgns.b <- mgns.b$residuals
shapiro.test(residuales_mgns.b)
## 
##  Shapiro-Wilk normality test
## 
## data:  residuales_mgns.b
## W = 0.96429, p-value = 5.831e-07
Moran.I(residuales_mgns.b, We)
## $observed
## [1] 0.05521322
## 
## $expected
## [1] -0.003205128
## 
## $sd
## [1] 0.004634968
## 
## $p.value
## [1] 0

El modelo no mejora cuando el NDVI es removido

Ninguno de los modelos utilizados explica la dependencia espacial de CEa a 150 posiblemente debido a que las variables medidas en campo no son las adecuadas con relacion a la CE