Defining an INLA stack with a formula interaction

[Stephen McNeill]

I am using INLA for a 2D SPDE example with several thousand (x,y) coordinates. It works, for a certain combination of predictors, which convinced me that I thought I understood how inla.stack() works. Now I have made a change to the formula I use that gives me an error in inla.stack() that has me confused. Clearly, I am doing something wrong, but I can't see what it is. I have distilled the example down to a simple reproducible example below.

>>>>> Beginning of pasted code

# First, define some toy data.

library("INLA")
library("sf")

# Build dummy data (content irrelevant)
N <- 100
set.seed(1234)
x1 <- rnorm(N,mean=2,sd=1)
x2 <- rnorm(N,mean=3,sd=1)
east <- sample(1:N,size=N)
north <- sample(1:N,size=N)
y <- 2.0*x1 + 3.0*x2 + rnorm(N,mean=0,sd=0.5)
df <- data.frame(x1=x1,x2=x2,y=y,east=east,north=north)
points_spdf <- as_Spatial(st_as_sf(df,coords=c("east","north")))
boundary_sp <- as(st_polygon(list(matrix(c(0,100,100,0,0,0,0,100,100,0),ncol=2,byrow=FALSE))),"Spatial")

# Mesh
MaxEdge <- 5
mesh <- inla.mesh.2d(points_spdf,
                     boundary=boundary_sp,
                     max.edge=c(1,5)*MaxEdge,
                     cutoff=MaxEdge/2,
                     offset=c(1,5))
plot(mesh)
points(df$east,df$north,pch=21,bg=2)

# Define projector matrices
A1 <- inla.spde.make.A(mesh,loc=points_spdf)

# Define the SPDE. Entirely fictitious
spde1 <- inla.spde2.pcmatern(mesh,
                             prior.range = c(10, 0.05),
                             prior.sigma = c(0.5, 0.05))

# Define the spatial field
w1.index <- inla.spde.make.index('w',n.spde=spde1$n.spde)

#--------------------------------------------------------------------------
# The following block works
#--------------------------------------------------------------------------

## Define the predictors
Xm <- model.matrix(~ x1 + x2,data=df)
Xm <- as.data.frame(Xm)
Xm_nointercept <- dplyr::select(Xm,-"(Intercept)") # Drop the intercept

## Build the stack
Nsub <- nrow(df)
StackFit <- inla.stack(
  tag = "Fit",
  data = list(y=df$y),
  A = list(1,1,A1),                  
  effects = list(  
    Intercept = rep(1,Nsub),    # Manual intercept
    Xm        = Xm_nointercept, # Covariates without the intercept (model matrix)
    w         = w1.index))      # Weights

##    Specify the model formula in terms of the response, covariates, and spatial correlated term
str <- paste0("y ~ -1 + Intercept + `",
              paste(colnames(Xm_nointercept),collapse="`+`"),
              "`",
              "+f(w, model = spde1)")
fm <- as.formula(str)

## Run the model in INLA - this works
m_inla <- inla(formula = fm,
               family = "Gamma",
               data=inla.stack.data(StackFit),
               control.compute = list(dic = TRUE, waic = TRUE),
               control.predictor = list(A = inla.stack.A(StackFit)),
               verbose = FALSE, debug = FALSE)
#--------------------------------------------------------------------------
# The following block does not work
#--------------------------------------------------------------------------

## Define the predictors
Xm <- model.matrix(~ x1*x2,data=df) # Note the interaction
Xm <- as.data.frame(Xm)
Xm_nointercept <- dplyr::select(Xm,-"(Intercept)")

## Build the stack
Nsub <- nrow(df)
StackFit <- inla.stack(
  tag = "Fit",
  data = list(y=df$y),
  A = list(1,1,A1),                  
  effects = list(  
    Intercept = rep(1,Nsub),    # Manual intercept
    Xm        = Xm_nointercept, # Covariates without the intercept (model matrix)
                                # Random effects (there are none)
    w         = w1.index))      # Weights

# Fails with the following message:
#
#   Error in match.names(clabs, names(xi)) :
#     names do not match previous names

> End of pasted code

In summary, if I exclude an interaction in my formula, I can generate a stack as expected, and the INLA model fits too. I have progressively added predictors and the model has improved with each addition. As I now suspect an interaction, I have added that interaction between two continuous explanatory variables. The inla.stack() call fails.

I have clearly done something wrong, but can't see it. It is clearly associated with the additional term in the model matric ("x1:x2" in the example above), but that looks ok to me.

I would be grateful for any pointers.

Stephen

[Elias T. Krainski]

Due to the nature of your objects, you can fix this using

effects = c(list(Intercept=rep(1,Nsub)), Xm_nointercept, w1.index),

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/eba55e2e-521b-406a-8fe5-22e1940ed8cfn%40googlegroups.com.

[tassiet...@gmail.com]

Hi Elias,

Ae you able to expand on the above explanation? I have a more complex model with multiple fixed effects and interactions. When I try your recommendation above I get an error saying :

Error in (function (data, A, effects, tag = "", compress = TRUE, remove.unused = TRUE) : length(A)=3 should be equal to length(effects)=22

So it seems that the specification for the effects above has created a list element for each of my covariate effects which are in my covariate data frame passed in the manner that   Xm_nointercept is above.

Why does the inla.stack function treat interaction effects differently and what is the expected format when the objects I want to pass are: (i) a user defined intercept, (ii) a covariate data frame that includes interactions, and (iii) a vector referencing the spde made using inla.spde.make.index..

Thanks in advance

Nick

[Finn Lindgren]

From what I recall, the middle entry also needs to be a list, with list(Xm_nointercept)

Each A matrix (you have three) should be accompanied by an effect list. The index variable is already a list, so it doesn’t need wrapping.  Versions where lists aren’t used take advantage of a simplified shortcut syntax, and are best to avoid.

Finn

On 20 Jul 2022, at 03:43, tassiet...@gmail.com <tassiet...@gmail.com> wrote:

Hi Elias,

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/2485059c-8ade-4232-abe2-ff12cb996bf9n%40googlegroups.com.

[tassiet...@gmail.com]

Thanks Finn,

I tried your suggestion, but am now back to the original error message " Error in match.names(clabs, names(xi)) : names do not match previous names"

I am used to creating the stack in the following kind of format:

  stack.1 <- inla.stack(data=list(z=dat$Count),
                        A=list(1, 1, A),
                        effects=list(
                          intercept = rep(1, nrow(dat)),
                          X = X,
                          i = i.index
                        ))

Where X is my covariate matrix created through specifying the model matrix for my data.

This creates a list with 3 elements, which usually works. I'm curious as to why this falls over when there is an interaction and what Elias means by " Due to the nature of your objects" above.

TIA 

Nick
  

[Finn Lindgren]

The X=X should then be list(X=X) instead.

Also check that the matrix column names are “safe”, I.e. don’t contain blanks or special characters.

Finn

On 20 Jul 2022, at 07:13, tassiet...@gmail.com <tassiet...@gmail.com> wrote:

Thanks Finn,

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/cd8e67f8-6b51-4b1e-80d7-257f21159dcbn%40googlegroups.com.

[tassiet...@gmail.com]

Thanks for the quick reply.

Looks like the issue is special characters in the covariate matrix. I substituted ":" with "_" and  I can now create the stack.

As ":" is the default syntax used for interactions when creating a model matrix this seems like it might be a common issue for others out there.

I note that I'll have to change the naming of my covariates that I pass in my formula to the call to inla as well.

Nick

[Finn Lindgren]

Hi,

Yes, this is a common issue with the R variable naming system; it
can't always propagate special names through all constructions.

However, in the inlabru package interface we've introduced a feature
to allow this to be done a bit easier. The following toy example
assumes that you have variables y and x1, and it then adds an x2
factor (just for code syntax illustration) and makes an interaction
model between x1 and x2:

mydata <- cbind(mydata, x2 = sample(x = factor(c("A", "B")), size =
nrow(mydata), replace = TRUE))
mycomp1 <- y ~ Intercept(1) + mix(~ -1 + x1:x2, model = "fixed")
fit1 <- bru(mycomp1,
family = "normal",
data = mydata
)

The "mix" component for technical reasons is treated as an iid "random
effect", but with precision fixed to the same precision as used for
other "fixed effects".
By giving the formula as input to the component, it feeds it to
MatrixModels::model.Matrix() to construct the model matrix for the
interaction component. You can also directly give it a pre-computed
matrix, but then later calls to predict() will be a bit more
complicated to get right.

Finn

> To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/7761439f-a407-4297-b97e-bec45666670cn%40googlegroups.com.



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