Aerial space-time interaction model

[Tim Meehan]

Hi all,

I would like to model aerial data with something like a type IV space-time interaction. I am experimenting with a few model statements, based on chapters 7 and 8 in Blangiardo and Cameletti (and Kranski).

I have seen these models implemented in INLA in two ways. First is the Knorr-Held type IV model in Chapter 7 (7.1.2). It looks like this (with a global fixed effect and exchangeable random effect added).

form1 <- count ~ 1 + firstyr +
  f(obser_route_idx1, model="iid") +

 # space and time stuff
  f(area_idx1, model="bym", graph=g1) +
  f(year_idx1, model="rw2") +
  f(year_idx2, model="iid") +
  f(area_idx3, model="besag", graph=g1, group=year_idx3,
    control.group=list(model="rw2")) 

It runs, but it takes a long time to run, has nonsensical output and, intuitively, it has a way more detailed specification than I need (my intuition is very often wrong). It also uses an RW2 when I'd rather a RW1 or AR1.

I have also run an apparently simpler model, more similar to that in Chapter 8 (8.3.1, also on some of Elias's teaching slides). It looks like this.

form2 <- count ~ 1 + firstyr +
  f(obser_route_idx1, model="iid") +

 # just the space-time interaction
  f(area_idx1, model="bym", graph=g1, group=year_idx1,
    control.group=list(model="rw1")) 

It runs really quickly and fits the data really well (but not too well). I understand that the second model is different model, and that it is deceptively simple, as bym adds both an exchangeable and structured term per year. 

My main question is, is the second model, with the iid and spatial terms grouped per year, a legit space-time model. Do you really need all the individual space and time terms in there along with the interaction as in the first model?

My secondary question is, in the second model, are both the iid and spatially-structured parts of the bym temporally structured according to the rw1?

Thanks for any insight.

Tim

[Esmail Abdul Fattah]

Hi Tim,

Here are some answers and hopes it helps!

1. It is not easy to figure out if the first model you showed for type4 is implemented correctly. You have five different ids: obser_route_idx1, area_idx1, year_idx1, year_idx2, area_idx3. I expect area_idx1 and obser_route_idx1 to be the same. 
2. You need to impose "the sum-to-zero" constraints for the interaction types 2,3 and 4, add one more constraint for RW2 and the default constraints for besag and RW1.
3. For example, if n is the number of years and m is the number of locations/areas/counties, then the constraints needed for type 4:

interact.constr <- matrix(0, n+m-1, mn)
for(i in 1:m){
  interact.constr[i, seq(1, n * m, by=m) + (i-1)] = 1
}
for(i in 1:(n-1)){
  interact.constr[i+m, 1:m + (i-1)*m] = 1
}

if(orderRW==2)
{
  v1 = rep(1:n,each=m)
  for(i in 1:(m-1)) interact.constr = rbind(interact.constr,v1*interact.constr[i,])
  constr.inter <- list(A = interact.constr, e = rep(0, n+m-1+m-1))
 
}else{
  constr.inter <- list(A = interact.constr, e = rep(0, n+m-1))
}

and this can be added through "extraconstr" in f function: extraconstr = constr.inter.

4. The speed you have for this type 4 when using inla() depends on the number of constraints which are necessary for inference.
5. You need to make sure the implementation is right before concluding with "nonsense" results.
6. How to use RW1 vs RW2 in inla:
- f(year_idx, model="rw1", constr=TRUE) 

- f(year_idx, model="rw2", constr=TRUE, extraconstr = list(A = matrix(1:n, 1, n), e = rep(0, 1) )

7. I can't answer your question about if the model is right or wrong or can be used/tested. This you need to decide according to your data/application and the criteria you choose to decide which model is better. 

8. I recommend you to read section 2 in this paper (shorturl.at/bdENT). It includes a summary about the different interpretations of the use of structured and unstructured time and space for interaction types 1,2,3 and 4. 

Best Regards,
Esmail  

[Tim Meehan]

Hi Esmail,

Thank you for your detailed response! 

1.  The area_idx1=area_idx2=area_idx3 were copies of the 1-85 area indices. Similarly, year_idx1=year_idx2=year_idx3 were copies of the 1-50 year indices.  The obser_route_idx1 was the index for a different iid random effect not related to the space-time components. I tried to use the same indices for the spacetime components as those in Blangiardo and Cameletti (2015, Ch 7).

2, 3, 4, 5, 6. Thanks for the info on constraints! These were not fully described in Blangiardo and Cameletti , and when I added them, models finished computing much more quickly and predictions were good, similar to an interaction-only model.

7, 8. The link you sent didn't work for me. But I am guessing you were linking to Knorr-Held 2000? I have read that paper and they cautioned against including an interaction effect that has components that are not given in the main effects. This made me wonder if there was something mathmatically/computationally incorrect about an interaction-only model. Your answer makes me think that this is not the case - that an interaction only model is OK to use if that is the right model for the job. Is that correct?

Thanks again,

Tim

[Esmail Abdul Fattah]

Hi Tim,

This is the paper with different link:

https://www.researchgate.net/publication/361458159_Approximate_Bayesian_Inference_for_the_Interaction_Types_1_2_3_and_4_with_Application_in_Disease_Mapping

There are different types of interactions, more than the four types mentioned in the paper.. 

Computationally, interaction-only model should be fine inla() but it does not mean this will make sense or will be relevant to your application. You need first to test the contribution you get from the other effects (structured and unstructured). Interaction is always based on main effects. 

Best Regards,
Esmail

[Finn Lindgren]

Hi all,

another perspective on the "main effects" and "interaction effects" is
that even though in some cases there are specific theoretical results
for special types of "main&interaction effect" models, like the
intrinsic stationary models in those "Type I,II,III,IV" models (plus a
collection of linear constraints that changes the model
interpretations), in many situations we don't view the model in that
way at all.

The prime example is spatial statistics; it's extremely common to have
a model with just "intercept + stationary random field", without also
having separate "easting" and "northing" components, and definitely
without having the spatial random field have kronecker product
precision structure (as would be required for some of the theory that
is used to claim some benefits to the special "interaction" models).
Technically, such f(easting,northing) models are can be seen as
"interaction effects between the two spatial coordinates", but I don't
recall seeing a f(northing)+f(easting) model except on global scales,
but then the philosophy is more that of having a "hierarchical model
with components acting on different length scales", rather than "main
effect" and "interaction". Linear effects northing+easting is
sometimes used to capture large scale structures, but they're treated
as a whole; most models are made to be invariant to the chosen
coordinate system.

What I wanted to get across here is that it's typically more important
to carefully consider what your model means, and what each model
component is intended to accomplish. The more I see the special "type
1,2,3,4... interaction models" applied to space and space-time
coordinates, the less I see them as being appropriate modelling tools
(partly due to the extreme difficulty in interpreting their stochastic
structure; intrinsic stationary random fields are already difficult,
and the complex constraints needed to make them identifiable makes it
even worse). Note: I am a collaborator on a paper-in-progress that
does involve these models, but I feel that their benefit is more of
the "it's possible to prove some theorems for them" type rather than
them being "generally useful". I prefer proper priors with problem
specific hierarchical structure.

Finn

> --
> 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/b41738b6-0001-484a-a8bd-acc92b805a07n%40googlegroups.com.



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

[Tim Meehan]

Thanks, Esmail, for the link to the paper. It was very interesting and I look forward to trying out INLAPLUS! Also, thanks for the different perspective, Finn. The connection with other spatial statistics methods was helpful. I had also been thinking about spatiotemporal GAMs as another example of a space-time 'interaction' model that often doesn't involve 'main effects'. In my case, I am using the model as a tool for prediction in space and time, not for inference about space or time.

[lyh...@gmail.com]

Hi,

Regarding the inclusion of spatial fixed/random main effects and interaction random effects in statistical models, I may provide some hints from an econometrics perspective.

You can read the following articles:

1. Baltagi, B. H., & Li, D. (2006). Prediction in the panel data model with spatial correlation: the case of liquor. Spatial Economic Analysis, 1(2), 175-185.

2. Pesaran, M. H., & Tosetti, E. (2011). Large panels with common factors and spatial correlation. Journal of Econometrics, 161(2), 182-202.

3. Sarafidis, V., & Wansbeek, T. (2012). Cross-sectional dependence in panel data analysis. Econometric Reviews, 31(5), 483-531.

4. Baltagi, B. H., Egger, P., & Pfaffermayr, M. (2013). A generalized spatial panel data model with random effects. Econometric reviews, 32(5-6), 650-685.

Best

[Tim Meehan]

Thanks! It is interesting to see how spatiotemporal structure is perceived and treated differently in the different fields.

[John Molitor]

Great conversation on interaction models (and highly relevant to my current work).

One thought I had relates to the reason behind including main effects along with interaction terms in simple models. If one has say sex (Male, Female) and age (Young, Middle, Old) then one could model the outcome simply with an interaction variables X.inter = {MY, MM, MO, FY, FM, FO}, where M=Male, Y=Young, etc. All combinations are covered here, so why bother with main effects?

One of the arguments for inclusion of main effects inclusion in this situation is to provide a kind of "smoothing" for covariate combinations with common elements. So we would expect that the parameter for Male Young (MY) to be similar to that of Male Midde-age (MM) since both categories contain Male. A model with a main effect for Sex would provide that smoothing. An interaction only model would treat all combinations as independent, despite presence of common elements in certain covariate combinations.

However, if we build some kind of multivariate structure for the interaction (Sex versus Age here) , as has been discussed above, then it seems that there is less of an argument for including main effects given that such smoothing may already be employed via the higher level Sex-by-Age Covariance structure. Of course there will still be some reasons to include main effects for interpretation reasons, but the interaction only model would seem reasonable in this situation. (Note that this is not the case when the interaction term is modeled with a simple iid N(0,sigma^2) - here main effects might still be necessary.)

John