Data stack with more than one A matrix

[pcstiles]

Hello,

I'm trying out some double-spde models to see how it would work with INLA (or if it's even possible to get these models to work). I essentially have two SPDE functions. One comes from a more "traditional" spatial mesh based on [x, y] locations and the other comes from a mesh of the point locations of two non-spatial variables, attempting to replicate a GAM smoothing function over two variables. My first question is if such a model is possible? My second question is, if it's possible, how would you define the data stack? I've tried the code below, but get the error that the arguments have differing numbers of rows, which I think comes from the A list. I've also included the code I'd use to run this model, although since I'm blocked on the data stack I don't know if this would generate any errors yet.

I appreciate any guidance you have here!

Best,

Pascale

stack.est = inla.stack(tag="est", data=list(y=dat2$y, offs=dat2$offset),

                       A=list(A1, A2),

                       effects=list(c(s_index1, s_index2, list(Intercept=1)), list(trap_type=dat2$trap_type2, yweek=dat2$yweek)))

form.inla = y~ -1 + Intercept + trap_type +

  f(s.field1, model=spde1) +

  f(s.field2, model=spde2) +

  f(yweek, model="rw2", constr=F)

m.inla = inla(form.inla,

              offset=log(inla.stack.data(stack.est)$offs),

              data=inla.stack.data(stack.est),

              family="nbinomial",

              control.predictor=list(A=inla.stack.A(stack.est), compute=T, link=1),

              control.compute=list(dic=T))

[Finn Lindgren]

Hi,

It’s not clear from the part of the code you included, but I suspect that your s_index1 and s_index2 are for the two different mesh/spde models, and that the part with week is something else? In that case, you need three A matrices and three corresponding index/effect lists.

This is a bit easier to setup with the inlabru interface, since that package handles all of the inla.stack() code internally instead.

Finn

On 14 Sep 2021, at 18:41, pcstiles <pascale...@gmail.com> wrote:

Hello,

--
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/8209d8d3-5468-44bf-b939-c76813e2e903n%40googlegroups.com.

[Pascale Stiles]

Hi Finn,

Thanks for your quick reply! You’re right, s_index1 and s_index2 are from two different spde models. The yweek part is for a separate (non-spde) function, which hasn’t needed an A matrix in the past when I’ve run it with only one spde model. The inlabru package might be interesting, but I’ve never used it before. Do you have a good starting point to recommend for that?

Best,

Pascale

[Finn Lindgren]

You must have had some A-matrix, even if it was something basic like “1” (I.e. identity); you probably need A=list(A1, A2,1), and separate you first effect list into two separate lists. The intercept can likely go in either of the three, but it may be clearest to put it with the other covariates.

For inlabru examples, got to https://inlabru-org.github.io/inlabru/articles/ where the last in the list is a basic 2D SPDE model example.

Finn

On 14 Sep 2021, at 19:20, Pascale Stiles <pascale...@gmail.com> wrote:



[Pascale Stiles]

Ah yes, I had an identity matrix in there! Thanks for the link. I’ll look into that, and I’ll test the stack as you suggested as well.

Best,

Pascale

[Finn Lindgren]

Note: where the inlabru examples use the keyword “coordinates”, use an explicit matrix construction of the relevant data columns instead, like

  cbind(xvar, yvar)

(The examples were constructed for models with only one spatial model components of this type, where some shortcuts could be taken.)

Finn

On 14 Sep 2021, at 19:40, Pascale Stiles <pascale...@gmail.com> wrote:



[Pascale Stiles]

Thanks for the tip. Looking at the vignette, it's not clear how you would implement a RW-type function (as in f(var, model="rw2") with INLA)? The ability to predict makes inlabru really interesting compared to inla, but I just want to make sure it can handle what I need the model to do.

Pascale

[Pascale Stiles]

Hi again,

I was testing a simple spatial-only model with inlabru just to see how it compares to INLA predictions at the same point locations and got this error (screenshot attached). Is this a memory issue? My code for both the inla and inlabru models and predictions is below for reference.

Thanks!

Pascale

mesh = inla.mesh.2d(loc=cbind(coords$x, coords$y),
                    loc.domain=borders[,c("long","lat")],
                    max.edge=c(7500, 20000),
                    cutoff=3000)

spde1 = inla.spde2.pcmatern(mesh=mesh, alpha=2,
                            prior.range=c(7000, 0.01),
                            prior.sigma=c(2.5, 0.01))

A1 = inla.spde.make.A(mesh=mesh,
                      loc=cbind(dat2$x, dat2$y))

s_index1 = inla.spde.make.index(name="spatial.field",
                                n.spde=spde1$n.spde)

stack1 = inla.stack(data=list(y=dat2$y, offs=dat2$offset),
                    A=list(A1, 1), # projector matrix
                    effects=list(c(s_index1, list(Intercept=1)), list(trap_type2=dat2$trap_type2)),
                    tag="est")

pred = centroids2@data
pred$offset = 1
pred$trap_type2 = 0

A2 = inla.spde.make.A(mesh=mesh,
                      loc=cbind(pred$x, pred$y))

stack2 = inla.stack(data=list(y=NA, offs=pred$offset),
                    A=list(A2, 1), # projector matrix
                    effects=list(c(s_index1, list(Intercept=1)), list(trap_type2=pred$trap_type2)),
                    tag="pred")

stack = inla.stack(stack1, stack2)

form = y ~ -1 + Intercept + trap_type2 +
  f(spatial.field, model=spde1)

m.inla = inla(form,
              offset=log(inla.stack.data(stack)$offs),
              data=inla.stack.data(stack),
              family="nbinomial",
              control.predictor=list(A=inla.stack.A(stack), compute=T, link=1))

i_pred = inla.stack.index(stack, "pred")$data
marginals = m.inla$marginals.fitted.values[i_pred]

pred$pred.inla = unlist(lapply(marginals, function(x) inla.qmarginal(0.5, x)))

comp = y ~ field(map=cbind(dat2$x, dat2$y), model=spde1) + Intercept + trap_type2

m.bru = bru(comp,
            family="nbinomial",
            data=dat2,
            formula=~field + Intercept + trap_type2,
            options=list(E=dat2$offset))

pred$pred.bru = predict(m.bru, pred[,c("x","y","trap_type2")], formula=~exp(field + Intercept + trap_type2))

[Finn Lindgren]

Hi, for a rw2 model, the inlabru version of f(var, model="rw2") is var(var, model="rw2") (unless you give the component a different name, which is optional)

Finn

--

Finn Lindgren
email: finn.l...@gmail.com

[Finn Lindgren]

I'm between many meetings right now, but I'm guessing it's because you need to setup the components to get their inputs from the data objects you provide, instead of hardcode it from an external object, so whether you have

cbind(coord$x,coord$y)

you probably meant

cbind(x,y)

and similarly for the other component. Both of the inputs need to be constructible from the variables in the one data object you give to the bru() call.

Finn

[Pascale Stiles]

Thanks, that worked! For the var(var, model=“rw2”) syntax, would you specify hyperpriors in the same way?

Pascale

[Finn Lindgren]

Yes. Most arguments are simply passed on to the inla f() function.

Finn

On 16 Sep 2021, at 06:24, Pascale Stiles <pascale...@gmail.com> wrote:



[Pascale Stiles]

That's what I suspected, but I thought I should double check. I do have one last question about inlabru for the moment. I know with "regular" inla you can implement functions from the splines package within the model formula - is this possible with inlabru? 

Thank you so much for all your guidance!

Pascale

[Finn Lindgren]

Hi,

I'd need a concrete code example to answer that; most things in inla work in inlabru, but the details in how to specify the model may be different.

Finn

[Pascale Stiles]

Hi Finn,

As a very simple example, the code would be something like m=inla(formula=y~x1+ns(x2, df=n), data=dat). The inlabru syntax doesn’t seem to allow you to use a function like ns() in this case. Does that make sense?

Pascale

[Finn Lindgren]

I wasn’t actually aware of this working in inla either; I assume that the ns() function returns a matrix that then gets used as a fixed matrix?

If that’s the case, there will be ways of specifying such models in inlabru, e.g. by calling ns directly in the predictor expression and multiplying it with a vector of latent coefficients.

Finn

[INLA help]

Not me either.  I guess ns() return a matrix and that is allowed

—

Haavard Rue

HelpDesk 

help@r-inla. org

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/0979B91F-EA37-42ED-98D3-4422356526C3%40gmail.com.

[Finn Lindgren]

Ok, that makes sense.

Currently, inlabru only supports matrix covariate input in the predictor expressions and not in the component definitions themselves, but I was thinking of how to add support for it just this week; when that’s done it will likely work like this:

 compname(ns(x1,df=n)),

just as for regular vector covariates. For precomputed named matrices, e.g. data=list(X=ns(x1,df=n)), the component definition ~ X would then also work, since those are converted to the more explicit ~ X(X, model=“linear”)

The main caveat is that I’m not sure if that equivalence actually holds in the inla f() definitions or not; I’ll discover that in testing, and solve it differently if needed; it’s a straightforward model, so it’s just the internal implementation details that may be different than expected.

Finn

On 17 Sep 2021, at 21:53, INLA help <he...@r-inla.org> wrote:



[Pascale Stiles]

When I run it like that, I'm getting this error:

data.bru = list(y=y, x=x, z=z, Intercept=rep(1, length(y)), X=ns(x, df=n))
comp = y ~ X(X, model="linear") + z + Intercept
result = bru(comp, data=data.bru, formula=~X + z + Intercept)
Error in Matrix::sparseMatrix(i = which(ok), j = rep(1, sum(ok)), x = input[ok],  :
  all(dims >= dims.min) is not TRUE

Do the Intercept and z terms also need to be coerced into matrices with the same dimensions?

FYI, I found that inla syntax in Chapter 9 of this book: https://becarioprecario.bitbucket.io/inla-gitbook/ch-smoothing.html

Pascale

[Finn Lindgren]

I was describing what it will likely look like when it’s implemented, not what is currently available.

Currently, you can do something like

components = ~ Xbeta(1:n,model=“factor_full”)

formula = y ~ ns(x1,df=n) %*% Xbeta

Finn

On 18 Sep 2021, at 09:23, Pascale Stiles <pascale...@gmail.com> wrote:



[Finn Lindgren]

Note: this likely won’t work with the current inlabru CRAN release, since list() data input etc was added to the devel version more recently.

Finn

On 18 Sep 2021, at 09:39, Finn Lindgren <finn.l...@gmail.com> wrote:

I was describing what it will likely look like when it’s implemented, not what is currently available.

[Pascale Stiles]

Ah ok, I misunderstood you. I've just tried that and got a non-conformable arguments error. The matrix multiplier works without issue outside the formula though. It's not a problem if it isn't possible with this version of inlabru since it does run with INLA.

Pascale

[Finn Lindgren]

Hi,

for some reason the plain "Xbeta" use in the predictor doesn't work as intended, but one can access the latent vector directly, which does work.

One important consideration is also that it's generally unsafe to have completely automated know selections, since if one does that in the predictor expression, one gets different values for estimation and for prediction, so it's much safer to be explicit (can be data dependent, but needs to be controlled to ensure the same values are used in all the calculations).

See the following example. I'll look into why the direct latent vector access is needed, and also look into other ways of specifying these types of models; they're essentially the same as the 1d-mesh based spde models but with independent priors; those are handled automatically via the main_mapper, and one can make user-defined mapper objects, which could be used to define the "x-to-basis function" mapping for these splines, just as it's done for the B-splines in the 1D mesh models.

library(INLA)
library(inlabru)

# Random locations on (0, 20)
data <- data.frame(x1 = runif(100, min = 0, max = 20))
# A linear predictor plus noise
data$y <- data$x1 + rnorm(100)

fit <- bru(
  components = ~ -1 + Xbeta(1:4, model = "factor_full"),
  formula =
    y ~
    splines::ns(x1, knots = c(5, 10, 15), Boundary.knots = c(0, 20)) %*%
    Xbeta_latent,
  allow_latent = TRUE,
  data = data
)

summary(fit)
fit$summary.random$Xbeta

[Finn Lindgren]

Found the issue with the Xbeta vs beta_latent. The code below is an alternative that sets allow_combine=TRUE, which does allow the use of Xbeta.
But personally I think the version with direct use of Xbeta_latent is clearer.

library(INLA)
library(inlabru)

# Random locations on (0, 20)
data <- data.frame(x1 = runif(100, min = 0, max = 20))
# A linear predictor plus noise
data$y <- data$x1 + rnorm(100)

fit <- bru(
  components = ~ -1 + Xbeta(1:4, model = "factor_full"),
  formula =
    y ~
    splines::ns(x1, knots = c(5, 10, 15), Boundary.knots = c(0, 20)) %*%

    Xbeta,
  allow_combine = TRUE,

  data = data
)

summary(fit)
fit$summary.random$Xbeta

[Finn Lindgren]

...and my understanding of how the splines::ns() function defines/uses the knots wasn't correct, so please investigate that before using, so you get the correct definitions!

The basis function behaviour I get doesn't match my expectations.

Finn

[Pascale Stiles]

Hi Finn,

Thank you so much for the update! I tried fitting a model as you described here and it works

Best,

Pascale