Se desea hallar el mejor modelo para la estimación de la CEa a dos profundidades: 75cm y 150cm en una zona de estudio, así pues se prueban cuatro modelos de regresión espacial los cuales tienen como formula general la descrita acontinuación.

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

La Conductividad Electrica aparente (CEa) tiene potencial para ser un parámetro de calidad del suelo, ya que puede identificar algunas propiedades geoeléctricas del suelo, además de zonas con riesgos de salinidad (Corwin et al., 2006).

Preparando datos

library(readxl)
datat <- read_excel("C:/Users/Sofia Hernandez/Documents/2020-2/Computacion estadistica/BD_MODELADO.xlsx")
head(datat)
## # A tibble: 6 x 8
##         X       Y CEa75 CEa150  NDVI   DEM SLOPE     Z
##     <dbl>   <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 843450. 955962.  7.24   18.0 0.863  199  6.39   193.
## 2 843455. 955962.  6.79   18.0 0.867  197. 1.98   193.
## 3 843494. 955951.  6.85   18.7 0.875  197  0.578  194.
## 4 843440. 955978.  7.14   18.3 0.846  197  1.18   194.
## 5 843469. 955973.  6.83   17.9 0.797  197  0.211  194.
## 6 843501. 955975.  6.70   18.4 0.758  198. 4.36   195.
library(spdep)
library(ape)
library(sp)
library(MVA)
library(Hmisc)
library(normtest)
library(nortest)
library(spatialreg)
library(dplyr)
res=datat[,3:4] # Variables respuesta
ex=datat[,5:8] # Variables explicativa
X=as.matrix(data.frame(ex))
options(digits = 9)

# Pasos matriz de peso
XYdata=as.matrix(datat[,1:2]) #Coordenadas
CE.d=as.matrix(dist(XYdata, diag=T, upper=T))
CE.d[1:5,1:5]
##             1           2          3          4          5
## 1  0.00000000  5.17585431 45.2932806 18.5048488 21.9357763
## 2  5.17585431  0.00000000 40.3796252 21.5015061 17.4254715
## 3 45.29328062 40.37962516  0.0000000 60.0542601 32.8990988
## 4 18.50484877 21.50150614 60.0542601  0.0000000 29.4799218
## 5 21.93577630 17.42547155 32.8990988 29.4799218  0.0000000
min(CE.d[CE.d!=0])
## [1] 5.17585431
max(CE.d)
## [1] 853.009894
dim(CE.d)
## [1] 313 313
CE.d.inv=as.matrix(1/CE.d)
diag(CE.d.inv)<-0
W=as.matrix(CE.d.inv) # Matriz de peso sin estandarización 
SUMAS=apply(W,1,sum)
We<-W/SUMAS # Matriz de peso Estandarizada
apply(We,1,FUN=sum)
##   1   2   3   4   5   6   7   8   9  10  11  12  13  14  15  16  17  18  19  20 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  21  22  23  24  25  26  27  28  29  30  31  32  33  34  35  36  37  38  39  40 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  41  42  43  44  45  46  47  48  49  50  51  52  53  54  55  56  57  58  59  60 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  61  62  63  64  65  66  67  68  69  70  71  72  73  74  75  76  77  78  79  80 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
##  81  82  83  84  85  86  87  88  89  90  91  92  93  94  95  96  97  98  99 100 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 
##   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1   1 
## 301 302 303 304 305 306 307 308 309 310 311 312 313 
##   1   1   1   1   1   1   1   1   1   1   1   1   1

Indice de Moran

Moran_CEa75=(Moran.I(datat$CEa75, We)); Moran_CEa75
## $observed
## [1] 0.268746829
## 
## $expected
## [1] -0.00320512821
## 
## $sd
## [1] 0.00466590594
## 
## $p.value
## [1] 0
Moran_CEa150=(Moran.I(datat$CEa150, We)); Moran_CEa150
## $observed
## [1] 0.160951018
## 
## $expected
## [1] -0.00320512821
## 
## $sd
## [1] 0.00465455026
## 
## $p.value
## [1] 0
library(dplyr)
bind_rows(Moran_CEa75, Moran_CEa150, id=NULL)
## # A tibble: 2 x 4
##   observed expected      sd p.value
##      <dbl>    <dbl>   <dbl>   <dbl>
## 1    0.269 -0.00321 0.00467       0
## 2    0.161 -0.00321 0.00465       0

Se observa dependencia espacial en las dos variables respuesta: CEa75 y CEa150.

Matrices de correlación

Se muestran todas las variables y su correlación.

# Exploracion de datos 
library(corrplot)
## corrplot 0.84 loaded
mcp=rcorr(as.matrix(datat[,3:8]), type="pearson")
mcs=rcorr(as.matrix(datat[,3:8]), 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")

Se evidencia correlación positiva en CEa75 para DEM y Z, pero negativa para SLOPE, lo que quiere decir que a mayor pendiente la CEa es menor. Por otro lado, CEa150 presenta correlación positiva con SLOPE, indicando que al aumentar la pendiente aumenta la CEa, asimismo para DEM y Z la correlación es negativa (cuando una aumenta el otro disminuye).

# Describiendo datos de otra forma
library(psych)
## 
## Attaching package: 'psych'
## The following object is masked from 'package:Hmisc':
## 
##     describe
## The following objects are masked from 'package:ggplot2':
## 
##     %+%, alpha
library(ggplot2)
options(digits = 5)
describe(datat[,3:7])
##        vars   n   mean   sd median trimmed  mad    min    max range  skew
## CEa75     1 313   9.77 1.53   9.65    9.76 1.40   6.17  13.60  7.44  0.10
## CEa150    2 313  18.50 0.74  18.44   18.47 0.64  16.80  21.51  4.71  0.57
## NDVI      3 313   0.83 0.03   0.84    0.84 0.02   0.70   0.88  0.17 -1.39
## DEM       4 313 205.11 4.55 204.83  205.13 5.68 196.17 213.33 17.17  0.06
## SLOPE     5 313   4.13 2.16   3.78    3.94 2.10   0.21  12.72 12.51  0.95
##        kurtosis   se
## CEa75     -0.28 0.09
## CEa150     1.22 0.04
## NDVI       2.48 0.00
## DEM       -1.13 0.26
## SLOPE      1.25 0.12

Se observa que a mayor profundidad se tienen valores mayores de CEa, lo que se puede deber a lo encontrado por Paggi et al. (2013) y Bosch Mayol et al. (2012) los cuales determinaron una estrecha asociación entre la CEa y la profundidad del suelo y el sodio soluble, respectivamente. Así tambien, se presenta una mayor dispersion en las variables DEM, SLOPE y CEa75, que en el NDVI.

Matriz de peso para los modelos

contnb=dnearneigh(coordinates(XYdata),0,380000,longlat = F) #380 mil grados
contnb
## Neighbour list object:
## Number of regions: 313 
## Number of nonzero links: 97656 
## Percentage nonzero weights: 99.681 
## Average number of links: 312
dlist <- nbdists(contnb, XYdata)
dlist <- lapply(dlist, function(x) 1/x) #Matriz de peso, lista para calcular el inverso de las distancias, lo mas distante me afecta menos
Wve=nb2listw(contnb,glist=dlist,style = "W") #Matriz de peso estilo W
Wve$weigths # pesos, todas las distancias valen
## NULL
dim(Wve) # No dá los datos, dá la info en forma de lista
## NULL

Modelos Espaciales: CEa 75cm

Se presenta el mapa de distribución de CEa75, donde a mayor tamaño de los puntos mayor valores de conductividad, como se observa en la parte media superior.

knitr::include_graphics("C:/Users/Sofia Hernandez/Documents/2020-2/Computacion estadistica/GEODA/ce75_puntos.png")

MODELO 1: modelo autoregresivo puro

En el cual no se tienen variables explicativas

\[Y = \lambda W Y + \epsilon\]

Dependencia espacial

map_CEa75=spautolm(CEa75~1,data=datat,listw=Wve)  #### MODELO
summary(map_CEa75)
## 
## Call: spautolm(formula = CEa75 ~ 1, data = datat, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -3.258254 -0.650679 -0.071829  0.824652  3.063002 
## 
## Coefficients: 
##             Estimate Std. Error z value Pr(>|z|)
## (Intercept)   5.6941     5.5177   1.032   0.3021
## 
## Lambda: 0.98811 LR test value: 162.5 p-value: < 2.22e-16 
## Numerical Hessian standard error of lambda: 0.011866 
## 
## Log likelihood: -494.82 
## ML residual variance (sigma squared): 1.347, (sigma: 1.1606)
## Number of observations: 313 
## Number of parameters estimated: 3 
## AIC: 995.65

Asume que \(\rho\) es igual a cero, ya que \(u\)=\(\epsilon\), lo que significa que es autoregresivo solo una vez, la conductividad electrica electrica aparente en poligono se puede explicar con la conductividad electrica electrica aparente de los “vecinos”.

Valores estimados del modelo 1

Para comparar entre los valores observados y determinar si el modelo es bueno.

CEa75E=map_CEa75$fit$fitted.values;CEa75E # Estimados
##       1       2       3       4       5       6       7       8       9      10 
##  8.6308  8.6776  8.9153  8.9056  8.8236  8.8560  8.9267  9.0756  9.0054  8.9769 
##      11      12      13      14      15      16      17      18      19      20 
##  8.8846  8.8911  9.0060  9.1184  9.0633  9.0784  9.0171  8.9640  9.0108  9.1408 
##      21      22      23      24      25      26      27      28      29      30 
##  9.1678  9.1535  9.1619  9.1763  9.0825  9.0567  9.0478  9.0672  9.0447  9.1880 
##      31      32      33      34      35      36      37      38      39      40 
##  9.1637  9.2192  9.2347  9.2775  9.2092  9.1145  9.1417  9.0571  9.0632  9.0851 
##      41      42      43      44      45      46      47      48      49      50 
##  9.3782  9.4087  9.3506  9.3033  9.1994  9.2065  9.2210  9.2981  9.3359  9.3437 
##      51      52      53      54      55      56      57      58      59      60 
##  9.1671  9.2176  9.1321  9.1396  9.2313  9.3898  9.3622  9.3639  9.3071  9.2695 
##      61      62      63      64      65      66      67      68      69      70 
##  9.2370  9.2178  9.2240  9.2513  9.3058  9.3910  9.4173  9.4125  9.3594  9.2655 
##      71      72      73      74      75      76      77      78      79      80 
##  9.3239  9.2164  9.2656  9.3774  9.3644  9.3646  9.3213  9.2619  9.2327  9.2311 
##      81      82      83      84      85      86      87      88      89      90 
##  9.2489  9.2862  9.3333  9.4118  9.4733  9.5075  9.4876  9.4922  9.4381  9.3610 
##      91      92      93      94      95      96      97      98      99     100 
##  9.4106  9.5162  9.3367  9.3264  9.2866  9.2545  9.2576  9.2849  9.3203  9.3817 
##     101     102     103     104     105     106     107     108     109     110 
##  9.4590  9.5114  9.5719  9.6277  9.5854  9.6143  9.6095  9.5499  9.5690  9.5446 
##     111     112     113     114     115     116     117     118     119     120 
##  9.6038  9.6738  9.3280  9.3182  9.2953  9.2850  9.3124  9.3656  9.4151  9.4837 
##     121     122     123     124     125     126     127     128     129     130 
##  9.5574  9.6224  9.6810  9.7327  9.7106  9.7839  9.7116  9.6820  9.6372  9.7157 
##     131     132     133     134     135     136     137     138     139     140 
##  9.8259  9.3429  9.3451  9.3261  9.3542  9.3945  9.4620  9.5138  9.5859  9.6685 
##     141     142     143     144     145     146     147     148     149     150 
##  9.7419  9.8130  9.8626  9.8637  9.9185  9.8194  9.8963  9.7944  9.9096 10.0180 
##     151     152     153     154     155     156     157     158     159     160 
##  9.3705  9.3963  9.3978  9.4401  9.4755  9.5495  9.5940  9.6980  9.7920  9.8731 
##     161     162     163     164     165     166     167     168     169     170 
##  9.9467 10.0081 10.0306 10.0494  9.9619 10.0184 10.0643 10.0020 10.1681 10.0806 
##     171     172     173     174     175     176     177     178     179     180 
##  9.4504  9.4740  9.5055  9.5540  9.6177  9.7048  9.8234  9.9248 10.0229 10.0915 
##     181     182     183     184     185     186     187     188     189     190 
## 10.1567 10.2052 10.1909 10.1217 10.2187 10.2377 10.2423 10.2441 10.2290  9.4994 
##     191     192     193     194     195     196     197     198     199     200 
##  9.5343  9.5627  9.6535  9.7414  9.8607 10.0037 10.0919 10.1894 10.2583 10.3270 
##     201     202     203     204     205     206     207     208     209     210 
## 10.3895 10.3665 10.4452 10.3516 10.3532 10.4099 10.1819 10.1901 10.1525  9.5697 
##     211     212     213     214     215     216     217     218     219     220 
##  9.6149  9.6843  9.7804  9.8942 10.0048 10.1651 10.2510 10.3339 10.4222 10.4946 
##     221     222     223     224     225     226     227     228     229     230 
## 10.5156 10.5303 10.5659 10.4488 10.3972 10.3374 10.0945 10.1034  9.6453  9.7259 
##     231     232     233     234     235     236     237     238     239     240 
##  9.8020  9.9227 10.0388 10.1558 10.3017 10.3613 10.4819 10.5475 10.5711 10.6251 
##     241     242     243     244     245     246     247     248     249     250 
## 10.6139 10.5648 10.4490 10.3797 10.2696 10.0658  9.7964  9.8558  9.9555 10.0398 
##     251     252     253     254     255     256     257     258     259     260 
## 10.1606 10.3221 10.3718 10.4630 10.5717 10.6140 10.6161 10.6422 10.6120 10.4919 
##     261     262     263     264     265     266     267     268     269     270 
## 10.3797 10.2952  9.9715 10.0433 10.1069 10.2565 10.3711 10.4114 10.5254 10.5896 
##     271     272     273     274     275     276     277     278     279     280 
## 10.6187 10.5942 10.6018 10.5574 10.4023 10.0749 10.1916 10.2864 10.3859 10.4369 
##     281     282     283     284     285     286     287     288     289     290 
## 10.5489 10.5654 10.5766 10.5609 10.5157 10.1271 10.1347 10.2299 10.3079 10.3587 
##     291     292     293     294     295     296     297     298     299     300 
## 10.4448 10.5030 10.4969 10.5031 10.4499 10.1627 10.1992 10.2353 10.3047 10.3435 
##     301     302     303     304     305     306     307     308     309     310 
## 10.4126 10.3965 10.4263 10.2179 10.1799 10.2282 10.2952 10.3118 10.1904 10.1663 
##     311     312     313 
## 10.1790 10.2662 10.1322
CEa75=datat$CEa75
DF1=data.frame(CEa75,CEa75E) # Comparación observado y estimado
head(DF1)
##    CEa75 CEa75E
## 1 7.2375 8.6308
## 2 6.7873 8.6776
## 3 6.8483 8.9153
## 4 7.1352 8.9056
## 5 6.8268 8.8236
## 6 6.7000 8.8560
colnames(DF1) <- c("CEa75","CEa75E")

Correlación entre lo observado y lo esperado

Se realiza la correlación y la comparación por medio de la libreria Metrics, dando

library(plotly)
cor<-plot_ly(x = ~CEa75, y = ~CEa75E, marker = list(size = 10, color = 'cyan',line = list(color = 'cyan3',width = 2)))## Modelo 1
cor <- cor %>% layout(title = "Modelo 1: Autoregresivo puro para CEa75", 
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE)) 
cor 
cor1<-cor(DF1$CEa75,DF1$CEa75E);cor1 # Correlación entre lo estimado y lo observado
## [1] 0.79772
library(Metrics)
comp1=rmse(DF1$CEa75,DF1$CEa75E);comp1 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 1.1606
rCEa75=map_CEa75$fit$residuals # Residuales modelo 1: diferencia entre lo observado y lo estimado
n_map=shapiro.test(rCEa75) # Normalidad de los residuales
n_map.=ifelse(n_map$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(rCEa75)), col="blue", pch=19)

Se realiza el Indice de Moran Montecarlo para los residuales, esto con el fin de encontrar la dependencia espacial, ya que cuando los puntos de los residuales son mas grandes (como se observa en el borde inferior), quiere decir que los lugares donde peor ajusto el modelo.

library(spdep)
McrCEa75=moran.mc(rCEa75, Wve, nsim=2000) 
McrCEa75_mc=ifelse(McrCEa75$p.value<0.05, "Dependencia espacial", "Independencia espacial")# Normalidad de los residuales

Se evidencia que los residuales tienen dependencia espacial.

library(plotly)
fig1.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF1$CEa75, opacity = 0.5))
fig1.1<-layout(fig1.1, title= "Modelo1: CEa75 observada") #Observado

fig1.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF1$CEa75E, opacity = 0.5))
fig1.2<-layout(fig1.2, title= "Modelo1: CEa75 estimada") 
# Estimado 
fig1.3<-subplot(fig1.1, fig1.2)%>% hide_legend()
fig1.3<-layout(fig1.3, title= "Modelo1: CEa75 observada y estimada");fig1.3
fig1.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF1$CEa75E-DF1$CEa75)*10, opacity = 0.5))
fig1.4<-layout(fig1.4, title= "Modelo1: Diferencia entre la CEa75 observada y estimada");fig1.4

MAP: Modelo Autorregresivo - S ARAR

El modelo es doblemente auto regresivo tanto para la respuesta como para los residuales, sin variables explicativas aún.

\[Y = \lambda W Y + u\\ u=\rho W u+ \epsilon\]

msarar_CEa75=sacsarlm(CEa75~1,data=datat,listw=Wve)  ## Modelo
summary(msarar_CEa75)
## 
## Call:sacsarlm(formula = CEa75 ~ 1, data = datat, 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.2577   -0.53   0.5961
## 
## Rho: 0.97863
## Approximate (numerical Hessian) standard error: 0.021305
##     z-value: 45.935, p-value: < 2.22e-16
## Lambda: 0.97863
## Approximate (numerical Hessian) standard error: 0.021299
##     z-value: 45.946, p-value: < 2.22e-16
## 
## LR test value: 254.56, p-value: < 2.22e-16
## 
## Log likelihood: -448.79 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)

El modelo ajusta tanto a lambda como a rho con un 2.1%.

CEa75E2=msarar_CEa75$fitted.values # Estimado
rCEa75.2=msarar_CEa75$residuals # Normalidad de residuales
n_sarar=shapiro.test(rCEa75.2)
n_sarar.=ifelse(n_sarar$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado

plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(rCEa75.2)), col="blue", pch=19) 

mc.sarar.CEa75 = moran.mc(rCEa75.2, Wve, nsim = 2000)
mc.sarar.CEa75mc=ifelse(mc.sarar.CEa75$p.value<0.05, "Dependencia espacial", "Independencia espacial")# Dependencia, se debe rompeer en dos partes en rho

Se realiza el Indice de Moran Montecarlo para los residuales: los puntos mas grandes observados en el centro y borde inferior indica un mal ajuste del modelo.

Correlación y comparación

library(Metrics)
comp2=rmse(DF1$CEa75,CEa75E2) # Comparar valores observados y estimados, mas pequeño mejor ajuste

cr2 <- plot_ly(x = ~CEa75, y = ~CEa75E2,marker = list(size = 10,
color = 'rgba(255, 182, 193, .9)',line =list(color ='rgba(152, 0, 0, .8)',
width = 2)))

cr2 <- cr2 %>% layout(title = "Modelo 2: ARAR para CEa75", yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) ## Modelo 4
cr2 ## Modelo 2
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
cor2<-cor(CEa75,CEa75E2);cor2
## [1] 0.81919
fig2.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF1$CEa75, opacity = 0.5))
fig2.1<-layout(fig2.1, title= "Modelo2: CEa75 observada") #Observado

fig2.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa75E2, opacity = 0.5))
fig2.2<-layout(fig2.2, title= "Modelo2: CEa75 estimada") 
# Estimado 
fig2.3<-subplot(fig2.1, fig2.2)%>% hide_legend()
fig2.3<-layout(fig2.3, title= "Modelo2: CEa75 observada y estimada");fig2.3
fig2.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa75E2-DF1$CEa75)*10, opacity = 0.5))
fig2.4<-layout(fig2.4, title= "Modelo1: Diferencia entre la CEa75 observada y estimada");fig2.4

Modelo espacial del error

En este modelo ya se establecen las variables explicativas, por lo que se deben escoger las que se ajustan mejor. Para ello, se corre el modelo con todas estas y se extraen las que tengan mayor significancia como se evidencia acontinuación.

# Modelos que permitan involucrar variables explicativas (X)
mser1=errorsarlm(formula=CEa75~CEa150+DEM+Z+SLOPE+NDVI,data=datat,listw=Wve)
summary(mser1)
## 
## Call:errorsarlm(formula = CEa75 ~ CEa150 + DEM + Z + SLOPE + NDVI, 
##     data = datat, listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.019160 -0.540466 -0.045367  0.513314  2.592838 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) -64.737579   5.752902 -11.2530 < 2.2e-16
## CEa150        0.859898   0.083054  10.3535 < 2.2e-16
## DEM           0.036792   0.020974   1.7542  0.079402
## Z             0.257034   0.028465   9.0299 < 2.2e-16
## SLOPE        -0.073067   0.024760  -2.9510  0.003168
## NDVI         -2.395368   1.907913  -1.2555  0.209301
## 
## 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.1 for error model
## ML residual variance (sigma squared): 0.76603, (sigma: 0.87523)
## Number of observations: 313 
## Number of parameters estimated: 8 
## AIC: 828.2, (AIC for lm: 925.56)
mser2=errorsarlm(formula=CEa75~CEa150+Z+SLOPE,data=datat,listw=Wve)
summary(mser2)
## 
## Call:errorsarlm(formula = CEa75 ~ CEa150 + Z + SLOPE, data = datat, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -2.150527 -0.558459 -0.045187  0.540349  2.578564 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -65.325177   5.620712 -11.622 < 2.2e-16
## CEa150        0.874324   0.083217  10.507 < 2.2e-16
## Z             0.286926   0.021256  13.498 < 2.2e-16
## SLOPE        -0.079881   0.024777  -3.224  0.001264
## 
## 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.65 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)

En este modelo, las variables explicativas que mejor se ajustaron fueron CEa150, Z y SLOPE, en donde solo se tiene a \(\lambda\), esto tambien se pudo evidenciar con la matriz de correlación realizada en la exploración de datos. Por otro lado, se hallo dependencia espacial por medio de Moran Montecarlo el cual indico el rechazo de la hipotesis nula para independencia.

### Modelo mser 2, presenta un menor AIC
rmser2=mser2$residuals
n_mser2=shapiro.test(rmser2)
n_mser2.=ifelse(n_mser2$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(rmser2)), col="blue", pch=19)

mc.mser2.CEa75 = moran.mc(rmser2, Wve, nsim = 2000)
mc.mser2.CEa75mc =ifelse(mc.mser2.CEa75$p.value<0.05, "Dependencia espacial", "Independencia espacial")
# Dependencia o independencia de los residuales

CEa75E3= mser2$fitted.values #Estimado
# correlacion
cr3<-plot_ly(x = ~CEa75, y = ~CEa75E3,marker = list(size = 10, color = 'cyan',line = list(color = 'cyan3',width = 2)))## Modelo 4

cr3 <- cr3 %>% layout(title = "Modelo 3: Modelo espacial del error para CEa75", yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) 
cr3## Modelo 3
comp3=rmse(CEa75,CEa75E3);comp3 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 0.8824
cor3<-cor(CEa75,CEa75E3);cor3
## [1] 0.82275
fig3.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF1$CEa75, opacity = 0.5))
fig3.1<-layout(fig3.1, title= "Modelo3: CEa75 observada") #Observado

fig3.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa75E3, opacity = 0.5))
fig3.2<-layout(fig3.2, title= "Modelo3: CEa75 estimada") 
# Estimado 
fig3.3<-subplot(fig3.1, fig3.2)%>% hide_legend()
fig3.3<-layout(fig3.3, title= "Modelo3: CEa75 observada y estimada");fig3.3
fig3.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa75E3-DF1$CEa75)*10, opacity = 0.5))
fig3.4<-layout(fig3.4, title= "Modelo3: Diferencia entre la CEa75 observada y estimada");fig3.4

Modelo 4: SACSAR

Apartir del modelo anterior se escogen las variables que mejor se ajustan, que son CEa150, Z (altura) y SLOPE (pendiente).

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

mlrho.CEa75=sacsarlm(formula=CEa75~CEa150+Z+SLOPE,data=datat,listw=Wve)
summary(mlrho.CEa75)
## 
## Call:sacsarlm(formula = CEa75 ~ CEa150 + Z + SLOPE, data = datat, 
##     listw = Wve)
## 
## Residuals:
##         Min          1Q      Median          3Q         Max 
## -2.11409261 -0.48405766 -0.00093026  0.51350443  2.20507300 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) -59.437373  14.726432 -4.0361 5.435e-05
## CEa150        0.857117   0.073143 11.7183 < 2.2e-16
## Z             0.213388   0.028874  7.3904 1.463e-13
## SLOPE        -0.065991   0.021538 -3.0639  0.002185
## 
## 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.020433
## 
## LR test value: 179.22, p-value: < 2.22e-16
## 
## Log likelihood: -367.79 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)

Se evidencia en el modelo que tanto \(\lambda\) como \(\rho\) son significativos, por lo que es necesario que ambos esten presentes.

### Modelo mser 1 menor, presenta un menor AIC
r.mlrho.CEa75=mlrho.CEa75$residuals
n_sacsar=shapiro.test(r.mlrho.CEa75) # Normalidad de los residuales
n_sacsar.=ifelse(n_sacsar$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(r.mlrho.CEa75)), col="blue", pch=19)

moran_mee = moran.mc(r.mlrho.CEa75, Wve, nsim = 2000)
moran_mee_mc=ifelse(moran_mee$p.value<0.05, "Dependencia espacial", "Independencia espacial") # Dependencia o independencia de los residuales. ¿se queda u o epsilon?
        
CEa75E4= mlrho.CEa75$fitted.values # Estimado

cr4 <- plot_ly(x = ~CEa75, y = ~CEa75E4,
               marker = list(size = 10,
                             color = 'rgba(255, 182, 193, .9)', 
                             line = list(color ='rgba(152, 0, 0, .8)',width = 2)))
cr4 <- cr4 %>% layout(title = 'Modelo 4: Modelo con lambda y rho para CEa75',
         yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) ## Modelo 4

cr4
comp4=rmse(DF1$CEa75,CEa75E4);comp4 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 0.76738
cor4<-cor(DF1$CEa75,CEa75E4);cor4
## [1] 0.86955

Se muestra dependencia espacial por medio del Moran Montecarlo ya que presenta un p.valor significativo, en la grafica los residuales son representados por puntos, donde al ser mas grandes indican un peor ajuste del modelo, es decir en el centro.

fig4.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF1$CEa75, opacity = 0.5))
fig4.1<-layout(fig4.1, title= "Modelo4: CEa75 observada") #Observado

fig4.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa75E4, opacity = 0.5))
fig4.2<-layout(fig4.2, title= "Modelo4: CEa75 estimada") 
# Estimado 
fig4.3<-subplot(fig4.1, fig4.2)%>% hide_legend()
fig4.3<-layout(fig4.3, title= "Modelo4: CEa75 observada y estimada");fig4.3
fig4.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa75E4-DF1$CEa75)*10, opacity = 0.5))
fig4.4<-layout(fig4.4, title= "Modelo4: Diferencia entre la CEa75 observada y estimada");fig4.4

### Comparación de modelos

AIC

A partir de los cuatro modelos planteados, se debe escoger el que se ajusta mejor, para esto se considera el AIC, el cual es una medida de la calidad relativa de un modelo estadıstico (Borrego, 2018).

Modelo=rep(c("MAP", "SARAR", "SEM", "SAC"))
AIC=c(995.65,905.59,829.3,749.59)
Correlación=c(cor1,cor2,cor3, cor4)
RMSE=c(comp1,comp2, comp3, comp4)
Normalidad_res= c(n_map.,n_sarar.,n_mser2.,n_sacsar.)
Moran_mc=c(McrCEa75_mc, mc.sarar.CEa75mc, mc.mser2.CEa75mc, moran_mee_mc )
Tab_CEA75=data.frame(Modelo,(round(AIC,3)),(round(Correlación,3)),(round(RMSE,3)),Normalidad_res, Moran_mc)
library(DT)
Tb_CEA75<-datatable(Tab_CEA75, class = 'cell-border stripe', editable = 'cell', rownames = FALSE, colnames = c("MODELO", "AIC", "CORRELACION", "RMSE", "NORMALIDAD_RES", "MORAN_MC"));Tb_CEA75

Se confirma atraves de la anterior tabla que el mejor modelo para predecir la Conductividad Electrica aparente a una profundidad de 75 cm es el Modelo 4, ya que presenta el AIC mas bajo, así como la comparación de la Raíz del error cuadrático medio (RMSE) que muestra la diferencias entre los valores predichos, en este modelo se aprecio el valor mas pequeño, indicando un mejor ajuste, finalmaente se obtuvo la correlacion mas alta, siendo de 0.87. En los cuatro modelos se evidencia normalidad y dependencia residual.

Con la diferencia entre lo observado y lo esperado de todos los modelos, se puede observar una mayor impresición de los residuales en el borde inferior y en el centro del area estudiada, sin embargo, varia en cada uno de los modelos.

La elección de variables explicativas tales como Z, SLOPE y CEa150 fueron las que mejor ajustaron en el modelo, en SLOPE se obtuvo que a mayor pendiente se obtienen menores valores de CEa, esto se puede deber a que la escorrentia del agua causa el lavado Segun investigaciones previas, reportaron que la CEa es influenciada por el contenido de partículas finas y Ar , lo que refleja la capacidad de retención de agua del suelo en zonas mas bajas de acumulación y por lo tanto, la variación espacial del rendimiento de grano de los cultivos (Herber, 2011). Asi tambien, las zonas bajas de CEa están asociadas a las zonas con lomas arenosas de baja capacidad de retención de hídrica (Urricariet et al., 2011) y MO (Simon, etal., 2013).

En el estudio de Simon et al., se obtuvo que la zona con una clase alta de CEa, estuvo asociada a zonas bajas del terreno, con mayor fertilidad potencial y rendimiento.

Se analizan los valores de las curvas de nivel del área muestreada mediante interpolación Akima

library(ggplot2)
library(akima)
## Warning: package 'akima' was built under R version 4.0.3
interpol = interp(x = XYdata[,1],
                  y = XYdata[,2],
                  z = CEa75,
                  nx = 500, ny = 500,
                  linear = F)
image(interpol)
contour(interpol, add = T)

Estimación de un valor en el modelo

General Spatial Model by Anselin (1988) or as an SAC model by LeSage and Kelly (2009).

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

\[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \\ E(Y-\lambda W Y) = E(\alpha 1_n+ X\beta + u)\\ E((I-\lambda W)Y) = \alpha + X\beta + u\\ (I-\lambda W) E(Y)) = \alpha + X\beta + u\\ (I-\lambda W) \hat{Y}) = \alpha + X\beta + u\\ \hat{Y} = (\alpha + X\beta + u)*((I-\lambda W)^-1)\] \[\hat{CEa75E} = (\alpha + X\beta + u)*((I-\lambda W)^-1)\] \((I-\lambda W)\) CONSTANTE

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

new_coor<-data.frame(843350, 956100, 0,0,0,0,0,0)
names(new_coor)<-c("X","Y","CEa75", "CEa150","NDVI","DEM", "SLOPE","Z")
new<-rbind(datat,new_coor)
XYdata1=as.matrix(new[,1:2])

### Matriz de peso
CE.d.est=as.matrix(dist(XYdata1, diag=T, upper=T))
CE.d.est[1:5,1:5]
##         1       2      3      4      5
## 1  0.0000  5.1759 45.293 18.505 21.936
## 2  5.1759  0.0000 40.380 21.502 17.425
## 3 45.2933 40.3796  0.000 60.054 32.899
## 4 18.5048 21.5015 60.054  0.000 29.480
## 5 21.9358 17.4255 32.899 29.480  0.000
min(CE.d.est[CE.d.est!=0])
## [1] 5.1759
max(CE.d.est)
## [1] 853.01
dim(CE.d.est)
## [1] 314 314
CE.d.est.inv=as.matrix(1/CE.d.est)
diag(CE.d.est.inv)<-0
W.e=as.matrix(CE.d.est.inv) # Matriz de peso sin estandarización 
#Estandarización
SUMAS=apply(W.e,1,sum)
We.e<-W.e/SUMAS # Matriz de peso Estandarizada

# Matriz Identidad
diag(CE.d.est.inv)<-314
Ident=as.matrix(CE.d.est.inv)
## Estimación 
Est1<-0.97212*We.e
Est2<-Ident-Est1
Est3<-solve(Est2)
Est4<-(-59.437373)*Est3
dim(Est4)
## [1] 314 314
Est4[314,314]
## [1] -0.18929

Modelos Espaciales: CEa 150cm

knitr::include_graphics("C:/Users/Sofia Hernandez/Documents/2020-2/Computacion estadistica/GEODA/ce150puntos.png")

MODELO 1: modelo autoregresivo puro (MAP)

En el cual no se tienen variables explicativas

\[Y = \lambda W Y + \epsilon\]

Dependencia espacial

map_CEa150=spautolm(CEa150~1,data=datat,listw=Wve)  #### MODELO
summary(map_CEa150)
## 
## Call: spautolm(formula = CEa150 ~ 1, data = datat, 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.79 
## ML residual variance (sigma squared): 0.40168, (sigma: 0.63378)
## Number of observations: 313 
## Number of parameters estimated: 3 
## AIC: 615.58

Asume que rho es igual a cero, ya que \(u\)=\(\epsilon\), lo que significa que es autoregresivo solo una vez, la conductividad electrica electrica aparente en poligono se puede explicar con la conductividad electrica electrica aparente de los “vecinos”.

Valores estimados del modelo 1

Para comparar entre los valores observados y determinar si el modelo es bueno.

CEa150E=map_CEa150$fit$fitted.values;CEa150E # Estimados
##      1      2      3      4      5      6      7      8      9     10     11 
## 18.528 18.523 18.573 18.592 18.587 18.584 18.606 18.662 18.650 18.640 18.613 
##     12     13     14     15     16     17     18     19     20     21     22 
## 18.610 18.621 18.692 18.681 18.694 18.662 18.643 18.650 18.712 18.727 18.702 
##     23     24     25     26     27     28     29     30     31     32     33 
## 18.721 18.720 18.692 18.666 18.636 18.595 18.584 18.732 18.719 18.727 18.732 
##     34     35     36     37     38     39     40     41     42     43     44 
## 18.725 18.712 18.658 18.650 18.597 18.579 18.538 18.928 18.948 18.805 18.761 
##     45     46     47     48     49     50     51     52     53     54     55 
## 18.704 18.715 18.709 18.734 18.723 18.692 18.626 18.607 18.570 18.535 18.511 
##     56     57     58     59     60     61     62     63     64     65     66 
## 18.860 18.908 18.846 18.787 18.727 18.705 18.705 18.699 18.700 18.688 18.709 
##     67     68     69     70     71     72     73     74     75     76     77 
## 18.684 18.645 18.619 18.599 18.579 18.532 18.496 18.478 18.843 18.820 18.785 
##     78     79     80     81     82     83     84     85     86     87     88 
## 18.738 18.705 18.681 18.693 18.687 18.682 18.661 18.664 18.646 18.607 18.585 
##     89     90     91     92     93     94     95     96     97     98     99 
## 18.546 18.478 18.451 18.469 18.795 18.785 18.737 18.720 18.705 18.683 18.672 
##    100    101    102    103    104    105    106    107    108    109    110 
## 18.653 18.639 18.621 18.620 18.621 18.589 18.595 18.550 18.500 18.460 18.421 
##    111    112    113    114    115    116    117    118    119    120    121 
## 18.430 18.435 18.777 18.752 18.718 18.710 18.688 18.662 18.651 18.606 18.598 
##    122    123    124    125    126    127    128    129    130    131    132 
## 18.588 18.579 18.586 18.564 18.539 18.492 18.453 18.411 18.423 18.448 18.746 
##    133    134    135    136    137    138    139    140    141    142    143 
## 18.737 18.698 18.681 18.650 18.616 18.581 18.558 18.552 18.544 18.527 18.538 
##    144    145    146    147    148    149    150    151    152    153    154 
## 18.505 18.454 18.424 18.407 18.418 18.451 18.488 18.717 18.683 18.664 18.630 
##    155    156    157    158    159    160    161    162    163    164    165 
## 18.588 18.550 18.510 18.506 18.504 18.484 18.452 18.468 18.421 18.382 18.389 
##    166    167    168    169    170    171    172    173    174    175    176 
## 18.397 18.423 18.470 18.513 18.494 18.627 18.607 18.564 18.517 18.473 18.439 
##    177    178    179    180    181    182    183    184    185    186    187 
## 18.443 18.428 18.403 18.376 18.388 18.354 18.351 18.383 18.401 18.434 18.489 
##    188    189    190    191    192    193    194    195    196    197    198 
## 18.555 18.593 18.577 18.527 18.482 18.443 18.393 18.376 18.371 18.348 18.330 
##    199    200    201    202    203    204    205    206    207    208    209 
## 18.321 18.313 18.324 18.346 18.372 18.412 18.469 18.495 18.560 18.667 18.698 
##    210    211    212    213    214    215    216    217    218    219    220 
## 18.499 18.447 18.393 18.353 18.324 18.305 18.299 18.267 18.283 18.280 18.278 
##    221    222    223    224    225    226    227    228    229    230    231 
## 18.335 18.378 18.392 18.446 18.492 18.547 18.599 18.698 18.410 18.366 18.315 
##    232    233    234    235    236    237    238    239    240    241    242 
## 18.277 18.265 18.250 18.246 18.228 18.264 18.272 18.306 18.361 18.407 18.418 
##    243    244    245    246    247    248    249    250    251    252    253 
## 18.515 18.545 18.596 18.602 18.336 18.306 18.255 18.226 18.231 18.223 18.206 
##    254    255    256    257    258    259    260    261    262    263    264 
## 18.230 18.270 18.300 18.324 18.398 18.435 18.460 18.552 18.547 18.267 18.229 
##    265    266    267    268    269    270    271    272    273    274    275 
## 18.188 18.204 18.197 18.201 18.255 18.284 18.329 18.372 18.432 18.466 18.481 
##    276    277    278    279    280    281    282    283    284    285    286 
## 18.242 18.181 18.189 18.194 18.227 18.279 18.316 18.368 18.418 18.461 18.215 
##    287    288    289    290    291    292    293    294    295    296    297 
## 18.204 18.173 18.186 18.221 18.264 18.320 18.360 18.422 18.431 18.204 18.210 
##    298    299    300    301    302    303    304    305    306    307    308 
## 18.184 18.216 18.275 18.323 18.361 18.406 18.247 18.246 18.274 18.348 18.361 
##    309    310    311    312    313 
## 18.276 18.316 18.324 18.374 18.339
CEa150=datat$CEa150
DF2=data.frame(CEa150,CEa150E) # Comparación observado y estimado
head(DF2)
##   CEa150 CEa150E
## 1 18.027  18.528
## 2 18.027  18.523
## 3 18.704  18.573
## 4 18.342  18.592
## 5 17.924  18.587
## 6 18.394  18.584
colnames(DF2) <- c("CEa150","CEa150E")

Correlación entre lo observado y lo esperado

Se realiza la correlación y la comparación por medio de la libreria Metrics, dando

library(plotly)
cor<-plot_ly(x = ~CEa150, y = ~CEa150E, marker = list(size = 10, color = 'cyan',line = list(color = 'cyan3',width = 2)))## Modelo 1
cor <- cor %>% layout(title = "Modelo 1: Autoregresivo puro para CEa150", 
yaxis = list(zeroline = FALSE),
xaxis = list(zeroline = FALSE)) 
cor 
corr1<-cor(DF2$CEa150,DF2$CEa150E);corr1 # Correlación entre lo estimado y lo observado
## [1] 0.66427
library(Metrics)
compp1=rmse(DF2$CEa150,DF2$CEa150E);compp1 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 0.63378
rCEa150=map_CEa150$fit$residuals # Residuales modelo 1: diferencia entre lo observado y lo estimado
n_map2=shapiro.test(rCEa150) # Normalidad de los residuales
n_map2.=ifelse(n_map2$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales
# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(rCEa150)), col="blue", pch=19)

Se realiza el Indice de Moran Montecarlo para los residuales, esto con el fin de encontrar si tienen dependencia espacial.

library(spdep)
McrCEa150=moran.mc(rCEa150, Wve, nsim=2000) 
McrCEa150mc=ifelse(McrCEa150$p.value<0.05, "Dependencia espacial", "Independencia espacial") # Dependencia o independencia de los residuales. ¿se queda u o epsilon?

Se evidencia que los residuales tienen dependencia espacial.

library(plotly)
figg1.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF2$CEa150, opacity = 0.5))
figg1.1<-layout(figg1.1, title= "Modelo1: CEa150 observada") #Observado

figg1.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF2$CEa150E, opacity = 0.5))
figg1.2<-layout(figg1.2, title= "Modelo1: CEa150 estimada") 
# Estimado 
figg1.3<-subplot(figg1.1, figg1.2)%>% hide_legend()
figg1.3<-layout(figg1.3, title= "Modelo1: CEa150 observada y estimada");figg1.3
figg1.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (DF2$CEa150E-DF2$CEa150)*10, opacity = 0.5))
figg1.4<-layout(figg1.4, title= "Modelo1: Diferencia entre la CEa150 observada y estimada");figg1.4

SARAR: Modelo Autorregresivo

El modelo es doblemente auto regresivo tanto para la respuesta como para los residuales mas no tiene variables explicativas

\[Y = \lambda W Y + u\\ u=\rho W u+ \epsilon\]

msarar_CEa150=sacsarlm(CEa150~1,data=datat,listw=Wve)  ## Modelo
summary(msarar_CEa150)
## 
## Call:sacsarlm(formula = CEa150 ~ 1, data = datat, 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.34 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 ajusta tanto a lambda como a rho

CEa150E2=msarar_CEa150$fitted.values # Estimado
rCEa150.2=msarar_CEa150$residuals # Normalidad de residuales
n_sarar2=shapiro.test(rCEa150.2)
n_sarar2.=ifelse(n_sarar2$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado

plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(rCEa150.2)), col="blue", pch=19) 
## Warning in plot.window(...): "title" não é um parâmetro gráfico
## Warning in plot.xy(xy, type, ...): "title" não é um parâmetro gráfico
## Warning in axis(side = side, at = at, labels = labels, ...): "title" não é um
## parâmetro gráfico

## Warning in axis(side = side, at = at, labels = labels, ...): "title" não é um
## parâmetro gráfico
## Warning in box(...): "title" não é um parâmetro gráfico
## Warning in title(...): "title" não é um parâmetro gráfico

mc.sarar.CEa150 = moran.mc(rCEa150.2, Wve, nsim = 2000)
mc.sarar.CEa150mc=ifelse(mc.sarar.CEa150$p.value<0.05, "Dependencia espacial", "Independencia espacial") ## Dependencia, se debe rompeer en dos partes en rho

Correlación y comparación

library(Metrics)
compp2=rmse(DF2$CEa150,CEa150E2) # Comparar valores observados y estimados, mas pequeño mejor ajuste

crr2 <- plot_ly(x = ~CEa150, y = ~CEa150E2,marker = list(size = 10,
color = 'rgba(255, 182, 193, .9)',line =list(color ='rgba(152, 0, 0, .8)',
width = 2)))

crr2 <- crr2 %>% layout(title = "Modelo 2: ARAR para CEa150", yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) ## Modelo 4
crr2 ## Modelo 2
## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter
## No scatter mode specifed:
##   Setting the mode to markers
##   Read more about this attribute -> https://plot.ly/r/reference/#scatter-mode
corr2<-cor(CEa150,CEa150E2);corr2
## [1] 0.68242
figg2.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF2$CEa150, opacity = 0.5))
figg2.1<-layout(figg2.1, title= "Modelo2: CEa150 observada") #Observado

figg2.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa150E2, opacity = 0.5))
figg2.2<-layout(figg2.2, title= "Modelo2: CEa150 estimada") 
# Estimado 
figg2.3<-subplot(figg2.1, figg2.2)%>% hide_legend()
figg2.3<-layout(figg2.3, title= "Modelo2: CEa150 observada y estimada");figg2.3
figg2.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa150E2-DF2$CEa150)*10, opacity = 0.5))
figg2.4<-layout(figg2.4, title= "Modelo1: Diferencia entre la CEa150 observada y estimada");figg2.4

Modelo espacial del error

Como se calculo en CEa75, se realiza el modelo con todas las variables y se dejan las que presenten una mayor significancia.

# Modelos que permitan involucrar variables explicativas (X)
mser1.2=errorsarlm(formula=CEa150~CEa75+DEM+Z+SLOPE+NDVI,data=datat,listw=Wve)
summary(mser1.2) # Que variables se pueden implementar
## 
## Call:errorsarlm(formula = CEa150 ~ CEa75 + DEM + Z + SLOPE + NDVI, 
##     data = datat, listw = Wve)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.46790 -0.32241 -0.02153  0.36060  1.99578 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##               Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 45.5088605  2.7607167 16.4844 < 2.2e-16
## CEa75        0.2959247  0.0286368 10.3337 < 2.2e-16
## DEM         -0.0093728  0.0123588 -0.7584 0.4482181
## Z           -0.1327355  0.0171932 -7.7202 1.155e-14
## SLOPE        0.0518339  0.0144655  3.5833 0.0003393
## NDVI        -1.1627276  1.1220178 -1.0363 0.3000703
## 
## 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.42 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)
mser2.2=errorsarlm(formula=CEa150~CEa75+SLOPE+Z,data=datat,listw=Wve)
summary(mser2.2) # Mejor ajuste con CEa75, SLOPE y Z
## 
## Call:errorsarlm(formula = CEa150 ~ CEa75 + SLOPE + Z, data = datat, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.470171 -0.318709 -0.014481  0.355224  2.014963 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error  z value  Pr(>|z|)
## (Intercept) 44.749718   2.699784  16.5753 < 2.2e-16
## CEa75        0.297488   0.028368  10.4868 < 2.2e-16
## SLOPE        0.052248   0.014422   3.6227 0.0002915
## Z           -0.143289   0.013322 -10.7558 < 2.2e-16
## 
## Lambda: 0.96928, LR test value: 50.495, p-value: 1.1945e-12
## Asymptotic standard error: 0.021655
##     z-value: 44.761, p-value: < 2.22e-16
## Wald statistic: 2003.5, p-value: < 2.22e-16
## 
## Log likelihood: -240.18 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)
### Modelo mser 1 menor, presenta un menor AIC
resm1.2=mser2.2$residuals
n_mser_2=shapiro.test(resm1.2) # Normalidad de los residuales
n_mser_2.=ifelse(n_mser_2$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales
# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(resm1.2)), col="blue", pch=19)

mc.mser1.CEa150 = moran.mc(resm1.2, Wve, nsim = 2000)
mc.mser1.CEa150mc =ifelse(mc.mser1.CEa150$p.value<0.05, "Dependencia espacial", "Independencia espacial") ## Dependencia, se debe rompeer en dos partes en rho

CEa150E3= mser2.2$fitted.values #Estimado
# correlacion
crr3<-plot_ly(x = ~CEa150, y = ~CEa150E3,marker = list(size = 10, color = 'cyan',line = list(color = 'cyan3',width = 2)))## Modelo 4

crr3 <- crr3 %>% layout(title = "Modelo 3: Modelo espacial del error para CEa150", yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) 
crr3## Modelo 3
compp3=rmse(CEa150,CEa150E3);compp3 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 0.51604
corr3<-cor(CEa75,CEa75E3);corr3
## [1] 0.82275

Muestra dependencia espacial por medio del Moran Montecarlo ya que presenta un p.valor pequeño. Asimismo, se muestra una significacia considerable en el NDVI, lo que indica que esta variable debe aparecer en el modelo junto con la “Z”.

figg3.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF2$CEa150, opacity = 0.5))
figg3.1<-layout(figg3.1, title= "Modelo3: CEa75 observada") #Observado

figg3.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa150E3, opacity = 0.5))
figg3.2<-layout(figg3.2, title= "Modelo3: CEa150 estimada") 
# Estimado 
figg3.3<-subplot(figg3.1, figg3.2)%>% hide_legend()
figg3.3<-layout(figg3.3, title= "Modelo3: CEa150 observada y estimada");figg3.3
figg3.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa150E3-DF2$CEa150)*10, opacity = 0.5))
figg3.4<-layout(figg3.4, title= "Modelo3: Diferencia entre la CEa75 observada y estimada");fig3.4

Modelo con SACSAR

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

Las variables explicativas que se van a implementar en el modelo son las utilizadas en el modelo anterior CEa75, Z y SLOPE.

mlrho.CEa150=sacsarlm(formula=CEa150~CEa75+Z+SLOPE,data=datat,listw=Wve)
summary(mlrho.CEa150)
## 
## Call:sacsarlm(formula = CEa150 ~ CEa75 + Z + SLOPE, data = datat, 
##     listw = Wve)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -1.443696 -0.278907 -0.020386  0.310970  1.959313 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##              Estimate Std. Error z value  Pr(>|z|)
## (Intercept) 21.463941  15.438935  1.3902  0.164454
## CEa75        0.290649   0.032187  9.0301 < 2.2e-16
## Z           -0.114901   0.017473 -6.5759 4.835e-11
## SLOPE        0.041795   0.013504  3.0949  0.001969
## 
## 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.46 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)

Se evidencia que en el modelo tanto lambda como rho son significativos, por lo que es necesario que ambos esten presentes.

### Modelo mser 1 menor, presenta un menor AIC
r.mlrho.CEa150=mlrho.CEa150$residuals
n_sac2=shapiro.test(r.mlrho.CEa150) # Normalidad de los residuales
n_sac2.=ifelse(n_sac2$p.value<0.05, "Datos no normales", "Datos normales")# Normalidad de los residuales

# Desajuste del modelo; donde se ven mas grande menos ajustado
plot(XYdata[,1], XYdata[,2], title = "Residuales del modelo", cex=(abs(r.mlrho.CEa150)), col="blue", pch=19)

moran_me2 = moran.mc(r.mlrho.CEa150, Wve, nsim = 2000)
moran_me2mc =ifelse(moran_me2$p.value<0.05, "Dependencia espacial", "Independencia espacial") # Dependencia o independencia de los residuales. ¿se queda u o epsilon?

CEa150E4= mlrho.CEa150$fitted.values # Estimado

crr4 <- plot_ly(x = ~CEa150, y = ~CEa150E4,
               marker = list(size = 10,
                             color = 'rgba(255, 182, 193, .9)', line = list(color ='rgba(152, 0, 0, .8)',width = 2)))
crr4 <- crr4 %>% layout(title = 'Modelo 4: Modelo con lambda y rho para CEa150',
         yaxis = list(zeroline = FALSE),
         xaxis = list(zeroline = FALSE)) ## Modelo 4

crr4
compp4=rmse(DF2$CEa150,CEa150E4);compp4 # Comparar valores observados y estimados, mas pequeño mejor ajuste
## [1] 0.48257
corr4<-cor(DF2$CEa150,CEa150E4);corr4
## [1] 0.76096

Muestra dependencia espacial por medio del Moran Montecarlo ya que presenta un p.valor pequeño.

figg4.1 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = DF2$CEa150, opacity = 0.5))
figg4.1<-layout(figg4.1, title= "Modelo4: CEa150 observada") #Observado

figg4.2 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = CEa150E4, opacity = 0.5))
figg4.2<-layout(figg4.2, title= "Modelo4: CEa150 estimada") 
# Estimado 
figg4.3<-subplot(figg4.1, figg4.2)%>% hide_legend()
figg4.3<-layout(figg4.3, title= "Modelo4: CEa150 observada y estimada");figg4.3
figg4.4 <- plot_ly(as.data.frame(XYdata), 
               x = ~X,
               y = ~Y,
               type = 'scatter', 
               mode = 'markers',
        marker = list(size = (CEa150E4-DF2$CEa150)*10, opacity = 0.5))
figg4.4<-layout(figg4.4, title= "Modelo4: Diferencia entre la CEa150 observada y estimada");figg4.4

Con la diferencia entre lo observado y lo esperado, se puede observar una mayor impresición en los bordes del area estudiada, esto puede deberse al efecto borde del lote, asimismo varia en cada uno de los modelos.

AIC
Modelo_=rep(c("MAP", "SARAR", "SEM", "SAC"))
AIC_=c(615.58,570.68,492.35,456.91)
Correlación_=c(cor1,cor2,cor3, cor4)
RMSE_=c(comp1,comp2, comp3, comp4)
Normalidad_res_= c(n_map2.,n_sarar2.,n_mser_2.,n_sac2.)
Moran_mc1=c(McrCEa150mc, mc.sarar.CEa150mc, mc.mser1.CEa150mc, moran_me2mc )
Tab_CEA150=data.frame(Modelo_,(round(AIC_,3)),(round(Correlación_,3)),(round(RMSE_,3)),Normalidad_res, Moran_mc1)
library(DT)
Tb_CEA150<-datatable(Tab_CEA150, class = 'cell-border stripe', editable = 'cell', rownames = FALSE, colnames = c("MODELO", "AIC", "CORRELACION", "RMSE", "NORMALIDAD_RES", "MORAN_MC1"));Tb_CEA150

Se confirma que el mejor modelo para predecir la CEa a una profundidad de 150 cm es el Modelo 4 (SACRAR), ya que tiene un AIC menor que los otros. De la misma manera, se obtuvo la comparación mas pequeña entre los valores observados y estimados en los modelos, lo cual indica un mejor ajuste, asimismo la correlación es la mejor. En todos los modelos se presenta normalidad y dependencia de los residuales.

Se analizan los valores de las curvas de nivel del área muestreada mediante interpolación Akima

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

Se concluye finalmente que el mejor modelo en la zona de estudio para predecir la CEa a 75 y 150 cm de profundidad es el modelo 4 SACRAR, que incluyte tanto a \(\lambda\) como a \(\rho\), en los cuales los residuales son normales y dependientes.

Referencias

  • Simón, M., Peralta, N. R., & Costa, J. L. (2013). Relación entre la conductividad eléctrica aparente con propiedades del suelo y nutrientes.

  • López-Lozano, R., Casterad Seral, M. A., & Herrero Isern, J. (2007). Relación entre desarrollo del cultivo, rasgos edáficos y rendimiento en una parcela de maíz mediante teledetección y SIG.

  • Munnaf, M. A., Haesaert, G., Van Meirvenne, M., & Mouazen, A. M. (2020). Map-based site-specific seeding of consumption potato production using high-resolution soil and crop data fusion. Computers and Electronics in Agriculture, 178, 105752.

  • Castro Franco, M., & Costa, J. L. CONDUCTIVIDAD ELÉCTRICA APARENTE Y SU RELACIÓN CON PROPIEDADES DEL SUELO, PARA LA DELIMITACIÓN DE ZONAS PARA MANEJO SITIO ESPECÍFICO.