´

options(digits = 12)
library(mvoutlier)
library(rworldmap)
library(rworldxtra)
library(ggmap)
library(sf)                                                                 
library(spdep)
library(corrplot)
library(ape)
library(sp)
library(MVA)
library(Hmisc)
library(normtest)
library(nortest)
library(corrplot)
library(psych) 
library(crayon)
library(pastecs)
library(readxl)
library(spatialreg)
XPABLO <- read_excel("XPABLO.xlsx")
library(clhs)
## Warning: package 'clhs' was built under R version 4.0.5
n = 0.80 *403
n = round(n,0)
dfCOS=data.frame(XPABLO)
X=as.matrix(data.frame(dfCOS[,4:15]))
#View(dfCOS)
###___matriz de coordenadas
df_michaell = data.frame(x=dfCOS$Long,
                     y=dfCOS$Lat,
                     COS=dfCOS$cos)
library(clhs)
res <- clhs(df_michaell, size = n, 
            iter = 100, progress = FALSE,
            simple = TRUE)

df_michaell = dfCOS[res,]
XY=as.matrix(df_michaell[,2:3])
COS.d=as.matrix(dist(XY, diag=T, upper=T))
##matriz inversa de distancias
COS.d.inv <-as.matrix(1/COS.d)
##asignando cero a la diag
diag(COS.d.inv) <- 0
##__la matriz de pesos basado en distancias
W=as.matrix(COS.d.inv)
##__verificando sumas por filas
SUMAS=apply(W,1,sum)
##__estandarizando la metriz
We=W/SUMAS
##verificando estandarizaci?n
##calculando ?ndice de Moran para variables
# Nueva construccion de matriz de pesos estandarizada
contnb=dnearneigh(coordinates(XY),0,380000,longlat = F)
dlist <- nbdists(contnb, XY)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W")
mod7=lagsarlm(formula=cos~z+Ca+Mg+K+Cu+Fe+Zn,data=df_michaell,listw=Wve,type="mixed")
corr.plot(x= df_michaell$cos, y= mod7$fitted.values, quan= 1/2, alpha = 0.025)

## $cor.cla
## [1] 0.715886361941
## 
## $cor.rob
## [1] 0.807971060766
shapiro.test(mod7$residuals)
## 
##  Shapiro-Wilk normality test
## 
## data:  mod7$residuals
## W = 0.9692815939, p-value = 2.40743182e-06

######################################3

## ver outliers
plot(y=dfCOS$cos, x=dfCOS$id)

boxplot(dfCOS$cos)

#ver outliers
(dfCOS$cos>2.4)
##   [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [13] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [25] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [37] FALSE FALSE FALSE FALSE FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE
##  [49] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [61] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [73] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [85] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
##  [97] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [109] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [121] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [133] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [145] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [157] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [169] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [193] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [205] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [217] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [229] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [241] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [253] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [265] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [277] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
## [289] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [301] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [313] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [325] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
## [337] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [349] FALSE  TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [361] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [373] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [385] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## [397] FALSE FALSE FALSE FALSE FALSE FALSE FALSE
table(dfCOS$cos>2.4)
## 
## FALSE  TRUE 
##   398     5

Remover valores fuera de la elipse de correlacion robusta

#Quitar outliers
dfCOS[c(42,120,288,336,338,350),c(1,15)]
##      id            cos
## 42   42 2.662200000000
## 120 120 0.964788427702
## 288 288 2.451335576011
## 336 336 2.665960492360
## 338 338 2.644800000000
## 350 350 2.601433527776
dfCOS2 <- dfCOS[-c(42,120,288,336,338,350),]
#dfCOS2
boxplot(dfCOS2$cos)

df_michaell = data.frame(x=dfCOS2$Long,
                     y=dfCOS2$Lat,
                     COS=dfCOS2$cos)
library(clhs)
res <- clhs(df_michaell, size = n, 
            iter = 100, progress = FALSE,
            simple = TRUE)

df_michaell = dfCOS2[res,]
XY=as.matrix(df_michaell[,2:3])
COS.d=as.matrix(dist(XY, diag=T, upper=T))
##matriz inversa de distancias
COS.d.inv <-as.matrix(1/COS.d)
##asignando cero a la diag
diag(COS.d.inv) <- 0
##__la matriz de pesos basado en distancias
W=as.matrix(COS.d.inv)
##__verificando sumas por filas
SUMAS=apply(W,1,sum)
##__estandarizando la metriz
We=W/SUMAS
##verificando estandarizaci?n
##calculando ?ndice de Moran para variables
# Nueva construccion de matriz de pesos estandarizada
contnb=dnearneigh(coordinates(XY),0,380000,longlat = F)
dlist <- nbdists(contnb, XY)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W")
#modelo autoregresivo puro
modelo.arp=spautolm(cos~1,data=df_michaell,listw=Wve)
summary(modelo.arp)
## 
## Call: spautolm(formula = cos ~ 1, data = df_michaell, listw = Wve)
## 
## Residuals:
##            Min             1Q         Median             3Q            Max 
## -0.89428092518 -0.29535528641  0.00953977075  0.27411708834  0.96330401274 
## 
## Coefficients: 
##                Estimate  Std. Error z value  Pr(>|z|)
## (Intercept) 1.204402807 0.167910458 7.17289 7.343e-13
## 
## Lambda: 0.876345188 LR test value: 14.2805881 p-value: 0.000157480581 
## Numerical Hessian standard error of lambda: 0.117876313 
## 
## Log likelihood: -140.9121339 
## ML residual variance (sigma squared): 0.138814045, (sigma: 0.372577569)
## Number of observations: 322 
## Number of parameters estimated: 3 
## AIC: 287.824268
res1 = modelo.arp$fit$residuals
Moran.I(res1, COS.d.inv)
## $observed
## [1] 0.0197105647257
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00540128604129
## 
## $p.value
## [1] 2.37883024552e-05

Este modelo auto regresivo puro no nos explica la dependencia espacial,puesto que la hipotesis nula de la prueba de Moran I se rechaza sobre los residuos de este modelo. Esto es debido a que tan solo se toman en cuenta la influencia de los vecinos respecto a la cantidad de carbono organico en el suelo

Segundo modelo

mod2 = sacsarlm(cos~1,data=df_michaell,listw=Wve)
summary(mod2)
## 
## Call:sacsarlm(formula = cos ~ 1, data = df_michaell, listw = Wve)
## 
## Residuals:
##          Min           1Q       Median           3Q          Max 
## -0.879832107 -0.298930680  0.010195262  0.269315080  0.963238250 
## 
## Type: sac 
## Coefficients: (numerical Hessian approximate standard errors) 
##                Estimate  Std. Error z value Pr(>|z|)
## (Intercept) 0.252979248 0.264005211 0.95824  0.33794
## 
## Rho: 0.778685824
## Approximate (numerical Hessian) standard error: 0.200594975
##     z-value: 3.88188101, p-value: 0.000103651596
## Lambda: 0.778685824
## Approximate (numerical Hessian) standard error: 0.201283825
##     z-value: 3.86859613, p-value: 0.000109463781
## 
## LR test value: 20.6042999, p-value: 3.35608632e-05
## 
## Log likelihood: -137.750277997 for sac model
## ML residual variance (sigma squared): 0.13558686, (sigma: 0.368221211)
## Number of observations: 322 
## Number of parameters estimated: 4 
## AIC: 283.500556, (AIC for lm: 300.104856)

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

res2 =mod2$residuals
Moran.I(res2, COS.d.inv)
## $observed
## [1] 0.01058332856
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00540114652675
## 
## $p.value
## [1] 0.0112050535561

Este segundo modelo no nos explica tampoco la dependencia espacial a pesar de tener el coeficiente Rho

Modelo 3

mser1=errorsarlm(formula=cos~z+Ca+Mg+K+Na+Fe+Cu+Zn+CE,data=df_michaell,listw=Wve)
summary(mser1)
## 
## Call:errorsarlm(formula = cos ~ z + Ca + Mg + K + Na + Fe + Cu + Zn + 
##     CE, data = df_michaell, listw = Wve)
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.7724767816 -0.1866900999 -0.0303596722  0.1667411310  0.8885842713 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept)  3.53840510e-01  1.89271824e-01  1.86948  0.0615556
## z            1.01247054e-03  1.54145486e-03  0.65683  0.5112916
## Ca           4.17243282e-02  5.36518760e-03  7.77686 7.3275e-15
## Mg          -9.84213783e-02  2.22957365e-02 -4.41436 1.0131e-05
## K            5.17211395e-01  1.37315959e-01  3.76658  0.0001655
## Na           1.48002575e-01  1.33999357e-01  1.10450  0.2693754
## Fe           4.20164208e-04  6.86392773e-05  6.12134 9.2793e-10
## Cu           6.44664029e-02  1.13510755e-02  5.67932 1.3523e-08
## Zn           6.27904899e-03  9.83547425e-03  0.63841  0.5232079
## CE          -3.48205502e-02  9.24249891e-02 -0.37674  0.7063639
## 
## Lambda: 0.890998139, LR test value: 14.1666008, p-value: 0.000167314244
## Asymptotic standard error: 0.0753457024
##     z-value: 11.8254673, p-value: < 2.220446e-16
## Wald statistic: 139.841677, p-value: < 2.220446e-16
## 
## Log likelihood: -33.317011536 for error model
## ML residual variance (sigma squared): 0.0710923041, (sigma: 0.266631401)
## Number of observations: 322 
## Number of parameters estimated: 12 
## AIC: 90.6340231, (AIC for lm: 102.800624)
res3 =mser1$residuals
Moran.I(res3, COS.d.inv)
## $observed
## [1] 0.017034563215
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539500935572
## 
## $p.value
## [1] 0.000187788745285

Este modelo no explica la dependencia espacial a pesar de que ya involucra las otras variables, el Zn es significativo ahora corremos el mismo modelo sin z,Zn,Na ni CE

mser1=errorsarlm(formula=cos~Ca+Mg+K+Cu+Fe,data=df_michaell,listw=Wve)
summary(mser1)
## 
## Call:errorsarlm(formula = cos ~ Ca + Mg + K + Cu + Fe, data = df_michaell, 
##     listw = Wve)
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.7892530297 -0.1885658579 -0.0298180831  0.1802472441  0.8601920233 
## 
## Type: error 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept)  4.39570211e-01  1.69684978e-01  2.59051 0.00958345
## Ca           4.13535712e-02  5.14077314e-03  8.04423 8.8818e-16
## Mg          -9.20475927e-02  2.11889625e-02 -4.34413 1.3983e-05
## K            5.27333699e-01  1.35970604e-01  3.87829 0.00010519
## Cu           6.80300462e-02  9.71381453e-03  7.00343 2.4976e-12
## Fe           4.44300831e-04  6.35744793e-05  6.98867 2.7751e-12
## 
## Lambda: 0.909026707, LR test value: 24.3878952, p-value: 7.87620375e-07
## Asymptotic standard error: 0.0632083148
##     z-value: 14.3814419, p-value: < 2.220446e-16
## Wald statistic: 206.825871, p-value: < 2.220446e-16
## 
## Log likelihood: -34.3526740894 for error model
## ML residual variance (sigma squared): 0.0714633262, (sigma: 0.267326254)
## Number of observations: 322 
## Number of parameters estimated: 8 
## AIC: 84.7053482, (AIC for lm: 107.093243)
res3 =mser1$residuals
Moran.I(res3, COS.d.inv)
## $observed
## [1] 0.0193507322175
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539523298208
## 
## $p.value
## [1] 3.12657411667e-05

Aunque quitando las variables que no eran significativas el modelo no mejora.

Como en el anterior modelo no fueron significativos el Na,Zn ni CE no lo pondremos en este modelo

mod4=lagsarlm(formula=cos~Ca+z+Mg+K+Cu+Fe,data=df_michaell,listw=Wve)
summary(mod4)
## 
## Call:lagsarlm(formula = cos ~ Ca + z + Mg + K + Cu + Fe, data = df_michaell, 
##     listw = Wve)
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.8083032760 -0.1899746911 -0.0270499048  0.1765175234  0.9040457945 
## 
## Type: lag 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept) -6.73029698e-01  1.83943931e-01 -3.65889 0.00025331
## Ca           3.66467981e-02  4.75498841e-03  7.70702 1.2879e-14
## z            2.56224684e-03  9.39340864e-04  2.72771 0.00637762
## Mg          -8.61007663e-02  2.03968566e-02 -4.22128 2.4292e-05
## K            5.37683777e-01  1.36094179e-01  3.95082 7.7883e-05
## Cu           6.73158097e-02  9.54826149e-03  7.05006 1.7883e-12
## Fe           4.32268235e-04  6.28704643e-05  6.87554 6.1757e-12
## 
## Rho: 0.788291576, LR test value: 9.07058463, p-value: 0.0025975401
## Asymptotic standard error: 0.132130627
##     z-value: 5.96600193, p-value: 2.43137399e-09
## Wald statistic: 35.593179, p-value: 2.4313741e-09
## 
## Log likelihood: -36.7830162893 for lag model
## ML residual variance (sigma squared): 0.0729741471, (sigma: 0.270137275)
## Number of observations: 322 
## Number of parameters estimated: 9 
## AIC: 91.5660326, (AIC for lm: 98.6366172)
## LM test for residual autocorrelation
## test value: 30.9233591, p-value: 2.68421489e-08
res4 = mod4$residuals
Moran.I(res4, COS.d.inv)
## $observed
## [1] 0.0232943435943
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539455295737
## 
## $p.value
## [1] 9.8003010307e-07

De igual manera este modelo no explica la dependencia espacial que tienen los datos

mod5=sacsarlm(formula=cos~z+Ca+Mg+K+Cu+Fe,data=df_michaell,listw=Wve)
summary(mod5)
## 
## Call:sacsarlm(formula = cos ~ z + Ca + Mg + K + Cu + Fe, data = df_michaell, 
##     listw = Wve)
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.7990787168 -0.1852818605 -0.0259186208  0.1713825797  0.8874972857 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept) -3.67939307e-01  5.80157184e-01 -0.63421 0.52594624
## z            1.19187249e-03  1.54130735e-03  0.77329 0.43935272
## Ca           3.94192537e-02  5.34382886e-03  7.37659 1.6231e-13
## Mg          -9.12968268e-02  2.10446885e-02 -4.33824 1.4363e-05
## K            5.22438378e-01  1.35172899e-01  3.86496 0.00011111
## Cu           6.80785267e-02  9.64289617e-03  7.05997 1.6653e-12
## Fe           4.47539720e-04  6.32109724e-05  7.08010 1.4406e-12
## 
## Rho: 0.60071083
## Asymptotic standard error: 0.448852365
##     z-value: 1.33832609, p-value: 0.180790161
## Lambda: 0.836403963
## Asymptotic standard error: 0.228087441
##     z-value: 3.66703208, p-value: 0.000245382016
## 
## LR test value: 17.616801, p-value: 0.000149472144
## 
## Log likelihood: -32.5099080918 for sac model
## ML residual variance (sigma squared): 0.0706591919, (sigma: 0.265817968)
## Number of observations: 322 
## Number of parameters estimated: 10 
## AIC: 85.0198162, (AIC for lm: 98.6366172)
res5 = mod5$residuals
Moran.I(res5, COS.d.inv)
## $observed
## [1] 0.0105547111964
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539468864391
## 
## $p.value
## [1] 0.0112778576741

Este modelo sigue sin explicar la dependencia, pero mejora al tener los coeficientes Lambda y Rho. Para este modelo no es significativo la z. Asi que vamos a volvel a correr el modelo sin esta variable

mod5=sacsarlm(formula=cos~Ca+Mg+K+Cu+Fe,data=df_michaell,listw=Wve)
summary(mod5)
## 
## Call:sacsarlm(formula = cos ~ Ca + Mg + K + Cu + Fe, data = df_michaell, 
##     listw = Wve)
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.8152020225 -0.1847059156 -0.0268585955  0.1802157428  0.8632258075 
## 
## Type: sac 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept) -2.73641066e-01  4.59720103e-01 -0.59523 0.55168702
## Ca           3.96634286e-02  5.26500461e-03  7.53341 4.9516e-14
## Mg          -9.24574128e-02  2.10363876e-02 -4.39512 1.1071e-05
## K            5.19530502e-01  1.35224291e-01  3.84199 0.00012204
## Cu           6.79062428e-02  9.65372416e-03  7.03420 2.0040e-12
## Fe           4.45710437e-04  6.32144714e-05  7.05077 1.7795e-12
## 
## Rho: 0.594694699
## Asymptotic standard error: 0.383652024
##     z-value: 1.55008878, p-value: 0.121120209
## Lambda: 0.873001232
## Asymptotic standard error: 0.15662158
##     z-value: 5.5739524, p-value: 2.49023804e-08
## 
## LR test value: 27.5034027, p-value: 1.06588902e-06
## 
## Log likelihood: -32.7949203294 for sac model
## ML residual variance (sigma squared): 0.0706673144, (sigma: 0.265833245)
## Number of observations: 322 
## Number of parameters estimated: 9 
## AIC: 83.5898407, (AIC for lm: 107.093243)
res5 = mod5$residuals
Moran.I(res5, COS.d.inv)
## $observed
## [1] 0.0126742944311
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539484145048
## 
## $p.value
## [1] 0.00342481837269

Al remover la variable z mejora el modelo pero no es suficiente para explicar la dependencia espacial de los datos

SAC Spatial Autocorrelation Model

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

mod6=sacsarlm(formula=cos~z+Ca+Mg+K+Cu+Zn+Fe,data=df_michaell,listw=Wve,type="mixed")
summary(mod6)
## 
## Call:
## sacsarlm(formula = cos ~ z + Ca + Mg + K + Cu + Zn + Fe, data = df_michaell, 
##     listw = Wve, type = "mixed")
## 
## Residuals:
##          Min           1Q       Median           3Q          Max 
## -0.733261962 -0.176199629 -0.027086356  0.161946639  0.838577039 
## 
## Type: sacmixed 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept)  4.33922492e-01  8.12443210e-01  0.53410 0.59327527
## z           -3.43265316e-03  3.07397229e-03 -1.11668 0.26412981
## Ca           4.17627760e-02  5.97404156e-03  6.99071 2.7351e-12
## Mg          -9.79590812e-02  2.23207205e-02 -4.38871 1.1403e-05
## K            4.86781567e-01  1.35877483e-01  3.58250 0.00034032
## Cu           6.24148409e-02  1.14495424e-02  5.45130 5.0004e-08
## Zn           1.16909311e-02  9.85107232e-03  1.18677 0.23531940
## Fe           4.03241465e-04  6.92033686e-05  5.82691 5.6465e-09
## lag.z        9.78414381e-03  8.82545604e-03  1.10863 0.26759088
## lag.Ca      -3.75010694e-02  4.90488147e-02 -0.76457 0.44452986
## lag.Mg      -1.26877225e-01  2.94715382e-01 -0.43051 0.66682641
## lag.K        2.03139675e+00  2.33479514e+00  0.87005 0.38427115
## lag.Cu       1.58934147e-01  1.93871555e-01  0.81979 0.41233529
## lag.Zn      -3.57709681e-01  1.53467872e-01 -2.33084 0.01976158
## lag.Fe      -1.28146336e-03  1.01950096e-03 -1.25695 0.20877116
## 
## Rho: 0.578268492
## Asymptotic standard error: 1.09228164
##     z-value: 0.529413361, p-value: 0.596518732
## Lambda: 0.314867038
## Asymptotic standard error: 1.55540876
##     z-value: 0.202433628, p-value: 0.839577741
## 
## LR test value: 35.4041016, p-value: 5.05543907e-05
## 
## Log likelihood: -23.612395789 for sacmixed model
## ML residual variance (sigma squared): 0.0674985622, (sigma: 0.259804854)
## Number of observations: 322 
## Number of parameters estimated: 18 
## AIC: 83.2247916, (AIC for lm: 100.628893)

GNS : General Nesting Spatial Model

\[Y = \lambda W Y +X \beta_1 +W X \beta_2+ u\\u =\rho W u +\epsilon\] \[Y = \lambda W Y +X \beta_1 +W X \beta_2+ u\\u =\rho W u +\epsilon\]

SDE Spatial Durbin Error Model

\[Y = X \beta_1 +W X \beta_2+ u\\u =\rho W u +\epsilon\]

mod7=lagsarlm(formula=cos~z+Ca+Mg+K+Cu+Fe+Zn,data=df_michaell,listw=Wve,type="mixed")
summary(mod7)
## 
## Call:
## lagsarlm(formula = cos ~ z + Ca + Mg + K + Cu + Fe + Zn, data = df_michaell, 
##     listw = Wve, type = "mixed")
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.7341936106 -0.1800899888 -0.0264019281  0.1612884498  0.8302475135 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept)  4.66774591e-01  7.16839104e-01  0.65116 0.51494533
## z           -3.61004375e-03  3.01014395e-03 -1.19929 0.23041414
## Ca           4.17661158e-02  5.95732569e-03  7.01088 2.3681e-12
## Mg          -9.72521785e-02  2.23003375e-02 -4.36102 1.2946e-05
## K            4.84971444e-01  1.35407696e-01  3.58156 0.00034154
## Cu           6.26711225e-02  1.14504061e-02  5.47327 4.4182e-08
## Fe           4.02305800e-04  6.93010678e-05  5.80519 6.4293e-09
## Zn           1.17856657e-02  9.78338020e-03  1.20466 0.22833384
## lag.z        9.58801267e-03  6.64867485e-03  1.44209 0.14927592
## lag.Ca      -3.97606976e-02  4.15021398e-02 -0.95804 0.33804274
## lag.Mg      -1.41974225e-01  2.67761388e-01 -0.53023 0.59595479
## lag.K        2.03206431e+00  1.95580105e+00  1.03899 0.29880782
## lag.Cu       1.54117315e-01  1.59182473e-01  0.96818 0.33295439
## lag.Fe      -1.42447206e-03  8.57574756e-04 -1.66105 0.09670399
## lag.Zn      -3.55749916e-01  1.29548734e-01 -2.74607 0.00603139
## 
## Rho: 0.666070899, LR test value: 3.75059679, p-value: 0.0527886606
## Asymptotic standard error: 0.206849779
##     z-value: 3.22007063, p-value: 0.00128159016
## Wald statistic: 10.3688549, p-value: 0.00128159016
## 
## Log likelihood: -23.7721932925 for mixed model
## ML residual variance (sigma squared): 0.0675215955, (sigma: 0.259849178)
## Number of observations: 322 
## Number of parameters estimated: 17 
## AIC: 81.5443866, (AIC for lm: 83.2949834)
## LM test for residual autocorrelation
## test value: 1.04428794, p-value: 0.306826304
res7 = mod7$residuals
Moran.I(res7, COS.d.inv)
## $observed
## [1] 0.00297771592969
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539541771784
## 
## $p.value
## [1] 0.258776342375

Este modelo sí explica la dependencia espacial, el coeficiente Rho es significativo y no se rechaza la hipotesis nula de la prueba de Moran I. Pero el modelo tiene variables que no son significativas, ahora vamos a correr el modelo sin z para intentar que mejore.

mod7.2=lagsarlm(formula=cos~Ca+Mg+K+Cu+Fe+Zn,data=df_michaell,listw=Wve,type="mixed")
summary(mod7.2)
## 
## Call:lagsarlm(formula = cos ~ Ca + Mg + K + Cu + Fe + Zn, data = df_michaell, 
##     listw = Wve, type = "mixed")
## 
## Residuals:
##           Min            1Q        Median            3Q           Max 
## -0.7182811285 -0.1872756782 -0.0178951865  0.1682243274  0.8581275307 
## 
## Type: mixed 
## Coefficients: (asymptotic standard errors) 
##                    Estimate      Std. Error  z value   Pr(>|z|)
## (Intercept)  1.16781804e+00  3.52191863e-01  3.31586 0.00091362
## Ca           4.13537829e-02  5.95929205e-03  6.93938 3.9382e-12
## Mg          -9.38202054e-02  2.21965034e-02 -4.22680 2.3704e-05
## K            4.72770216e-01  1.35443555e-01  3.49053 0.00048206
## Cu           6.28737118e-02  1.14704185e-02  5.48138 4.2202e-08
## Fe           4.00466312e-04  6.89262422e-05  5.81007 6.2447e-09
## Zn           1.17241055e-02  9.80010029e-03  1.19633 0.23156974
## lag.Ca      -4.29547244e-02  4.00037203e-02 -1.07377 0.28292657
## lag.Mg      -1.96238925e-02  2.35747123e-01 -0.08324 0.93365969
## lag.K        1.33671771e+00  1.64542803e+00  0.81238 0.41657192
## lag.Cu       1.12034144e-01  1.49333794e-01  0.75023 0.45311840
## lag.Fe      -1.36852575e-03  8.27734903e-04 -1.65334 0.09826205
## lag.Zn      -4.06415616e-01  1.23244447e-01 -3.29764 0.00097502
## 
## Rho: 0.744779829, LR test value: 6.47093674, p-value: 0.010965272
## Asymptotic standard error: 0.163120502
##     z-value: 4.56582601, p-value: 4.97531318e-06
## Wald statistic: 20.8467671, p-value: 4.97531318e-06
## 
## Log likelihood: -24.7889430019 for mixed model
## ML residual variance (sigma squared): 0.067823763, (sigma: 0.260429958)
## Number of observations: 322 
## Number of parameters estimated: 15 
## AIC: 79.577886, (AIC for lm: 84.0488227)
## LM test for residual autocorrelation
## test value: 0.897001057, p-value: 0.343587113
res7.2 = mod7.2$residuals
Moran.I(res7.2, COS.d.inv)
## $observed
## [1] 0.00334319752884
## 
## $expected
## [1] -0.00311526479751
## 
## $sd
## [1] 0.00539545976014
## 
## $p.value
## [1] 0.231299547964

El modelo no mejora al retirar la elevacion

shapiro.test(res7) # Son normales 
## 
##  Shapiro-Wilk normality test
## 
## data:  res7
## W = 0.9916485445, p-value = 0.0661522561
#plot(df_michaell$cos, mod7$fitted.values)
corr.plot(x= df_michaell$cos, y= mod7$fitted.values, quan= 1/2, alpha = 0.025)

## $cor.cla
## [1] 0.735010080331
## 
## $cor.rob
## [1] 0.807146701434
library(ggplot2)
library(plotly)
## Warning: package 'plotly' was built under R version 4.0.5
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:crayon':
## 
##     style
## The following object is masked from 'package:Hmisc':
## 
##     subplot
## The following object is masked from 'package:ggmap':
## 
##     wind
## 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
library(viridis)
## Warning: package 'viridis' was built under R version 4.0.5
## Loading required package: viridisLite
## Warning: package 'viridisLite' was built under R version 4.0.5
dataq=df_michaell
dataq$ajustados =mod7$fitted.values 
q<-ggplot(dataq,aes(x=Long,y=Lat,color=ajustados))+
  geom_point(size = 3) +
  ggtitle("Predicciones del modelo 7")+
  labs(x="x",y="y")+
  scale_colour_gradient(low = "yellow", high = "red", na.value = NA)
ggplotly(q)