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

cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br


[Esteartigofoipublicadopelaprimeiravezem[Thisarticlewasfirstpublishedon Programação R – DataScience +, e gentilmente contribuiu para os R-blogueiros]. (Você pode relatar um problema sobre o conteúdo desta página aqui)


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.

Você está interessado em postar comentários? Publique no DataScience + através do seu editor RStudio.

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.

Leia Também  foreach 1.5.0 agora disponível no CRAN

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.

Leia Também  Por que R? Webinar - Redes Neurais para Modelar Interações Moleculares com Tensorflow
cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br

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 ))
Leia Também  Screenager: tempos de triagem no bioRxiv

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

Para Deixe um comentário para o autor, siga o link e comente no blog: Programação R - DataScience +.

R-bloggers.com oferece atualizações diárias por email sobre notícias e tutoriais do R sobre o aprendizado do R e muitos outros tópicos. Clique aqui se você deseja publicar ou encontrar um emprego em ciência da informação / dados.


Deseja compartilhar seu conteúdo com R-blogueiros? clique aqui se você tiver um blog ou aqui se não tiver.



cupom com desconto - o melhor site de cupom de desconto cupomcomdesconto.com.br