Asignación final computación estadística

Nicolle Salamanca

library(corrplot)

package 㤼㸱corrplot㤼㸲 was built under R version 3.6.3corrplot 0.84 loaded

library(psych)

package 㤼㸱psych㤼㸲 was built under R version 3.6.3

library(Metrics)

package 㤼㸱Metrics㤼㸲 was built under R version 3.6.3

library(ggcorrplot)

package 㤼㸱ggcorrplot㤼㸲 was built under R version 3.6.3Loading required package: ggplot2
package 㤼㸱ggplot2㤼㸲 was built under R version 3.6.3
Attaching package: 㤼㸱ggplot2㤼㸲

The following objects are masked from 㤼㸱package:psych㤼㸲:

    %+%, alpha

library(ggpubr)

package 㤼㸱ggpubr㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'data.table':
  method           from
  print.data.table     

library(plotly)

package 㤼㸱plotly㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'htmlwidgets':
  method           from         
  print.htmlwidget tools:rstudio

Attaching package: 㤼㸱plotly㤼㸲

The following object is masked from 㤼㸱package:ggplot2㤼㸲:

    last_plot

The following object is masked from 㤼㸱package:stats㤼㸲:

    filter

The following object is masked from 㤼㸱package:graphics㤼㸲:

    layout

library(spdep)

package 㤼㸱spdep㤼㸲 was built under R version 3.6.3Loading required package: sp
package 㤼㸱sp㤼㸲 was built under R version 3.6.3Loading required package: spData
package 㤼㸱spData㤼㸲 was built under R version 3.6.3To 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
package 㤼㸱sf㤼㸲 was built under R version 3.6.3Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

library(ape)

package 㤼㸱ape㤼㸲 was built under R version 3.6.3Registered S3 method overwritten by 'ape':
  method   from 
  plot.mst spdep

Attaching package: 㤼㸱ape㤼㸲

The following object is masked from 㤼㸱package:ggpubr㤼㸲:

    rotate

library(sp)
library(MVA)

package 㤼㸱MVA㤼㸲 was built under R version 3.6.3Loading required package: HSAUR2
package 㤼㸱HSAUR2㤼㸲 was built under R version 3.6.3Loading required package: tools

library(Hmisc)

package 㤼㸱Hmisc㤼㸲 was built under R version 3.6.3Loading required package: lattice
Loading required package: survival
package 㤼㸱survival㤼㸲 was built under R version 3.6.3Loading required package: Formula
package 㤼㸱Formula㤼㸲 was built under R version 3.6.3
Attaching package: 㤼㸱Hmisc㤼㸲

The following object is masked from 㤼㸱package:ape㤼㸲:

    zoom

The following object is masked from 㤼㸱package:plotly㤼㸲:

    subplot

The following object is masked from 㤼㸱package:psych㤼㸲:

    describe

The following objects are masked from 㤼㸱package:base㤼㸲:

    format.pval, units

library(normtest)
library(nortest)
library(readxl)

dm<- BD_MODELADO

dmm<- as.matrix(dm[,c(3:8)])
r75 <- dmm[,1]  #Variable respuesta 75 cm
r150 <- dmm[,2] #Variable respuesta 150 
e75 <- dmm[,2:6] #Variable explicativa 75 cm
e150 <- dmm[,c(1,3:6)] #Variable explicativa 150 cm
me75 <- as.matrix(e75) 
me150 <- as.matrix(e150)
xydata <- as.matrix(dm[,1:2]) #coordenadas para elaborar la matriz de pesos
mpesos <- as.matrix(dist(xydata, diag = T, upper = T)) #matriz de pesos "W"
mpesosinv <- as.matrix(1/mpesos)
diag(mpesosinv) <- 0
W <- as.matrix(mpesosinv)
suma <- apply(W, 1, sum)
We <- W/suma #Matriz de pesos estandarizada
apply(We, 1, sum) #Si la matriz de pesos estandarizada es correcta, la suma de cada fila debe ser 1

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

moran75 <- Moran.I(r75, We);moran75 #Indices de Moran

$observed
[1] 0.2687468

$expected
[1] -0.003205128

$sd
[1] 0.004665906

$p.value
[1] 0

moran150 <- Moran.I(r150, We);moran150

$observed
[1] 0.160951

$expected
[1] -0.003205128

$sd
[1] 0.00465455

$p.value
[1] 0

moranNDVI <- Moran.I(dmm[,3], We);moranNDVI

$observed
[1] 0.09750403

$expected
[1] -0.003205128

$sd
[1] 0.004644979

$p.value
[1] 0

moranDEM <- Moran.I(dmm[,4], We);moranDEM

$observed
[1] 0.3096708

$expected
[1] -0.003205128

$sd
[1] 0.004672384

$p.value
[1] 0

moranSLOPE <- Moran.I(dmm[,5], We);moranSLOPE

$observed
[1] 0.06993324

$expected
[1] -0.003205128

$sd
[1] 0.004654307

$p.value
[1] 0

moranz <- Moran.I(dmm[,6], We);moranz

$observed
[1] 0.3505031

$expected
[1] -0.003205128

$sd
[1] 0.004667935

$p.value
[1] 0

#Se observa dependencia espacial en los datos

#Matrices de correlación de Pearson y Spearman
mcp <- rcorr(as.matrix(dmm[,1:6]), type = 'p') 
mcs <- rcorr(as.matrix(dmm[,1:6]), type = 's')
mcorp <- mcp$r
mcors <- mcs$r

cp<-ggcorrplot(mcorp,hc.order = TRUE,type = "upper",outline.color = "white",ggtheme = ggplot2::theme_gray,lab = TRUE,lab_size = 2,tl.cex = 10)+labs(title = "Pearson")
cs<-ggcorrplot(mcors,hc.order = TRUE,type = "upper",outline.color = "white",ggtheme = ggplot2::theme_gray,lab = TRUE,lab_size = 2,tl.cex = 10)+labs(title = "Spearman")
ggarrange(cp,cs) + labs(title = "Spearman")

#La variable DEM presenta la mayor correlaciión positiva con la variable CE 75 y Z

desc <- describe(dmm[,1:6]);desc

dmm[, 1:6] 

 6  Variables      313  Observations
-------------------------------------------------------------------------------
Avg_CEa_07 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      313        1    9.769    1.728    7.193    7.950 
     .25      .50      .75      .90      .95 
   8.808    9.645   10.835   11.830   12.446 

lowest :  6.165389  6.180742  6.283077  6.358915  6.561038
highest: 13.227629 13.262604 13.318076 13.424340 13.601094
-------------------------------------------------------------------------------
Avg_CEa_15 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      313        1     18.5   0.8114    17.41    17.64 
     .25      .50      .75      .90      .95 
   18.03    18.44    18.88    19.45    19.68 

lowest : 16.80343 16.80468 16.88126 16.97415 17.01939
highest: 20.57460 20.68170 20.93098 21.20233 21.51339
-------------------------------------------------------------------------------
NDVI 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      313        1   0.8346  0.03075   0.7753   0.7987 
     .25      .50      .75      .90      .95 
  0.8233   0.8404   0.8552   0.8652   0.8702 

lowest : 0.704767 0.715250 0.737705 0.739356 0.750449
highest: 0.874883 0.874916 0.875065 0.876872 0.877252
-------------------------------------------------------------------------------
DEM 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      143        1    205.1    5.245    197.7    199.1 
     .25      .50      .75      .90      .95 
   201.5    204.8    209.3    211.5    212.3 

lowest : 196.1667 196.7500 196.8333 196.8889 197.0000
highest: 212.7500 212.8333 213.1111 213.1667 213.3333
-------------------------------------------------------------------------------
SLOPE 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      299        1    4.128    2.374    1.361    1.667 
     .25      .50      .75      .90      .95 
   2.519    3.781    5.385    6.911    7.796 

lowest :  0.210996  0.281328  0.397830  0.455108  0.577682
highest: 11.049718 11.253322 11.296110 11.876602 12.718485
-------------------------------------------------------------------------------
Avg_z 
       n  missing distinct     Info     Mean      Gmd      .05      .10 
     313        0      313        1    202.5    4.174    196.5    197.8 
     .25      .50      .75      .90      .95 
   199.8    202.5    205.5    206.8    207.6 

lowest : 193.0512 193.2986 193.5659 193.9931 194.4116
highest: 208.9608 210.9083 211.9958 212.1226 212.3547
-------------------------------------------------------------------------------

par(mfrow = c(1,2))
boxplot(dmm[,1], main = 'Boxplot CE 75', ylab= "CE 75 cm " , col= "dark blue")
hist(dmm[,1], main = 'Histograma CE 75', xlab= "CE 75 cm ", col= "dark blue")

par(mfrow = c(1,2))
boxplot(dmm[,2], main = 'Boxplot CE 150', ylab= "CE 150 cm " , col= "dark red")
hist(dmm[,2], main = 'Histograma CE 150', xlab= "CE 150 cm ", col= "dark red")

#Evaluación de la normalidad
cvm75<-cvm.test(BD_MODELADO$Avg_CEa_07);cvm75 

    Cramer-von Mises normality test

data:  BD_MODELADO$Avg_CEa_07
W = 0.13062, p-value = 0.04289

cvm150<-cvm.test(BD_MODELADO$Avg_CEa_15);cvm150

    Cramer-von Mises normality test

data:  BD_MODELADO$Avg_CEa_15
W = 0.19809, p-value = 0.005641

sh75<-sf.test(BD_MODELADO$Avg_CEa_07); sh75

    Shapiro-Francia normality test

data:  BD_MODELADO$Avg_CEa_07
W = 0.99381, p-value = 0.2013

sh150<-sf.test(BD_MODELADO$Avg_CEa_15); sh150

    Shapiro-Francia normality test

data:  BD_MODELADO$Avg_CEa_15
W = 0.97701, p-value = 0.0001414

#Modelado
#Modelo clásico no espacial
modc75 <- lm(BD_MODELADO$Avg_CEa_07~BD_MODELADO$SLOPE+BD_MODELADO$Avg_z+ BD_MODELADO$Avg_CEa_15) 
summary(modc75)

Call:
lm(formula = BD_MODELADO$Avg_CEa_07 ~ BD_MODELADO$SLOPE + BD_MODELADO$Avg_z + 
    BD_MODELADO$Avg_CEa_15)

Residuals:
     Min       1Q   Median       3Q      Max 
-2.41298 -0.71115 -0.05541  0.64834  3.01718 

Coefficients:
                        Estimate Std. Error t value Pr(>|t|)    
(Intercept)            -73.98404    4.74765 -15.583  < 2e-16 ***
BD_MODELADO$SLOPE       -0.12525    0.02799  -4.475 1.07e-05 ***
BD_MODELADO$Avg_z        0.33680    0.01834  18.369  < 2e-16 ***
BD_MODELADO$Avg_CEa_15   0.86887    0.09270   9.373  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.05 on 309 degrees of freedom
Multiple R-squared:  0.5315,    Adjusted R-squared:  0.527 
F-statistic: 116.9 on 3 and 309 DF,  p-value: < 2.2e-16

modc150 <- lm(BD_MODELADO$Avg_CEa_15~BD_MODELADO$SLOPE+BD_MODELADO$Avg_z+ BD_MODELADO$Avg_CEa_07) 
summary(modc150)

Call:
lm(formula = BD_MODELADO$Avg_CEa_15 ~ BD_MODELADO$SLOPE + BD_MODELADO$Avg_z + 
    BD_MODELADO$Avg_CEa_07)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.50608 -0.39071  0.01281  0.36509  2.04703 

Coefficients:
                       Estimate Std. Error t value Pr(>|t|)    
(Intercept)            47.55651    2.11792  22.454  < 2e-16 ***
BD_MODELADO$SLOPE       0.07595    0.01503   5.053 7.45e-07 ***
BD_MODELADO$Avg_z      -0.15733    0.01123 -14.009  < 2e-16 ***
BD_MODELADO$Avg_CEa_07  0.25480    0.02718   9.373  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.5686 on 309 degrees of freedom
Multiple R-squared:  0.4107,    Adjusted R-squared:  0.405 
F-statistic: 71.78 on 3 and 309 DF,  p-value: < 2.2e-16

res75mc <- modc75$residuals
shapiro.test(res75mc)

    Shapiro-Wilk normality test

data:  res75mc
W = 0.99092, p-value = 0.05044

res150mc <- modc150$residuals
shapiro.test(res150mc)

    Shapiro-Wilk normality test

data:  res150mc
W = 0.99144, p-value = 0.06673

cvm.test(res75mc)

    Cramer-von Mises normality test

data:  res75mc
W = 0.1071, p-value = 0.08941

cvm.test(res150mc)

    Cramer-von Mises normality test

data:  res150mc
W = 0.035689, p-value = 0.7586

moranres75mc<-(Moran.I(res75mc,We));moranres75mc 

$observed
[1] 0.1580117

$expected
[1] -0.003205128

$sd
[1] 0.004665226

$p.value
[1] 0

moranres150mc<-(Moran.I(res150mc,We));moranres150mc

$observed
[1] 0.08362486

$expected
[1] -0.003205128

$sd
[1] 0.004659635

$p.value
[1] 0

cla75e<-modc75$fitted.values
cla75res<-modc75$residuals
dfclas75<-data.frame(cla75e,cla75res,BD_MODELADO$Avg_CEa_07 ,index<-c(1:313))

#Comparación CE 75 estimada y observada por el modelo 
plot(BD_MODELADO$Avg_CEa_07, cla75e, cex=0.5, pch =18, col = 'dark green', xlab= "CE 75 observada", ylab="CE 75 estimada")

cla150e<-modc150$fitted.values
cla150res<-modc150$residuals
dfclas150<-data.frame(cla150e,cla150res,BD_MODELADO$Avg_CEa_15 ,index<-c(1:313))

#Comparación CE 150 estimada y observada por el modelo 
plot(BD_MODELADO$Avg_CEa_15, cla150e, cex=0.5, pch =18, col = 'dark green', xlab= "CE 150 observada", ylab="CE 150 estimada")

#Como se muestra dependencia especial, no es aplicable un modelo clásico

#Modelo autoregresivo puro
contnb <- dnearneigh(coordinates(xydata),0,854, longlat = F)
contnb

Neighbour list object:
Number of regions: 313 
Number of nonzero links: 97656 
Percentage nonzero weights: 99.68051 
Average number of links: 312 

dlist <- nbdists(contnb, xydata)
dlist <- lapply(dlist, function(x) 1/x)
wve <- nb2listw(contnb, glist = dlist, style = 'W')

#Modelo 75 cm
map75 <- spautolm(BD_MODELADO$Avg_CEa_07~1, data = BD_MODELADO, listw = wve) 

Function spautolm moved to the spatialreg packageRegistered 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

summary(map75)

Call: spatialreg::spautolm(formula = formula, data = data, listw = listw, 
    na.action = na.action, family = family, method = method, 
    verbose = verbose, trs = trs, interval = interval, zero.policy = zero.policy, 
    tol.solve = tol.solve, llprof = llprof, control = control)

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.011871 

Log likelihood: -494.8231 
ML residual variance (sigma squared): 1.347, (sigma: 1.1606)
Number of observations: 313 
Number of parameters estimated: 3 
AIC: 995.65

resmap75<-map75$fit$residuals
moranmcmap75<-moran.mc(resmap75,wve,nsim = 2000);moranmcmap75

    Monte-Carlo simulation of Moran I

data:  resmap75 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.16722, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(resmap75)

    Shapiro-Wilk normality test

data:  resmap75
W = 0.99351, p-value = 0.1971

ma75e <- as.data.frame(map75$fit['fitted.values'])
ma75e1 <- map75$fit$fitted.values
dfma75 <- data.frame(r75, ma75e1)
colnames(dfma75) <-  c('CE75_obs','CE75_est')
resmap75 <- map75$fit$residuals
cor(dfma75$CE75_obs, dfma75$CE75_est) #correlacción entre lo estimado y lo observado

[1] 0.7977199

#Comparación CE 75 estimada y observada por el modelo MAP
plot(dfma75$CE75_obs, dfma75$CE75_est, cex=0.5, pch =18, col = 'dark green', xlab= "CE 75 observada", ylab="CE 75 estimada")

datamap<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfma75$CE75_obs, dfma75$CE75_est)


colnames(datamap) <- c('x', 'y', 'CE75_obs', 'CE75_est')
MAP<-ggplot(data = datamap, aes(x = x, y = y)) +
  geom_point(cex = datamap$CE75_obs*0.2) +
  geom_point(color = datamap$CE75_est)
MAP

#Modelo 150 cm
map150 <- spautolm(BD_MODELADO$Avg_CEa_15~1, data = BD_MODELADO, listw = wve) 

Function spautolm moved to the spatialreg package

summary(map150)

Call: spatialreg::spautolm(formula = formula, data = data, listw = listw, 
    na.action = na.action, family = family, method = method, 
    verbose = verbose, trs = trs, interval = interval, zero.policy = zero.policy, 
    tol.solve = tol.solve, llprof = llprof, control = control)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.453255 -0.397645 -0.042934  0.322283  2.953512 

Coefficients: 
            Estimate Std. Error z value  Pr(>|z|)
(Intercept)  19.3151     1.5515   12.45 < 2.2e-16

Lambda: 0.97691 LR test value: 86.774 p-value: < 2.22e-16 
Numerical Hessian standard error of lambda: 0.023 

Log likelihood: -304.7918 
ML residual variance (sigma squared): 0.40168, (sigma: 0.63378)
Number of observations: 313 
Number of parameters estimated: 3 
AIC: 615.58

resmap150<-map150$fit$residuals
moranmcmap150<-moran.mc(resmap150,wve,nsim = 2000);moranmcmap150

    Monte-Carlo simulation of Moran I

data:  resmap150 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.094365, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(resmap150)

    Shapiro-Wilk normality test

data:  resmap150
W = 0.95729, p-value = 6.37e-08

ma150e <- as.data.frame(map150$fit['fitted.values'])
ma150e1 <- map150$fit$fitted.values
dfma150 <- data.frame(r150, ma150e1)
colnames(dfma150) <-  c('CE150_obs','CE150_est')
resmap150 <- map150$fit$residuals
cor(dfma150$CE150_obs, dfma150$CE150_est) #correlacción entre lo estimado y lo observado

[1] 0.6642671

#Comparación CE 150 estimada y observada por el modelo MAP
plot(dfma150$CE150_obs, dfma150$CE150_est, cex=0.5, pch =18, col = 'dark green', xlab= "CE 150 observada", ylab="CE 150 estimada")

datamap150<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfma150$CE150_obs, dfma150$CE150_est)


colnames(datamap150) <- c('x', 'y', 'CE150_obs', 'CE150_est')
MAP150<-ggplot(data = datamap150, aes(x = x, y = y)) +
  geom_point(cex = datamap150$CE150_obs*0.2) +
  geom_point(color = datamap150$CE150_est)
MAP150

#Modelo de autocorrelación espacial SAC 75 cm
sac<-sacsarlm(formula=BD_MODELADO$Avg_CEa_07~ BD_MODELADO$DEM+ BD_MODELADO$Avg_CEa_15,data=BD_MODELADO,listw=wve)

Function sacsarlm moved to the spatialreg package

summary(sac)

Call:spatialreg::sacsarlm(formula = formula, data = data, listw = listw, 
    listw2 = listw2, na.action = na.action, Durbin = Durbin, 
    type = type, method = method, quiet = quiet, zero.policy = zero.policy, 
    tol.solve = tol.solve, llprof = llprof, interval1 = interval1, 
    interval2 = interval2, trs1 = trs1, trs2 = trs2, control = control)

Residuals:
      Min        1Q    Median        3Q       Max 
-2.237973 -0.565571  0.025546  0.503591  2.576355 

Type: sac 
Coefficients: (asymptotic standard errors) 
                         Estimate Std. Error z value  Pr(>|z|)
(Intercept)            -39.956408  66.055529 -0.6049    0.5453
BD_MODELADO$DEM          0.121747   0.022010  5.5313 3.178e-08
BD_MODELADO$Avg_CEa_15   0.708038   0.077897  9.0894 < 2.2e-16

Rho: 0.98007
Asymptotic standard error: 0.65942
    z-value: 1.4863, p-value: 0.13721
Lambda: 0.97946
Asymptotic standard error: 0.67979
    z-value: 1.4408, p-value: 0.14963

LR test value: 233.02, p-value: < 2.22e-16

Log likelihood: -395.063 for sac model
ML residual variance (sigma squared): 0.69855, (sigma: 0.83579)
Number of observations: 313 
Number of parameters estimated: 6 
AIC: 802.13, (AIC for lm: 1031.1)

sac75e<-sac$fitted.values
sac75re<-sac$residuals
moran75sac <- moran.mc(sac75re, wve, nsim = 2000);moran75sac

    Monte-Carlo simulation of Moran I

data:  sac75re 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.11394, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(sac75re)

    Shapiro-Wilk normality test

data:  sac75re
W = 0.99707, p-value = 0.8468

dfsac75<-data.frame(BD_MODELADO$Avg_CEa_07,sac75e)
colnames(dfsac75)<- c("CE75_Obsac","CE75_Estsac")
rmse(dfsac75$CE75_Obsac,dfsac75$CE75_Estsac)

[1] 0.835792

#Comparación CE 150 estimada y observada por el modelo 
plot(dfsac75$CE75_Obsac, dfsac75$CE75_Estsac, cex=0.5, pch =18, col = 'dark green', xlab= "CE 75 observada", ylab="CE 75 estimada")

cor(dfsac75$CE75_Obsac,dfsac75$CE75_Estsac)

[1] 0.8544081

datasac75<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfsac75$CE75_Obsac, dfsac75$CE75_Estsac)


colnames(datasac75) <- c('x', 'y', 'CE75_Obsac', 'CE75_Estsac')
SAC75<-ggplot(data = datasac75, aes(x = x, y = y)) +
  geom_point(cex = datasac75$CE75_Obsac*0.2) +
  geom_point(color = datasac75$CE75_Estsac)
SAC75

#Modelo de autocorrelación espacial SAC 150 cm
sac150<-sacsarlm(formula=BD_MODELADO$Avg_CEa_15 ~ BD_MODELADO$DEM+ BD_MODELADO$Avg_CEa_07,data=BD_MODELADO,listw=wve)

Function sacsarlm moved to the spatialreg package

summary(sac150)

Call:spatialreg::sacsarlm(formula = formula, data = data, listw = listw, 
    listw2 = listw2, na.action = na.action, Durbin = Durbin, 
    type = type, method = method, quiet = quiet, zero.policy = zero.policy, 
    tol.solve = tol.solve, llprof = llprof, interval1 = interval1, 
    interval2 = interval2, trs1 = trs1, trs2 = trs2, control = control)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.577747 -0.324598 -0.064869  0.311256  2.152924 

Type: sac 
Coefficients: (asymptotic standard errors) 
                        Estimate Std. Error z value  Pr(>|z|)
(Intercept)             9.856517  21.001115  0.4693    0.6388
BD_MODELADO$DEM        -0.052041   0.010629 -4.8960 9.780e-07
BD_MODELADO$Avg_CEa_07  0.229127   0.034447  6.6516 2.899e-11

Rho: 0.96203
Asymptotic standard error: 0.53967
    z-value: 1.7826, p-value: 0.074647
Lambda: 0.97103
Asymptotic standard error: 0.41252
    z-value: 2.3539, p-value: 0.018576

LR test value: 125.32, p-value: < 2.22e-16

Log likelihood: -248.1623 for sac model
ML residual variance (sigma squared): 0.27499, (sigma: 0.52439)
Number of observations: 313 
Number of parameters estimated: 6 
AIC: 508.32, (AIC for lm: 629.64)

sac150e<-sac$fitted.values
sac150re<-sac$residuals
moran150sac <- moran.mc(sac150re, wve, nsim = 2000);moran150sac

    Monte-Carlo simulation of Moran I

data:  sac150re 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.11394, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(sac150re)

    Shapiro-Wilk normality test

data:  sac150re
W = 0.99707, p-value = 0.8468

dfsac150<-data.frame(BD_MODELADO$Avg_CEa_15 ,sac150e)
colnames(dfsac150)<- c("CE150_Obsac","CE150_Estsac")
rmse(dfsac150$CE150_Obsac,dfsac150$CE150_Estsac)

[1] 8.821918

#Comparación CE 150 estimada y observada por el modelo 
plot(dfsac150$CE150_Obsac, dfsac150$CE150_Estsac, cex=0.5, pch =18, col = 'dark green', xlab= "CE 150 observada", ylab="CE 150 estimada")

cor(dfsac150$CE150_Obsac,dfsac150$CE150_Estsac)

[1] 0.05522331

datasac150<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfsac150$CE150_Obsac, dfsac150$CE150_Estsac)


colnames(datasac150) <- c('x', 'y', 'CE150_Obsac', 'CE150_Estsac')
SAC150<-ggplot(data = datasac150, aes(x = x, y = y)) +
  geom_point(cex = datasac150$CE150_Obsac*0.2) +
  geom_point(color = datasac150$CE150_Estsac)
SAC150

#Modelo espacial del error SEM 75 CM
ese75 <- errorsarlm(BD_MODELADO$Avg_CEa_07~BD_MODELADO$SLOPE+ BD_MODELADO$Avg_z+BD_MODELADO$Avg_CEa_15+ BD_MODELADO$DEM,listw=wve)

Function errorsarlm moved to the spatialreg package

summary(ese75) 

Call:spatialreg::errorsarlm(formula = formula, data = data, listw = listw, 
    na.action = na.action, Durbin = Durbin, etype = etype, method = method, 
    quiet = quiet, zero.policy = zero.policy, interval = interval, 
    tol.solve = tol.solve, trs = trs, control = control)

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
BD_MODELADO$SLOPE       -0.074849   0.024782  -3.0203  0.002525
BD_MODELADO$Avg_z        0.251732   0.028220   8.9203 < 2.2e-16
BD_MODELADO$Avg_CEa_15   0.871288   0.082765  10.5273 < 2.2e-16
BD_MODELADO$DEM          0.039380   0.020925   1.8819  0.059845

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.8867 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)

ese75e<-ese75$fitted.values
ese75res <- ese75$residuals

moran75ese <- moran.mc(ese75res, wve, nsim = 2000);moran75ese

    Monte-Carlo simulation of Moran I

data:  ese75res 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.12895, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(ese75res)

    Shapiro-Wilk normality test

data:  ese75res
W = 0.99235, p-value = 0.1078

dfese75<-data.frame(BD_MODELADO$Avg_CEa_07 ,ese75e)
colnames(dfese75)<- c("CE75_Obsese","CE75_Estese")
rmse(dfese75$CE75_Obsese,dfese75$CE75_Estese)

[1] 0.8774362

cvm.test(ese75res)

    Cramer-von Mises normality test

data:  ese75res
W = 0.084219, p-value = 0.182

#Comparación CE 75 estimada y observada por el modelo 
plot(dfese75$CE75_Obsese, dfese75$CE75_Estese, cex=0.5, pch =18, col = 'dark green', xlab= "CE 75 observada", ylab="CE 75 estimada")

cor(dfese75$CE75_Obsese,dfese75$CE75_Estese)

[1] 0.8246131

dataese75<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfese75$CE75_Obsese, dfese75$CE75_Estese)


colnames(dataese75) <- c('x', 'y', 'CE75_Obsese', 'CE75_Estese')
ESEM75<-ggplot(data = dataese75, aes(x = x, y = y)) +
  geom_point(cex = dataese75$CE75_Obsese*0.2) +
  geom_point(color = dataese75$CE75_Estese)
ESEM75

#Modelo espacial del error SEM 150 CM
ese150 <- errorsarlm(BD_MODELADO$Avg_CEa_15 ~BD_MODELADO$SLOPE+ BD_MODELADO$Avg_z+BD_MODELADO$Avg_CEa_07+ BD_MODELADO$DEM,listw=wve)

Function errorsarlm moved to the spatialreg package

summary(ese150) 

Call:spatialreg::errorsarlm(formula = formula, data = data, listw = listw, 
    na.action = na.action, Durbin = Durbin, etype = etype, method = method, 
    quiet = quiet, zero.policy = zero.policy, interval = interval, 
    tol.solve = tol.solve, trs = trs, control = control)

Residuals:
      Min        1Q    Median        3Q       Max 
-1.477728 -0.316994 -0.014091  0.367283  2.009859 

Type: error 
Coefficients: (asymptotic standard errors) 
                        Estimate Std. Error z value  Pr(>|z|)
(Intercept)            45.035550   2.729375 16.5003 < 2.2e-16
BD_MODELADO$SLOPE       0.051310   0.014481  3.5432 0.0003953
BD_MODELADO$Avg_z      -0.136386   0.016857 -8.0910 6.661e-16
BD_MODELADO$Avg_CEa_07  0.299387   0.028492 10.5079 < 2.2e-16
BD_MODELADO$DEM        -0.008240   0.012332 -0.6682 0.5040078

Lambda: 0.96901, LR test value: 50.337, p-value: 1.2949e-12
Asymptotic standard error: 0.021843
    z-value: 44.363, p-value: < 2.22e-16
Wald statistic: 1968.1, p-value: < 2.22e-16

Log likelihood: -239.9531 for error model
ML residual variance (sigma squared): 0.26594, (sigma: 0.51569)
Number of observations: 313 
Number of parameters estimated: 7 
AIC: 493.91, (AIC for lm: 542.24)

ese150e<-ese150$fitted.values
ese150res <- ese150$residuals

moran150ese <- moran.mc(ese150res, wve, nsim = 2000);moran150ese

    Monte-Carlo simulation of Moran I

data:  ese150res 
weights: wve  
number of simulations + 1: 2001 

statistic = 0.07571, observed rank = 2001, p-value = 0.0004998
alternative hypothesis: greater

shapiro.test(ese150res)

    Shapiro-Wilk normality test

data:  ese150res
W = 0.98823, p-value = 0.01217

dfese150<-data.frame(BD_MODELADO$Avg_CEa_15 ,ese150e)
colnames(dfese150)<- c("CE150_Obsese","CE150_Estese")
rmse(dfese150$CE150_Obsese,dfese150$CE150_Estese)

[1] 0.5156905

cvm.test(ese150res)

    Cramer-von Mises normality test

data:  ese150res
W = 0.068, p-value = 0.2961

#Comparación CE 150 estimada y observada por el modelo 
plot(dfese150$CE150_Obsese, dfese150$CE150_Estese, cex=0.5, pch =18, col = 'dark green', xlab= "CE 150 observada", ylab="CE 150 estimada")

cor(dfese150$CE150_Obsese,dfese150$CE150_Estese)

[1] 0.7184014

dataese150<-data.frame(x = dm$Avg_X_MCB , y = dm$Avg_Y_MCE, dfese150$CE150_Obsese, dfese150$CE150_Estese)


colnames(dataese150) <- c('x', 'y', 'CE150_Obsese', 'CE150_Estese')
ESEM150<-ggplot(data = dataese150, aes(x = x, y = y)) +
  geom_point(cex = dataese150$CE150_Obsese*0.2) +
  geom_point(color = dataese150$CE150_Estese)
ESEM150

tabla<- modelos
tabla

NA

De acuerdo a los valores de AIC, el modelo más adecuado para 75 cm es SAC, y para 150 cm es SEM.

#Variables relacionadas

ce75vsaltitud <- plot_ly(x = BD_MODELADO$Avg_X_MCB, y = BD_MODELADO$Avg_Y_MCE, z = BD_MODELADO$Avg_z, size = I(90))%>%
            layout(title = 'CE 75 vs Altitud (z)',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Altitud")
    )
  )%>%
add_markers(color = BD_MODELADO$Avg_CEa_07)
ce75vsaltitud

`arrange_()` is deprecated as of dplyr 0.7.0.
Please use `arrange()` instead.
See vignette('programming') for more help
This warning is displayed once every 8 hours.
Call `lifecycle::last_warnings()` to see where this warning was generated.

ce150vsaltitud <- plot_ly(x = BD_MODELADO$Avg_X_MCB, y = BD_MODELADO$Avg_Y_MCE, z = BD_MODELADO$Avg_z, size = I(90))%>%
            layout(title = 'CE 150 vs Altitud (z)',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Altitud")
    )
  )%>%
add_markers(color = BD_MODELADO$Avg_CEa_15)
ce150vsaltitud

ce75vsslope <- plot_ly(x = BD_MODELADO$Avg_X_MCB, y = BD_MODELADO$Avg_Y_MCE,  z = BD_MODELADO$SLOPE, size = I(90))%>%
            layout(title = 'CE 75 vs Slope',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Slope")
    )
  )%>%
add_markers(color = BD_MODELADO$Avg_CEa_07)
ce75vsslope

ce150vsslope <- plot_ly(x = BD_MODELADO$Avg_X_MCB, y = BD_MODELADO$Avg_Y_MCE,  z = BD_MODELADO$SLOPE, size = I(90))%>%
            layout(title = 'CE 150 vs Slope',
                  scene = list(
                              xaxis = list(title = "Longitud"),
                              yaxis = list(title = "Latitud"),
                              zaxis = list(title = "Slope")
    )
  )%>%
add_markers(color = BD_MODELADO$Avg_CEa_15)
ce150vsslope

#Coordenada de interés
plot(x=BD_MODELADO$Avg_X_MCB,y=BD_MODELADO$Avg_Y_MCE,pch=18,col="blue", xlab = 'Longitud',ylab = 'Latitud')
points(x=843473,y=956093,col="red",pch=18)

X<- BD_MODELADO$Avg_X_MCB
Y<- BD_MODELADO$Avg_Y_MCE
CE75<-BD_MODELADO$Avg_CEa_07
CE150<- BD_MODELADO$Avg_CEa_15
NDVI<- BD_MODELADO$NDVI
DEM<- BD_MODELADO$DEM
SLOPE<- BD_MODELADO$SLOPE
Z<- BD_MODELADO$Avg_z

#Nueva matriz de datos
nuevop <- data.frame(843473,956093,0,0,0,0,0,0) 
names(nuevop) <- c("X","Y","CE75","CE150","NDVI","DEM","SLOPE","Z")
nuevoW<- rbind(dm,nuevop)
nw<-as.matrix(dist(cbind(nuevoW$X, nuevoW$Y)))
invnw<- 1/nw
invnw[is.infinite(invnw)] <- 0
mident<-diag(314) #Nueva matriz de identidad