make_unique_inputs fails on some matrix inputs [was: Question regarding the new interaction feature]

[joendijk]

Hi,

Since inlabru 2.5.0 it is now possible to add a fixed effect interaction. However, I am not entirely sure how to write this in R code. It would be nice to have a few more examples. Here's the model formula that you've seen before in my previous post. It includes your feedback on how write the interacting wind direction and windspeed:

form3c = resp ~ Intercept(1) + depth + temp + sal + windspeed + wind(fwinddir, model = "rw1", mapper = bru_mapper_index(8), cyclic = TRUE, replicate = windspeed) + ship(fveskey, model = "iid", n = n_veskey) + trend(coordinates, group = cyear, model = inla.spde2.matern(mesh1c), control.group = list(model="ar1")) + season(coordinates, group = week, group_mapper = bru_mapper_index(26), model = inla.spde2.matern(mesh1c), control.group = list(model="ar1"))

I would now like to add another interaction between wave height (using windspeed as proxy) and depth (to see how deep the effect of rough weather goes and whether this actually affects shallow areas directly). Since depth and windspeed are already in the form (and windspeed also as a component in "wind"), can I leave these as written, or should it all go between brackets as in the example:

~ name(~ -1 + a + b + a:b, model = "fixed")

Or can I simply put it like:

column(wind:depth, model="fixed")

Looking forward to any feedback. Many thanks again!

[finnlindgren]

I think you need a component

~ somename(~ -1 + windspeed:depth, model="fixed")

The "fixed" model simply calls MatrixModels::model.Matrix() on the provided formula (or you supply code that constructs a matrix instead), so you can check what it (hopefully) does by calling model.Matrix(theformula, data = thedata).

Note: Outside of formulas, ":" means something completely different.

[finnlindgren]

Putting it into context, you can either have

~ ... + windspeed + depth + somename(~ -1 + windspeed:depth, model="fixed")

or

~ ... + somename(~ -1 + windspeed + depth + windspeed:depth, model="fixed")

or

~ ... + somename(~ -1 + windspeed*depth, model="fixed")

Either of those should give the same model, but with slightly different internal representations (the last two give the same model matrix for the "fixed" component).

[joendijk]

Hi, thanks for the answers!

I did stumble upon some issues while testing the different forms.

~ ... + windspeed + depth + somename(~ -1 + windspeed:depth, model="fixed")

gives me the following error. Error in duplicated.default(temp, fromLast = fromLast, ...) : duplicated() applies only to vectors

While using: ~ ... + somename(~ -1 + windspeed + depth + windspeed:depth, model="fixed")

and ~ ... + somename(~ -1 + windspeed*depth, model="fixed") do work. (The model is still running, but so far it goes on without issue).

Indicating that the interacting covariates must be between the parentheses and not outside.

[finnlindgren]

The error looks like a bug; the component specification is correct and should be allowed. I'm guessing either there's a misspelled variable similar to duplicated in the code, making it pick up some function object instead, or there's a missing check on what method to use to detect duplicates. Probably the latter. Can you provide the traceback() output?

[joendijk]

Yes sure, no problem.

traceback() 12: duplicated.default(temp, fromLast = fromLast, ...) 11: unique.matrix(do.call(rbind, inp)) 10: make_uniqueinputs(inp, allow_list = TRUE) 9: add_mapper(component$main, label = component$label, lhoods = lh, env = component$env, require_indexed = FALSE) 8: add_mappers.component(components[[k]], lhoods) 7: add_mappers(components[[k]], lhoods) 6: add_mappers.component_list(object, lhoods = lhoods) 5: add_mappers(object, lhoods = lhoods) 4: component_list.list(components, lhoods) 3: component_list(components, lhoods) 2: bru_model(components, lhoods) 1: bru(form4, family = "nbinomial", data = alldat, options = list(verbose = TRUE))

[finnlindgren]

It seems that unique.matrix doesn't work for all matrices, but then I'm not sure why the 3rd version does work; the only difference I can see is that the failing matrix only has one column, and it doesn't matter if it's a MatrixModels object or a Matrix object:

library(MatrixModels)
M1 <- model.Matrix(~-1+a*b,data=data.frame(a=c(1,1,2,2),b=c(1,2,1,2)))
M2 <- model.Matrix(~-1+a:b,data=data.frame(a=c(1,1,2,2),b=c(1,2,1,2)))
M1b <- INLA::inla.as.sparse(M1)
M2b <- INLA::inla.as.sparse(M2)
str(M1)
#> Formal class 'ddenseModelMatrix' [package "MatrixModels"] with 6 slots
#>   ..@ x        : num [1:12] 1 1 2 2 1 2 1 2 1 2 ...
#>   ..@ Dim      : int [1:2] 4 3
#>   ..@ Dimnames :List of 2
#>   .. ..$ : chr [1:4] "1" "2" "3" "4"
#>   .. ..$ : chr [1:3] "a" "b" "a:b"
#>   ..@ factors  : list()
#>   ..@ assign   : int [1:3] 1 2 3
#>   ..@ contrasts: list()
str(M2)
#> Formal class 'ddenseModelMatrix' [package "MatrixModels"] with 6 slots
#>   ..@ x        : num [1:4] 1 2 2 4
#>   ..@ Dim      : int [1:2] 4 1
#>   ..@ Dimnames :List of 2
#>   .. ..$ : chr [1:4] "1" "2" "3" "4"
#>   .. ..$ : chr "a:b"
#>   ..@ factors  : list()
#>   ..@ assign   : int 1
#>   ..@ contrasts: list()
str(M1b)
#> Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
#>   ..@ i       : int [1:12] 0 1 2 3 0 1 2 3 0 1 ...
#>   ..@ j       : int [1:12] 0 0 0 0 1 1 1 1 2 2 ...
#>   ..@ Dim     : int [1:2] 4 3
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:4] "1" "2" "3" "4"
#>   .. ..$ : chr [1:3] "a" "b" "a:b"
#>   ..@ x       : num [1:12] 1 1 2 2 1 2 1 2 1 2 ...
#>   ..@ factors : list()
str(M2b)
#> Formal class 'dgTMatrix' [package "Matrix"] with 6 slots
#>   ..@ i       : int [1:4] 0 1 2 3
#>   ..@ j       : int [1:4] 0 0 0 0
#>   ..@ Dim     : int [1:2] 4 1
#>   ..@ Dimnames:List of 2
#>   .. ..$ : chr [1:4] "1" "2" "3" "4"
#>   .. ..$ : chr "a:b"
#>   ..@ x       : num [1:4] 1 2 2 4
#>   ..@ factors : list()
unique.matrix(M1)
#> 4 x 3 Matrix of class "dgeMatrix"
#>   a b a:b
#> 1 1 1   1
#> 2 1 2   2
#> 3 2 1   2
#> 4 2 2   4
unique.matrix(M2)
#> Error in duplicated.default(temp, fromLast = fromLast, ...): duplicated() applies only to vectors
unique.matrix(M1b)
#> 4 x 3 sparse Matrix of class "dgTMatrix"
#>   a b a:b
#> 1 1 1   1
#> 2 1 2   2
#> 3 2 1   2
#> 4 2 2   4
unique.matrix(M2b)
#> Error in duplicated.default(temp, fromLast = fromLast, ...): duplicated() applies only to vectors

I'll look into it. As a workaround, it may be possible to work around the issue by specifying a mapper explicitly (requires knowning how many columns the model matrix will have):

label(~-1+windspeed:depth, model="fixed", mapper = bru_mapper_matrix(c("windspeed:depth")))

I think that will bypass the code that calls unique.matrix for the component.

[finnlindgren]

I believe 1bfc320 fixes the label(~-1+a:b,model="fixed") bug. The unique.matrix method fails on single-column Matrix and ModelMatrix package matrices. By adding an extra constant column before calling the function, and then removing it again, it works.