allow group argument with spde models

[ThierryO]

The current implementation yields an error when adding the group argument. Below is a minimal example. The goal is to fit a model with different seasonal pattern for each year. The patterns has a common range and variance.

load packages

library(dplyr)
library(tidyr)
library(purrr)
library(ggplot2)
library(INLA)
library(inlabru)

set parameters

n_year <- 10
n_sample <- 20
doy <- 250:330
sd_shift <- 5
sd_season <- 10
size <- 0.5

generate dataset

data.frame(
  year = seq_len(n_year),
  intercept = rnorm(n_year, mean = 4, sd = 0.2),
  shift = rnorm(n_year, mean = 0, sd = sd_shift)
) %>%
  mutate(
    doy = map(
      year,
      function(x){sample(doy, n_sample)}
    )
  ) %>%
  unnest() %>%
  mutate(
    eta = intercept + 20 * dnorm(doy, mean = mean(doy) + shift, sd = sd_season),
    count = rpois(n(), lambda = exp(eta))
  ) -> dataset
ggplot(dataset, aes(x = doy, y = count)) +
  geom_point() +
  geom_line(aes(y = exp(eta))) +
  facet_wrap(~year)

fit model with common seasonal pattern

mesh_doy <- inla.mesh.1d(doy, boundary = "free")
spde_doy <- inla.spde2.pcmatern(
  mesh_doy,
  prior.range = c(10, 0.01),
  prior.sigma = c(3, 0.01)
)
model <- bru(
  count ~ 
    season(map = doy, model = spde_doy) +
    trend(map = year, model = "iid", n = n_year), 
  family = "poisson", 
  data = dataset
)

fit model with different seasonal pattern per year this doesn't work at the moment

model <- bru(
  count ~ 
    season(map = doy, model = spde_doy, group = year), 
  family = "poisson", 
  data = dataset
)

[finnlindgren]

Can you include the exact error message as well please, and which branch generates that error?

For a group model one needs to specify a group model with control.group as well, so I'd expect an error when that's not there (but a different one when it is there).

Independent replicates should also give error messages in the devel branch, but not in the backend branch I believe, with the "replicate" argument instead of "group".

[ThierryO]

here is the updated code, the error and the session info (devel branch)

model <- bru(
  count ~ 
    season(
      map = doy, 
      model = spde_doy, 
      group = year, 
      control.group = list(model = "ar1")
    ), 
  family = "poisson", 
  data = dataset
)

Error in INLA::f(xxx, ..., group = group, model = model) : object 'year' not found

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-23

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)

[ThierryO]

On same code yields a different error on the backend branch

Error in (function (data, model, stackmaker, n = 10, result = NULL, family, : INLA returned message: object 'season.group' not found

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-23

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.7.999 2018-08-23 Github (fbachl/inlabru@83c9d6d) 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]

Thanks!

[finnlindgren]

Sorry for the long wait; I believe this was fixed in the devel branch quite a while ago (2.1.13.999). Needed a small tweak in unnest() for new tidyverse versions, but otherwise it seems ok:

library(dplyr)
#> 
#> Attaching package: 'dplyr'
#> The following objects are masked from 'package:stats':
#> 
#>     filter, lag
#> The following objects are masked from 'package:base':
#> 
#>     intersect, setdiff, setequal, union
library(tidyr)
library(purrr)
library(ggplot2)

# set parameters

n_year <- 10
n_sample <- 20
doy <- 250:330
sd_shift <- 5
sd_season <- 10
size <- 0.5

# generate dataset

data.frame(
  year = seq_len(n_year),
  intercept = rnorm(n_year, mean = 4, sd = 0.2),
  shift = rnorm(n_year, mean = 0, sd = sd_shift)
) %>%
  mutate(
    doy = map(
      year,
      function(x){sample(doy, n_sample)}
    )
  ) %>%
  unnest(cols = c(doy)) %>%
  mutate(
    eta = intercept + 20 * dnorm(doy, mean = mean(doy) + shift, sd = sd_season),
    count = rpois(n(), lambda = exp(eta))
  ) -> dataset
ggplot(dataset, aes(x = doy, y = count)) +
  geom_point() +
  geom_line(aes(y = exp(eta))) +
  facet_wrap(~year)

library(INLA)
#> Loading required package: Matrix
#> 
#> Attaching package: 'Matrix'
#> The following objects are masked from 'package:tidyr':
#> 
#>     expand, pack, unpack
#> Loading required package: sp
#> Loading required package: parallel
#> Loading required package: foreach
#> 
#> Attaching package: 'foreach'
#> The following objects are masked from 'package:purrr':
#> 
#>     accumulate, when
#> This is INLA_20.04.18 built 2020-04-24 13:46:03 UTC.
#> See www.r-inla.org/contact-us for how to get help.
#> To enable PARDISO sparse library; see inla.pardiso()
library(inlabru)
#> 
#> Attaching package: 'inlabru'
#> The following object is masked from 'package:purrr':
#> 
#>     map

# fit model with common seasonal pattern

mesh_doy <- inla.mesh.1d(doy, boundary = "free")
spde_doy <- inla.spde2.pcmatern(
  mesh_doy,
  prior.range = c(10, 0.01),
  prior.sigma = c(3, 0.01)
)
model <- bru(
  count ~
    season(map = doy, model = spde_doy) +
    trend(map = year, model = "iid", n = n_year),
  family = "poisson",
  data = dataset
)

model <- bru(
  count ~
    season(
      map = doy,
      model = spde_doy,
      group = year,
      control.group = list(model = "ar1")
    ),
  family = "poisson",
  data = dataset
)

INLA:::summary.inla(model)
#> 
#> Call:
#>    "(function (formula, family = \"gaussian\", contrasts = NULL, data, 
#>    quantiles = c(0.025, 0.5, 0.975), E = NULL, offset = NULL, scale = 
#>    NULL, weights = NULL, Ntrials = NULL, strata = NULL, link.covariates = 
#>    NULL, verbose = FALSE, lincomb = NULL, selection = NULL, 
#>    control.compute = list(), control.predictor = list(), control.family = 
#>    list(), control.inla = list(), control.results = list(), control.fixed 
#>    = list(), control.mode = list(), control.expert = list(), 
#>    control.hazard = list(), control.lincomb = list(), control.update = 
#>    list(), only.hyperparam = FALSE, inla.call = 
#>    inla.getOption(\"inla.call\"), inla.arg = inla.getOption(\"inla.arg\"), 
#>    num.threads = inla.getOption(\"num.threads\"), blas.num.threads = 
#>    inla.getOption(\"blas.num.threads\"), keep = inla.getOption(\"keep\"), 
#>    working.directory = inla.getOption(\"working.directory\"), silent = 
#>    inla.getOption(\"silent\"), debug = inla.getOption(\"debug\"), 
#>    .parent.frame = parent.frame()) { old.options = options() 
#>    on.exit(options(old.options)) options(OutDec = \".\") my.time.used = 
#>    numeric(4) my.time.used[1] = Sys.time() all.hyper = list() data.orig = 
#>    data if (missing(formula)) { stop(\"Usage: inla(formula, family, data, 
#>    ...); see ?inla\\n\") } if (missing(data)) { stop(\"Missing 
#>    data.frame/list `data'. Leaving `data' empty might lead 
#>    to\\n\\t\\tuncontrolled behaviour, therefore is it required.\") } if 
#>    (!is.data.frame(data) && !is.list(data)) { stop(\"\\n\\tArgument `data' 
#>    must be a data.frame or a list.\") } if (!missing(weights) && 
#>    !is.null(weights)) { if 
#>    (!inla.getOption(\"enable.inla.argument.weights\")) { 
#>    stop(paste(\"Argument 'weights' must be enabled before use due to the 
#>    risk of mis-interpreting the results.\\n\", \"\\tUse 
#>    'inla.setOption(\\\"enable.inla.argument.weights\\\", TRUE)' to enable 
#>    it; see ?inla\")) } } family.alias = list(list(from = \"normal\", to = 
#>    \"gaussian\")) for (i in seq_along(family.alias <<<the rest is not 
#>    shown>>>" 
#> Time used:
#>     Pre = 0.705, Running = 0.621, Post = 0.235, Total = 1.56 
#> Fixed effects:
#>            mean    sd 0.025quant 0.5quant 0.975quant  mode kld
#> Intercept 4.204 0.086      4.032    4.204      4.375 4.204   0
#> 
#> Random effects:
#>   Name     Model
#>     season SPDE2 model
#> 
#> Model hyperparameters:
#>                   mean    sd 0.025quant 0.5quant 0.975quant  mode
#> Range for season 4.921 0.849      3.535    4.822      6.852 4.610
#> Stdev for season 0.306 0.044      0.231    0.302      0.404 0.293
#> 
#> Expected number of effective parameters(stdev): 56.20(3.28)
#> Number of equivalent replicates : 3.56 
#> 
#> Deviance Information Criterion (DIC) ...............: 1678.19
#> Deviance Information Criterion (DIC, saturated) ....: 472.62
#> Effective number of parameters .....................: 57.36
#> 
#> Watanabe-Akaike information criterion (WAIC) ...: 1757.47
#> Effective number of parameters .................: 110.84
#> 
#> Marginal log-Likelihood:  -901.58 
#> Posterior marginals for the linear predictor and
#>  the fitted values are computed

^{Created on 2020-04-30 by the reprex package (v0.3.0)}