Unidade III – Modelo Linear Misto de Variação Espacial

Temas: Definição de dado geoestatístico. Especificação do modelo linear misto de variação espacial. Generalização de modelos do variograma. A função de covariância Whittle-Matérn. Estimativa de parâmetros. O método da máxima verossimilhança (MV) e da máxima verossimilhança restrita (MVR).

1 Definição e propriedades

1.1 Os dados

Dados geoestatísticos são dados observados de uma variável ambiental, \(y(\boldsymbol{s}_i)\), em um conjunto finito de locais, \(\boldsymbol{s}_i\), \(i = 1, 2, \ldots, n\), no interior de um determinado domínio, \(G \in R^d\), \(d \in (1, 2, 3)\), juntamente com dados de covariáveis ambientais espacialmente exaustivas, \(\boldsymbol{x}_j(\boldsymbol{s}_i)\), \(j = 1, 2, \ldots, p\).

1.2 O modelo

O modelo usado para os dados geoestatísticos \(y(\boldsymbol{s}_i)\) é o modelo linear misto de variação espacial. Segundo esse modelo, os dados geoestatísticos são uma realização de um campo aleatório \(Y(\boldsymbol{s}_i)\) que podem ser descritos como a combinação aditiva de efeitos fixos, efeitos aleatórios e erro aleatório independente. Esse modelo é denotado por

\[Y(\boldsymbol{s}_i) = Z(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i) = \boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta} + B(\boldsymbol{s}_i) + \varepsilon(\boldsymbol{s}_i)\]

\(Z(\boldsymbol{s}_i)\) é o chamado sinal e possui dois componentes. O primeiro, \(\boldsymbol{x}(\boldsymbol{s}_i)^\text{T}\boldsymbol{\beta}\) representa os efeitos fixos, ou seja, a tendência espacial de origem desterminística, especificamente, a relação de dependência (causa e efeito) entre as covariáveis ambientais e a variável ambiental. Aqui \(\boldsymbol{\beta}\) são os coeficientes desconhecidos da regressão linear e \(\text{T}\) denota transposição.

O segundo componente do sinal, \(B(\boldsymbol{s}_i)\), representa os efeitos aleatórios, especificamente, um campo aleatório Gaussiano estacionário não-observável, também chamado de variável latente. Esse campo aleatório Gaussiano é inteiramente descrito por sua função de média e função de covariância.

Por fim, \(\varepsilon(\boldsymbol{s}_i)\) é um erro (ou ruído) independente e identicamente distribuído (iid), descrito por uma distribuição Gaussina de probabilidade, cujo parâmetro desconhecido de escala é \(\tau\). Tradicionalmente, \(\tau^2\) é chamado de efeito pepita ou ainda de variância residual. Note que alguns autores – como Papritz (2015) – denotam esse erro aleatório independente simplesmente por \(\varepsilon_i\) a fim de destacar sua independência.

2 Função de covariância Whittle-Matérn

Pela 1ª Lei da Geografia, sabemos que todos os indivíduos de uma população são parecidos entre si, mas que os indivíduos mais próximos se parecem mais uns com os outros do que os indivíduos mais distantes. Assim, o comportamento empírico mais comumumente observado para uma estrutura de covariância espacial estacionária é aquele em que a correlação entre o sinal em dois locais, \(Z(\boldsymbol{s}_i)\) e \(Z(\boldsymbol{s}_j)\), diminui à medida que a distância de separação entre esses dois locais, \(h = ||\boldsymbol{s}_i - \boldsymbol{s}_j||\), aumenta (Diggle e Ribeiro Jr (2007) denota a distãncia de separação por \(u\)).

Quando a correlação espacial de \(Z\) reduz rapidamente com o aumento de \(h\), então diz-se que o processo espacial é áspero ou rugoso. Nesse caso, observa-se que \(Z\) muda rapidamente, ou mesmo quase abruptamente, no espaço. Do contrário, quando a correlação espacial se mantém mesmo após grandes distâncias de separação, o que é evidenciado pela pequena taxa de mudança em \(Z\) no espaço, diz-se que o processo é liso ou suave.

Uma das funções de covariância mais utilizadas para modelar esse tipo de comportamento empírico é a função de Whittle-Matérn (Matérn 1960). O principal diferencial dessa função em relação às demais é que a mesma possui um parâmetro extra \(\nu\) chamado de parâmetro de suavidade (alguns autores como Diggle e Ribeiro Jr (2007) denotam esse parâmetro por \(\kappa\)). Quando \(\nu\) é pequeno, ou seja, \(\nu \to 0\), significa que o processo espacial é mais rugoso. Do contrário, quando \(\nu\) é grande, ou seja, \(\nu \to \infty\), então o processo espacial é mais suave.

A função de covariância de Whittle-Matérn é dada por

\[\rho(h) = \frac{1}{2^{\nu-1} \varGamma(\nu)} \left(\frac{h}{\alpha}\right)^\nu K_\nu\left(\frac{h}{\alpha}\right)\]

onde \(K_\nu(\cdot)\) denota a função de Bessel modificada de segundo tipo de ordem \(\nu\) e \(\varGamma(\cdot)\) é a função gama. O parâmetro \(\alpha > 0\), também conhecido como alcance ou parâmetro de distância, é aquele que determina a taxa com que a autocorrelação espacial diminui à medida que aumenta a distância de separação \(h\). Assim, combinados os parâmetros \(\nu\) e \(\alpha\), a função de covariância de Whittle-Matérn consegue adorar uma grande variedade de formas, o que faz com que ela seja considerada um generalização de um grande número de modelos teóricos do variograma (2.1). Por exemplo, se \(\nu \to \infty\), então a função de pode ser considerada equivalente à função de covariância Gaussiana. Quando \(\nu \to 1\), ela corresponde à função elementar de covariância de Whittle. Quando \(\nu = 0.5\), então é o mesmo que a função exponencial. Se \(\alpha \to \infty\) e \(\nu > 0\) ou \(\nu \to 0\), então ela aproxima-se das funções potência e logarítmica, respectivamente.

h <- seq(0, 1, 0.01)
alpha <- 0.20
nu <- c(0.1, 0.2, 0.5, 1, 2, 4, 8)
gamma <- 1 - sapply(nu, function (x) geoR::matern(u = h, phi = alpha, kappa = x))
plot(gamma, type = 'n', xlim = range(h), xlab = 'Distância de separação', ylab = 'Semivariância')
for (i in 1:ncol(gamma)) {
  lines(x = h, y = gamma[, i])
  text(x = h[25], y = gamma[25, i], labels = nu[i], pos = 1)
}

Figure 2.1: Curvas da função de covariância de Whittle-Matérn para valores do parâmetro de suavidade entre 0.1 e 8, e valores constantes dos parâmetros de distância, 0.20, e variância, 1.

3 Método da máxima verossimilhança (restrita)

Todo modelo estocástico é especificado por um conjunto de parâmetros que descreve uma função de distribuição de probabilidade. Usar um modelo estocástico para modelar um conjunto de dados consiste, na prática, em estimar esses parâmetros a partir dos dados existentes. O conjunto de parâmetros estimados deve ser tal que descreva, dentre todos os conjuntos de parâmetros possíveis, a função de distribuição de probabilidade com a maior chance de ter gerado exatamente aqueles dados. Como uma função dos parâmetros de um modelo estocástico, estimados a partir de um conjunto de dados, é chamada de função de verossimilhança, ou simplesmente verossimilhança, o método de estimação é chamado método da máxima verossimilhança. Em outras palavras, o método da máxima verossimilhança procura pelo conjunto de parâmetros do modelo estocástico que melhor explica os dados observados.

No contexto do modelo linear misto de variação espacial, nossa intenção é, a partir do conjunto finito \(n\) de observações espaciais da variável ambiental, \(\boldsymbol{y}^\text{T} = \{y(\boldsymbol{s}_1), y(\boldsymbol{s}_2), \ldots, y(\boldsymbol{s}_n)\}\) estimar parâmetros que geralmente são de interesse científico direto, como os coeficientes da regressão linear, \(\boldsymbol{\beta}\), estimar os parâmetros que definem a estrutura de covariância de um modelo para o sinal \(Z(\boldsymbol{s})\), sendo eles \(\boldsymbol{\theta}^\text{T} = (\sigma^2, \sigma_n^2, \ldots, \alpha)\), além de \(\tau^2\), e estimar o vetor não-observável de efeitos aleatórios \(\boldsymbol{B}\). No vetor de parâmetros \(\boldsymbol{\theta}^\text{T}\), \(\ldots\) denota parâmetros adicionais do modelo, tal como o parâmetro de suavidade, \(\nu\), da função de covariância Whittle-Matérn.

Um dos problemas do método da MV(R) é que a solução encontrada pelos algorítmos otimizadores pode ser multimodal. Isso significa que mais de um conjunto de valores dos parâmetros do modelo podem retornar valores muito similares da log-verossimilhança. A maioria desses conjuntos consistem em máximos locais, quando o objetivo é encontrar o máximo global. Para evitar esse tipo de problema é preciso fornecer uma boa estimativa inicial dos parâmetros para o otimizador, o que requer um bom conhecimento da estrutura dos dados e do processo estocástico subjascente.

Contudo, mesmo fornecendo uma boa estimativa inicial, não há garantia de que os valores retornados pelo máximo global façam qualquer sentido, significando que não sejam exatamente ótimos. Mais uma vez, identificar isso requer um bom conhecimento da estrutura dos dados e do processo estocástico subjascente. A solução mais comumente adotada aqui é a avaliação de uma perfil da da função logarítmica da verossimilhança.

4 Exercício

# 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)

Vamos descarregar os arquivos do mapa pedológico para um diretório local. O nome desse diretório é definido abaixo pelo objeto data_folder. Altere o caminho para esse diretório conforme necessário. Caso você já tenha descarregado os arquivos do mapa pedológico, então data_folder deve ser o caminho para o diretório onde esses arquivos se encontram.

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])
  }
}

Agora você já pode carregar o mapa pedológico para o seu ambiente de trabalho. Para isso, use a função raster::shapefile. Note que se o sistema operacional do seu computador for o MS Windows, pode ser necessário incluir o o argumento encoding = 'UTF-8' na função raster::shapefile para garantir que os carácteres especiais usados nos nomes das unidades de mapeamento seja reconhecidos. Verifique se a estrutura de pedologia25 está conforme o esperado para um objeto do tipo SpatialPolygonsDataFrame. Note que a função sp::spTransform é usada para projetar as coordenadas original no plano cartesiano (UTM).

pedologia25 <- 
  glue('{data_folder}pedologia25.shp') %>% 
  raster::shapefile(stringsAsFactors = TRUE) %>% 
  sp::spTransform(wgs84utm22s)
col_soil_um <- terrain.colors(nlevels(pedologia25$um))
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 slot

geologia50 <- 
  glue('{data_folder}geologia50.shp') %>% 
  raster::shapefile(stringsAsFactors = TRUE) %>% 
  sp::spTransform(wgs84utm22s)
col_geo_um <- topo.colors(nlevels(geologia50$geo))
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 slot

pontos400_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), ]

sp::spplot(
  pedologia25, col.regions = col_soil_um, alpha.regions = 0.3, colorkey = FALSE) +
  latticeExtra::as.layer(
    lattice::xyplot(coord_y ~ coord_x, data = as.data.frame(pontos400in@coords), col = 'red', pch = 17))

sp::spplot(
  geologia50, col.regions = col_geo_um, alpha.regions = 0.5, colorkey = FALSE) +
  latticeExtra::as.layer(
    lattice::xyplot(coord_y ~ coord_x, data = as.data.frame(pontos400in@coords), col = 'red', pch = 17))

4.1 Análise exploratória

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-16

A Figura ?? mostra os resíduos da regressão linear ajustada aos dados do conteúdo de areia.

op <- par(mfrow = c(2, 2))
plot(lm_fit, which = 1:4)

par(op)

Conforme mostra a Figura 4.1, os resíduos da regressão linear possuem dependência espacial. Partindo da primeira classe de distância de separação, a semivariância aumenta de maneira relativamente rápida e atinge certa estabilidade a partir de 1000-1500 m. Avaliada a evolução da semivariância nas diferentes direções, principalmente nas primeiras classes de distância de separação, não há evidência clara da existência de estruturas de autocorrelação espacial dependentes da direção. Assim, é possível assumir que o processo espacial é isotrópico (do contrário seria chamado anisotrópico).

limites <- 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 4.1: Variograma direcional dos resíduos do modelo linear do conteúdo de areia na camada superficial do solo.

Assumida a isotropia do processo espacial, podemos computar o semivariograma amostral omnidirecional (ou independente da direção). Em seguida, vamos ajustar ao variograma amostral um modelo exponencial do variograma. Para isso usamos o método dos quadrados mínimos não-lineares ponderados, sendo a ponderação definida conforme o método de Cressie (1993, sec. 2.6.2). O processo de estimativa dos parâmetros do modelo exponencial do variograma é conduzido via otimização usando a função stats::optim(method = "BFGS"). Os valores produzidos pelo otimizador a cada 10 iterações são impressos no console.

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+36

O ajuste do modelo exponencial ao variograma amostral dos resíduos do conteúdo de areia é mostrado na Figura 4.2. Note como a curva ajustada passa pelo centro de massa dos doze pontos do variograma amostral, uma característica do método dos quadrados mínimos.

plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância')
lines(vario_fit, col = "red", lty = 'dashed')

Figure 4.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.

Os parâmetros ajustados do modelo exponencial do variograma são mostrados abaixo. Note que o otimizador precisou de 74 iterações para convergir, ou seja, encontrar os valores ótimos – ótimo condicional aos dados – para o modelo exponencial. Um diferencial do pacote georob – comparado ao pacote gestat – é que o mesmo retorna, além dos parâmetros estimados, os limites inferior (Lower) e superior (Upper) de confiança dessas estimativas. Note como o intervalo de confiança da estimativa da variância residual (nugget), \(\tau^2\).

summary(vario_fit)

## 
## Call:georob::fit.variogram.model(sv = vario, variogram.model = "RMexp", 
##     param = c(variance = 20000, nugget = 2000, scale = 500), 
##     weighting.method = "cressie", method = "BFGS")
## 
## Convergence in 74 function and 22 Jacobian/gradient evaluations
## 
## Residual Sum of Squares: 32.43287 
## 
## Residuals (epsilon):
##     Min      1Q  Median      3Q     Max 
## -1270.3  -695.5  -172.5   507.2  1979.9 
## 
## Variogram:  RMexp 
##                Estimate   Lower   Upper
## variance        27503.6 26653.9 28380.5
## snugget(fixed)      0.0      NA      NA
## nugget           2622.5  1918.5  3584.9
## scale             571.8   514.5   635.5

4.2 (RE)ML Gaussiano

reml_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 4.3: 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.

prof_reml_fit_scale <- georob::profilelogLik(object = reml_fit, values = data.frame(scale = seq(50, 350, by = 10)))

plot(loglik ~ scale, prof_reml_fit_scale, type = "l")
abline(v = summary(reml_fit)$param.aniso[[1]]['scale', ], lty = c("dashed", rep('dotted', 2)), col = 'red')
abline(h = reml_fit$loglik - 0.5 * qchisq(0.95, 1), lty = "dotted")

Figure 4.4: Perfil da função logarítmica da verossimilhança restrita para o parâmetro alcance (scale). A linhas verticais representam a estimativa e intervalo de confiança do parâmetro alcance. A interseção entre a linha horizontal e o perfil indica a região de 95% de confiança para o parâmetro alcance segundo o teste da razão de verossimilhança.

op <- par(mfrow = c(1, 2))
plot(variance ~ scale, prof_reml_fit_scale, ylim = c(0, max(variance)), type = "l")
plot(nugget ~ scale, prof_reml_fit_scale, ylim = c(0, max(nugget)), type = "l")

Figure 4.5: Relação das estimativas da variância parcial (variance) e da variância residual (nugget) com as estimativas do alcance (scale).

par(op)

Mesmo que nós não tenhamos observações replicadas de localizações idênticas, é possível fazer suposições sobre a natureza da variação espacialmente não-correlacionada (erro de medida versus variação espacial de pequena escala) No pacote georob, a variãncia devida aos erros de medida é modelada pelo parâmetro nugget. Enquanto isso, a variância devida à variação espacial em pequena escala processo do sinal é modelada pelo parâmetro snugget.

Se nós soubermos a variância do erro de medida a priori, então é possível fazer o seguinte. Primeiro, procedemos na análise dos dados conforme feito até aqui. A estimativa do nugget será, nesse caso, a soma do nugget e snugget, o último tendo sido mantido fixo (snugget = 0) durante a modelagem. Nesse caso, nugget é 2622.5268441. Vamos supor que sabemos a priori que a variância do erro de medida é igual à 1/4 desse valor, ou seja, 655.631711. Assim, a variância do componente espacialmente não-correlacionado do sinal será igual a 1966.8951331.

Feitas essas suposições, precisamos ajustar o modelo novamente, agora mantendo ambos nugget e snugget fixos.

reml_fit_error <- 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']] * 0.25,
            snugget = vario_fit$variogram.object[[1]]$param[['nugget']] * 0.75,
            scale = vario_fit$variogram.object[[1]]$param[['scale']]),
  fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE),
  tuning.psi = 1000, control = georob::control.georob(initial.fixef = 'lm'))
summary(reml_fit_error)

## 
## 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"]] * 
##             0.25, snugget = vario_fit$variogram.object[[1]]$param[["nugget"]] * 
##             0.75, scale = vario_fit$variogram.object[[1]]$param[["scale"]]), 
##     fit.param = georob::default.fit.param(nugget = FALSE, snugget = FALSE), 
##     tuning.psi = 1000, control = georob::control.georob(initial.fixef = "lm"))
## 
## Tuning constant:  1000 
## 
## Convergence in 9 function and 6 Jacobian/gradient evaluations
## 
## Estimating equations (gradient)
## 
##                                   xi         scale
##   Gradient           : -1.318966e-02 -8.548980e-03
## 
## Maximized restricted log-likelihood: -2328.991 
## 
## Predicted latent variable (B):
##     Min      1Q  Median      3Q     Max 
## -346.68  -74.40   11.67  126.13  415.07 
## 
## Residuals (epsilon):
##      Min       1Q   Median       3Q      Max 
## -63.5387  -9.8631   0.1225  10.1568  68.1847 
## 
## Standardized residuals:
##      Min       1Q   Median       3Q      Max 
## -3.53802 -0.50768  0.00637  0.47715  3.07858 
## 
## 
## Gaussian REML estimates
## 
## Variogram:  RMexp 
##                Estimate   Lower Upper
## variance        25060.7 19314.5 32516
## snugget(fixed)      0.0      NA    NA
## nugget(fixed)    2622.5      NA    NA
## scale             229.0   154.8   339
## 
## 
## Fixed effects coefficients:
##                                            Estimate Std. Error t value
## (Intercept)                                 495.914     41.917  11.831
## umArgissolo Vermelho                        -67.548     49.422  -1.367
## umCambissolo - Neossolo                     -57.250     44.644  -1.282
## umNeossolo Flúvico                          137.354     76.834   1.788
## umNeossolo Litólico                          -8.472     34.267  -0.247
## umNeossolo Litólico - Neossolo Regolí­tico  -65.331     67.357  -0.970
## umNeossolo Regolítico                       -92.496     64.447  -1.435
## umNeossolo Regolí­tico                      106.673     76.067   1.402
## umPlanossolo Háplico                        -72.355     49.016  -1.476
## geobotucatu                                 142.596     25.925   5.500
## geocaturrita                                194.218     43.344   4.481
## georiolito                                 -153.613     32.202  -4.770
##                                            Pr(>|t|)    
## (Intercept)                                 < 2e-16 ***
## umArgissolo Vermelho                         0.1725    
## umCambissolo - Neossolo                      0.2005    
## umNeossolo Flúvico                           0.0747 .  
## umNeossolo Litólico                          0.8049    
## umNeossolo Litólico - Neossolo Regolí­tico   0.3327    
## umNeossolo Regolítico                        0.1521    
## umNeossolo Regolí­tico                       0.1617    
## umPlanossolo Háplico                         0.1408    
## geobotucatu                                7.16e-08 ***
## geocaturrita                               9.97e-06 ***
## georiolito                                 2.67e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error (sqrt(nugget)): 51.21 
## 
## Robustness weights: 
##  All 376 weights are ~= 1.

plot(vario, type = "b", xlab = 'Distância de separação (m)', ylab = 'Semivariância', lty = 'dashed')
lines(vario_fit, col = "magenta")
lines(reml_fit, col = "navyblue")
lines(reml_fit_error, col = "orange")

Figure 4.6: 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.

Referências

Diggle, Peter J., e Paulo Justiniano Ribeiro Jr. 2007. Model-based Geostatistics. 1 ed. New York: Springer. http://www.springer.com/earth+sciences+and+geography/book/978-0-387-32907-9.

Matérn, Bertil. 1960. «Spatial variation: stochastic models and their application to some problems in forest surveys and other sampling investigations». Tese de doutoramento, Stockholm: Statens Skogsforskningsinstitut. http://pub.epsilon.slu.se/10033/1/medd_statens_skogsforskningsinst_049_05.pdf.

Papritz, Andreas. 2015. georob: Robust Geostatistical Analysis of Spatial Data. https://CRAN.R-project.org/package=georob.