inlabru-org / inlabru

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

how to pass parameters to control.fixed? #46

Open ThierryO opened 5 years ago

ThierryO commented 5 years ago

Create some example data

set.seed(20180827)
n <- 1000
dataset <- data.frame(
  Area = runif(n, min = 1, max = 4),
  X = rnorm(n, mean = 2, sd = 1)
)
dataset$logArea <- log(dataset$Area)
dataset$eta <- 2 + 0.7 * dataset$X
dataset$lambda <- exp(dataset$eta) * dataset$Area
dataset$Y <- rpois(n, dataset$lambda)

Model it with INLA. First using the defaults, then fixed the parameter of logArea to 1 by setting the prior to mean = 1 and var = 0.

library(INLA)
inla(
  Y ~ X + logArea,
  family = "poisson",
  data = dataset
)$summary.fixed
inla(
  Y ~ X + logArea,
  family = "poisson",
  control.fixed = list(
    mean = list(logArea = 1, default = 0),
    prec = list(logArea = 1e20, default = 0.001)
  ),
  data = dataset
)$summary.fixed

bru() seems to ignore the options when passed as a list. Passing the output of bru.options() yields an error.

library(inlabru)
m <- bru(
  Y ~ X + logArea,
  family = "poisson",
  options = list(
    control.fixed = list(
      mean = list(logArea = 1, default = 0),
      prec = list(logArea = 1e20, default = 0.001)
    )
  ),
  data = dataset
)
summary(m)
bru(
  Y ~ X + logArea,
  family = "poisson",
  options = bru.options(
    control.fixed = list(
      mean = list(logArea = 1, default = 0),
      prec = list(logArea = 1e20, default = 0.001)
    )
  ),
  data = dataset
)

Error in (function (data, model, stackmaker, n = 10, result = NULL, family, : INLA returned message: unused argument (inla.options = list(list(TRUE, TRUE, TRUE), list("auto"), list(list(1, 0), list(1e+20, 0.001)))) In addition: Warning message: In (function (formula, family = "gaussian", contrasts = NULL, data, :

Error in (function (data, model, stackmaker, n = 10, result = NULL, family, : INLA returned message: unused argument (inla.options = list(list(TRUE, TRUE, TRUE), list("auto"), list(list(1, 0), list(1e+20, 0.001))))

SessionInfo

Session info --------------------------------------------------------------------------------------------------
 setting  value                       
 version  R version 3.5.1 (2018-07-02)
 system   x86_64, linux-gnu           
 ui       RStudio (1.1.453)           
 language nl:en                       
 collate  nl_NL.UTF-8                 
 tz       Europe/Brussels             
 date     2018-08-27                  

Packages ------------------------------------------------------------------------------------------------------
 package      * version date       source                         
 assertthat     0.2.0   2017-04-11 CRAN (R 3.5.1)                 
 cli            1.0.0   2017-11-05 CRAN (R 3.5.1)                 
 colorspace     1.3-2   2016-12-14 CRAN (R 3.5.1)                 
 crayon         1.3.4   2017-09-16 CRAN (R 3.5.1)                 
 dichromat      2.0-0   2013-01-24 CRAN (R 3.5.1)                 
 digest         0.6.15  2018-01-28 CRAN (R 3.5.1)                 
 fansi          0.2.3   2018-05-06 CRAN (R 3.5.1)                 
 ggplot2      * 3.0.0   2018-07-03 CRAN (R 3.5.1)                 
 glue           1.3.0   2018-07-17 CRAN (R 3.5.1)                 
 graphics     * 3.5.1   2018-07-03 local                          
 grDevices    * 3.5.1   2018-07-03 local                          
 grid           3.5.1   2018-07-03 local                          
 gtable         0.2.0   2016-02-26 CRAN (R 3.5.1)                 
 inlabru      * 2.1.9   2018-08-23 Github (fbachl/inlabru@d77d689)
 labeling       0.3     2014-08-23 CRAN (R 3.5.1)                 
 lattice        0.20-35 2017-03-25 CRAN (R 3.5.0)                 
 lazyeval       0.2.1   2017-10-29 CRAN (R 3.5.1)                 
 magrittr       1.5     2014-11-22 CRAN (R 3.5.1)                 
 MASS           7.3-50  2018-04-30 CRAN (R 3.5.0)                 
 Matrix       * 1.2-14  2018-04-09 CRAN (R 3.5.0)                 
 methods      * 3.5.1   2018-07-03 local                          
 mgcv           1.8-24  2018-06-18 CRAN (R 3.5.0)                 
 munsell        0.5.0   2018-06-12 CRAN (R 3.5.1)                 
 nlme           3.1-137 2018-04-07 CRAN (R 3.5.0)                 
 pillar         1.3.0   2018-07-14 CRAN (R 3.5.1)                 
 plyr           1.8.4   2016-06-08 CRAN (R 3.5.1)                 
 R6             2.2.2   2017-06-17 CRAN (R 3.5.1)                 
 RColorBrewer   1.1-2   2014-12-07 CRAN (R 3.5.1)                 
 Rcpp           0.12.17 2018-05-18 CRAN (R 3.5.1)                 
 reshape2       1.4.3   2017-12-11 CRAN (R 3.5.1)                 
 rgdal          1.3-4   2018-08-03 cran (@1.3-4)                  
 rgeos          0.3-28  2018-06-08 CRAN (R 3.5.1)                 
 rlang          0.2.1   2018-05-30 CRAN (R 3.5.1)                 
 scales         0.5.0   2017-08-24 CRAN (R 3.5.1)                 
 sp           * 1.3-1   2018-06-05 CRAN (R 3.5.1)                 
 stats        * 3.5.1   2018-07-03 local                          
 stringi        1.2.3   2018-06-12 CRAN (R 3.5.1)                 
 stringr      * 1.3.1   2018-05-10 CRAN (R 3.5.1)                 
 tibble       * 1.4.2   2018-01-22 CRAN (R 3.5.1)                 
 tools          3.5.1   2018-07-03 local                          
 utf8           1.1.4   2018-05-24 CRAN (R 3.5.1)                 
 utils        * 3.5.1   2018-07-03 local                          
 viridisLite    0.3.0   2018-02-01 CRAN (R 3.5.1)                 
 withr          2.1.2   2018-03-15 CRAN (R 3.5.1)    
finnlindgren commented 5 years ago

I have an email thread with another user from a few weeks back that looked at this in some detail, including providing some solutions. I’ll try to extract the relevant bits and put here, awaiting a full solution in the package.

finnlindgren commented 5 years ago

Even if inlabru would recognise control.fixed and pass it on to inla(), this would not have the desired effect. The reason is that internally, fixed effect terms are re-interpreted as f(..., model="linear") terms before the inla() call, and control.fixed doesn't apply to such model components.

Until we come up with a way to automate it by using the information in control.fixed in the re-interpretation step, the following workaround can be used: Replace

Y ~ X + logArea

by

Y ~ X(map=X, model="linear", mean.linear=0, prec.linear=0.001) +
  logArea(map=logArea, model="linear", mean.linear=1, prec.linear=1e20)

Full example attached: code46.zip

Note that this also applies to the intercept term, so to control the prior for that you need to include it explicitly, with + Intercept(map=Inter, model="linear", prec.linear=...) and have Inter=1 in the data data.frame.

Dev notes: The sensible thing seems to be if we can parse control.fixed, but I'm not sure if the point in the processing where the inla formula is constructed is too early for that to work. (If it is, then it needs to be adjusted in the code before calling inla)

ThierryO commented 5 years ago

The workaround works fine.

I would suggest to update the documentation that mentions control.fixed (?bru and ?bru.options) rather than putting a lot of effort in the parsing.

This issue can be closed.

finnlindgren commented 5 years ago

Ah, I didn’t realise that’s in the documentation! Will fix it!

finnlindgren commented 5 years ago

I've updated the documentation in commit 82764e7 and 48cd2cb, which also enables only partially overriding control.* parameters. I'll keep this open as a wishlist enhancement item.

finnlindgren commented 5 years ago

Updated code (removing no longer needed control.fixed from the correct bru() call, and more precise comment on the output): code46.zip