another knmodels question: ar1 specification possible?

69 views
Skip to first unread message

pav

unread,
Jun 24, 2021, 9:56:14 AM6/24/21
to R-inla discussion group

A working form of the help-page example for inla.knmodels() is:

tgraph <- sparseMatrix(i=c(2:10, 1:9), j=c(1:9, 2:10), x=1)
res <- inla.knmodels(y ~ f(time, model='bym2', graph=tgraph) +
                       f(space, model='bym2', graph=graph),
                     data=dat, family='poisson', E=dat$E, progress=TRUE,
                     control.st=list(time=time, space=space,
                                     spacetime=spacetime, graph=graph,
type=c(4, '4c', '4d')),
                     control.compute=list(dic=TRUE, waic=TRUE, cpo=TRUE),
                     control.inla=list(strategy='gaussian'))

But what if one wants to specify, say, an ar1 process for the
time terms? The obvious suggestion is to replace

    f(time, model='bym2', graph=tgraph)

with

   f(time, model='ar1')

But this produces an error:

Error in which(r.size == n) : 
  'list' object cannot be coerced to type 'integer'

as do attempts with specifying rw1, rw2, or iid. Can anyone suggest a fix? 

Elias T. Krainski

unread,
Jun 24, 2021, 10:19:02 AM6/24/21
to pav, R-inla discussion group
Hi,

The issue with AA being singular happens at the current INLA version when using PARDISO.

For now, please try the following code:

### doing the type 4 (step by step)
R.t <- crossprod(diff(Diagonal(m)))
R.s <- Diagonal(n, colSums(graph))-graph
st.C <- kronecker(R.t, R.s)

str(head(dat,2), , 30)
kronecker(matrix(1,1,3), diag(4)) ### each area sum-to-zero constraint type
kronecker(diag(3), matrix(1,1,4)) ### each time sum-to-zero constraint type

f4 <- y ~ f(time, model='ar1') + ##bym2', graph=tgraph) +
    f(space, model='bym2', graph=graph) +
    f(spacetime, model='generic0', Cmatrix=st.C,
      extraconstr=list(A=rbind(kronecker(matrix(1, 1, m), diag(n))[-1,],
                               kronecker(diag(m), matrix(1, 1, n))),
                       e=rep(0, n+m-1)))

res4 <- inla(f4, 'poisson', E=dat$expected,
             data=dat[c('space', 'time', 'spacetime', 'y')],
             verbose=TRUE)
###             control.inla=list(strategy='gaussian')) ## use this if using PARDISO (will be solved in the next INLA versions)

best regards,
Elias

--
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/7b1c4968-76f4-4f3f-be1c-626c722f43aan%40googlegroups.com.
Reply all
Reply to author
Forward
0 new messages