Unidade IV – Predição Espacial e Validação Estatística
Temas: O melhor preditor linear não enviesado empírico. Krigagem simples, ordinária e universal. Funções de covariância e os pesos da krigagem. Variância do erro de predição. Suporte de predição. Validação estatística de predições espaciais. Validação cruzada.
1 Predição espacial
Predição significa dizer algo sobre um indivíduo antes de observá-lo. Segundo a terminologia geoestatística clássica, a predição espacial é chamada de krigagem.
A krigagem é um método de interpolação linear que usa os dados das observações e os parâmetros estimados empiricamente da função de covariância para fazer predições espaciais. Assim, de maneira geral, o valor predito do sinal \(Z\) num determinado ponto \(\boldsymbol{s}_0\) será o resultado da soma linear ou média ponderada dos dados em sua vizinhança. O principal detalhe da krigagem é a determinação dos pesos usados para calcular essa média ponderada.
Como na maioria dos problemas de predição linear, o objetivo da krigagem é fazer uma predição que minimize o erro de predição, usando como medida a média do erro quadrado de predição \(\text{E}[\hat{Z} - Z]^2\). Assim, os pesos da krigagem são escolhidos de maneira que a média do erro quadrado de predição seja mínima. Nesse sentido, podemos dizer que a krigagem é um método de predição linear ótimo e que produz valores não-enviezados, ou seja, o melhor preditor linear não enviesado empírico.
1.1 Krigagem simples e krigagem ordinária
Segundo a terminologia geoestatística clássica, a krigagem simples consiste na predição espacial dos efeitos aleatórios, mais especificamente, de uma superfície \(\hat{B}(\boldsymbol{s})\). Conforme vimos anteriormente, os efeitos aleatórios são uma variável aleatória não-observável, que apresenta correlação espacial, e média igual à 0 que, assim como a variância, é constante em toda a área de estudo. Contudo, o termo krigagem simples também é usado para os casos em que nós simplesmente conhecemos a média – ou podemos assumir que a conhecemos – e, assim, decidimos modelar separadamente os efeitos fixos e os efeitos aleatórios.
Contudo, em geral, nós não conhecemos a média espacial. Nesse caso é preciso estimá-la a partir dos dados que possuímos em mãos. Se não tivermos segurança de que a média é constante em toda a área de estudo, então podemos simplesmente pressupor que seja constante em curtos intervalos de distância (hipótese intrínseca). Nesse caso, estaremos então fazendo a predição espacial tanto dos efeitos aleatórios como da média espacial (efeito fixo). Essa prática é conhecida como krigagem ordinária na terminologia geoestatística clássica.
A predição espacial do sinal \(Z\) em um ponto \(\boldsymbol{s}_0\) usando krigagem ordinária é feita usando a seguinte equação:
\[\hat{Z}(\boldsymbol{s}_0) = \sum_{i=1}^n \lambda_i y(\boldsymbol{s}_i)\]
onde \(\boldsymbol{\lambda}_i\) são os pesos da predição, também conhecidos na terminologia geoestatística clássica como pesos da krigagem. Esses são atribuídos a cada uma das \(n\) observações que possuímos em mãos. Na krigagem ordinária, a principal propriedade dos pesos é que precisam somar, obrigatoriamente, 1, ou seja,
\[\sum_{i=1}^n \lambda_i = 1\]
Essa imposição garante que as predições sejam não-enviesadas.
1.2 Krigagem universal
2 Validação estatística
O procedimento da validação cruzada consiste na partição aleatória do conjunto completo de dados em \(k\) subconjuntos com aproximadamente o mesmo número de observações. O número de subconjuntos é variável, geralmente entre 2 e \(n\).
A cada passo da validação cruzada, um dos subconjuntos é deixado separado para ser utilizado como conjunto de dados de validação. Os demais \(k - 1\) subconjuntos são utilizados para constituir o conjunto de calibração do modelo linear misto. Calibrado o modelo linear misto, faz-se a predição dos valores da variável de interesse nas observações do conjunto deixado separado, o conjunto de dados de validação. Esse procedimento é repetido até que cada subconjunto \(k\) seja, em algum momento, deixado separado para constituir o conjunto de dados de validação enquanto os outros \(k - 1\) segmentos são utilizados para estimar os parâmetros do modelo linear misto. A partir das predições realizadas para cada um dos \(k\) subconjuntos são calculados os erros para avaliar a qualidade das predições.
A validação externa difere da validação cruzada pelo fato de que os dados usados para validação não são nunca usados para estimar os parâmetros do modelo. Assim, a validação cruzada é uma medida obtida na fase inicial de trabalho, usando os dados que temos em mãos. Enquanto isso, a validação externa é sempre uma fase posterior, realizada usando dados obtidos no campo depois de já termos, por exemplo, produzido o mapa da variável de interesse com a qual estamos trabalhando.
Como a validação cruzada reusa os dados já usados para estimar os parâmetros do modelo, ela costuma ser otimista em relação à qualidade do modelo. Por outro lado, como a validação externa usa dados coletados somente depois de termos realizado as predições espaciais, sua avaliação costuma ser mais rígida, geralmente mais rigorosa do que a validação cruzada.
# Pacotes
library(magrittr)
library(dplyr)
library(glue)
library(lattice)
library(latticeExtra)
library(georob)
library(sp)
# Sistemas de referência de coordenadas (Fonte: http://spatialreference.org/ref/epsg/)
wgs84utm22s <- sp::CRS('+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs')
sirgas2000 <- sp::CRS('+proj=longlat +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +no_defs')
# Rampas de cores
col_soil_var <- topo.colors(100)data_folder <- '../data/'
ext <- c('dbf', 'prj', 'shp', 'shx')
files <- c(glue('pedologia25.{ext}'), glue('geologia50.{ext}'))
download <- !all(files %in% list.files(data_folder))
if (download) {
url <- 'https://github.com/samuel-rosa/UFSM-SOL-843/tree/master/data/'
url <- glue('{url}{files}')
destfile <- glue('{data_folder}{files}')
for (i in 1:length(files)) {
download.file(url = url[i], destfile = destfile[i])
}
}pedologia25 <-
glue('{data_folder}pedologia25.shp') %>%
raster::shapefile(stringsAsFactors = TRUE) %>%
sp::spTransform(wgs84utm22s)
str(pedologia25, 2)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 25 obs. of 1 variable:
## ..@ polygons :List of 25
## ..@ plotOrder : int [1:25] 12 8 11 5 6 18 24 10 22 23 ...
## ..@ bbox : num [1:2, 1:2] 226655 6715750 231890 6721725
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slotgeologia50 <-
glue('{data_folder}geologia50.shp') %>%
raster::shapefile(stringsAsFactors = TRUE) %>%
sp::spTransform(wgs84utm22s)
str(geologia50, 2)## Formal class 'SpatialPolygonsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 9 obs. of 1 variable:
## ..@ polygons :List of 9
## ..@ plotOrder : int [1:9] 2 9 7 5 8 6 4 3 1
## ..@ bbox : num [1:2, 1:2] 226655 6715750 231890 6721725
## .. ..- attr(*, "dimnames")=List of 2
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slotpontos400_o <- febr::observations('ctb0003', which.cols = 'all', progress = FALSE)
pontos400_l <- febr::layers('ctb0003', which.cols = 'all', missing.data = 'keep', progress = FALSE)
id <- c('dataset_id', 'observacao_id')
pontos400 <-
merge(pontos400_o, pontos400_l, by.x = id, by.y = id) %>%
select(observacao_id, coord_x, coord_y, taxon_sibcs_2009, ca_kcl_aas, argila_, areia_)
rm(pontos400_l, pontos400_o)
sp::coordinates(pontos400) <- ~ coord_x + coord_y
sp::proj4string(pontos400) <- sirgas2000
pontos400 <- sp::spTransform(pontos400, wgs84utm22s)
pontos400$um <- sp::over(x = pontos400, y = pedologia25) %>% unlist()
pontos400$geo <- sp::over(x = pontos400, y = geologia50) %>% unlist()
pontos400in <- pontos400[!is.na(pontos400$um) & !is.na(pontos400$geo), ]lm_fit <- lm(areia_ ~ um + geo, pontos400in)
summary(lm_fit)##
## Call:
## lm(formula = areia_ ~ um + geo, data = pontos400in)
##
## Residuals:
## Min 1Q Median 3Q Max
## -482.00 -84.08 -4.46 97.10 387.72
##
## Coefficients:
## Estimate Std. Error t value
## (Intercept) 548.26 38.67 14.178
## umArgissolo Vermelho -80.77 48.05 -1.681
## umCambissolo - Neossolo -82.96 44.12 -1.881
## umNeossolo Flúvico 110.64 113.63 0.974
## umNeossolo Litólico -17.99 34.93 -0.515
## umNeossolo Litólico - Neossolo Regolítico -57.66 59.85 -0.963
## umNeossolo Regolítico -152.48 58.98 -2.585
## umNeossolo Regolítico 106.39 81.38 1.307
## umPlanossolo Háplico -87.66 49.93 -1.756
## geobotucatu 176.73 27.31 6.472
## geocaturrita 204.60 37.39 5.471
## georiolito -229.83 27.84 -8.255
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## umArgissolo Vermelho 0.0936 .
## umCambissolo - Neossolo 0.0608 .
## umNeossolo Flúvico 0.3309
## umNeossolo Litólico 0.6070
## umNeossolo Litólico - Neossolo Regolítico 0.3360
## umNeossolo Regolítico 0.0101 *
## umNeossolo Regolítico 0.1919
## umPlanossolo Háplico 0.0800 .
## geobotucatu 3.12e-10 ***
## geocaturrita 8.32e-08 ***
## georiolito 2.83e-15 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 158.6 on 364 degrees of freedom
## Multiple R-squared: 0.6388, Adjusted R-squared: 0.6278
## F-statistic: 58.51 on 11 and 364 DF, p-value: < 2.2e-16limites <- seq(0, 1500, length.out = 15)
residuals(lm_fit) %>%
georob::sample.variogram(
locations = pontos400in@coords, lag.dist.def = limites,
xy.angle.def = c(0, 22.5, 67.5, 112.5, 157.5, 180)) %>%
plot(type = "b", ylab = 'Semivariância', xlab = 'Distância de separação (km)')
Figure 2.1: Variograma direcional dos resíduos do modelo linear do conteúdo de areia na camada superficial do solo.
vario <-
residuals(lm_fit) %>%
georob::sample.variogram(
locations = pontos400in@coords, lag.dist.def = limites)
vario_fit <-
georob::fit.variogram.model(
vario, variogram.model = 'RMexp', param = c(variance = 20000, nugget = 2000, scale = 500),
weighting.method = "cressie", method = "BFGS")##
##
## variance snugget nugget scale
## Variogram parameters Inf 0e+00 Inf 0e+00
##
##
## variance snugget nugget scale
## Variogram parameters Inf 0.000000e+00 3.517747e+158 0.000000e+00
##
##
## variance snugget nugget scale
## Variogram parameters 2.069175e+241 0.000000e+00 2.239118e+34 6.482568e-82
##
##
## variance snugget nugget scale
## Variogram parameters 5.058054e+51 0.000000e+00 3.242197e+09 8.346904e-15
##
##
## variance snugget nugget scale
## Variogram parameters 6.048124e+13 0.000000e+00 3.491338e+04 2.205368e-01
##
##
## variance snugget nugget scale
## Variogram parameters 1.013092e+50 0.000000e+00 3.654724e-45 4.045095e+106
##
##
## variance snugget nugget scale
## Variogram parameters 3.181483e+13 0.000000e+00 5.772228e-07 2.877972e+23
##
##
## variance snugget nugget scale
## Variogram parameters 1.475465e+48 0.000000e+00 2.289240e+08 9.043636e-09
##
##
## variance snugget nugget scale
## Variogram parameters 1.551143e+13 0.000000e+00 8.506122e+03 3.355399e+00
##
##
## variance snugget nugget scale
## Variogram parameters 2.501513e+04 0.000000e+00 1.333848e+12 2.440327e+03
##
##
## variance snugget nugget scale
## Variogram parameters 5.930840e-73 0.000000e+00 2.112996e-03 4.443735e+36reml_fit <- georob::georob(
areia_ ~ um + geo, pontos400in, locations = ~ coord_x + coord_y, variogram.model = 'RMexp',
param = c(variance = vario_fit$variogram.object[[1]]$param[['variance']],
nugget = vario_fit$variogram.object[[1]]$param[['nugget']],
scale = vario_fit$variogram.object[[1]]$param[['scale']]),
tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(reml_fit)##
## Call:georob::georob(formula = areia_ ~ um + geo, data = pontos400in,
## locations = ~coord_x + coord_y, variogram.model = "RMexp",
## param = c(variance = vario_fit$variogram.object[[1]]$param[["variance"]],
## nugget = vario_fit$variogram.object[[1]]$param[["nugget"]],
## scale = vario_fit$variogram.object[[1]]$param[["scale"]]),
## tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
##
## Tuning constant: 1000
##
## Convergence in 13 function and 8 Jacobian/gradient evaluations
##
## Estimating equations (gradient)
##
## eta scale
## Gradient : -2.498554e-02 -1.454839e-01
##
## Maximized restricted log-likelihood: -2328.15
##
## Predicted latent variable (B):
## Min 1Q Median 3Q Max
## -401.79 -82.80 10.03 128.73 427.66
##
## Residuals (epsilon):
## Min 1Q Median 3Q Max
## -9.332793 -1.331201 -0.004369 1.248425 8.968548
##
## Standardized residuals:
## Min 1Q Median 3Q Max
## -3.440757 -0.530451 -0.002582 0.479414 3.028140
##
##
## Gaussian REML estimates
##
## Variogram: RMexp
## Estimate Lower Upper
## variance 2.710e+04 2.139e+04 34316.7
## snugget(fixed) 0.000e+00 NA NA
## nugget 3.277e+02 1.511e-01 710529.6
## scale 1.966e+02 1.342e+02 287.8
##
##
## Fixed effects coefficients:
## Estimate Std. Error t value
## (Intercept) 501.17 40.93 12.244
## umArgissolo Vermelho -73.01 49.39 -1.478
## umCambissolo - Neossolo -61.94 44.39 -1.395
## umNeossolo Flúvico 122.26 70.76 1.728
## umNeossolo Litólico -12.93 34.16 -0.378
## umNeossolo Litólico - Neossolo Regolítico -69.27 66.78 -1.037
## umNeossolo Regolítico -93.44 62.93 -1.485
## umNeossolo Regolítico 107.93 76.75 1.406
## umPlanossolo Háplico -75.13 48.81 -1.539
## geobotucatu 146.13 26.06 5.608
## geocaturrita 198.27 42.72 4.641
## georiolito -158.14 31.62 -5.002
## Pr(>|t|)
## (Intercept) < 2e-16 ***
## umArgissolo Vermelho 0.1402
## umCambissolo - Neossolo 0.1638
## umNeossolo Flúvico 0.0849 .
## umNeossolo Litólico 0.7053
## umNeossolo Litólico - Neossolo Regolítico 0.3003
## umNeossolo Regolítico 0.1385
## umNeossolo Regolítico 0.1605
## umPlanossolo Háplico 0.1246
## geobotucatu 4.06e-08 ***
## geocaturrita 4.85e-06 ***
## georiolito 8.85e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error (sqrt(nugget)): 18.1
##
## Robustness weights:
## All 376 weights are ~= 1.plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância')
lines(vario_fit, col = "red", lty = 'dashed')
lines(reml_fit, col = "blue", lty = 'dashed')
Figure 2.2: Variograma amostral (preto) dos resíduos do modelo linear do conteúdo de areia na camada superficial do solo e a função exponencial (vermelho) a ele ajustada.
grid <- sp::spsample(pedologia25, 10000, type = 'regular')
grid <-
sp::SpatialPointsDataFrame(
coords = grid@coords,
data = data.frame(
um = sp::over(grid, pedologia25) %>% unlist(),
geo = sp::over(grid, geologia50) %>% unlist()),
proj4string = grid@proj4string)
str(grid)## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
## ..@ data :'data.frame': 10005 obs. of 2 variables:
## .. ..$ um : Factor w/ 9 levels "Argissolo Bruno-Acinzentado",..: 1 1 1 1 1 9 9 9 9 1 ...
## .. ..$ geo: Factor w/ 4 levels "basalto","botucatu",..: 3 3 3 3 3 3 3 3 3 3 ...
## ..@ coords.nrs : num(0)
## ..@ coords : num [1:10005, 1:2] 229605 229649 229694 229738 229782 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:10005] "um1" "um2" "um3" "um4" ...
## .. .. ..$ : chr [1:2] "x1" "x2"
## ..@ bbox : num [1:2, 1:2] 226693 6715777 231856 6721689
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "x1" "x2"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr "+proj=utm +zone=22 +south +ellps=WGS84 +datum=WGS84 +units=m +no_defs +towgs84=0,0,0"colnames(grid@coords) <- colnames(pontos400in@coords)
pred_ponto <- predict(reml_fit, newdata = grid, control = georob::control.predict.georob(extended.output = TRUE))
sp::gridded(pred_ponto) <- TRUE
str(pred_ponto)## Formal class 'SpatialPixelsDataFrame' [package "sp"] with 7 slots
## ..@ data :'data.frame': 10005 obs. of 8 variables:
## .. ..$ pred : num [1:10005] 737 747 756 763 766 ...
## .. ..$ se : num [1:10005] 151 151 151 149 149 ...
## .. ..$ lower : num [1:10005] 441 451 461 470 474 ...
## .. ..$ upper : num [1:10005] 1033 1044 1051 1055 1057 ...
## .. ..$ trend : num [1:10005] 699 699 699 699 699 ...
## .. ..$ var.pred : num [1:10005] 5682 5599 5861 6289 6528 ...
## .. ..$ cov.pred.target: num [1:10005] 5000 4894 5138 5549 5764 ...
## .. ..$ var.target : num [1:10005] 27096 27096 27096 27096 27096 ...
## .. ..- attr(*, "variogram.object")=List of 1
## .. .. ..$ :List of 9
## .. .. .. ..$ variogram.model: chr "RMexp"
## .. .. .. ..$ param : Named num [1:4] 27096 0 328 197
## .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
## .. .. .. ..$ fit.param : Named logi [1:4] TRUE FALSE TRUE TRUE
## .. .. .. .. ..- attr(*, "names")= chr [1:4] "variance" "snugget" "nugget" "scale"
## .. .. .. ..$ isotropic : logi TRUE
## .. .. .. ..$ aniso : Named num [1:5] 1 1 90 90 0
## .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
## .. .. .. ..$ fit.aniso : Named logi [1:5] FALSE FALSE FALSE FALSE FALSE
## .. .. .. .. ..- attr(*, "names")= chr [1:5] "f1" "f2" "omega" "phi" ...
## .. .. .. ..$ sincos :List of 6
## .. .. .. .. ..$ co: num 6.12e-17
## .. .. .. .. ..$ so: num 1
## .. .. .. .. ..$ cp: num 6.12e-17
## .. .. .. .. ..$ sp: num 1
## .. .. .. .. ..$ cz: num 1
## .. .. .. .. ..$ sz: num 0
## .. .. .. ..$ rotmat : num [1:2, 1:2] 1.00 -6.12e-17 6.12e-17 1.00
## .. .. .. ..$ sclmat : Named num [1:2] 1 1
## .. .. .. .. ..- attr(*, "names")= chr [1:2] "" "f1"
## .. ..- attr(*, "psi.func")= chr "logistic"
## .. ..- attr(*, "tuning.psi")= num 1000
## .. ..- attr(*, "type")= chr "signal"
## ..@ coords.nrs : num(0)
## ..@ grid :Formal class 'GridTopology' [package "sp"] with 3 slots
## .. .. ..@ cellcentre.offset: Named num [1:2] 226693 6715777
## .. .. .. ..- attr(*, "names")= chr [1:2] "coord_x" "coord_y"
## .. .. ..@ cellsize : Named num [1:2] 44.1 44.1
## .. .. .. ..- attr(*, "names")= chr [1:2] "coord_x" "coord_y"
## .. .. ..@ cells.dim : Named int [1:2] 118 135
## .. .. .. ..- attr(*, "names")= chr [1:2] "coord_x" "coord_y"
## ..@ grid.index : int [1:10005] 15879 15880 15881 15882 15883 15884 15885 15886 15887 15756 ...
## ..@ coords : num [1:10005, 1:2] 229605 229649 229694 229738 229782 ...
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:10005] "um1" "um2" "um3" "um4" ...
## .. .. ..$ : chr [1:2] "coord_x" "coord_y"
## ..@ bbox : num [1:2, 1:2] 226671 6715754 231878 6721712
## .. ..- attr(*, "dimnames")=List of 2
## .. .. ..$ : chr [1:2] "coord_x" "coord_y"
## .. .. ..$ : chr [1:2] "min" "max"
## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
## .. .. ..@ projargs: chr NAat <- pred_ponto@data[, c("pred", "lower", "upper")] %>% range()
at <- seq(at[1], at[2], length.out = 20)
sp::spplot(pred_ponto, zcol = c("lower", "pred", "upper"), at = at, main = "prediction")
validacao <- georob::cv(reml_fit)## Warning in FUN(X[[i]], ...): for factor 'um' some levels are missing in some calibration set(s),
## respective factors are treated as offset terms in model## for factor ' um ' some levels are missing in some calibration set(s),
## respective factors are treated as offset terms in modelsummary(validacao)##
## Statistics of cross-validation prediction errors
## me mede rmse made qne msse medsse
## 14.6560 9.9873 153.0066 133.5758 145.5232 0.9352 0.3223
## crps
## 86.20421 - sum((validacao$pred$data - validacao$pred$pred)^2) / sum((validacao$pred$data - mean(validacao$pred$data))^2)## [1] 0.6527338plot(validacao)