inlabru-org / inlabru

inlabru
https://inlabru-org.github.io/inlabru/
88 stars 21 forks source link

NAs in factor error message when there are no NAs #27

Closed ASeatonSpatial closed 6 years ago

ASeatonSpatial commented 6 years ago

Following error message encountered while working through 06_2d_lgcp.r from 5-day workshop.

Error in (function (data, model, stackmaker, n = 10, result = NULL, family,  : 
  INLA returned message: With control.fixed = list(expand.factor.strategy='model.matrix'), then NA's in factor are not allowd. Please use strategy 'inla' instead.

However, there are no NAs:

> sum(is.na(gcov$vegetation$vegetation))
[1] 0
> cov <- over(pixels(mesh), gcov$vegetation)
> sum(is.na(cov))
[1] 0

I attempted model fit with only the vegetation factor covariate using:

cmp1 = coordinates ~ veg(map = gcov$vegetation, model = "factor") - Intercept
fit1 = lgcp(cmp1, nests, samplers = boundary, domain = list(coordinates = mesh))

Attempted to change the strategy using lgcp(..., options = ) but had no success.

fbachl commented 6 years ago

As far as I can tell this problem only occurs on some (not all) Windows machines. Running bru/lgcp using the parameter options = list(control.fixed = list(expand.factor.strategy='inla')) solved the issue on at least one machine. Please let me know if that does not work for you. It is hard for me to properly address this issue until I'm actually sitting in front of a computer where this happens.

finnlindgren commented 6 years ago

Almost every time inla.stack is used in combination with factors, this error message should be expected, as inla.stack by design adds NAs in required places (this is the whole purpose of inla.stack). The problem is that factors are ambiguous, and INLA can't always know which factor handling strategy. (Is there an intercept in the model? Drop the first level from each factor? Analyse the factors for rank deficiency and drop levels based on that?) In regular inla, one should therefore nearly always use expand.factor.strategy="inla", and leave it to the user to make sure it's not a rank-deficient model. This is clearly not ideal. Without expand.factor.strategy="inla", and the automatic formula intercept-handling, model.matrix() is used by R-INLA, but with special Intercept terms (both in R-INLA and inlabru I think) this doesn't work, as it doesn't know that Intercept is, well, an intercept, and doesn't notice rank-deficiencies because of that!

inlabru should probably always set expand.factor.strategy="inla", but then handle factors in a more user-transparent manner (consistent with basic R-INLA when appropriate), also taking into account the special meaning of Intercept.

Having a well defined intended factor handling behaviour would help in debugging this issue. I know we discussed this, but don't remember the outcome.

fbachl commented 6 years ago

I think the outcome was to set the default factor strategy to "inla". My plan was to deal with that issue once I'm streamlining the way inlabru handles different latent models. Since that's exactly what I'm currently doing I'll also try to deal with the factor issue.

ASeatonSpatial commented 6 years ago

options = list(control.fixed = list(expand.factor.strategy='inla')) solved the problem. I'm not familiar with the internals of INLA so the error was hard for me to understand why it was occurring. The suggested solution works though, i was just missing the outer list() when setting lgcp options. ]

Following the logic of the lgcp help file I thought this was what i should be doing - options = bru.options(control.fixed = list(expand.factor.strategy = "inla")))

finnlindgren commented 6 years ago

"Solved" in the limited sense "removed the error message" or also the more desirable "will produce the correct result for all models"? With the extra interpretation layer of inlabru we need to check this a bit more before calling it fixed.

ASeatonSpatial commented 6 years ago

Also should mention I'm on ubuntu 17.10 so problem is not restricted to some Windows machines

ASeatonSpatial commented 6 years ago

'solved' as in 'removed the error message'. Of course I'm trusting that expand.factor.strategy = 'inla' is doing something sensible.

finnlindgren commented 6 years ago

With the 'inla' factor strategy and the specific model cmp1 = coordinates ~ veg(map = gcov$vegetation, model = "factor") - Intercept then yes, the result should be correct.

The potential problems occur with + Intercept and when there is more than one factor component in the model.

finnlindgren commented 6 years ago

As of aa2209e the default is to set expand.factor.strategy="inla"