Posterior eval error when using experimental mode

[Connor Gascoigne]

Hi All,

I am attempting to fit a binomial model on survey data and then generate samples from the posterior using inla.posterior.sample() and inla.posterior.sample.eval().

Due to the data set I am using being large, I can only fit the model when I use inla.mode = 'experimental' instead of inla.mode = 'classic'. This seems to cause an issue when I get down to using the inla.posterior.sample.eval() function as the model fit using inla.mode = 'experimental' returns the following error "Error in a.sample$latent[.contents$start[i]:(.contents$start[i] + .contents$length[i] -  : subscript out of bounds" whereas the model fit using inla.mode = 'classic' does not return this error.

To highlight this, below I have included an example of the code I am running with some fake data. Any variables labelled with 1 (e.g., fit1, samp1, eval1a) are relating to the model fit using inla.mode = 'experimental' and any variables labelled with 2 are relating to the model fit using inla.model = 'classic'. The aforementioned error is returned when I attempt to run eval1a but not eval2a.

Looking at the error, it says there are issues relating to subscript in the latent part of the sample where the subscript is found from the model contents. Both model contents (fit1$misc$configs$contents and fit2$misc$configs$contents) are exactly the same and when looking at the latent sample, the only difference was the experimental mode labels the latent column V1 rather than sample:1 like the classic mode. I changed this to see if it did anything but unfortnely the error still persisted and I am at a bit of a loss as what to do now...

If anyone has any advice as how to run inla.mode = 'experimental' with inla.posterior.sample.eval(), it would be greatly appreciated! 

Thank you,

Connor

P.s. Also, on a really quick second note... at the bottom eval1c and eval2c are more similar to what I am doing with the actual data. For the function 'my.fun' used within eval1c and eval2c, is this the correct way to handle factors and random effects that are put into the model?

###############################################################

# packages ----

library(tidyverse)
library(INLA)

# data ----

# setting the seed
set.seed(1234)

# length of data
n = 100

dat <-
  data.frame(
    # random id
    id = 1:n,
    # age: continuous, adults between 16 and 100
    age = sample(16:100, n, replace = TRUE),
    # sex: factor 1 for male and 2 for female
    sex = sample(1:2, n, replace = TRUE) %>% factor(., levels = 1:2),
    # marital status: factor 1-9 for different classifications
    marStat = sample(1:9, n, replace = TRUE) %>% factor(., levels = 1:9),
    # observed outcome: binary 0 or 1
    y = sample(0:1, n, replace = TRUE))

# formula ----

formula <- y ~ 1 + age + sex + marStat + f(id, model = 'iid')

# model fit ----

verbose = TRUE

## using inla.mode = 'experimental' ----

fit1<-
  INLA::inla(formula,
             family = 'binomial',
             data = dat,
             control.compute = list(config = TRUE),
             control.predictor = list(compute = FALSE, link = 1),
             control.inla = list(strategy = 'adaptive', int.strategy = 'auto'),
             inla.mode = 'experimental',
             lincomb = NULL,
             verbose = verbose)

## using inla.mode = 'classic' ----

fit2<-
  INLA::inla(formula,
             family = 'binomial',
             data = dat,
             control.compute = list(config = TRUE),
             control.predictor = list(compute = FALSE, link = 1),
             control.inla = list(strategy = 'adaptive', int.strategy = 'auto'),
             inla.mode = 'classic',
             lincomb = NULL,
             verbose = verbose)

# prediction ----

## number of samples ----

nsim = 100

## stop predictor being sampled for compuational reasons ----

# contents tags
cs1 <- fit1$misc$configs$contents$tag
cs1 <- cs1[cs1 != "Predictor"]
select1 <- list()

cs2 <- fit2$misc$configs$contents$tag
cs2 <- cs2[cs2 != "Predictor"]
select2 <- list()


for (i in 1:length(cs1)) {
  select1[[i]] <- 0
  names(select1)[i] <- cs1[i]
 
  select2[[i]] <- 0
  names(select2)[i] <- cs2[i]
}

## samples from posterior ----

samp1 <- INLA::inla.posterior.sample(n = nsim, result = fit1, intern = TRUE, selection = select1, verbose = verbose)
samp2 <- INLA::inla.posterior.sample(n = nsim, result = fit2, intern = TRUE, selection = select2, verbose = verbose)

## evaluating the posterior samples ----

### simple example only wishing to return the intercept ----

eval1a <- INLA::inla.posterior.sample.eval(samp1, fun = function(x) Intercept)
eval2a <- INLA::inla.posterior.sample.eval(samp2, fun = function(x) Intercept)

### looking into error ----

#### difference in the content ----

all.equal(fit1$misc$configs$contents, fit2$misc$configs$contents)

#### difference in the samples ----

test1a <- samp1[[1]]$latent
test2a <- samp2[[1]]$latent

### changing experimental column name from V1 sample:1 like classic ----

for(i in 1:nsim){
  colnames(samp1[[i]]$latent) <- paste0('sample:', i)
}

test1b <- samp1[[1]]$latent
test2b <- samp2[[1]]$latent

### re-attempt eval ----


eval1b <- INLA::inla.posterior.sample.eval(samp1, fun = function(x) Intercept)
eval2b <- INLA::inla.posterior.sample.eval(samp2, fun = function(x) Intercept)


### similar example to what I am attempting to do ----

newdata <-
  dat  %>%
  # expand factors
  ## need to recreate model matrix for factor terms
  ## drop the intercept columns after
  cbind(model.matrix(~ sex, data = .),
        model.matrix(~ marStat, data = .)) %>%
  dplyr::select(-contains('(Intercept)'))

my.fun <- function(newdata = NA){
 
  linearPredictor <-
    Intercept +
    age * newdata$age +
    sex2 * newdata$sex2 +
    marStat2 * newdata$marStat2 +
    marStat3 * newdata$marStat3 +
    marStat4 * newdata$marStat4 +
    marStat5 * newdata$marStat5 +
    marStat6 * newdata$marStat6 +
    marStat7 * newdata$marStat7 +
    marStat8 * newdata$marStat8 +
    marStat9 * newdata$marStat9 +
    id[newdata$id]
 
  return(linearPredictor)
 
}

eval1c <- INLA::inla.posterior.sample.eval(samp1, fun = my.fun, newdata = newdata)
eval2c <- INLA::inla.posterior.sample.eval(samp2, fun = my.fun, newdata = newdata)

[Helpdesk]

if you remove the 'selection' argument, it runs fine here.

that argument is a little tricky as the internal indexing change, and it
does not speed up the calculations either...

##, selection = select1)

> --
> 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/e6e3a3f3-05bc-447e-b9d2-95f378dc98d9n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Connor Gascoigne]

Hi Håvard,

Thank you so much for getting back to me. I did keep the 'selection' argument in as I thought it sped things up, but it all works fine now after removing it.

Again, thank you!

Connor