Problems running full ST model

119 views
Skip to first unread message

simon brewer

unread,
Mar 25, 2021, 2:48:03 PM3/25/21
to R-inla discussion group

Hi all

I am trying to fit a fairly large spatial-temporal model, but am having problems using the full dataset. The data have about 3000 locations and 365 time steps, and fits well, with no problems to about 60% of the data (all spatial locations but 60% of the time steps).  

I found another answer here (https://groups.google.com/g/r-inla-discussion-group/c/uf2ZGh4jmWc) that suggested adding a small non-negative diagonal value in control.inla, and then using the output of that model as starting values in a second model with diagonal = 0. If I do this, then the first model will run on the full dataset, but I cannot get the second model to run. This crashes with the following error:

inla.mkl: smtp-pardiso.c:757: GMRFLib_pardiso_solve_core: Assertion `store->pstore->done_with_chol == GMRFLib_TRUE' failed.

Segmentation fault

Error in inla.inlaprogram.has.crashed() : 

  The inla-program exited with an error. Unless you interupted it yourself, please rerun with verbose=TRUE and check the output carefully.

  If this does not help, please contact the developers at <he...@r-inla.org>.

 

Do you have any suggestions about how to fix this? Or where I may be going wrong? I’ve copied the code below. 


Thank you!
Simon Brewer

## Priors

### BYM2

prec_u <- .5/.31

prec_alpha <- .01

phi_u <- .5

phi_alpha <- 2/3

prior_bym2 <- list(

  prec = list(

    prior = "pc.prec",

    param = c(prec_u, prec_alpha)),

  phi = list(

    prior = "pc",

    param = c(phi_u, phi_alpha))

)


### RW1 From INLA Germany example 

rw1_u = .2/.31

rw1_alpha = .01

prior_rw1 <- list(theta = list(prior="pc.prec", param=c(rw1_u, rw1_alpha)))


## Model formula

f8 <- TESTSNEW7POS_CDC ~ scale(Pop_m) + scale(gini) + scale(Uninsured) + scale(Production) + RUCC + scale(PCP) + scale(RPL_THEMES) +

  f(date1, model = "rw1", hyper = prior_rw1, 

scale.model = TRUE) +

f(struct, model = "bym2", graph = g, hyper = prior_bym2, 

scale.model = TRUE)


## Model w/ diagonal value

res1 <- inla(f8, data = dat2,

control.predictor = list(compute = TRUE),

               control.compute = list(dic = TRUE, waic = TRUE), 

               control.fixed = list(prec.intercept = 1),

               control.inla = list(strategy = "gaussian",

                                   int.strategy = "eb",

                                   diagonal = 5e-2),

               verbose = TRUE, num.threads = 32)


## Model w/ diagonal value

res1 <- inla(f8, data = dat2,

control.predictor = list(compute = TRUE),

               control.compute = list(dic = TRUE, waic = TRUE), 

               control.fixed = list(prec.intercept = 1),

               control.inla = list(strategy = "gaussian",

                                   int.strategy = "eb",

                                   diagonal = 5e-2),

               verbose = TRUE, num.threads = 32)


simon brewer

unread,
Mar 25, 2021, 7:53:26 PM3/25/21
to R-inla discussion group
Apologies - I copied the first model code twice in the previous email. Here's the code with the second model included:

## Model w/ restart

 res1 <- inla(f8, data = dat2,

            control.predictor = list(compute = TRUE), 

control.compute = list(dic = TRUE, waic = TRUE),

control.fixed = list(prec.intercept = 1),

control.mode = list(result = res1, restart = TRUE),

control.inla = list(diagonal = 0),

verbose = TRUE, num.threads = 32,

)


Helpdesk

unread,
Mar 27, 2021, 8:02:41 AM3/27/21
to simon brewer, R-inla discussion group
On Thu, 2021-03-25 at 11:48 -0700, simon brewer wrote:
> control.predictor = list(compute = TRUE),
>                control.compute = list(dic = TRUE, waic = TRUE), 


if you remove those, the RAM requirement would go down. 

also set num.threads to

num.threads="A:-1"

where A=number of hyperparamters.


I guess you're using the PARDISO library ?

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

Helpdesk

unread,
Mar 27, 2021, 8:28:24 AM3/27/21
to simon brewer, R-inla discussion group

it runs fine here


maybe its a RAM issue, it seems to peak at about 6Gb

see if its run when setting

control.inla=list(int.strategy="eb")

also, please get the license for the PARDISO library, as it will run
more smooth and faster

with PARDISO 20 seconds
without PARDISO 127 seconds


let me know how it goes.

PS use an academic email address when requesting the license....

H


On Thu, 2021-03-25 at 11:48 -0700, simon brewer wrote:
> --
> 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/7560ccc4-ba8b-45b7-9dcd-c42216c3a8d6n%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org
Screenshot from 2021-03-27 15-23-31.png

Helpdesk

unread,
Mar 27, 2021, 8:30:28 AM3/27/21
to simon brewer, R-inla discussion group

I can rerun here, if you attach data (or some other data that give the
same issues)
> --
> 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/9a66fce1-c4c0-484d-b246-c3ad76a56878n%40googlegroups.com

simon brewer

unread,
Mar 27, 2021, 11:07:56 AM3/27/21
to R-inla discussion group
Thank you for the quick reply. I'll try your suggestions today, and get back to you. If I'm still running into the same problem I'll send on the data. 

Best wishes
Simon

Reply all
Reply to author
Forward
0 new messages