options(digits = 12)
library(rworldmap)## Loading required package: sp## ### Welcome to rworldmap ##### For a short introduction type : vignette('rworldmap')library(rworldxtra)
library(ggmap)## Loading required package: ggplot2## 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) ## Linking to GEOS 3.9.0, GDAL 3.2.1, PROJ 7.2.1library(spdep)## Loading required package: spData## To access larger datasets in this package, install the spDataLarge
## package with: `install.packages('spDataLarge',
## repos='https://nowosad.github.io/drat/', type='source')`library(ape)## Registered S3 method overwritten by 'ape':
## method from
## plot.mst spdeplibrary(sp)
library(MVA)## Loading required package: HSAUR2## Loading required package: toolslibrary(Hmisc)## 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, unitslibrary(normtest)
library(nortest)
library(corrplot)## corrplot 0.84 loadedlibrary(psych) ##
## Attaching package: 'psych'## The following object is masked from 'package:Hmisc':
##
## describe## The following objects are masked from 'package:ggplot2':
##
## %+%, alphalibrary(crayon)##
## Attaching package: 'crayon'## The following object is masked from 'package:psych':
##
## %+%## The following object is masked from 'package:ggplot2':
##
## %+%library(pastecs)
library(readxl)
library(spatialreg)## Loading required package: Matrix## Registered S3 methods overwritten by 'spatialreg':
## method from
## residuals.stsls spdep
## deviance.stsls spdep
## coef.stsls spdep
## print.stsls spdep
## summary.stsls spdep
## print.summary.stsls spdep
## residuals.gmsar spdep
## deviance.gmsar spdep
## coef.gmsar spdep
## fitted.gmsar spdep
## print.gmsar spdep
## summary.gmsar spdep
## print.summary.gmsar spdep
## print.lagmess spdep
## summary.lagmess spdep
## print.summary.lagmess spdep
## residuals.lagmess spdep
## deviance.lagmess spdep
## coef.lagmess spdep
## fitted.lagmess spdep
## logLik.lagmess spdep
## fitted.SFResult spdep
## print.SFResult spdep
## fitted.ME_res spdep
## print.ME_res spdep
## print.lagImpact spdep
## plot.lagImpact spdep
## summary.lagImpact spdep
## HPDinterval.lagImpact spdep
## print.summary.lagImpact spdep
## print.sarlm spdep
## summary.sarlm spdep
## residuals.sarlm spdep
## deviance.sarlm spdep
## coef.sarlm spdep
## vcov.sarlm spdep
## fitted.sarlm spdep
## logLik.sarlm spdep
## anova.sarlm spdep
## predict.sarlm spdep
## print.summary.sarlm spdep
## print.sarlm.pred spdep
## as.data.frame.sarlm.pred spdep
## residuals.spautolm spdep
## deviance.spautolm spdep
## coef.spautolm spdep
## fitted.spautolm spdep
## print.spautolm spdep
## summary.spautolm spdep
## logLik.spautolm spdep
## print.summary.spautolm spdep
## print.WXImpact spdep
## summary.WXImpact spdep
## print.summary.WXImpact spdep
## predict.SLX spdep##
## Attaching package: 'spatialreg'## The following objects are masked from 'package:spdep':
##
## anova.sarlm, as.spam.listw, as_dgRMatrix_listw, as_dsCMatrix_I,
## as_dsCMatrix_IrW, as_dsTMatrix_listw, bptest.sarlm, can.be.simmed,
## cheb_setup, coef.gmsar, coef.sarlm, coef.spautolm, coef.stsls,
## create_WX, deviance.gmsar, deviance.sarlm, deviance.spautolm,
## deviance.stsls, do_ldet, eigen_pre_setup, eigen_setup, eigenw,
## errorsarlm, fitted.gmsar, fitted.ME_res, fitted.sarlm,
## fitted.SFResult, fitted.spautolm, get.ClusterOption,
## get.coresOption, get.mcOption, get.VerboseOption,
## get.ZeroPolicyOption, GMargminImage, GMerrorsar, griffith_sone,
## gstsls, Hausman.test, HPDinterval.lagImpact, impacts, intImpacts,
## Jacobian_W, jacobianSetup, l_max, lagmess, lagsarlm, lextrB,
## lextrS, lextrW, lmSLX, logLik.sarlm, logLik.spautolm, LR.sarlm,
## LR1.sarlm, LR1.spautolm, LU_prepermutate_setup, LU_setup,
## Matrix_J_setup, Matrix_setup, mcdet_setup, MCMCsamp, ME, mom_calc,
## mom_calc_int2, moments_setup, powerWeights, predict.sarlm,
## predict.SLX, print.gmsar, print.ME_res, print.sarlm,
## print.sarlm.pred, print.SFResult, print.spautolm, print.stsls,
## print.summary.gmsar, print.summary.sarlm, print.summary.spautolm,
## print.summary.stsls, residuals.gmsar, residuals.sarlm,
## residuals.spautolm, residuals.stsls, sacsarlm, SE_classic_setup,
## SE_interp_setup, SE_whichMin_setup, set.ClusterOption,
## set.coresOption, set.mcOption, set.VerboseOption,
## set.ZeroPolicyOption, similar.listw, spam_setup, spam_update_setup,
## SpatialFiltering, spautolm, spBreg_err, spBreg_lag, spBreg_sac,
## stsls, subgraph_eigenw, summary.gmsar, summary.sarlm,
## summary.spautolm, summary.stsls, trW, vcov.sarlm, Wald1.sarlmsetwd("~/university/Computacion/Clase 20 (15-07)")
XPABLO <- read_excel("XPABLO.xlsx")
View(XPABLO)dfCOS=data.frame(XPABLO)
X=as.matrix(data.frame(dfCOS[,4:15]))
###___matriz de coordenadaslibrary(clhs)## Warning: package 'clhs' was built under R version 4.0.5n = 0.80 *403
n = round(n,0)
df_michaell = data.frame(x=dfCOS$Long,
y=dfCOS$Lat,
Ca=dfCOS$Ca)res <- clhs(df_michaell, size = n,
iter = 100, progress = FALSE,
simple = TRUE)
df_michaell = dfCOS[res,]XY=as.matrix(df_michaell[,2:3])
Ca.d=as.matrix(dist(XY, diag=T, upper=T))
##matriz inversa de distancias
Ca.d.inv <-as.matrix(1/Ca.d)
##asignando cero a la diag
diag(Ca.d.inv) <- 0
##__la matriz de pesos basado en distancias
W=as.matrix(Ca.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(Ca~1,data=df_michaell,listw=Wve)
summary(modelo.arp)##
## Call: spautolm(formula = Ca ~ 1, data = df_michaell, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.467193908 -2.296903697 -0.181648921 1.895044703 11.113958891
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 7.18493196 3.95492681 1.8167 0.069262
##
## Lambda: 0.950990984 LR test value: 39.83807 p-value: 2.75914291e-10
## Numerical Hessian standard error of lambda: 0.0484384499
##
## Log likelihood: -861.204844172
## ML residual variance (sigma squared): 12.0971968, (sigma: 3.47810247)
## Number of observations: 322
## Number of parameters estimated: 3
## AIC: 1728.40969library(spatialreg)\[Y = \lambda W Y + \epsilon\]
res1 = modelo.arp$fit$residuals
Moran.I(res1, Ca.d.inv)## $observed
## [1] 0.0491621205315
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.00539458106464
##
## $p.value
## [1] 0mod2 = sacsarlm(Ca~1,data=df_michaell,listw=Wve)
summary(mod2)##
## Call:sacsarlm(formula = Ca ~ 1, data = df_michaell, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -7.534924402 -2.325997483 -0.168860323 1.876297531 10.903581223
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.36960423e-01 1.05839214e+06 0 1
##
## Rho: 0.917794045
## Asymptotic standard error: 635264.808
## z-value: 1.44474246e-06, p-value: 0.999998847
## Lambda: 0.917794097
## Asymptotic standard error: 635264.415
## z-value: 1.44474344e-06, p-value: 0.999998847
##
## LR test value: 61.817004, p-value: 3.77475828e-14
##
## Log likelihood: -850.215377195 for sac model
## ML residual variance (sigma squared): 11.1712676, (sigma: 3.34234463)
## Number of observations: 322
## Number of parameters estimated: 4
## AIC: 1708.43075, (AIC for lm: 1766.24776)\[Y = \lambda W Y + u \\ u= \rho Wu + \epsilon\]
res2 =mod2$residuals
Moran.I(res2, Ca.d.inv)## $observed
## [1] 0.0309369107046
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.00539428182023
##
## $p.value
## [1] 2.74309686077e-10mser1=errorsarlm(formula=Ca~Mg+K+Zn+cos,data=df_michaell,listw=Wve)
summary(mser1)##
## Call:errorsarlm(formula = Ca ~ Mg + K + Zn + cos, data = df_michaell,
## listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.757581733 -1.753384347 -0.209381016 1.174501683 11.727160508
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.6957706483 3.6784663825 0.73285 0.463649
## Mg 1.4453841667 0.1905596967 7.58494 3.3307e-14
## K 2.8285097628 1.4143996239 1.99980 0.045522
## Zn -0.4296151722 0.0724732308 -5.92792 3.0680e-09
## cos 3.3521438355 0.4163884678 8.05052 8.8818e-16
##
## Lambda: 0.959805969, LR test value: 49.1031912, p-value: 2.42850184e-12
## Asymptotic standard error: 0.0282511819
## z-value: 33.9740111, p-value: < 2.220446e-16
## Wald statistic: 1154.23343, p-value: < 2.220446e-16
##
## Log likelihood: -772.197763426 for error model
## ML residual variance (sigma squared): 6.95067662, (sigma: 2.63641359)
## Number of observations: 322
## Number of parameters estimated: 7
## AIC: 1558.39553, (AIC for lm: 1605.49872)res3 =mser1$residuals
Moran.I(res3, Ca.d.inv)## $observed
## [1] 0.0549041992838
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.00536764874266
##
## $p.value
## [1] 0\[Y = \lambda W Y +X \beta + \epsilon \]
mod4=lagsarlm(formula=Ca~Mg+Cu+cos,data=df_michaell,listw=Wve)
summary(mod4)##
## Call:lagsarlm(formula = Ca ~ Mg + Cu + cos, data = df_michaell, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.33192157 -1.75115582 -0.31417196 1.29701837 12.54151516
##
## Type: lag
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -5.5631585469 0.5915042412 -9.40510 < 2.22e-16
## Mg 1.6275863688 0.1733954564 9.38656 < 2.22e-16
## Cu -0.4059121610 0.0941979995 -4.30914 1.6389e-05
## cos 3.7194793488 0.4076707045 9.12373 < 2.22e-16
##
## Rho: 0.953796921, LR test value: 45.6194363, p-value: 1.43609569e-11
## Asymptotic standard error: 0.032361338
## z-value: 29.4733462, p-value: < 2.220446e-16
## Wald statistic: 868.678139, p-value: < 2.220446e-16
##
## Log likelihood: -783.806669139 for lag model
## ML residual variance (sigma squared): 7.47716217, (sigma: 2.73444001)
## Number of observations: 322
## Number of parameters estimated: 6
## AIC: 1579.61334, (AIC for lm: 1623.23277)
## LM test for residual autocorrelation
## test value: 118.523538, p-value: < 2.220446e-16res4 = mod4$residuals
Moran.I(res4, Ca.d.inv)## $observed
## [1] 0.0546673115997
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.00536867488249
##
## $p.value
## [1] 0mod5=sacsarlm(formula=Ca~Mg+Zn+cos,data=df_michaell,listw=Wve)
summary(mod5)##
## Call:sacsarlm(formula = Ca ~ Mg + Zn + cos, data = df_michaell, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.254434758 -1.582975670 -0.275360939 1.197206645 11.531896956
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -4.7588703493 6.0236526574 -0.79003 0.42951
## Mg 1.6076304792 0.1679678849 9.57106 < 2.22e-16
## Zn -0.4001732024 0.0703379613 -5.68929 1.2757e-08
## cos 3.3972727906 0.3852743011 8.81780 < 2.22e-16
##
## Rho: 0.911741174
## Asymptotic standard error: 0.458236792
## z-value: 1.98967256, p-value: 0.0466270169
## Lambda: 0.933192287
## Asymptotic standard error: 0.3515788
## z-value: 2.65429055, p-value: 0.00794753618
##
## LR test value: 71.9434819, p-value: 2.22044605e-16
##
## Log likelihood: -762.855564035 for sac model
## ML residual variance (sigma squared): 6.48717723, (sigma: 2.54699376)
## Number of observations: 322
## Number of parameters estimated: 7
## AIC: 1539.71113, (AIC for lm: 1607.65461)res5 = mod5$residuals
Moran.I(res5, Ca.d.inv)## $observed
## [1] 0.0338208974686
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.0053666386981
##
## $p.value
## [1] 5.8790750046e-12\[Y = \lambda W Y +X \beta + u\\u =\rho W u +\epsilon \]
mod6=sacsarlm(formula=Ca~Mg+Zn+cos,data=df_michaell,listw=Wve,type="mixed")
summary(mod6)##
## Call:sacsarlm(formula = Ca ~ Mg + Zn + cos, data = df_michaell, listw = Wve,
## type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -5.660583693 -1.635916678 -0.243605414 1.166362652 11.750428060
##
## Type: sacmixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -11.79157323430 7.85168954125 -1.50179 0.1331519
## Mg 1.82181085937 0.16820320790 10.83101 < 2.22e-16
## Zn -0.35215846279 0.07095691500 -4.96299 6.9416e-07
## cos 3.04016480954 0.39553268810 7.68625 1.5099e-14
## lag.Mg -3.76311902078 1.42595379122 -2.63902 0.0083146
## lag.Zn 0.00238222341 1.25730432886 0.00189 0.9984882
## lag.cos 11.17278956417 7.67291162147 1.45613 0.1453555
##
## Rho: 0.90484024
## Asymptotic standard error: 0.938938615
## z-value: 0.963684127, p-value: 0.335204315
## Lambda: 0.91862483
## Asymptotic standard error: 0.807730263
## z-value: 1.13729158, p-value: 0.255416412
##
## LR test value: 88.5608752, p-value: < 2.220446e-16
##
## Log likelihood: -754.546867391 for sacmixed model
## ML residual variance (sigma squared): 6.17218327, (sigma: 2.48438791)
## Number of observations: 322
## Number of parameters estimated: 10
## AIC: 1529.09373, (AIC for lm: 1607.65461)Los coeficientes Lambda y Rho no son significativos
mod7=lagsarlm(formula=Ca~Mg+K+Cu+cos,data=df_michaell,listw=Wve,type="mixed")
summary(mod7)##
## Call:lagsarlm(formula = Ca ~ Mg + K + Cu + cos, data = df_michaell,
## listw = Wve, type = "mixed")
##
## Residuals:
## Min 1Q Median 3Q Max
## -6.442045210 -1.561357009 -0.282300217 1.189838713 12.719521817
##
## Type: mixed
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -18.786824947 5.609297784 -3.34923 0.00081037
## Mg 1.724782910 0.213340340 8.08465 6.6613e-16
## K 2.444688412 1.413895398 1.72904 0.08380109
## Cu -0.317196374 0.100585139 -3.15351 0.00161319
## cos 2.784891827 0.454683078 6.12491 9.0736e-10
## lag.Mg -5.331507830 1.686250465 -3.16175 0.00156822
## lag.K 2.989369835 15.420149737 0.19386 0.84628451
## lag.Cu 1.354248310 0.771523091 1.75529 0.07920937
## lag.cos 12.939972613 4.336070708 2.98426 0.00284263
##
## Rho: 0.95171472, LR test value: 41.213329, p-value: 1.36488598e-10
## Asymptotic standard error: 0.0338935678
## z-value: 28.0795084, p-value: < 2.220446e-16
## Wald statistic: 788.458791, p-value: < 2.220446e-16
##
## Log likelihood: -770.780373662 for mixed model
## ML residual variance (sigma squared): 6.89801355, (sigma: 2.62640696)
## Number of observations: 322
## Number of parameters estimated: 11
## AIC: 1563.56075, (AIC for lm: 1602.77408)
## LM test for residual autocorrelation
## test value: 84.5851018, p-value: < 2.220446e-16res7 = mod7$residuals
Moran.I(res7, Ca.d.inv)## $observed
## [1] 0.045909720117
##
## $expected
## [1] -0.00311526479751
##
## $sd
## [1] 0.00536864001007
##
## $p.value
## [1] 0El modelo no explica la dependencia espacial
Ninguno de los modelos explican la dependencia del Ca, es decir que las distancias y las otras variables no son suficientes para explicarla. Posiblemente sea porque las variables medidas no son las necesarias para explicar el calcio.