Is it possible to interact a fixed effect with structured random effect

104 views
Skip to first unread message

Sebastian Rodriguez

unread,
May 18, 2022, 2:25:07 PM5/18/22
to R-inla discussion group
Hello everyone!

Hopefully this is a fairly basic question. I'm wondering if its possible in R-INLA to interact a fixed effect term with a structured random effect term.

My fixed effect term is simply an indicator variable, specifying the before and after period of a treatment policy. Lets call this treatment.
My structured random effect is time with a seasonal variation structure. Lets call this time.
Additionally, my model has an SPDE structured spatial random effect, which has an i.i.d. grouped structure. Therefore, all my data, along with covariates, is in an inla.stack object.

I want my formula to look something like what follows: 
f ~ 0 + b0 + treatment : f(time, model = "seasonal",  season.length = 12) + f(spatial.field, model = spde, group = spatial.field.group, control.group = list(model = 'iid'))

However, my INLA doesnt take this formula, outputting the following error:
Error in treatment:f(time, model = "seasonal", season.length = 12) :
  NA/NaN argument
In addition: Warning messages:
1: In treatment:f(time, model = "seasonal", season.length = 12) :
  numerical expression has 169632 elements: only the first used
2: In treatment:f(time, model = "seasonal", season.length = 12) :
  numerical expression has 50 elements: only the first used

If I change my formula to instead use treatment * f(time, model = "seasonal",  season.length = 12), I get the following error:
Error in model.frame.default(new.fix.formula, data.same.len, na.action = inla.na.action) :
  invalid type (list) for variable 'f(time, model = "seasonal", season.length = 12)'
In addition: Warning message:
In rt[i] <- ind :
  number of items to replace is not a multiple of replacement length


I'm wondering if the issue is arising because what I am trying to do doesnt make too much sense to begin with or if its due to the structure of the inla.stack object and how the time random effect is being treated.

Thanks in advance,
Sebastian

Finn Lindgren

unread,
May 18, 2022, 2:48:42 PM5/18/22
to Sebastian Rodriguez, R-inla discussion group
It depends a bit on what you want the “interaction” to mean in mathematical terms, but it seems to me that if your “treatment” variable is discrete, you just need to construct replicates, with f(…, replicate=treatment), assuming treatment is encoded as an integer  (1 and 2, and optionally further with 3,4,5), etc.
For continuous variables, you can use f(time, treatment, …) to define what _mathematically_ would be treatment_i*f(time_i) for observation “i “.

Finn

On 18 May 2022, at 19:25, Sebastian Rodriguez <sebastianro...@u.northwestern.edu> wrote:

Hello everyone!
--
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/df1e4cf1-a473-42e4-b04f-0365cc931923n%40googlegroups.com.

Sebastian Rodriguez

unread,
May 18, 2022, 3:14:41 PM5/18/22
to R-inla discussion group
Hi Finn,

Thanks for the quick reply! My treatment variable is encoded as a boolean variable (0 before treatment, 1 after treatment). My time variable is just an enumeration of the months in my data (in my case 1,...,48). Given your description, I would want to do the replicates option, where my new structured temporal random effect would now look like: f(time, replicate = treatment, model = "seasonal", season.length = 12).

Thank a bunch,
Sebastian

Finn Lindgren

unread,
May 18, 2022, 3:53:16 PM5/18/22
to Sebastian Rodriguez, R-inla discussion group
Yes, but you probably need to convert the 0/1 to 1/2 to work with the replicate feature, as it expects a 1-based index.

Finn

On 18 May 2022, at 20:14, Sebastian Rodriguez <sebastianro...@u.northwestern.edu> wrote:

Hi Finn,

Sebastian Rodriguez

unread,
May 18, 2022, 3:55:29 PM5/18/22
to R-inla discussion group
Hi Finn,

Thanks for the clarification. I was wondering why I was getting some strange messaged with my 0/1 encoding.

-Sebastian
Reply all
Reply to author
Forward
0 new messages