Modelado de CE75
Daniel Antonio Rodriguez Parra
library(spdep)## Loading required package: sp## 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')`## Loading required package: sf## Linking to GEOS 3.5.1, GDAL 2.2.2, PROJ 4.9.2library(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## Loading required package: ggplot2##
## 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(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_dgRMatrix_listw, as_dsCMatrix_I, as_dsCMatrix_IrW,
## as_dsTMatrix_listw, as.spam.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.sarlmlibrary(readxl)Análisis exploratorio de los datos
library(readxl)
BD_MODELADO <- read_excel("/cloud/project/COMPUTACION Y ESTADISTICA/clase 15/BD_MODELADO.xlsx")
#View(BD_MODELADO)
data1=data.frame(BD_MODELADO )
head(data1,10)## Avg_X_MCB Avg_Y_MCE Avg_CEa_07 Avg_CEa_15 NDVI DEM SLOPE
## 1 843449.6 955962.0 7.237480 18.02656 0.863030 199.0000 6.385167
## 2 843454.8 955962.4 6.787250 18.02737 0.866502 197.1667 1.981082
## 3 843493.6 955951.3 6.848250 18.70444 0.874883 197.0000 0.577682
## 4 843439.6 955977.6 7.135162 18.34237 0.845838 197.0000 1.175075
## 5 843468.7 955972.8 6.826763 17.92409 0.797179 197.0000 0.210996
## 6 843500.6 955975.3 6.699966 18.39441 0.758272 197.6667 4.357386
## 7 843525.7 955982.6 6.180742 17.84332 0.763436 199.7500 6.628445
## 8 843413.4 956014.0 8.539024 18.75812 0.823320 197.1667 1.462050
## 9 843433.6 956013.4 8.869958 18.85396 0.759923 197.3333 1.663344
## 10 843468.8 956010.1 7.231308 18.34269 0.757382 197.6667 3.541936
## Avg_z
## 1 193.0512
## 2 193.2986
## 3 193.5659
## 4 194.4116
## 5 193.9931
## 6 195.3814
## 7 196.6780
## 8 194.9936
## 9 196.1356
## 10 197.8522ce75=data1[,3] #variable respuesta# X = matriz de variables explicativas
# CE150, NDVI, slope, y z(altura)
X=as.matrix(data1[,c(4:5,7:8)])
head (X)## Avg_CEa_15 NDVI SLOPE Avg_z
## [1,] 18.02656 0.863030 6.385167 193.0512
## [2,] 18.02737 0.866502 1.981082 193.2986
## [3,] 18.70444 0.874883 0.577682 193.5659
## [4,] 18.34237 0.845838 1.175075 194.4116
## [5,] 17.92409 0.797179 0.210996 193.9931
## [6,] 18.39441 0.758272 4.357386 195.3814# Coordenadas
options(digits = 12)
xydata=as.matrix(data1[,1:2]) #coordenadas geograficas
distancias = as.matrix(dist(xydata))
#distancias[1:5,1:5]
distancias.inv= as.matrix(1/distancias)
diag(distancias.inv)= 0
min(distancias[distancias!=0])## [1] 5.17585430653w= as.matrix(distancias.inv) #matrix de pesos
sumas= apply(w, 1, sum)
we= w/sumas # estandarizado
apply(we,1,sum) # todos los datos dan 1 y significa que quedó bien estandarizado## 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 1max(distancias)## [1] 853.009893533dim(distancias)## [1] 313 313Indices de Morán.
moran.ce75= (Moran.I(data1[, 3], we));moran.ce75## $observed
## [1] 0.268746828627
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00466590593614
##
## $p.value
## [1] 0moran.ce150= (Moran.I(data1[, 4], we));moran.ce150## $observed
## [1] 0.160951018381
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00465455026174
##
## $p.value
## [1] 0moran.NDVI= (Moran.I(data1[, 5], we));moran.NDVI## $observed
## [1] 0.0975040302526
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00464497872749
##
## $p.value
## [1] 0moran.DEM= (Moran.I(data1[, 6], we));moran.DEM## $observed
## [1] 0.309670845207
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00467238374652
##
## $p.value
## [1] 0moran.slope= (Moran.I(data1[, 7], we));moran.slope## $observed
## [1] 0.0699332392626
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00465430682803
##
## $p.value
## [1] 0moran.z= (Moran.I(data1[, 8], we));moran.z## $observed
## [1] 0.350503097268
##
## $expected
## [1] -0.00320512820513
##
## $sd
## [1] 0.00466793519805
##
## $p.value
## [1] 0matriz de correlacion (mat)
matp=(rcorr(as.matrix(data1[, 3:8]), type= "pearson"))
mats=(rcorr(as.matrix(data1[, 3:8]), type= "spearman"))
mcorp=matp$r
mcors=mats$r
library (corrplot)## corrplot 0.84 loadedpar(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")
library (psych)##
## Attaching package: 'psych'## The following object is masked from 'package:Hmisc':
##
## describe## The following objects are masked from 'package:ggplot2':
##
## %+%, alphaoptions(digits = 5)
describe(data1[,3:8])## vars n mean sd median trimmed mad min max range skew
## Avg_CEa_07 1 313 9.77 1.53 9.65 9.76 1.40 6.17 13.60 7.44 0.10
## Avg_CEa_15 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
## Avg_z 6 313 202.47 3.66 202.49 202.57 4.34 193.05 212.35 19.30 -0.14
## kurtosis se
## Avg_CEa_07 -0.28 0.09
## Avg_CEa_15 1.22 0.04
## NDVI 2.48 0.00
## DEM -1.13 0.26
## SLOPE 1.25 0.12
## Avg_z -0.54 0.21boxplot(data1$Avg_CEa_07)
hist(data1$Avg_CEa_07)
shapir =c()
for (j in 3:8) {
shapir[j] = sf.test(data1[,j])$p.value
}
shapir## [1] NA NA 2.0125e-01 1.4139e-04 5.1932e-12 1.6549e-06 4.5506e-08
## [8] 5.7707e-04Modelos Espaciales
contnb = dnearneigh(x = coordinates(xydata),
d1 = 0, d2 = 380000, longlat = F)
contnb## Neighbour list object:
## Number of regions: 313
## Number of nonzero links: 97656
## Percentage nonzero weights: 99.681
## Average number of links: 312dlist <- nbdists(contnb, xydata)
dlist <- lapply(dlist, function(x) 1/x)
Wve=nb2listw(contnb,glist=dlist,style = "W") # hasta aqui se estandarizó la matrix de pesos.Primer modelo
Modelo Autorregresivo Puro (map.)
\[Y = \lambda W Y + u\]
library(spatialreg)
# map. = Modelo Autorregresivo Puro (map.) libro ARBIA 2014
map.=spautolm(ce75~1,data=data1,listw=Wve) #### modelo
summary(map.)##
## Call: spautolm(formula = ce75 ~ 1, data = data1, 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.011876
##
## 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.65ce75E1=map.$fit$fitted.values ;ce75E1 # valores 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.1322DF_map.=data.frame(ce75,ce75E1) #dataframe
colnames(DF_map.) <- c("ce75_o","ce75_e");DF_map.## ce75_o ce75_e
## 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
## 7 6.1807 8.9267
## 8 8.5390 9.0756
## 9 8.8700 9.0054
## 10 7.2313 8.9769
## 11 7.3722 8.8846
## 12 7.5568 8.8911
## 13 6.6135 9.0060
## 14 8.7076 9.1184
## 15 8.6195 9.0633
## 16 9.4434 9.0784
## 17 7.9488 9.0171
## 18 7.6172 8.9640
## 19 6.9522 9.0108
## 20 8.9010 9.1408
## 21 8.3623 9.1678
## 22 9.2462 9.1535
## 23 9.5656 9.1619
## 24 9.5142 9.1763
## 25 7.7654 9.0825
## 26 7.7404 9.0567
## 27 8.0054 9.0478
## 28 6.5610 9.0672
## 29 6.2831 9.0447
## 30 8.3191 9.1880
## 31 9.0395 9.1637
## 32 8.9674 9.2192
## 33 10.1804 9.2347
## 34 10.3069 9.2775
## 35 10.3879 9.2092
## 36 8.0793 9.1145
## 37 7.4166 9.1417
## 38 7.7941 9.0571
## 39 6.3589 9.0632
## 40 7.2514 9.0851
## 41 9.2399 9.3782
## 42 8.8082 9.4087
## 43 9.6902 9.3506
## 44 10.1558 9.3033
## 45 8.8046 9.1994
## 46 8.4342 9.2065
## 47 9.1565 9.2210
## 48 9.2740 9.2981
## 49 10.2077 9.3359
## 50 10.9029 9.3437
## 51 7.7615 9.1671
## 52 8.0063 9.2176
## 53 7.5334 9.1321
## 54 6.5999 9.1396
## 55 6.1654 9.2313
## 56 9.8012 9.3898
## 57 9.8834 9.3622
## 58 9.3561 9.3639
## 59 9.4172 9.3071
## 60 8.7556 9.2695
## 61 8.6986 9.2370
## 62 8.7034 9.2178
## 63 8.2734 9.2240
## 64 8.3495 9.2513
## 65 9.4067 9.3058
## 66 9.4376 9.3910
## 67 10.1062 9.4173
## 68 10.6228 9.4125
## 69 9.7106 9.3594
## 70 8.7927 9.2655
## 71 8.1573 9.3239
## 72 7.5864 9.2164
## 73 7.5712 9.2656
## 74 6.8278 9.3774
## 75 10.4136 9.3644
## 76 9.1669 9.3646
## 77 8.5040 9.3213
## 78 8.7885 9.2619
## 79 8.6623 9.2327
## 80 8.7192 9.2311
## 81 8.3419 9.2489
## 82 8.2930 9.2862
## 83 8.8941 9.3333
## 84 9.6408 9.4118
## 85 9.4923 9.4733
## 86 9.7830 9.5075
## 87 11.1631 9.4876
## 88 9.3082 9.4922
## 89 8.1564 9.4381
## 90 8.2873 9.3610
## 91 8.9510 9.4106
## 92 7.0400 9.5162
## 93 9.4233 9.3367
## 94 9.0178 9.3264
## 95 8.6484 9.2866
## 96 8.5786 9.2545
## 97 8.4992 9.2576
## 98 8.4041 9.2849
## 99 8.7421 9.3203
## 100 9.3693 9.3817
## 101 9.5602 9.4590
## 102 9.7545 9.5114
## 103 9.5505 9.5719
## 104 9.4888 9.6277
## 105 11.0769 9.5854
## 106 9.9988 9.6143
## 107 9.7593 9.6095
## 108 8.1344 9.5499
## 109 8.2830 9.5690
## 110 9.0620 9.5446
## 111 8.1718 9.6038
## 112 7.5538 9.6738
## 113 9.1050 9.3280
## 114 9.5407 9.3182
## 115 8.6032 9.2953
## 116 8.7890 9.2850
## 117 8.8159 9.3124
## 118 8.6719 9.3656
## 119 9.6392 9.4151
## 120 9.7246 9.4837
## 121 9.8616 9.5574
## 122 10.1280 9.6224
## 123 10.3889 9.6810
## 124 9.9523 9.7327
## 125 11.1444 9.7106
## 126 10.9161 9.7839
## 127 9.7106 9.7116
## 128 8.0851 9.6820
## 129 9.4378 9.6372
## 130 8.4613 9.7157
## 131 7.9542 9.8259
## 132 9.1584 9.3429
## 133 9.3207 9.3451
## 134 8.9145 9.3261
## 135 8.7622 9.3542
## 136 9.0318 9.3945
## 137 8.9581 9.4620
## 138 9.3907 9.5138
## 139 9.7992 9.5859
## 140 10.4915 9.6685
## 141 10.6666 9.7419
## 142 10.2208 9.8130
## 143 9.5765 9.8626
## 144 11.3255 9.8637
## 145 11.4924 9.9185
## 146 7.9689 9.8194
## 147 8.1211 9.8963
## 148 9.4615 9.7944
## 149 9.1197 9.9096
## 150 9.4900 10.0180
## 151 9.0331 9.3705
## 152 9.0388 9.3963
## 153 8.9306 9.3978
## 154 8.8756 9.4401
## 155 9.4281 9.4755
## 156 8.7883 9.5495
## 157 9.1774 9.5940
## 158 9.7940 9.6980
## 159 10.4948 9.7920
## 160 10.6376 9.8731
## 161 10.6512 9.9467
## 162 9.8736 10.0081
## 163 11.1766 10.0306
## 164 11.4201 10.0494
## 165 10.3732 9.9619
## 166 9.6119 10.0184
## 167 8.9063 10.0643
## 168 10.2531 10.0020
## 169 10.4814 10.1681
## 170 13.0589 10.0806
## 171 9.3501 9.4504
## 172 9.2424 9.4740
## 173 9.0993 9.5055
## 174 9.0631 9.5540
## 175 9.2049 9.6177
## 176 8.9118 9.7048
## 177 10.2908 9.8234
## 178 10.6832 9.9248
## 179 11.1029 10.0229
## 180 11.5557 10.0915
## 181 10.2627 10.1567
## 182 10.8338 10.2052
## 183 11.7262 10.1909
## 184 10.0572 10.1217
## 185 10.2259 10.2187
## 186 9.1905 10.2377
## 187 10.4669 10.2423
## 188 11.0885 10.2441
## 189 10.5376 10.2290
## 190 9.5880 9.4994
## 191 8.8905 9.5343
## 192 9.1787 9.5627
## 193 8.9396 9.6535
## 194 9.1215 9.7414
## 195 9.1614 9.8607
## 196 9.6241 10.0037
## 197 11.1248 10.0919
## 198 11.8900 10.1894
## 199 11.2165 10.2583
## 200 11.0746 10.3270
## 201 11.5755 10.3895
## 202 11.6253 10.3665
## 203 12.5340 10.4452
## 204 12.6649 10.3516
## 205 10.7292 10.3532
## 206 9.7107 10.4099
## 207 13.2276 10.1819
## 208 9.1261 10.1901
## 209 6.8942 10.1525
## 210 8.5036 9.5697
## 211 8.4513 9.6149
## 212 8.9304 9.6843
## 213 9.9790 9.7804
## 214 10.0391 9.8942
## 215 10.8348 10.0048
## 216 11.1278 10.1651
## 217 11.7919 10.2510
## 218 12.8631 10.3339
## 219 12.3663 10.4222
## 220 12.0968 10.4946
## 221 11.7383 10.5156
## 222 12.7405 10.5303
## 223 13.6011 10.5659
## 224 11.4573 10.4488
## 225 11.1662 10.3972
## 226 11.2972 10.3374
## 227 11.0391 10.0945
## 228 7.9633 10.1034
## 229 9.0486 9.6453
## 230 8.8815 9.7259
## 231 10.8351 9.8020
## 232 10.5065 9.9227
## 233 9.8400 10.0388
## 234 10.9167 10.1558
## 235 11.0109 10.3017
## 236 13.4243 10.3613
## 237 12.6752 10.4819
## 238 12.4478 10.5475
## 239 12.4449 10.5711
## 240 11.6324 10.6251
## 241 12.9714 10.6139
## 242 13.2626 10.5648
## 243 10.6916 10.4490
## 244 10.1106 10.3797
## 245 10.2149 10.2696
## 246 9.3301 10.0658
## 247 8.8087 9.7964
## 248 9.2843 9.8558
## 249 9.9659 9.9555
## 250 10.3226 10.0398
## 251 9.7532 10.1606
## 252 10.1070 10.3221
## 253 12.0642 10.3718
## 254 13.3181 10.4630
## 255 11.8400 10.5717
## 256 11.6916 10.6140
## 257 12.4710 10.6161
## 258 11.8600 10.6422
## 259 11.9740 10.6120
## 260 12.4945 10.4919
## 261 9.8456 10.3797
## 262 10.7921 10.2952
## 263 9.8996 9.9715
## 264 10.4534 10.0433
## 265 10.9725 10.1069
## 266 9.7149 10.2565
## 267 11.0642 10.3711
## 268 12.0266 10.4114
## 269 12.1360 10.5254
## 270 12.2982 10.5896
## 271 12.1225 10.6187
## 272 12.7160 10.5942
## 273 11.4330 10.6018
## 274 11.0975 10.5574
## 275 12.0781 10.4023
## 276 9.7875 10.0749
## 277 9.9779 10.1916
## 278 10.4012 10.2864
## 279 10.1850 10.3859
## 280 11.6950 10.4369
## 281 12.2476 10.5489
## 282 12.0259 10.5654
## 283 11.2564 10.5766
## 284 11.5159 10.5609
## 285 11.6114 10.5157
## 286 9.5776 10.1271
## 287 9.8945 10.1347
## 288 10.3465 10.2299
## 289 10.3866 10.3079
## 290 9.9725 10.3587
## 291 11.3632 10.4448
## 292 11.5113 10.5030
## 293 11.3215 10.4969
## 294 11.4038 10.5031
## 295 12.2376 10.4499
## 296 10.0131 10.1627
## 297 9.6453 10.1992
## 298 10.9614 10.2353
## 299 9.8268 10.3047
## 300 9.6928 10.3435
## 301 11.4779 10.4126
## 302 11.5633 10.3965
## 303 10.7877 10.4263
## 304 10.0407 10.2179
## 305 11.5077 10.1799
## 306 9.8912 10.2282
## 307 9.9730 10.2952
## 308 10.5907 10.3118
## 309 9.5653 10.1904
## 310 9.0025 10.1663
## 311 9.7625 10.1790
## 312 9.2256 10.2662
## 313 9.3946 10.1322plot(DF_map.$ce75_o, DF_map.$ce75_e)
plot (xydata[,1], xydata[,2], cex=0.09*abs(DF_map.$ce75_o), pch=19,col= "brown",
main= "Ce75 OBSERVADOS VS ESPERADOS (map.)" )
points(xydata[,1], xydata[,2], cex=1.55, pch=1, col= floor(abs(DF_map.$ce75_e)+3))
cor(DF_map.$ce75_o, DF_map.$ce75_e) #correlacion## [1] 0.79772La correlacion entre ce75 observado y ce75 esperado es de 0,79 lo cual se considera que hay buena correlación.
Análisis de los residuales del primer modelo
res_map. = map.$fit$residuals # Residuales "u" del modelo
plot(xydata[,1], xydata[,2], cex=0.4*abs(res_map.), pch=19, col="green", main="RESIDUALES DEL MODELO") # puntos mas grandes mayor desajuste
library(spdep)
I.Moran_res_map. = moran.mc(res_map., Wve, nsim = 2000);I.Moran_res_map.##
## Monte-Carlo simulation of Moran I
##
## data: res_map.
## weights: Wve
## number of simulations + 1: 2001
##
## statistic = 0.167, observed rank = 2001, p-value = 5e-04
## alternative hypothesis: greatershapiro.test(res_map.)##
## Shapiro-Wilk normality test
##
## data: res_map.
## W = 0.994, p-value = 0.2skewness.norm.test(res_map.)##
## Skewness test for normality
##
## data: res_map.
## T = 0.0399, p-value = 0.77Con el Moran Montecarlo se obtiene que los residuales presentan dependencia.
Con la prueba shapiro se evidencia que los residuales presentan normalidad y por lo tanto el modelo se puede hacer.
Con la prueba de skewness.norm.test se comprueba la simetria de los datos.
Despeje de la ecuacion para extrapolar un punto de CE75
\[Y = \lambda W Y + \epsilon\] \[Y - \lambda W Y = \epsilon\\ ( I - \lambda W) Y = \epsilon\\ Y = ( I - \lambda W)^{-1}\epsilon\]
\[E(Y) = E((I-\lambda W)^{-1}\epsilon)\\ E(Y) = E((I-\lambda W)^{-1})E(\epsilon)\\ \require{cancel} E(Y) = E((I-\lambda W)^{-1})\cancelto{0}{E(\epsilon)}\\ \] ### Curvas de nivel para Ce75 estimada
library(ggplot2)
library(akima) # Interpolador Lineal.
interpol = interp(x = xydata[,1],
y = xydata[,2],
z = DF_map.$ce75_e,
nx = 500, ny = 500,
linear = F,
)
image(interpol)
contour(interpol, add = T)
### contorno convexo
cc = chull(xydata)
cc = c(cc, cc[1])
plot(xydata, pch = 16, col="lightblue", main= "CONTORNO CONVEXO")
lines(xydata[cc,], type = 'l')
points(843628,956183,col="orange ",pch=19,cex=0.7)
Segundo modelo
Modelo Autorregresivo - S ARAR (doble autogresivo)
\[Y = \lambda W Y + u\\ u=\rho W u+ \epsilon\]
library(spatialreg)
map._sarar=sacsarlm(ce75~1,data=data1,listw=Wve) #### modelo## Warning in sacsarlm(ce75 ~ 1, data = data1, listw = Wve): inversion of asymptotic covariance matrix failed for tol.solve = 2.22044604925031e-16
## reciprocal condition number = 2.88892e-20 - using numerical Hessian.summary(map._sarar)##
## Call:sacsarlm(formula = ce75 ~ 1, data = data1, 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.2488 -0.5315 0.5951
##
## Rho: 0.97863
## Approximate (numerical Hessian) standard error: 0.021304
## z-value: 45.936, p-value: < 2.22e-16
## Lambda: 0.97863
## Approximate (numerical Hessian) standard error: 0.02128
## z-value: 45.988, 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)ce75E2=map._sarar$fitted.values #valores estimados
res2=map._sarar$residuals
moran_msarar = moran.mc(res2, Wve, nsim = 2000)
moran_msarar##
## Monte-Carlo simulation of Moran I
##
## data: res2
## weights: Wve
## number of simulations + 1: 2001
##
## statistic = 0.105, observed rank = 2001, p-value = 5e-04
## alternative hypothesis: greatershapiro.test(res2)##
## Shapiro-Wilk normality test
##
## data: res2
## W = 0.994, p-value = 0.27Este segundo modelo calcula rho y lambda y arroja un mejor AIC que nuestro primer modelo. Por otro lado los residuales son normales, pero no son independientes, por lo tanto a pesar de que el AIC del segundo modelo es mejor, ambos modelos no son buenos.
DF_sarar.=data.frame(ce75,ce75E2) #dataframe
colnames(DF_sarar.) <- c("ce75_o","ce75_e")
plot(DF_sarar.$ce75_o, DF_sarar.$ce75_e)
plot (xydata[,1], xydata[,2], cex=0.09*abs(DF_sarar.$ce75_o), pch=19,col= "brown",
main= "Ce75 OBSERVADOS VS ESPERADOS (sarar)" )
points(xydata[,1], xydata[,2], cex=1.55, pch=1, col= floor(abs(DF_sarar.$ce75_e)+3))
Modelos que permitan involucrar variables explicativas (X).
TERCER MODELO
Modelo espacial del error (meer)
library(spatialreg)
meer1= spatialreg::errorsarlm(formula = ce75~ Avg_CEa_15+NDVI+DEM+SLOPE+Avg_z, data=data1,listw=Wve)
summary(meer1)##
## Call:spatialreg::errorsarlm(formula = ce75 ~ Avg_CEa_15 + NDVI + DEM +
## SLOPE + Avg_z, data = data1, 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
## Avg_CEa_15 0.859898 0.083054 10.3535 < 2.2e-16
## NDVI -2.395368 1.907913 -1.2555 0.209301
## DEM 0.036792 0.020974 1.7542 0.079402
## SLOPE -0.073067 0.024760 -2.9510 0.003168
## Avg_z 0.257034 0.028465 9.0299 < 2.2e-16
##
## 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)meer2= spatialreg::errorsarlm(formula = ce75~ Avg_CEa_15+DEM+SLOPE+Avg_z, data=data1,listw=Wve)
summary(meer2)##
## Call:spatialreg::errorsarlm(formula = ce75 ~ Avg_CEa_15 + DEM + SLOPE +
## Avg_z, data = data1, listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.068942 -0.573110 -0.041672 0.535538 2.620533
##
## Type: error
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -66.334322 5.621356 -11.8004 < 2.2e-16
## Avg_CEa_15 0.871288 0.082765 10.5273 < 2.2e-16
## DEM 0.039380 0.020925 1.8819 0.059845
## SLOPE -0.074849 0.024782 -3.0203 0.002525
## Avg_z 0.251732 0.028220 8.9203 < 2.2e-16
##
## Lambda: 0.98246, LR test value: 98.998, p-value: < 2.22e-16
## Asymptotic standard error: 0.012369
## z-value: 79.427, p-value: < 2.22e-16
## Wald statistic: 6308.6, p-value: < 2.22e-16
##
## Log likelihood: -406.89 for error model
## ML residual variance (sigma squared): 0.76989, (sigma: 0.87744)
## Number of observations: 313
## Number of parameters estimated: 7
## AIC: 827.77, (AIC for lm: 924.77)meer3= spatialreg::errorsarlm(formula = ce75~ Avg_CEa_15+SLOPE+Avg_z, data=data1,listw=Wve)
summary(meer3)##
## Call:spatialreg::errorsarlm(formula = ce75 ~ Avg_CEa_15 + SLOPE +
## Avg_z, data = data1, 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
## Avg_CEa_15 0.874324 0.083217 10.507 < 2.2e-16
## SLOPE -0.079881 0.024777 -3.224 0.001264
## Avg_z 0.286926 0.021256 13.498 < 2.2e-16
##
## 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)ce75E3=meer2$fitted.values ;ce75E3 # valores estimados## 1 2 3 4 5 6 7 8 9 10
## 6.3532 6.7756 7.3939 7.3357 6.8736 7.2830 7.0122 7.8515 8.1919 7.9993
## 11 12 13 14 15 16 17 18 19 20
## 7.6302 8.1429 7.7701 7.7603 8.0725 8.5316 8.8586 8.8486 8.7354 8.5649
## 21 22 23 24 25 26 27 28 29 30
## 8.5080 9.0925 9.0784 9.4300 9.0792 9.3877 8.9749 8.3689 8.3520 8.5051
## 31 32 33 34 35 36 37 38 39 40
## 9.3293 9.0010 9.6560 10.3564 9.7567 9.6052 8.6408 8.6504 7.8907 8.1190
## 41 42 43 44 45 46 47 48 49 50
## 9.1649 8.9475 8.7870 8.7676 8.5070 8.7471 9.1622 8.8809 10.0379 10.1420
## 51 52 53 54 55 56 57 58 59 60
## 8.2969 8.6488 8.5750 8.3688 8.1238 9.4405 8.8773 8.5014 7.9879 7.9175
## 61 62 63 64 65 66 67 68 69 70
## 8.0564 8.2047 8.3752 8.8247 9.4331 8.8741 9.3841 10.0038 9.7229 8.3482
## 71 72 73 74 75 76 77 78 79 80
## 8.5531 8.2909 8.3447 8.0076 9.6205 9.3041 8.4874 8.1497 8.4294 8.8485
## 81 82 83 84 85 86 87 88 89 90
## 8.9191 8.9489 9.1848 9.2346 8.8159 9.4021 10.9351 10.2658 8.6519 8.4828
## 91 92 93 94 95 96 97 98 99 100
## 8.7076 7.4685 9.6090 9.0594 8.8344 8.3204 8.3815 9.2387 9.7174 9.7465
## 101 102 103 104 105 106 107 108 109 110
## 9.2965 9.5506 9.2984 9.6589 11.3372 11.0704 10.8214 8.8952 8.9115 9.6076
## 111 112 113 114 115 116 117 118 119 120
## 8.5637 8.2955 9.2810 9.3237 9.1736 8.8181 9.3431 9.6646 9.5618 9.7107
## 121 122 123 124 125 126 127 128 129 130
## 9.5056 9.7497 10.1168 10.1546 11.1583 10.8729 10.2898 8.8806 9.3207 9.0695
## 131 132 133 134 135 136 137 138 139 140
## 8.4394 9.6502 9.2622 9.5210 9.5152 9.4385 9.4358 9.5719 9.9303 9.9843
## 141 142 143 144 145 146 147 148 149 150
## 10.6455 10.5166 10.3898 10.9336 10.9366 8.6919 8.9915 9.4702 9.0948 9.1036
## 151 152 153 154 155 156 157 158 159 160
## 9.7642 9.7229 9.4785 9.7882 9.5009 9.3508 9.6227 9.9667 10.7250 10.5330
## 161 162 163 164 165 166 167 168 169 170
## 10.7792 10.1142 10.4690 10.7656 9.7834 9.8710 9.1241 9.2484 9.7441 11.1493
## 171 172 173 174 175 176 177 178 179 180
## 10.1785 9.7858 10.1951 9.5586 9.4422 9.4849 10.9345 10.7785 10.2968 10.5936
## 181 182 183 184 185 186 187 188 189 190
## 9.8921 10.2048 10.8858 10.4220 10.0822 9.3887 9.2444 9.6301 10.3046 9.9582
## 191 192 193 194 195 196 197 198 199 200
## 9.9911 9.9405 9.0755 9.6818 9.9213 11.2390 10.6223 10.0835 10.4443 10.1427
## 201 202 203 204 205 206 207 208 209 210
## 9.9971 11.3874 11.0027 10.0444 9.0156 9.1751 11.8651 9.5032 8.1381 9.8152
## 211 212 213 214 215 216 217 218 219 220
## 10.0437 10.2247 10.2100 10.1339 10.6796 11.3393 11.4148 10.5773 10.5025 10.5143
## 221 222 223 224 225 226 227 228 229 230
## 9.9023 10.8281 11.2410 9.7372 9.9216 10.2024 11.4924 9.5610 10.1095 9.9093
## 231 232 233 234 235 236 237 238 239 240
## 10.1908 10.1147 9.5827 10.6926 10.6708 11.6418 10.9434 10.8073 11.1414 11.8199
## 241 242 243 244 245 246 247 248 249 250
## 11.4020 11.5493 9.9570 9.9562 9.7423 11.2172 9.5892 9.6998 9.9913 10.1681
## 251 252 253 254 255 256 257 258 259 260
## 9.6035 10.2710 10.8523 11.4702 10.9743 10.6657 10.7546 12.3622 11.9538 11.5415
## 261 262 263 264 265 266 267 268 269 270
## 9.8961 10.2375 9.8532 10.2680 11.4023 9.8920 10.2530 10.7561 11.1563 11.2824
## 271 272 273 274 275 276 277 278 279 280
## 10.7296 10.6027 11.9470 12.4094 11.7311 10.4373 10.8366 10.1619 10.5037 10.6382
## 281 282 283 284 285 286 287 288 289 290
## 11.4428 11.2023 10.8375 10.6842 11.9037 10.8233 10.7984 10.8727 10.4037 10.0450
## 291 292 293 294 295 296 297 298 299 300
## 10.6968 11.0776 11.2510 11.0079 11.7251 10.7181 10.8037 11.1407 10.4984 9.7902
## 301 302 303 304 305 306 307 308 309 310
## 10.6029 11.0347 10.6894 10.8660 11.0810 11.1877 9.5084 10.8270 11.2705 10.9890
## 311 312 313
## 11.7630 10.5324 11.1356DF_meer=data.frame(ce75,ce75E3) #dataframe
colnames(DF_meer) <- c("ce75_o","ce75_e")
plot(DF_meer$ce75_o, DF_meer$ce75_e)
plot (xydata[,1], xydata[,2], cex=0.09*abs(DF_meer$ce75_o), pch=19,col= "orange",
main= "Ce75 OBSERVADOS VS ESPERADOS (meer2)" )
points(xydata[,1], xydata[,2], cex=1.55, pch=1, col= floor(abs(DF_meer$ce75_e)+3))
reem3=meer2$residuals
shapiro.test(reem3)##
## Shapiro-Wilk normality test
##
## data: reem3
## W = 0.992, p-value = 0.11moran_mee = moran.mc(reem3, Wve, nsim = 2000)
moran_mee##
## Monte-Carlo simulation of Moran I
##
## data: reem3
## weights: Wve
## number of simulations + 1: 2001
##
## statistic = 0.129, observed rank = 2001, p-value = 5e-04
## alternative hypothesis: greaterDel tercer modelo la simulacion dos dio mejor AIC y a esta simulacion le hacemos la prueba de normalidad y Moran Montecarlo obteniendo que los residuales son normales y presentan dependencia espacial, entonces la fórmula del modelo 3 se escribe.
Cuarto Modelo
Modelo con lambda y rho \[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \\ u=\rho W u + \epsilon\]
mlrho=sacsarlm(formula=ce75~ Avg_CEa_15+DEM+SLOPE+Avg_z,data=data1,listw=Wve)
summary(mlrho)##
## Call:sacsarlm(formula = ce75 ~ Avg_CEa_15 + DEM + SLOPE + Avg_z, data = data1,
## listw = Wve)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.090770 -0.476848 -0.031738 0.518306 2.235748
##
## Type: sac
## Coefficients: (asymptotic standard errors)
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -60.138673 15.556909 -3.8657 0.0001108
## Avg_CEa_15 0.854968 0.072736 11.7543 < 2.2e-16
## DEM 0.028548 0.018830 1.5160 0.1295090
## SLOPE -0.062334 0.021580 -2.8885 0.0038705
## Avg_z 0.187865 0.031034 6.0535 1.417e-09
##
## Rho: 0.97458
## Asymptotic standard error: 0.38031
## z-value: 2.5626, p-value: 0.01039
## Lambda: 0.9722
## Asymptotic standard error: 0.41632
## z-value: 2.3352, p-value: 0.019532
##
## LR test value: 179.63, p-value: < 2.22e-16
##
## Log likelihood: -366.57 for sac model
## ML residual variance (sigma squared): 0.58432, (sigma: 0.76441)
## Number of observations: 313
## Number of parameters estimated: 8
## AIC: 749.14, (AIC for lm: 924.77)ce75E4=mlrho$fitted.values ;ce75E4 # valores estimados## 1 2 3 4 5 6 7 8 9 10
## 6.1279 6.5454 7.3261 7.2249 6.7108 7.0973 6.8288 7.8528 8.0567 7.7674
## 11 12 13 14 15 16 17 18 19 20
## 7.3735 7.7735 7.5168 7.7710 7.9830 8.3830 8.5028 8.4498 8.3819 8.4946
## 21 22 23 24 25 26 27 28 29 30
## 8.4325 8.9820 8.9236 9.2259 8.7382 9.0409 8.6222 8.1162 8.0251 8.4362
## 31 32 33 34 35 36 37 38 39 40
## 9.1095 8.8510 9.5153 10.1973 9.5166 9.2860 8.3656 8.2796 7.6480 7.9943
## 41 42 43 44 45 46 47 48 49 50
## 9.3501 9.1547 8.9525 8.8529 8.5018 8.6866 8.9541 8.7487 9.8879 9.9566
## 51 52 53 54 55 56 57 58 59 60
## 8.0768 8.4212 8.2643 8.1647 7.9886 9.6201 8.9970 8.6504 8.1080 7.9581
## 61 62 63 64 65 66 67 68 69 70
## 8.0198 8.1274 8.3259 8.7074 9.2727 8.7969 9.2813 9.8812 9.4929 8.1887
## 71 72 73 74 75 76 77 78 79 80
## 8.4122 8.0270 8.0993 7.9103 9.5697 9.3214 8.5282 8.1385 8.3142 8.7200
## 81 82 83 84 85 86 87 88 89 90
## 8.7218 8.8607 9.0849 9.1607 8.7827 9.3257 10.7694 10.0157 8.5548 8.3191
## 91 92 93 94 95 96 97 98 99 100
## 8.5448 7.5050 9.4952 8.9972 8.8122 8.2552 8.2708 9.0519 9.5118 9.5925
## 101 102 103 104 105 106 107 108 109 110
## 9.2546 9.4828 9.2330 9.5531 11.1596 10.8450 10.5335 8.8002 8.8240 9.4442
## 111 112 113 114 115 116 117 118 119 120
## 8.5799 8.3924 9.1763 9.2902 9.0801 8.6752 9.1497 9.4628 9.3720 9.5926
## 121 122 123 124 125 126 127 128 129 130
## 9.4591 9.6786 10.0201 10.0506 11.0776 10.8457 10.1001 8.8730 9.2275 9.0920
## 131 132 133 134 135 136 137 138 139 140
## 8.6791 9.5535 9.1657 9.3587 9.3257 9.2435 9.2613 9.4265 9.7811 9.9565
## 141 142 143 144 145 146 147 148 149 150
## 10.5348 10.4617 10.2980 10.9100 10.9572 8.7661 9.1068 9.4980 9.2774 9.4778
## 151 152 153 154 155 156 157 158 159 160
## 9.6016 9.6049 9.3339 9.6034 9.3226 9.2299 9.4734 9.8566 10.5966 10.5077
## 161 162 163 164 165 166 167 168 169 170
## 10.7672 10.1143 10.5348 10.8880 9.8784 9.9517 9.3419 9.4631 10.1322 11.4968
## 171 172 173 174 175 176 177 178 179 180
## 9.9966 9.6195 9.9751 9.4130 9.3089 9.3705 10.7146 10.6921 10.3837 10.6335
## 181 182 183 184 185 186 187 188 189 190
## 9.9888 10.3623 11.0410 10.5721 10.3372 9.7466 9.6305 10.0688 10.8097 9.7392
## 191 192 193 194 195 196 197 198 199 200
## 9.7839 9.7485 8.9526 9.5480 9.8050 11.0798 10.6064 10.2563 10.5655 10.3542
## 201 202 203 204 205 206 207 208 209 210
## 10.2235 11.5976 11.3464 10.4292 9.4268 9.6680 12.1889 10.0098 8.7098 9.5916
## 211 212 213 214 215 216 217 218 219 220
## 9.8016 10.0388 10.0462 10.0208 10.5702 11.2650 11.4211 10.7757 10.7507 10.8214
## 221 222 223 224 225 226 227 228 229 230
## 10.2172 11.2543 11.6369 10.2041 10.3368 10.5830 11.7355 9.9328 9.8498 9.7367
## 231 232 233 234 235 236 237 238 239 240
## 10.0477 10.0404 9.5915 10.6782 10.7531 11.7348 11.2149 11.1520 11.3492 11.8576
## 241 242 243 244 245 246 247 248 249 250
## 11.7939 11.9204 10.3648 10.3651 10.0821 11.4209 9.5034 9.6436 9.9370 10.1165
## 251 252 253 254 255 256 257 258 259 260
## 9.6725 10.3769 10.9829 11.6721 11.3086 11.0769 11.2118 12.3208 12.0862 11.8367
## 261 262 263 264 265 266 267 268 269 270
## 10.3094 10.5625 9.8088 10.2186 11.3079 10.0126 10.4177 10.9205 11.3857 11.6064
## 271 272 273 274 275 276 277 278 279 280
## 11.1379 11.0012 11.9642 12.3350 11.9435 10.3894 10.8097 10.2303 10.6127 10.8255
## 281 282 283 284 285 286 287 288 289 290
## 11.6732 11.4876 11.2010 11.0737 12.0017 10.7151 10.6671 10.8368 10.4422 10.1237
## 291 292 293 294 295 296 297 298 299 300
## 10.8883 11.2880 11.4860 11.2977 11.9540 10.5797 10.6962 11.0691 10.5547 9.9604
## 301 302 303 304 305 306 307 308 309 310
## 10.8038 11.2307 10.9808 10.7611 10.9403 11.0873 9.6829 11.0031 11.1715 10.8688
## 311 312 313
## 11.6250 10.6711 10.9587DF_mlrho=data.frame(ce75,ce75E4) #dataframe
colnames(DF_mlrho) <- c("ce75_o","ce75_e")
plot(DF_mlrho$ce75_o, DF_mlrho$ce75_e)
plot (xydata[,1], xydata[,2], cex=0.09*abs(DF_mlrho$ce75_o), pch=19,col= "orange",
main= "Ce75 OBSERVADOS VS ESPERADOS (mlrho)" )
points(xydata[,1], xydata[,2], cex=1.55, pch=1, col= floor(abs(DF_mlrho$ce75_e)+3))
En el cuarto modelo, lambda y rho son significativos siendo bueno el modelo y el que mejor tiene AIC y por lo tanto es el que sirve posiblemente para explicar la ce75 es: \[Y=\lambda W Y+ \alpha 1_n+ X\beta + u \\ u=\rho W u + \epsilon\]
De los cuatro modelos trabajados, el que mejor explica la conductividad electrica aparente (CEa) es el modelo de lambda y rho (mlrho) con un AIC de 749.14. Este modelo utiliza como variables explicativas a la conductividad electrica 1.5m (Avg_CEa_15), el modelo digital de elevacion (DEM), la pendiente (SLOPE) y y la altura sobre el nivel del mar (Avg_z,). La CEa esta relacionada con la altitud y la pendiende, en donde se reporta que los sitios con mayor altitud y menor pendiente presentan mayor CEa en el suelo (Coitiño et al 2015). por su parte los DEM como su nombre lo indica es un modelo que se basa en datos obtenidos de un lugar para hacer una representacion de la realidad, y toma en cuenta datos de elevacion, pendiente, por lo tanto eso explica la relacion que tiene con la CEa. por su parte. La CEa juega un papel muy importante en la agriculturav debido a que puede ayudar a estimar la variabilidad espacial del suelo y por lo tanto su productividad (Corwin et al., 2006 e Paggi et al., 2013)
Coitiño-López, Javier, Barbazán, Mónica, & Ernst, Oswaldo. (2015). Conductividad eléctrica aparente para delimitar zonas de manejo en un suelo agrícola con reducida variabilidad en propiedades físico-químicas. Agrociencia Uruguay, 19(1), 102-111. Recuperado en 12 de diciembre de 2020, de http://www.scielo.edu.uy/scielo.php?script=sci_arttext&pid=S2301-15482015000100012&lng=es&tlng=es.
Corwin, DL; SM Lesch; JD Oster & SR Kaffka. 2006. Monitoring management-induced spatio-temporal changes in soil quality through soil sampling directed by apparent electrical conductivity. Geoderma 131: 369-387
Identificacion de series de suelos mediante el uso de sensores de conductividad eléctrica aparente en el sudoeste Bonraense. Paggi M., Peralta M., Calandroni M., Cabria F., Costa J. y Aparicio V. CIENC SUELO (ARGENTINA) 31(2): 175-188, 2013
.