Help with Spatially Varying Coefficient (SVC) model using SPDE

[Tim Meehan]

Hi all,

I am trying to construct an SVC model in INLA using Matern/SPDEs. I have successfully built a similar model using the INLA Besag approach, but am having trouble feeding variables to INLA through stacks for the SPDE approach.

I am attaching data and some markdown output because this seems the best way to describe what I'm trying to accomplish.

Thanks for any help you can offer. I hope to write a tutorial on this type of model once I get this stuff figured out. If you think this would be simpler using inlabru, I am happy to learn that syntax instead.

Best, 

Tim

[Elias T. Krainski]

Hi, 

You just need to multiply the data projector matrix with the covariate (A * covariate) as it multiplies each line of A with each element of the covariate. You can check formula 8.1 in https://becarioprecario.bitbucket.io/spde-gitbook/ch-stapp.html#sec:dynamic

best regards,

Elias

--
You received this message because you are subscribed to the Google Groups "R-inla discussion group" group.
To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.
To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e4ff87ca-4469-4670-aec9-ea205eb9d0c2n%40googlegroups.com.

[Tim Meehan]

Hi Elias,

Thanks for your initial response. After reading your message, I started to think of the model differently. Instead of calculating an SVC trend as part of the model, and then also trying to get the annual residual around the trend for annual abundance indices, I began thinking it is better to just use a spacetime model to get the abundance index across space per year, and then compute a trend from those estimates after modeling using the posteriors. 

So the model will be something like:

abundance index (per site and year) ~ zero centered IID site effect + 

                                                                    SVC for count effort (per site) + 

                                                                    space time intercept (per site and year)

I am using code from your book, Chapter 7, to try to get a model working. But I have not figured out the SVC part. I tried using the weights argument to inla.spde.make.A() function to multiply A by the spatially varying covariate. But I get errors when trying to run inla() with this A matrix. Really, it seems like I am coding the covariate wrong. It seems to me that all I need is a spatially varying coefficient, but the A matrix has a spatial and temporal dimension. It is a weird model because it has an IID intercept, a spatial covariate effect, and a spatiotemporal intercept.  For example:

## setup
library(INLA); library(splancs); library(raster)
setwd(getwd())
source("R/spde-book-functions.R"); load('data/prmesh1.RData')

## get data
data(PRborder); data(PRprec)
coords <- as.matrix(PRprec[sample(1:nrow(PRprec)), 1:2])

## make spacetime field
k <- 12 # number time steps
rho <- 0.7 # temporal cor
params <- c(variance = 1, kappa = 1)
x.k <- book.rspde(coords, range = sqrt(8) / params[2],
                  sigma = sqrt(params[1]), n = k, mesh = prmesh1,
                  return.attributes = TRUE)
x <- x.k
for (j in 2:k) {x[, j] <- rho * x[, j - 1] + sqrt(1 - rho^2) * x.k[, j]}

## make covariate
n <- nrow(coords)
ccov <- factor(sample(LETTERS[1:3], n * k, replace = TRUE))
beta <- -1:1

## make response
sd.y <- 0.01
y <- beta[unclass(ccov)] + x + rnorm(n * k, 0, sd.y)

## package data
isel <- sample(1:(n * k), n * k / 2)
dat <- data.frame(y = as.vector(y), w = ccov,
                  time = rep(1:k, each = n),
                  xcoo = rep(coords[, 1], k),
                  ycoo = rep(coords[, 2], k))[isel, ]

##########################################################
# here i change the data to be more like the data i will be modeling.

# i make y like an INTEGER COUNT and change the covariate to CONTINUOUS. i also 

# add a site ID column for the IID RANDOM SITE EFFECT in the model

dat$y <- round((dat$y^2), 0); summary(dat$y); hist(dat$y) # make it an integer count
dat$w <- as.numeric(dat$w); summary(dat$w); hist(dat$w) # make it a continuous cov
dat$logw <- log(dat$w)
dat$siteID <- as.integer(factor(paste(dat$ycoo, dat$xcoo)))
summary(dat)

##########################################################

## spde
spde <- inla.spde2.pcmatern(mesh = prmesh1,
                            prior.range = c(0.5, 0.01), # P(range < 0.05) = 0.01
                            prior.sigma = c(1, 0.01)) # P(sigma > 1) = 0.01

###########################################################

# here i make two sets of indices because i hope to have the spacetime

# effect (intercept) and also a spatially varying coefficient for the covariate (logw)
iset <- inla.spde.make.index('i', n.spde = spde$n.spde,

                             n.group = k) # space time intercept
iset2 <- inla.spde.make.index('i2', n.spde = spde$n.spde,
                             n.group = k) # SVC logw

##########################################################


#########################################################

# here I make two sets of A matrices for the same reason as above.

# the second one for the SVC term has a weights argument so that 

# the covariate is multiplied by the data projector matrix
A <- inla.spde.make.A(mesh = prmesh1, # space time intercept
                      loc = cbind(dat$xcoo, dat$ycoo), group = dat$time)
A2 <- inla.spde.make.A(mesh = prmesh1, # SVC logw
                      loc = cbind(dat$xcoo, dat$ycoo), group = dat$time,
                      weights = dat$logw)

#########################################################  

## data stack
sdat <- inla.stack(
  data = list(y = dat$y),
  A = list(A, A2, 1, 1),
  effects = list(iset, iset2,
                 logw = dat$logw,
                 siteID = dat$siteID),
  tag = 'stdata')

## priors
h.spec <- list(theta = list(prior = 'pccor1', param = c(0, 0.9)))
prec.prior <- list(prior = 'pc.prec', param = c(1, 0.01))

###########################################################

# this is the part that I can't get to work. when i try to get the spatially
# varying coefficient term in there, and then run the inla() function, i get an error:
# "There are one or more NA's in 'group' where 'idx' in f(idx,...) is not NA: # idx = 'i2'"
form <- y ~ 0 + 
  f(siteID, model="iid", constr=T) + # random site effect
  f(i2, model = spde, # the problem term that is supposed to code the SVC
    group = i.group, control.group = list(model = 'ar1', hyper = h.spec)) +
  f(i, model = spde, # space time intercept
    group = i.group, control.group = list(model = 'ar1', hyper = h.spec))

############################################################  

inla.setOption(inla.mode="experimental")
res <- inla(form,  data = inla.stack.data(sdat), family="nbinomial",
            control.predictor = list(compute = TRUE,
                                     A = inla.stack.A(sdat)),
            control.family = list(hyper = list(theta = prec.prior)),
            control.compute = list(waic = T, cpo = T,
                                   config = T),
            verbose = T,
            num.threads = 4,
            control.inla = list(strategy = 'adaptive', int.strategy = 'eb'))
summary(res)

## make plotting grid
stepsize <- 4 * 1 / 111
nxy <- round(c(diff(range(coords[, 1])),
               diff(range(coords[, 2]))) / stepsize)
projgrid <- inla.mesh.projector(
  prmesh1, xlim = range(coords[, 1]),
  ylim = range(coords[, 2]), dims = nxy)

## extract predictions
xmean <- list()
for (j in 1:k){
  xmean[[j]] <- inla.mesh.project(
    projgrid, res$summary.random$i$mean[iset$i.group == j])
}

## mask prediction grid
xy.in <- inout(projgrid$lattice$loc,
               cbind(PRborder[, 1], PRborder[, 2]))

## make a prediction raster stack and plot
year_stk <- stack()
for (j in 1:k) {
  xmean[[j]][!xy.in] <- NA
  ras <- flip(raster(t(xmean[[j]])))
  year_stk <- addLayer(year_stk, ras)
}
plot(year_stk)

## regression across stack layers to get trend per pixel
year_stk <- year_stk[[c(6:9)]]
yr <- 1:nlayers(year_stk)
X <- cbind(1, yr)
invXtX <- solve(t(X) %*% X) %*% t(X)
get_trend <- function(y) (invXtX %*% y)[2]
trend_ras <- calc(year_stk, get_trend)
plot(year_stk)
plot(trend_ras)

[Tim Meehan]

Ahh Haaa!

I got it the SVC term to run by removing the group argument from the  inla.spde.make.index() and  inla.spde.make.A() functions. Now I get the correct number of random effect estimates for all three terms!

I'll keep exploring and if I get it all worked out, I'll post a clean example.

[Elias T. Krainski]

Hi Tim,

The example in the book works for spacetime varying regression coefficient. As you want the coefficient to vary only in space, you have to simplify things: the index set, the A matrix and the formula.

i2 <- 1:spde$n.spde ### just spatial index

A2 <- inla.spde.make.A(mesh = prmesh1, # SVC logw
                      loc = cbind(dat$xcoo, dat$ycoo),

                      weights = dat$logw)

### A2 = diag(covariate)  %*% A, each line (linked to one data point) of A is multiplied by the corresponding covariate value.

f(i2, model = spde) ## the distribution for the random effect indexed by 'i2' is spatially structured 

Then this term in the linear predictor reads as

...  + u_i * covariate_i + ...

where u_i = A_i * u_0, A_i being the usual projector matrix and u_0 the random field at the mesh nodes. Thus, the 'trick' is to pre-multiply A by the covariate so you don't need to have the covariate in the formula, as  f(index, covariate, ...) as in the other cases in INLA when you have other distributions for the coefficients (besag, ar, and such).

best regards,

Elias

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/fee131ad-a997-4b3e-8ae2-65024c500b8bn%40googlegroups.com.

[Tim Meehan]

Thanks, Elias! I am finally beginning to understand how index sets, A matrices, and stacks relate to one another. Best, Tim

[luanh...@gmail.com]

This is a very helpful thread.

I am fitting a similar model, but need to explore the temporally varying effect of a covariate. In this case, how should I revise the code? I used the f(index, covariate, ...) approach for besag, ar1 before for areal data, but not for geostatistical data using SPDE.

Thanks,

[Finn Lindgren]

There’s now an SVC vignette both for inla and for inlabru.

For spde models in particular, I strongly suggest looking at the inlabru version, since that one is easier to expand to space-time:

https://inlabru-org.github.io/inlabru/articles/web/svc.html

Expanding to an ar1 “group” model works the same as in plain inla.

Finn

On 7 Dec 2022, at 04:04, luanh...@gmail.com <luanh...@gmail.com> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/e423dddf-a90a-4325-a94e-823ac14c4896n%40googlegroups.com.