library(rworldmap)
## Warning: package 'rworldmap' was built under R version 4.0.5
## Loading required package: sp
## Warning: package 'sp' was built under R version 4.0.5
## ### Welcome to rworldmap ###
## For a short introduction type : vignette('rworldmap')
library(rworldxtra)
## Warning: package 'rworldxtra' was built under R version 4.0.5
library(ggmap)
## Warning: package 'ggmap' was built under R version 4.0.5
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.0.5
## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.
## Please cite ggmap if you use it! See citation("ggmap") for details.
library(sf)
## Warning: package 'sf' was built under R version 4.0.5
## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1
library(spdep)
## Warning: package 'spdep' was built under R version 4.0.5
## Loading required package: spData
## Warning: package 'spData' was built under R version 4.0.5
## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`
library(ape)
## Warning: package 'ape' was built under R version 4.0.5
## Registered S3 method overwritten by 'ape':
## method from
## plot.mst spdep
library(sp)
library(MVA)
## Warning: package 'MVA' was built under R version 4.0.5
## Loading required package: HSAUR2
## Warning: package 'HSAUR2' was built under R version 4.0.5
## Loading required package: tools
library(Hmisc)
## Warning: package 'Hmisc' was built under R version 4.0.5
## Loading required package: lattice
## Loading required package: survival
## Loading required package: Formula
##
## Attaching package: 'Hmisc'
## The following object is masked from 'package:ape':
##
## zoom
## The following objects are masked from 'package:base':
##
## format.pval, units
library(normtest)
library(nortest)
library(corrplot)
## corrplot 0.90 loaded
library(psych)
## Warning: package 'psych' was built under R version 4.0.5
##
## Attaching package: 'psych'
## The following object is masked from 'package:Hmisc':
##
## describe
## The following objects are masked from 'package:ggplot2':
##
## %+%, alpha
library(crayon)
## Warning: package 'crayon' was built under R version 4.0.4
##
## Attaching package: 'crayon'
## The following object is masked from 'package:psych':
##
## %+%
## The following object is masked from 'package:ggplot2':
##
## %+%
library(pastecs)
## Warning: package 'pastecs' was built under R version 4.0.5
library(readxl)
## Warning: package 'readxl' was built under R version 4.0.4
library(clhs)
## Warning: package 'clhs' was built under R version 4.0.5
library(spatialreg)
## Warning: package 'spatialreg' was built under R version 4.0.5
## Loading required package: Matrix
##
## Attaching package: 'spatialreg'
## The following objects are masked from 'package:spdep':
##
## as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, can.be.simmed, cheb_setup,
## create_WX, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
## errorsarlm, get.ClusterOption, get.coresOption, get.mcOption,
## get.VerboseOption, get.ZeroPolicyOption, GMargminImage, GMerrorsar,
## griffith_sone, gstsls, Hausman.test, impacts, intImpacts,
## Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
## lextrS, lextrW, lmSLX, LU_prepermutate_setup, LU_setup,
## Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
## mom_calc_int2, moments_setup, powerWeights, sacsarlm,
## SE_classic_setup, SE_interp_setup, SE_whichMin_setup,
## set.ClusterOption, set.coresOption, set.mcOption,
## set.VerboseOption, set.ZeroPolicyOption, similar.listw, spam_setup,
## spam_update_setup, SpatialFiltering, spautolm, spBreg_err,
## spBreg_lag, spBreg_sac, stsls, subgraph_eigenw, trW
data = read_excel("~/R/Computacion/XPABLO.xlsx")
Matriz de coordenadas
df=data.frame(data)
X=as.matrix(data.frame(df[,4:15])); head(X)
## z MO Ca Mg K Na CICE CE
## [1,] 120 2.089051 7.828356 1.5564959 0.1751623 0.2912804 9.851295 0.1299864
## [2,] 119 1.649450 3.952001 0.7713688 0.4958336 0.1364539 5.355657 0.1257592
## [3,] 111 1.647495 5.877168 1.2313661 0.2734062 0.1347315 7.516672 0.2874496
## [4,] 114 2.476012 5.617296 1.1285410 0.2171382 0.1628644 7.125839 0.4147940
## [5,] 115 3.005717 11.439130 2.3604502 0.5010103 0.2915674 14.592158 0.2694840
## [6,] 109 1.929737 7.486470 1.5635561 0.2435279 0.1154978 9.409052 0.4095100
## Fe Cu Zn cos
## [1,] 133.28 3.47 1.686 1.2116495
## [2,] 29.73 1.46 1.402 0.9566813
## [3,] 237.27 4.33 2.589 0.9555473
## [4,] 331.43 4.23 4.159 1.4360867
## [5,] 281.18 3.82 2.946 1.7433161
## [6,] 258.28 3.40 3.430 1.1192475
Muestreo espacial
n = 0.80 *403
n = round(n,0)
df_ = data.frame(x=df$Long,
y=df$Lat,
K=df$K)
res <- clhs(df_, size = n,
iter = 100, progress = FALSE,
simple = TRUE)
df_ = df_[res,]; head(df)
## id Long Lat z MO Ca Mg K Na
## 1 1 -72.58261 8.078605 120 2.089051 7.828356 1.5564959 0.1751623 0.2912804
## 2 2 -72.57632 8.078605 119 1.649450 3.952001 0.7713688 0.4958336 0.1364539
## 3 3 -72.58261 8.084898 111 1.647495 5.877168 1.2313661 0.2734062 0.1347315
## 4 4 -72.57632 8.084898 114 2.476012 5.617296 1.1285410 0.2171382 0.1628644
## 5 5 -72.58261 8.091191 115 3.005717 11.439130 2.3604502 0.5010103 0.2915674
## 6 6 -72.57632 8.091191 109 1.929737 7.486470 1.5635561 0.2435279 0.1154978
## CICE CE Fe Cu Zn cos mod1 mod2 mod3 mod4
## 1 9.851295 0.1299864 133.28 3.47 1.686 1.2116495 1.25 1.128 1.194 1.23743
## 2 5.355657 0.1257592 29.73 1.46 1.402 0.9566813 1.27 1.108 1.200 1.15116
## 3 7.516672 0.2874496 237.27 4.33 2.589 0.9555473 1.28 1.238 1.291 1.33531
## 4 7.125839 0.4147940 331.43 4.23 4.159 1.4360867 1.25 1.218 1.270 1.35233
## 5 14.592158 0.2694840 281.18 3.82 2.946 1.7433161 1.25 1.484 1.523 1.55687
## 6 9.409052 0.4095100 258.28 3.40 3.430 1.1192475 1.28 1.228 1.267 1.31463
Primer mapa de puntos : Zona de arroz
#Rango de coordenadas
mapa.puntos = getMap(resolution = "high")
plot(mapa.puntos, xlim = c(-72.6, -72.3), ylim = c(8.0, 8.4), asp = 1)
points(df_$x, df_$y, col = "purple", cex = .6, pch = 16)
Segundo mapa de puntos : Lagunas de arroz
ggplot(df_, aes(x, y)) + geom_point(colour="darkgreen")+ ggtitle("Lagunas de arróz")+labs(y ="Latitud",x ="Longitud")
Comparando población con muestras
plot(df$Long, df$Lat)
points(df_$x, df_$y, pch=16, col = 10* df_$K)
XY=as.matrix(df[,2:3])
k.d=as.matrix(dist(XY, diag=T, upper=T))
MatrĆz inversa de distancias
k.d.inv <-as.matrix(1/k.d)
Asignando 0 a la diagonal
diag(k.d.inv) <- 0
Matriz de peso basada en distancias
w=as.matrix(k.d.inv)
Verificando sumas por filas
sumas = apply(w, 1, sum); head(sumas)
## 1 2 3 4 5 6
## 3644.500 3775.069 4007.650 4218.188 4183.214 4563.151
Nueva construcción de matriz de pesos estandarizados
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"); Wve
## Characteristics of weights list object:
## Neighbour list object:
## Number of regions: 403
## Number of nonzero links: 162006
## Percentage nonzero weights: 99.75186
## Average number of links: 402
##
## Weights style: W
## Weights constants summary:
## n nn S0 S1 S2
## W 403 162409 403 4.967671 1617.208
Primer modelo \[Y = \lambda W Y + \epsilon\]
Modelo autoregresivo puro
modelo.arp=spautolm(K~1,data=df,listw=Wve)
summary(modelo.arp)
##
## Call: spautolm(formula = K ~ 1, data = df, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.261475 -0.106784 -0.022284 0.075295 0.551268
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20872 0.10283 2.0297 0.04239
##
## Lambda: 0.93328 LR test value: 31.281 p-value: 2.2325e-08
## Numerical Hessian standard error of lambda: 0.065028
##
## Log likelihood: 224.5228
## ML residual variance (sigma squared): 0.018969, (sigma: 0.13773)
## Number of observations: 403
## Number of parameters estimated: 3
## AIC: -443.05
res1 = modelo.arp$fit$residuals
Moran.I(res1, k.d.inv)
## $observed
## [1] 0.02593183
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004258952
##
## $p.value
## [1] 2.50866e-11
shapiro.test(res1)
##
## Shapiro-Wilk normality test
##
## data: res1
## W = 0.9606, p-value = 6.342e-09
ad.test(res1)
##
## Anderson-Darling normality test
##
## data: res1
## A = 4.0015, p-value = 5.648e-10
sf.test(res1)
##
## Shapiro-Francia normality test
##
## data: res1
## W = 0.96077, p-value = 5.08e-08
cvm.test(res1)
##
## Cramer-von Mises normality test
##
## data: res1
## W = 0.65484, p-value = 1.342e-07
Segundo modelo \[Y = \lambda W Y + u \\ u= \rho Wu + \epsilon\]
mod2 = sacsarlm(K~1,data=df,listw=Wve)
## Warning in sacsarlm(K ~ 1, data = df, listw = Wve): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## nĆŗmero de condición recĆproco = 5.20016e-19 - using numerical Hessian.
summary(mod2)
##
## Call:sacsarlm(formula = K ~ 1, data = df, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.282232 -0.099132 -0.021571 0.074928 0.550390
##
## Type: sac
## Coefficients: (numerical Hessian approximate standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.015328 0.064853 0.2364 0.8132
##
## Rho: 0.85709
## Approximate (numerical Hessian) standard error: 0.13773
## z-value: 6.2229, p-value: 4.8818e-10
## Lambda: 0.85709
## Approximate (numerical Hessian) standard error: 0.13556
## z-value: 6.3224, p-value: 2.5759e-10
##
## LR test value: 43.362, p-value: 3.8378e-10
##
## Log likelihood: 230.5632 for sac model
## ML residual variance (sigma squared): 0.018325, (sigma: 0.13537)
## Number of observations: 403
## Number of parameters estimated: 4
## AIC: -453.13, (AIC for lm: -413.76)
res2 =mod2$residuals
Moran.I(res2, k.d.inv)
## $observed
## [1] 0.01220681
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004258454
##
## $p.value
## [1] 0.0005592702
Modelo 3 - Spatial lag model \[Y = \lambda W Y +X \beta + \epsilon \]
mser1=errorsarlm(formula=K~Ca+Mg+Cu+Zn,data=df,listw=Wve)
summary(mser1)
##
## Call:errorsarlm(formula = K ~ Ca + Mg + Cu + Zn, data = df, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.260048 -0.071225 -0.017940 0.049051 0.569193
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.0299800 0.0269814 1.1111 0.26651
## Ca 0.0082847 0.0017826 4.6475 3.361e-06
## Mg 0.0636701 0.0065734 9.6860 < 2.2e-16
## Cu -0.0015472 0.0040780 -0.3794 0.70440
## Zn 0.0113919 0.0032678 3.4861 0.00049
##
## Lambda: 0.75653, LR test value: 6.0434, p-value: 0.013959
## Asymptotic standard error: 0.15953
## z-value: 4.7423, p-value: 2.1129e-06
## Wald statistic: 22.49, p-value: 2.1129e-06
##
## Log likelihood: 330.6833 for error model
## ML residual variance (sigma squared): 0.01128, (sigma: 0.10621)
## Number of observations: 403
## Number of parameters estimated: 7
## AIC: -647.37, (AIC for lm: -643.32)
res3 =mser1$residuals
Moran.I(res3, k.d.inv)
## $observed
## [1] 0.007194041
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004241071
##
## $p.value
## [1] 0.02244095
Spatial error model
mod4=lagsarlm(formula=K~Ca+z+Mg+Cu,data=df,listw=Wve)
summary(mod4)
##
## Call:lagsarlm(formula = K ~ Ca + z + Mg + Cu, data = df, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.264163 -0.073460 -0.018022 0.056891 0.566211
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.04298810 0.07748581 -0.5548 0.57904
## Ca 0.00734479 0.00165234 4.4451 8.786e-06
## z -0.00013919 0.00037474 -0.3714 0.71031
## Mg 0.05962408 0.00654979 9.1032 < 2.2e-16
## Cu 0.00627996 0.00312743 2.0080 0.04464
##
## Rho: 0.40225, LR test value: 1.7994, p-value: 0.17978
## Asymptotic standard error: 0.2419
## z-value: 1.6629, p-value: 0.096331
## Wald statistic: 2.7653, p-value: 0.096331
##
## Log likelihood: 323.4067 for lag model
## ML residual variance (sigma squared): 0.011748, (sigma: 0.10839)
## Number of observations: 403
## Number of parameters estimated: 7
## AIC: -632.81, (AIC for lm: -633.01)
## LM test for residual autocorrelation
## test value: 6.9985, p-value: 0.0081577
res4 = mod4$residuals
Moran.I(res4, k.d.inv)
## $observed
## [1] 0.01000845
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004241978
##
## $p.value
## [1] 0.003221228
mod5=sacsarlm(formula=K~Ca+Mg+Cu+Zn,data=df,listw=Wve)
summary(mod5)
##
## Call:sacsarlm(formula = K ~ Ca + Mg + Cu + Zn, data = df, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.257121 -0.072217 -0.016213 0.048832 0.566822
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.0599480 0.1001759 -0.5984 0.5495549
## Ca 0.0083417 0.0017648 4.7266 2.284e-06
## Mg 0.0617844 0.0068616 9.0044 < 2.2e-16
## Cu -0.0021386 0.0041372 -0.5169 0.6052173
## Zn 0.0116345 0.0032709 3.5570 0.0003752
##
## Rho: 0.35886
## Asymptotic standard error: 0.39149
## z-value: 0.91667, p-value: 0.35932
## Lambda: 0.66204
## Asymptotic standard error: 0.31633
## z-value: 2.0929, p-value: 0.036362
##
## LR test value: 7.2986, p-value: 0.026009
##
## Log likelihood: 331.3109 for sac model
## ML residual variance (sigma squared): 0.011255, (sigma: 0.10609)
## Number of observations: 403
## Number of parameters estimated: 8
## AIC: -646.62, (AIC for lm: -643.32)
res5 = mod5$residuals
Moran.I(res5, k.d.inv)
## $observed
## [1] 0.004782086
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004241265
##
## $p.value
## [1] 0.08652356
\[Y = \lambda W Y +X \beta + u\\u =\rho W u +\epsilon \]
mod6=sacsarlm(formula=K~Ca+Mg+Cu+Zn,data=df,listw=Wve,type="mixed")
summary(mod6)
##
## Call:sacsarlm(formula = K ~ Ca + Mg + Cu + Zn, data = df, listw = Wve,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.270981 -0.065464 -0.018122 0.051767 0.571218
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1208978 0.2233859 0.5412 0.5883655
## Ca 0.0082244 0.0020977 3.9207 8.829e-05
## Mg 0.0620351 0.0073759 8.4105 < 2.2e-16
## Cu -0.0009405 0.0041933 -0.2243 0.8225348
## Zn 0.0116238 0.0033125 3.5091 0.0004496
## lag.Ca -0.0021975 0.0215387 -0.1020 0.9187361
## lag.Mg 0.0718716 0.2075396 0.3463 0.7291151
## lag.Cu -0.0625863 0.0742318 -0.8431 0.3991614
## lag.Zn -0.0132292 0.0565608 -0.2339 0.8150673
##
## Rho: 0.52073
## Asymptotic standard error: 2.0831
## z-value: 0.24998, p-value: 0.8026
## Lambda: 0.61822
## Asymptotic standard error: 1.7581
## z-value: 0.35163, p-value: 0.72511
##
## LR test value: 11.375, p-value: 0.077449
##
## Log likelihood: 333.3492 for sacmixed model
## ML residual variance (sigma squared): 0.011135, (sigma: 0.10552)
## Number of observations: 403
## Number of parameters estimated: 12
## AIC: -642.7, (AIC for lm: -643.32)
res6 = mod6$residuals
Moran.I(res6, k.d.inv)
## $observed
## [1] 0.003403076
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004239908
##
## $p.value
## [1] 0.164732
mod7=lagsarlm(formula=K~Ca+Mg+cos+Cu,data=df,listw=Wve,type="mixed")
summary(mod7)
##
## Call:lagsarlm(formula = K ~ Ca + Mg + cos + Cu, data = df, listw = Wve,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.276927 -0.069153 -0.013564 0.053039 0.580889
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.09751193 0.21600449 0.4514 0.6517
## Ca 0.00257805 0.00217411 1.1858 0.2357
## Mg 0.06869492 0.00734669 9.3505 < 2.2e-16
## cos 0.08234599 0.01690740 4.8704 1.114e-06
## Cu 0.00050213 0.00351144 0.1430 0.8863
## lag.Ca 0.00710751 0.01306043 0.5442 0.5863
## lag.Mg 0.00791296 0.05722090 0.1383 0.8900
## lag.cos -0.14803000 0.17367523 -0.8523 0.3940
## lag.Cu -0.03630500 0.03234406 -1.1225 0.2617
##
## Rho: 0.7184, LR test value: 4.8584, p-value: 0.027511
## Asymptotic standard error: 0.18094
## z-value: 3.9703, p-value: 7.1789e-05
## Wald statistic: 15.763, p-value: 7.1789e-05
##
## Log likelihood: 337.7557 for mixed model
## ML residual variance (sigma squared): 0.0109, (sigma: 0.1044)
## Number of observations: 403
## Number of parameters estimated: 11
## AIC: -653.51, (AIC for lm: -650.65)
## LM test for residual autocorrelation
## test value: 5.2685, p-value: 0.021715
res7 = mod7$residuals
Moran.I(res7, k.d.inv)
## $observed
## [1] 0.005774499
##
## $expected
## [1] -0.002487562
##
## $sd
## [1] 0.004239806
##
## $p.value
## [1] 0.05133266
shapiro.test(res7)
##
## Shapiro-Wilk normality test
##
## data: res7
## W = 0.92941, p-value = 7.3e-13
ad.test(res7)
##
## Anderson-Darling normality test
##
## data: res7
## A = 5.0875, p-value = 1.362e-12
sf.test(res7)
##
## Shapiro-Francia normality test
##
## data: res7
## W = 0.9266, p-value = 1.2e-11
cvm.test(res7)
##
## Cramer-von Mises normality test
##
## data: res7
## W = 0.84651, p-value = 7.615e-09
hist(res1);hist(res1);hist(res2);hist(res3);hist(res4);hist(res5);hist(res6);hist(res7)
outliers::grubbs.test(res7) # El valor mas alto es un outlier
##
## Grubbs test for one outlier
##
## data: res7
## G.310 = 5.55708, U = 0.92299, p-value = 2.958e-06
## alternative hypothesis: highest value 0.580888607775061 is an outlier
which.max(res7) # Maximo
## 310
## 310
plot(df$K,mod7$fitted.values,xlab="Observados",ylab="Estimados",col="purple")
library(mvoutlier)
## Warning: package 'mvoutlier' was built under R version 4.0.5
## Loading required package: sgeostat
## Registered S3 method overwritten by 'GGally':
## method from
## +.gg ggplot2
## sROC 0.1-2 loaded
##
## Attaching package: 'mvoutlier'
## The following object is masked _by_ '.GlobalEnv':
##
## X
corr.plot(df$K,mod7$fitted.values,
quan=1/2, alpha=0.025,
xlab="Observados",
ylab="Estimados")
## $cor.cla
## [1] 0.6892811
##
## $cor.rob
## [1] 0.7337593
data = read_excel("~/R/Computacion/XPABLO.xlsx")
df=data.frame(data)
df_ka = df
atipicosk = which(df_ka$K>0.6)
df_ka = df_ka[-atipicosk,]
XY=as.matrix(df_ka[,2:3])
k.d=as.matrix(dist(XY, diag=T, upper=T))
k.d.inv <-as.matrix(1/k.d)
diag(k.d.inv) <- 0
w=as.matrix(k.d.inv)
sumas = apply(w, 1, sum)
We = w/sumas
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")
mod7b=lagsarlm(formula = K ~ Ca + Mg + cos + Cu, data = df_ka, listw = Wve, type = "mixed")
summary(mod7b)
##
## Call:lagsarlm(formula = K ~ Ca + Mg + cos + Cu, data = df_ka, listw = Wve,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.270264 -0.062007 -0.013858 0.049083 0.381029
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.1410944 0.1889618 0.7467 0.45526
## Ca 0.0050643 0.0019703 2.5704 0.01016
## Mg 0.0590258 0.0070580 8.3629 < 2.2e-16
## cos 0.0673137 0.0151678 4.4379 9.082e-06
## Cu 0.0034673 0.0031565 1.0985 0.27200
## lag.Ca -0.0017326 0.0113559 -0.1526 0.87874
## lag.Mg 0.0195636 0.0573612 0.3411 0.73306
## lag.cos -0.1316332 0.1524876 -0.8632 0.38801
## lag.Cu -0.0430414 0.0297409 -1.4472 0.14784
##
## Rho: 0.79805, LR test value: 7.5141, p-value: 0.0061217
## Asymptotic standard error: 0.13472
## z-value: 5.9236, p-value: 3.1501e-09
## Wald statistic: 35.089, p-value: 3.1501e-09
##
## Log likelihood: 377.0249 for mixed model
## ML residual variance (sigma squared): 0.0085353, (sigma: 0.092387)
## Number of observations: 393
## Number of parameters estimated: 11
## AIC: -732.05, (AIC for lm: -726.54)
## LM test for residual autocorrelation
## test value: 10.627, p-value: 0.0011147
res7b = mod7b$residuals
Moran.I(res7b, k.d.inv) # No hay depedencia espacial
## $observed
## [1] 0.009954865
##
## $expected
## [1] -0.00255102
##
## $sd
## [1] 0.004373581
##
## $p.value
## [1] 0.004244226
shapiro.test(res7b)
##
## Shapiro-Wilk normality test
##
## data: res7b
## W = 0.96802, p-value = 1.42e-07
ad.test(res7b)
##
## Anderson-Darling normality test
##
## data: res7b
## A = 3.0519, p-value = 1.138e-07
sf.test(res7b)
##
## Shapiro-Francia normality test
##
## data: res7b
## W = 0.96666, p-value = 4.412e-07
cvm.test(res7b)
##
## Cramer-von Mises normality test
##
## data: res7b
## W = 0.51666, p-value = 1.909e-06
library(mvoutlier)
color.plot(cbind(df_ka$K,mod7b$fitted.values))
## $outliers
## 1 2 3 4 5 6 7 8 9 10 11 12 13
## FALSE TRUE FALSE FALSE FALSE FALSE FALSE TRUE TRUE FALSE FALSE FALSE FALSE
## 14 15 16 17 18 19 20 21 22 23 24 25 26
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 27 28 29 30 31 32 33 34 35 36 37 38 39
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 40 41 42 43 44 45 46 47 48 49 50 51 52
## FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 53 54 55 56 57 58 59 60 61 62 63 64 65
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 66 67 68 69 70 71 72 73 74 75 76 77 78
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 79 80 81 82 83 84 85 86 87 88 89 90 91
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 92 93 94 95 96 97 98 99 100 101 102 103 104
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 105 106 107 108 109 110 111 112 113 114 115 116 117
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE
## 118 119 120 121 122 123 124 125 126 127 128 129 130
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 131 132 133 134 135 136 138 139 140 141 142 143 144
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 145 146 147 148 149 150 151 152 153 154 155 156 157
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 158 159 160 162 163 164 165 166 167 168 169 170 171
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 172 173 174 175 176 177 178 179 180 181 182 183 184
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 185 186 187 188 189 191 192 193 194 195 196 197 198
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 199 200 201 202 203 204 205 206 207 208 209 210 211
## FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 212 213 214 215 216 217 218 219 220 221 222 223 224
## FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 225 226 228 229 230 231 232 233 234 235 236 237 238
## FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 239 240 241 242 243 244 245 246 247 248 249 250 251
## FALSE FALSE FALSE FALSE TRUE TRUE TRUE FALSE FALSE FALSE FALSE FALSE FALSE
## 252 253 254 255 257 258 259 260 261 262 263 264 265
## FALSE FALSE TRUE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 266 267 268 269 270 271 273 274 275 276 277 278 279
## FALSE FALSE TRUE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE TRUE
## 280 282 283 284 285 286 287 288 289 290 291 292 293
## TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE
## 294 295 296 297 298 299 300 301 302 303 304 305 306
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 307 308 309 311 312 313 314 315 316 317 318 319 320
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 321 322 323 324 325 326 327 328 329 330 331 332 333
## FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 334 335 336 337 338 339 340 341 342 343 344 345 346
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 347 348 349 350 351 352 353 354 355 356 357 358 359
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 360 361 362 363 364 365 366 367 368 369 370 371 372
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 374 376 377 378 379 380 381 382 383 384 385 386 387
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE TRUE FALSE FALSE FALSE FALSE TRUE
## 388 389 390 391 392 393 394 395 396 397 398 399 400
## FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
## 401 402 403
## FALSE FALSE FALSE
##
## $md
## 1 2 3 4 5 6 7
## 0.62084834 4.54859726 1.29233179 0.34145257 2.11535600 0.55774336 1.84128593
## 8 9 10 11 12 13 14
## 3.17622513 4.39017143 1.07911203 0.81677413 1.08675449 1.94567704 0.56125526
## 15 16 17 18 19 20 21
## 1.39473251 2.45069654 1.29055590 1.24594564 1.25942131 1.05246984 0.54185558
## 22 23 24 25 26 27 28
## 0.36563379 0.92720357 2.08696211 1.78409994 1.75896806 0.83863036 1.07517159
## 29 30 31 32 33 34 35
## 0.42231854 1.43597660 0.32676251 1.70826651 0.61350732 0.95496252 0.32178734
## 36 37 38 39 40 41 42
## 1.40435318 1.42053936 1.34962557 1.54686497 0.73093579 0.71797572 2.13455880
## 43 44 45 46 47 48 49
## 0.93511484 0.50590312 0.83428493 2.98272035 0.19132302 0.77432342 1.41730560
## 50 51 52 53 54 55 56
## 1.23442508 0.78709827 2.10977717 1.22559351 1.25869075 0.82044113 0.08433578
## 57 58 59 60 61 62 63
## 1.19170548 1.21593420 1.03974168 1.87872865 0.24333172 1.17048598 1.39080028
## 64 65 66 67 68 69 70
## 1.73861452 1.35106383 1.88191657 1.48997294 1.23767736 0.49218306 1.11409660
## 71 72 73 74 75 76 77
## 0.22882592 0.90909996 1.59389425 0.95454863 2.07104363 1.44140661 0.82053399
## 78 79 80 81 82 83 84
## 1.80755955 0.95668590 1.90629429 2.60910100 0.31259210 1.59938939 1.54754663
## 85 86 87 88 89 90 91
## 1.12135639 1.69565908 1.87690514 1.84055333 0.79683086 1.04192663 0.70782098
## 92 93 94 95 96 97 98
## 2.09134283 0.43247913 1.22243888 1.42576411 1.56666485 1.28205585 1.64165688
## 99 100 101 102 103 104 105
## 1.20138152 1.45143578 0.55393515 0.32058105 0.47945200 2.37760730 0.76919691
## 106 107 108 109 110 111 112
## 0.47041565 1.89642948 1.46845260 1.26106677 0.51313759 1.31283764 1.23255764
## 113 114 115 116 117 118 119
## 2.81877511 1.15481956 2.10039873 0.71915501 1.93265658 1.34269423 0.80226306
## 120 121 122 123 124 125 126
## 0.19021976 0.28924244 1.62557644 1.73893673 1.46552708 1.72255644 1.37182504
## 127 128 129 130 131 132 133
## 0.23653803 1.15480834 1.08900589 1.27092080 2.59491637 1.17151604 1.64177561
## 134 135 136 138 139 140 141
## 0.50897886 1.17742169 1.49816600 0.54878306 1.85642000 0.67308124 1.41150444
## 142 143 144 145 146 147 148
## 0.93965830 1.08039670 1.03376027 0.29598710 0.57782043 0.99950188 2.22689329
## 149 150 151 152 153 154 155
## 1.86382067 1.12543712 0.73699468 0.52669394 1.16308684 1.53070195 0.31158496
## 156 157 158 159 160 162 163
## 0.73309365 1.42014151 1.52196585 1.40643760 1.14995492 1.53924162 0.96379519
## 164 165 166 167 168 169 170
## 0.72470149 0.86483357 0.52853709 1.18222165 0.57085744 1.46811367 0.55853449
## 171 172 173 174 175 176 177
## 1.24629508 0.56491377 0.50118231 0.42292247 0.85980987 0.89885321 0.83835675
## 178 179 180 181 182 183 184
## 0.49077532 1.45731370 0.91453815 0.67764157 0.92193205 0.59388223 1.55445070
## 185 186 187 188 189 191 192
## 1.66674569 0.79817271 1.78502507 2.10733326 0.85345344 2.48534635 1.10509939
## 193 194 195 196 197 198 199
## 1.87218589 2.64742106 0.39199292 0.84276265 0.68537528 1.28510680 0.48982099
## 200 201 202 203 204 205 206
## 1.31854889 3.13198497 1.45201787 0.51669586 2.70860940 1.06716214 1.79898903
## 207 208 209 210 211 212 213
## 2.10965979 0.92272275 0.62594643 1.11717817 0.67974168 0.44280443 1.47617417
## 214 215 216 217 218 219 220
## 0.81695352 0.37012470 3.08575229 1.99288570 0.82863720 0.95055650 1.70905486
## 221 222 223 224 225 226 228
## 2.27866082 1.08965189 2.44230495 0.18356721 0.73679118 0.75393974 1.09869805
## 229 230 231 232 233 234 235
## 2.84388369 0.78641430 2.70765157 0.62696050 1.57123847 0.49444884 1.32539954
## 236 237 238 239 240 241 242
## 0.29528048 1.47453838 0.62684202 0.78336392 0.43175764 2.06735838 1.67743922
## 243 244 245 246 247 248 249
## 3.35566608 2.73753485 2.87510133 1.08551401 1.37022281 1.11683947 0.13401558
## 250 251 252 253 254 255 257
## 1.39127316 0.90751894 1.36426659 1.91767333 3.44962900 1.74009933 3.31334235
## 258 259 260 261 262 263 264
## 2.13143394 0.42167430 1.60920380 1.37477131 0.46450320 1.04427283 0.67842058
## 265 266 267 268 269 270 271
## 0.73393905 1.79678655 0.37036992 3.09594321 1.88981842 1.15720456 3.04425623
## 273 274 275 276 277 278 279
## 0.28301299 1.03582421 1.86038626 1.33837086 1.54179952 0.79490183 3.22957519
## 280 282 283 284 285 286 287
## 3.87538387 1.18846974 0.70584863 2.34294104 0.82084399 1.07409269 0.94639819
## 288 289 290 291 292 293 294
## 1.56054917 1.00816486 1.81704511 3.42384637 2.71288127 1.69367234 2.02626297
## 295 296 297 298 299 300 301
## 1.59707859 2.35286749 1.20390028 0.50389629 1.36699660 1.05588807 0.33221581
## 302 303 304 305 306 307 308
## 0.48576542 1.46283726 0.46246562 0.26649456 1.44150525 0.86685465 1.09564689
## 309 311 312 313 314 315 316
## 0.91784195 1.22168314 1.51338633 1.20449123 0.38593946 0.37628116 1.00881626
## 317 318 319 320 321 322 323
## 0.80245500 2.49004325 1.71591305 1.02813291 1.18767637 0.95432891 2.01702598
## 324 325 326 327 328 329 330
## 0.42845675 1.36133109 4.27211765 2.26026909 0.96557881 1.07462369 0.75846547
## 331 332 333 334 335 336 337
## 1.59155955 0.84114352 1.42787653 0.64787494 1.38875257 2.22374523 0.16892966
## 338 339 340 341 342 343 344
## 2.07484163 1.83552494 0.80455474 1.84076611 0.96537270 0.84357798 0.10853720
## 345 346 347 348 349 350 351
## 1.10226047 0.90713211 0.98270770 1.70691585 0.54414205 1.96386761 1.68580024
## 352 353 354 355 356 357 358
## 0.90132557 0.80807963 0.61512554 0.73276203 0.67868197 0.72372687 0.76524311
## 359 360 361 362 363 364 365
## 1.47592743 1.35619893 2.24024454 0.93930568 0.99575813 0.78788883 0.64187309
## 366 367 368 369 370 371 372
## 1.15204012 2.23880014 0.30888568 1.54057442 0.13339473 1.44349517 1.34106625
## 374 376 377 378 379 380 381
## 0.90482950 1.20638847 1.17610415 1.27073488 2.36123950 0.98070690 0.89393795
## 382 383 384 385 386 387 388
## 4.41076841 0.89020930 1.24689790 0.59596538 1.01516515 3.02551977 1.44300926
## 389 390 391 392 393 394 395
## 1.11405480 1.38042222 1.26772404 0.74931616 0.90564675 1.17417205 0.29358795
## 396 397 398 399 400 401 402
## 1.19876095 0.72339538 1.68534452 1.02572488 0.24212420 2.24782265 0.41615523
## 403
## 2.25979032
##
## $euclidean
## 1 2 3 4 5 6 7 8
## 2.6700125 4.2220862 2.8681865 2.7214101 5.1691401 2.8570800 3.5776450 3.0137952
## 9 10 11 12 13 14 15 16
## 5.5988880 2.4703354 3.7743492 4.4412405 3.0149862 3.6834799 1.2860612 2.1662896
## 17 18 19 20 21 22 23 24
## 1.4848436 1.4520555 3.0233803 2.1653050 3.5867010 2.6236334 2.3120114 1.4419269
## 25 26 27 28 29 30 31 32
## 2.2069504 1.0289041 4.0375257 2.4132054 2.8823397 1.8669189 3.2102309 1.7061344
## 33 34 35 36 37 38 39 40
## 3.7953241 4.2195465 3.1223179 1.2624946 2.6725961 1.4530263 3.9662090 2.5194485
## 41 42 43 44 45 46 47 48
## 3.3718497 5.1463488 1.8467063 2.4434268 2.0097438 2.5127881 2.9320239 2.0877439
## 49 50 51 52 53 54 55 56
## 4.4769427 4.4176885 2.8279248 1.2181003 2.6924835 1.5162768 2.5562794 3.1229840
## 57 58 59 60 61 62 63 64
## 3.9619690 4.3075351 2.0382867 1.2986187 2.7901349 1.7340253 1.4808160 1.7286892
## 65 66 67 68 69 70 71 72
## 1.4675196 2.5435166 3.5894429 4.4795078 3.3240918 3.6829857 3.1912897 2.9945963
## 73 74 75 76 77 78 79 80
## 1.0766357 1.8416452 3.1960117 1.7282311 3.0959692 4.1638843 3.2349240 4.8851293
## 81 82 83 84 85 86 87 88
## 4.3540746 2.7504527 1.4640715 1.1975282 2.9748238 0.8827726 0.6838664 4.4533376
## 89 90 91 92 93 94 95 96
## 2.4558729 3.2452153 3.9544970 4.5889585 3.5419423 3.2119753 1.4460134 1.2303394
## 97 98 99 100 101 102 103 104
## 2.0235326 2.4418797 1.5598907 1.6642299 2.9759816 3.0177956 3.5477586 5.0497588
## 105 106 107 108 109 110 111 112
## 2.7517226 2.5041310 1.7203565 1.3427832 3.5337473 3.5485648 4.4861473 1.4911018
## 113 114 115 116 117 118 119 120
## 5.0758178 1.5738519 3.8984496 3.3574039 2.7085400 4.0043973 2.1252429 2.8022552
## 121 122 123 124 125 126 127 128
## 2.8092202 1.9772547 2.2272417 1.1953057 1.6324226 1.6146494 3.3375536 1.8705335
## 129 130 131 132 133 134 135 136
## 1.6608818 1.6446760 1.5305639 1.5436182 2.9838359 2.8428069 3.7950531 3.3547039
## 138 139 140 141 142 143 144 145
## 3.1396744 3.6660159 3.3957879 1.5709590 2.3126543 2.5037314 4.3721427 2.9789715
## 146 147 148 149 150 151 152 153
## 2.5552244 3.7005658 5.8565011 5.0380077 4.1418165 2.1260603 2.3719586 3.0542845
## 154 155 156 157 158 159 160 162
## 1.5367178 3.4445052 2.9439375 4.7814616 3.6461745 2.6337235 3.8111507 5.0168384
## 163 164 165 166 167 168 169 170
## 1.8209177 2.1658768 3.7814902 3.7073692 1.6796705 2.3133081 2.0297423 3.6556097
## 171 172 173 174 175 176 177 178
## 4.5471982 3.5399516 3.6399176 2.6705091 4.1249788 4.1432742 3.8807069 3.5387616
## 179 180 181 182 183 184 185 186
## 3.8746846 2.6809191 2.1918790 3.9689156 3.5516395 3.8547240 5.1716128 3.6114987
## 187 188 189 191 192 193 194 195
## 5.1815395 4.7818833 2.6605062 6.0470509 3.9348776 5.2933262 5.2062391 3.0672058
## 196 197 198 199 200 201 202 203
## 1.9653025 3.5041712 2.7907128 3.6474942 4.2671669 6.9602829 4.6850553 3.2808292
## 204 205 206 207 208 209 210 211
## 6.0524396 3.3344433 5.3251012 5.7097417 4.0693945 3.8241490 1.6463664 2.6437776
## 212 213 214 215 216 217 218 219
## 3.1151041 4.9225972 3.3234344 3.3228445 5.4610745 4.0678326 3.8956092 4.1132945
## 220 221 222 223 224 225 226 228
## 1.2884138 5.1247719 4.4347089 5.9913677 2.8828339 3.7189586 2.9149462 3.7695194
## 229 230 231 232 233 234 235 236
## 6.2254240 3.7984534 6.3117315 3.7986216 1.1454978 3.3640058 2.2878150 2.9316659
## 237 238 239 240 241 242 243 244
## 4.3771203 3.3269953 2.2276125 3.3996055 5.1813813 4.4656102 6.2464869 5.2768863
## 245 246 247 248 249 250 251 252
## 6.6413043 4.1482154 4.7483540 2.7071066 3.1532636 4.1634169 2.1937236 1.8783687
## 253 254 255 257 258 259 260 261
## 3.7997195 6.6133191 4.3533652 5.6561721 1.8662155 3.5749565 4.4823694 4.7206973
## 262 263 264 265 266 267 268 269
## 2.7793853 2.4630753 2.7164165 3.4678332 3.1870333 2.8821787 6.6288121 5.3015117
## 270 271 273 274 275 276 277 278
## 4.3649337 6.7982442 2.9136650 2.1725578 2.3857218 1.3594141 4.3457263 3.9869700
## 279 280 282 283 284 285 286 287
## 6.4708716 6.9209710 4.4581378 2.2904866 3.9564929 3.9921056 1.6686249 2.5096320
## 288 289 290 291 292 293 294 295
## 4.6255885 2.3988527 4.4267003 6.7466300 6.3125391 1.6920361 4.7012084 1.0999479
## 296 297 298 299 300 301 302 303
## 5.5405313 1.5253908 3.2086893 1.3474801 2.1025754 2.7185476 2.7061572 1.3444750
## 304 305 306 307 308 309 311 312
## 2.4552986 2.8443905 3.4998392 1.9461729 3.0186526 3.1651973 3.7380482 1.1859319
## 313 314 315 316 317 318 319 320
## 3.2005785 3.5083730 3.4896104 2.8243778 3.5493033 3.1023099 1.6845828 4.2661733
## 321 322 323 324 325 326 327 328
## 1.6859578 4.2713634 2.8500625 3.0960034 2.0084102 3.8209576 5.2359999 3.9379084
## 329 330 331 332 333 334 335 336
## 3.5245563 2.7800005 1.4305702 1.9861289 1.3164349 3.4546507 3.9502320 5.8929983
## 337 338 339 340 341 342 343 344
## 3.0972428 5.2248674 4.9521871 2.8251926 1.2175192 2.5605186 3.2891590 3.1563647
## 345 346 347 348 349 350 351 352
## 4.3445486 2.9506154 1.8807165 5.2353861 2.7488935 4.2252858 5.1193289 3.8463463
## 353 354 355 356 357 358 359 360
## 2.2177172 3.1816246 3.9606536 3.0203007 2.4520886 3.9764175 3.7613877 1.4962766
## 361 362 363 364 365 366 367 368
## 1.4471883 3.4395445 3.9336795 2.6938592 3.3840455 4.0724212 1.7393814 3.1897974
## 369 370 371 372 374 376 377 378
## 1.5325522 3.1209418 4.7873195 2.4446095 2.7846394 4.2395631 4.1242012 1.5817403
## 379 380 381 382 383 384 385 386
## 0.7537195 2.4580688 1.9001681 5.0086114 1.9646687 2.5546639 2.7906276 3.4077181
## 387 388 389 390 391 392 393 394
## 5.0100468 4.5922988 4.3303897 1.9833110 1.7246773 3.9307135 2.1707786 1.5544813
## 395 396 397 398 399 400 401 402
## 3.1046795 3.9598388 2.2143539 4.7767389 2.3791770 2.9722903 4.3190013 3.1816523
## 403
## 2.1498471
corr.plot(df_ka$K,mod7b$fitted.values,
quan=1/2, alpha=0.025,
xlab="Observados",
ylab="Estimados")
## $cor.cla
## [1] 0.7037577
##
## $cor.rob
## [1] 0.7335219
library(akima)
## Warning: package 'akima' was built under R version 4.0.5
dataw=df_ka
dataw$ajustados =mod7b$fitted.values
akima.li <- interp(dataw$Long, dataw$Lat, dataw$ajustados,
nx = 50, ny = 50,
linear = F)
image(akima.li)
contour(akima.li,add=T)