Regressão espacial na parte R da parte 2: INLA

Categoria

• Modelagem Avançada

Tag

• Visualização de Dados
• Programação R
• espacial

Dez meses após a parte 1 da regressão espacial em R (oh meu Deus, para onde foram esses meses?), Aqui está uma (esperançosamente aguardada) segunda parte desta vez usando o INLA, um pacote que é útil em muitas situações.

O que será isso

Existem muitos tipos diferentes de dados espaciais e todos vêm com modelos específicos. Aqui vamos nos concentrar nos chamados modelos geoestatísticos ou de referência de ponto. Basicamente, você coletou algumas informações em locais diferentes e deseja levar em consideração o fato de que os locais mais próximos têm mais probabilidade de mostrar valores semelhantes do que os locais mais separados.

Na parte 1, vimos como ajustar a regressão espacial da seguinte forma:

[ y_i sim mathcal{N}(mu_i, sigma) ]

onde (i ) indexa as diferentes linhas do seu conjunto de dados, (y ) é a variável de resposta, ( mu ) é um vetor de valores esperados e ( sigma ) é o desvio padrão residual.

[ mu_i = X_i * beta + u(s_i) ]

os valores esperados dependem de algumas covariáveis ​​que estão na matriz de design X, multiplicadas por algum coeficiente para estimar (os ( beta ) ‘) mais um efeito aleatório espacial (u ) que depende da localização ( s_i ) do ponto de dados.

A parte crucial desses modelos é estimar o efeito aleatório espacial, na maioria das vezes é modelado como uma distribuição normal multivariada:

[ u(s_i) sim mathcal{MVN}(0, Sigma) ]

e a matriz de variância-covariância ( ( Sigma )) é preenchida usando a função de correlação Matern:

[ cov(i, j) = delta * Matern(d_{ij}, kappa) ]

a covariância entre dois locais no conjunto de dados depende da distância ( (d )), do intervalo da função Materna ( ( kappa )) e da variação espacial ( ( delta )). Os parâmetros a serem estimados nesses tipos de modelos são, portanto: (i) o coeficiente dos efeitos das covariáveis ​​( ( beta )), a variação no efeito espacial ( ( delta )), a faixa do espaço efeito ( ( kappa )) e a variação residual ( ( sigma )).

Como na primeira parte, ajustaremos esse modelo a um conjunto de dados de exemplo do pacote geoR:

library(geoR)

# the ca20 dataset
# load the example dataset,
# calcium content in soil samples in Brazil
data(ca20)
# put this in a data frame
dat <- data.frame(x = ca20$coords[,1], 
                  y = ca20$coords[,2], 
                  calcium = ca20$data, 
                  elevation = ca20$covariate[,1], 
                  region = factor(ca20$covariate[,2]))

O modelo que queremos ajustar é:

[ calcium_i = intercept + elevation_i + region_i + u(s_i) ]

EM LOS ANGELES

O INLA é um pacote que permite ajustar uma ampla gama de modelos, usa a aproximação de Laplace para ajustar modelos bayesianos muito, muito mais rapidamente do que algoritmos como o MCMC. O INLA permite o ajuste de modelos geoestatísticos por meio da equação diferencial parcial estocástica (SPDE), um bom local para obter informações mais detalhadas sobre esses dois gitbooks: spde-gitbook e inla-gitbook.

A instalação de um modelo espacial no INLA requer um conjunto de etapas específicas:

1. Crie uma malha para aproximar o efeito espacial
2. Crie uma matriz de projeção para vincular as observações à malha
3. Defina a equação diferencial parcial estocástica
4. opcionalmente, especifique um conjunto de dados para derivar previsões de modelo
5. Coloque tudo junto em uma pilha de objetos
6. Ajuste o modelo

Vamos seguir estas etapas:

A malha

Para estimar o efeito aleatório espacial, o INLA usa uma malha, que pode ser facilmente definida da seguinte forma:

library(INLA)
# meshes in 2D space can be created as follow:
mesh <- inla.mesh.2d(loc = dat[,c("x", "y")],
                     max.edge = c(50, 5000))

inla.mesh.2d precisa da localização das amostras, além de algumas informações sobre a precisão da malha. Uma malha mais precisa fornecerá uma melhor estimativa do efeito espacial (a previsão será mais suave), mas isso tem o custo de tempos computacionais mais longos. Aqui nós especificamos a malha dizendo que a distância máxima entre dois nós está entre 50 e 5000 metros. Você pode explorar interativamente como definir sua malha usando a função meshbuilder do INLA.

A matriz de projeção

A matriz de projeção faz a ligação entre os dados observados e o efeito espacial estimado pelo modelo. É fácil criar:

Amat <- inla.spde.make.A(mesh, loc = as.matrix(dat[,c("x", "y")]))

Defina o SPDE

O efeito espacial no INLA está sendo estimado usando algumas máquinas matemáticas complexas denominadas equação diferencial parcial estocástica. A idéia básica é que podemos estimar um efeito espacial contínuo usando um conjunto de ponto discreto (os nós definidos na malha) e função de base, semelhante a splines de regressão. Isso facilita muito a estimativa dos campos espaciais. Agora, a parte complicada de definir o SPDE é definir os anteriores para o alcance do efeito espacial (até onde você precisa ir para que dois pontos sejam independentes, o parâmetro ( kappa )) e o variação no efeito espacial (quão variável é o campo de um ponto para o próximo, o parâmetro ( delta )).

spde <- inla.spde2.pcmatern(mesh, 
                            prior.range = c(500, 0.5),
                            prior.sigma = c(2, 0.05))

A configuração desses prévios pode ser complicada, pois terá um grande impacto no modelo ajustado. Felizmente, como o ajuste de modelos no INLA é relativamente rápido, é fácil fazer análises de sensibilidade e entender o que funciona. O anterior para o intervalo corresponde à seguinte fórmula:

[ P(kappa < kappa_0) = p ]

onde ( kappa ) é o intervalo e ( kappa_0 ) e (p ) são os dois valores a serem passados ​​para a função. Acima, dissemos: a probabilidade de o alcance do efeito espacial estar abaixo de 500 metros é de 0,5. Para o parâmetro de variação, o prior é definido de maneira semelhante:

[ P(delta > delta_0) = p ]

Assim, nas linhas acima, dissemos: a probabilidade de a variação no efeito espacial ser maior que 2 é 0,05.

A definição desses antecedentes pode ser muito confusa no começo; recomendo que você experimente alguns valores, verifique os modelos e leia em torno para obter alguma inspiração. Eu recomendo definir priors bastante fortes para o ( delta ), se o prior for muito vago (p próximo a 0,5), o modelo poderá ter algum problema, especialmente para modelos mais complexos.

Quando tivermos tudo isso, podemos colocar tudo em uma pilha:

# create the data stack
dat_stack <- inla.stack(data = list(calcium = dat$calcium), # the response variable
                        A = list(Amat, 1, 1, 1), # the projection matrix
                        effects = list(i = 1:spde$n.spde, # the spatial effect
                                       Intercept = rep(1, nrow(dat)), # the intercept
                                       elevation = dat$elevation,
                                       region = factor(dat$region)))

A parte principal aqui é que, no argumento (A ), especificamos a projeção dos diferentes efeitos. O efeito espacial é nomeado Eu (mas podemos nomear o que quisermos) e é indexado pelo número de nós da malha. Lembre-se, quanto mais fina a malha, mais fina será a estimativa do efeito espacial. Esse efeito espacial Eu está vinculado aos dados por meio da matriz de projeção (Amat ). O outro efeito está todos diretamente vinculado aos dados, portanto não há necessidade de matrizes de projeção.

Defina as previsões

No INLA, geralmente é mais fácil derivar previsões de modelo passando os novos dados que queremos usar para fazer previsões diretamente para o ajuste do modelo. Em outras palavras, precisamos definir esses novos dados antes ajuste do modelo.

Primeiro, derivaremos novos dados para prever apenas o efeito de elevação e região:

# a newdata to get the predictions
modmat <- expand.grid(elevation = seq(min(dat$elevation), 
                                      max(dat$elevation), 
                                      length.out = 10),
                      region = unique(dat$region))

# the stack for these predictions
pred_stack_fixef <- inla.stack(data = list(calcium = NA),
                               A = list(1, 1, 1),
                               effects = list(Intercept = rep(1, nrow(modmat)),
                                              elevation = modmat$elevation,
                                              region = factor(modmat$region)),
                               tag = "prd_fixef")

O principal a notar aqui é que definimos calcium=NA no argumento de dados, o modelo estimará esses valores com base nos efeitos e nos parâmetros dos modelos. Nesta pilha, também usamos uma tag prd_fixef isso nos permitirá, mais tarde, capturar mais facilmente os valores previstos.

Como ajustamos um modelo espacial, também queremos fazer previsões no espaço. Poderíamos fazer isso com base apenas no campo espacial, mas é mais interessante derivar essas previsões espaciais, considerando também as outras covariáveis ​​(elevação e região). A derivação dessas pilhas de previsão é um pouco mais envolvida, pois precisaremos de valores de elevação e região não apenas nos locais observados, mas no espaço. Isso exigirá algumas etapas não relacionadas ao INLA que podem parecer assustadoras, mas o objetivo final é basicamente definir rasters com informações de elevação e região a partir dos dados que temos:

library(raster)
library(fields) # for Tps

## first we define an empty raster to hold the coordinates of the predictions
r <- raster(xmn = min(dat$x), xmx = max(dat$x),
            ymn = min(dat$y), ymx = max(dat$y),
            resolution = 25)

## the we use thin-plate spline to derive elevation across the data
elev_m <- Tps(dat[,c("x","y")], dat$elevation)

## put this into a raster
elev <- interpolate(r, elev_m)

## for the region info we create a SpatialPolygons 
## based on the coordinates given in the ca20 object
pp <- SpatialPolygons(list(Polygons(list(Polygon(ca20[[5]])), ID = "reg1"),
                           Polygons(list(Polygon(ca20[[6]])), ID = "reg2"),
                           Polygons(list(Polygon(ca20[[7]])), ID = "reg3")))
# turn the SpatialPolygon into a raster object
region <- rasterize(pp, r)

# the new data frame with coordinates from the raster
# plus elevation and region information
newdat <- as.data.frame(xyFromCell(r, cell = 1:ncell(r)))
newdat$elevation <- values(elev)
newdat$region <- factor(values(region))
# remove NAs
newdat <- na.omit(newdat)

# create a new projection matrix for the points
Apred <- inla.spde.make.A(mesh,
                          loc = as.matrix(newdat[,c("x", "y")]))

# put this in a new stack
pred_stack_alleff <- inla.stack(data = list(calcium = NA),
                               A = list(Apred, 1, 1, 1),
                               effects = list(i = 1:spde$n.spde,
                                              Intercept = rep(1, nrow(newdat)),
                                              elevation = newdat$elevation,
                                              region = factor(newdat$region)),
                               tag = "prd_alleff")

Mais uma vez, a ideia era obter informações de elevação e região de uma varredura da região de interesse, uma vez que tivemos isso no newdat No objeto, definimos uma nova matriz de projeção, pois, desta vez, também usaremos o campo espacial para derivar previsões. Em seguida, colocamos todas essas informações em uma nova pilha com uma nova tag.

Ajuste o modelo

Estamos finamente (quase) prontos para ajustar o modelo, apenas precisamos juntar os dados observados e as duas pilhas de previsão:

# put all the stacks together
all_stack <- inla.stack(dat_stack, pred_stack_fixef,
                      pred_stack_alleff)

Agora podemos ajustar o modelo:

# fit the model
m_inla <- inla(calcium ~ -1 + Intercept + elevation + region + f(i, model = spde),
               data = inla.stack.data(all_stack),
               control.predictor = list(A = inla.stack.A(all_stack), compute = TRUE),
               quantiles = NULL)

Isso deve levar cerca de 30 segundos para ser executado. Para ajustar o modelo, o primeiro argumento é uma fórmula que descreve o modelo, precisamos usar -1 para remover a interceptação interna e montá-la separadamente. Então usamos o f() para especificar um efeito aleatório indexado por Eu seguindo o modelo SPDE definido acima. O resto é apenas passar os dados e a matriz de projeção. Observe que definimos compute=TRUE para que o modelo estime os valores de cálcio que foram dados como NAs.

Podemos obter o resumo do modelo:

summary(m_inla)
## 
## Call:
##    c("inla(formula = calcium ~ -1 + Intercept + elevation + region + 
##    ", " f(i, model = spde), data = inla.stack.data(all_stack), 
##    quantiles = NULL, ", " control.predictor = list(A = 
##    inla.stack.A(all_stack), compute = TRUE))" ) 
## Time used:
##     Pre = 3.33, Running = 10.2, Post = 0.179, Total = 13.7 
## Fixed effects:
##             mean     sd   mode kld
## Intercept 29.608 17.259 29.628   0
## elevation  1.336  1.711  1.311   0
## region1    2.353 16.212  2.350   0
## region2   10.665 16.107 10.668   0
## region3   16.587 16.211 16.589   0
## 
## Random effects:
##   Name     Model
##     i SPDE2 model
## 
## Model hyperparameters:
##                                            mean     sd    mode
## Precision for the Gaussian observations   0.031  0.009   0.027
## Range for i                             202.389 47.472 187.915
## Stdev for i                               7.450  1.045   7.152
## 
## Expected number of effective parameters(stdev): 71.96(22.09)
## Number of equivalent replicates : 2.47 
## 
## Marginal log-Likelihood:  -667.01 
## Posterior marginals for the linear predictor and
##  the fitted values are computed

Alguns elementos importantes:

• Efeitos fixos: são as estimativas do modelo para os coeficientes de interceptação, elevação e região
• Efeitos aleatórios: a definição do efeito aleatório espacial
• Hiperparâmetros do modelo: o desvio residual (dado como precisão), a faixa do efeito espacial ( ( kappa )) e o desvio no efeito espacial ( ( delta ))

Usando uma função útil de Greg Albery, podemos derivar a correlação prevista entre dois pontos dados à distância:

# devtools::install_github("https://github.com/gfalbery/ggregplot")
library(ggregplot)
library(dplyr)

INLARange(ModelList = m_inla, MaxRange = 1000, MeshList = mesh, WNames = "i") +
  labs(x = "Distance in meters") +
  theme(legend.position = "none")

E agora podemos observar primeiro as previsões do modelo apenas a partir dos efeitos fixos, calculando a média das variações espaciais:

## first we create an index to easily find these 
## prediction within the fitted model
id_fixef <- inla.stack.index(all_stack, "prd_fixef")$data

## add to modmat the prediction and their sd
modmat$calcium <- m_inla$summary.fitted.values[id_fixef, "mean"]
modmat$sd <- m_inla$summary.fitted.values[id_fixef, "sd"]

## a plot with the original data
ggplot(dat, aes(x = elevation, y = calcium)) +
  geom_ribbon(data = modmat, aes(ymin = calcium - 2 * sd,
                                 ymax = calcium + 2 * sd,
                                 fill = region),
              alpha = 0.2) +
  geom_line(data = modmat, aes(color = region)) +
  geom_point(aes(color = region)) 

Usando a tag que definimos anteriormente, podemos pegar facilmente as linhas relevantes no objeto de modelo que são armazenadas dentro do summary.fitted.values quadro de dados no objeto. Aí obtemos a média e o desvio padrão que, em seguida, plotamos junto com os dados originais.

E agora vem o mapa legal:

# again get the correct indices
id_alleff <- inla.stack.index(all_stack, "prd_alleff")$data

# now add the model predictions
newdat$pred <- m_inla$summary.fitted.values[id_alleff, "mean"]
newdat$sd <- m_inla$summary.fitted.values[id_alleff, "sd"]
# get lower and upper confidence interval
newdat$lower_ci <- with(newdat, pred - 2 * sd)
newdat$upper_ci <- with(newdat, pred + 2 * sd)

# some data wraggling
nn <- pivot_longer(newdat, cols = c("pred", "lower_ci", "upper_ci"))

ggplot(nn, aes(x=x, y=y, fill=value)) +
  geom_raster() +
  facet_wrap(~name) +
  scale_fill_continuous(type = "viridis")

Lá vemos uma das vantagens da análise bayesiana: podemos obter estimativas de variação e incerteza em torno de quaisquer quantidades estimadas do modelo.

Voila! Isso é tudo por hoje, da próxima vez continuaremos ao longo das linhas bayesianas com o Processo Gaussiano em greta, fiquem animados!

Algumas referências úteis para leitura adicional:

• Alguns tutoriais bastante úteis e muito profundos sobre o INLA (não apenas modelos espaciais): flutterbys.com.au
• O site oficial: r-inla.org
• Outro tutorial para: dados espaciais da rede

Post relacionado

• Notícias manchetes análise de texto
• Algoritmo Genético no Aprendizado de Máquina usando Python
• Pacote poderoso para aprendizado de máquina, ajuste de hiperparâmetro (pesquisa em grade e aleatória), aplicativo brilhante
• Analisando HTML e Aplicando Aprendizado de Máquina Não Supervisionado. Parte 3: Análise de Componentes Principais (PCA) usando Python
• Analisando HTML e Aplicando Aprendizado de Máquina Não Supervisionado. Parte 2: Clustering Aplicado Usando Python

var vglnk = {key: '949efb41171ac6ec1bf7f206d57e90b8'}; (função (d, t) {var s = d.createElement