bym pc prior problem

[Patrick Brown]

Hi folks.  The bym2 model isn't working for me now, I think there's a problem with the pc.bym.phi function.  The prior for logit-phi is monotone decreasing from -3, regardless of what I set lambda to.  Please help!

p

library(INLA)

priorTable = INLA:::inla.pc.bym.phi(

graph = system.file(

"demodata/germany.graph", 

package="INLA"),

lambda=-log(0.9)/0.9, 

return.as.table = TRUE, 

debug=TRUE)

priorMat = matrix(scan(

text=gsub("[[:alpha:]]+:", "", priorTable),

quiet=TRUE), ncol=2)

plot(priorMat[,1], 

exp(priorMat[,2]), 

xlab='logit phi', 

ylab='dens')

[Håvard Rue]

Patrick,

are you happy with this one?

phi = 1/(1+exp(-seq(-10, 10, len=1000)))
for(graph in system.file("demodata/germany.graph", package="INLA")) {
    prior = INLA:::inla.pc.bym.phi(graph = graph)
    dev.new()
    plot(phi, prior(phi), type="l",  lwd=2, main = graph)

--
Håvard Rue
Department of Mathematical Sciences
Norwegian University of Science and Technology
N-7491 Trondheim, Norway
Voice: +47-7359-3533 URL : http://www.math.ntnu.no/~hrue
Mobile: +47-9260-0021 Email: havar...@math.ntnu.no

R-INLA: www.r-inla.org

[Patrick Brown]

yes, that looks good when returning the function.  the table is different.

library(INLA)

phi = 1/(1+exp(-seq(-10, 10, len=1000)))

priorFun = INLA:::inla.pc.bym.phi(

graph = system.file(

"demodata/germany.graph", 

package="INLA"),

u=0.2, alpha=0.8)

dev.new()

plot(phi, priorFun(phi), type='l')

priorTable = INLA:::inla.pc.bym.phi(

graph = system.file(

"demodata/germany.graph", 

package="INLA"),

u=0.2, alpha=0.8,

return.as.table = TRUE)

priorMat = matrix(scan(

text=gsub("[[:alpha:]]+:", "", priorTable),

quiet=TRUE), ncol=2)

dev.new()

plot(priorMat[,1], 

exp(priorMat[,2]), 

xlab='logit phi', 

ylab='dens')

[Håvard Rue]

the function return

prior(phi)



the table returns log-prior of log(phi/(1-phi))

doing this involves adding a log(jacobian)

?

[Patrick Brown]

Done some more digging.

The 'stable' version of INLA produces a table with small values of logit phi, but the 'testing' inla table only includes values of logit phi above -2.77.  The prior mode of logit phi is much smaller than -2.77 and the optimizer was trying logit phi of about -20 when searching for the posterior mode.

The table needs to go further to the left.  I'll try a few more things.

p

[Patrick Brown]

Priors seem to be wrong in the 'testing' version of INLA, though even the 'stable' version is very insensitive to 'u'.  Attached are two sets of plots from testing and stable INLA.

# install.packages("INLA", repos="https://www.math.ntnu.no/inla/R/stable", lib = '/tmp')

#library('INLA', lib.loc = '/tmp', verbose=TRUE)

library('INLA', verbose=TRUE)

phi = 1/(1+exp(-seq(-5, 5, len=1000))) 

graph = system.file("demodata/germany.graph", package="INLA")

ualpha = expand.grid(u=c(0.1, 0.9), alpha = c(0.01, 0.5))

pdf(paste("/tmp/inla", inla.version("version"), ".pdf", sep=""))

par(mfrow=c(2,2))

for(D in 1:nrow(ualpha)) {

  prior = INLA:::inla.pc.bym.phi(

    graph = graph, u=ualpha[D,'u'], alpha=ualpha[D,'alpha']) 

  priorLogit = inla.tmarginal(

    function(x) log(x/(1-x)),

    cbind(phi, exp(prior(phi)))

  )

  plot(priorLogit, type="l",  lwd=2, 

    main = paste('u=',ualpha[D,'u'], 'alpha=',ualpha[D,'alpha']))

}

dev.off()

Probably related to whatever's causing aa < bb below

phi = 0.05

graph = system.file(

  "demodata/germany.graph", 

  package="INLA")

adjust.for.con.comp = TRUE

Q = INLA:::inla.pc.bym.Q(graph)

n = dim(Q)[1]

g = INLA:::inla.read.graph(Q)

nc = g$cc$n

rankdef = nc

log.q1.det = INLA:::inla.sparse.det.bym(Q, 

  adjust.for.con.comp = adjust.for.con.comp)

constr = list(A = matrix(0, nc, n), e = rep(0, nc))

for(k in 1:nc) {

  constr$A[k, g$cc$nodes[[k]]] = 1

}

Qinv.d = diag(inla.qinv(

    Q + Diagonal(n) * max(diag(Q)) * eps, 

    constr))

f = mean(Qinv.d) - 1.0

aa = n*phi*f

if (phi <= 0.5) {

  eps = sqrt(.Machine$double.eps)

  rdef = 0

} else {

  eps = 0

  rdef = 1

}

bb = n * log((1.0 - phi)/phi) +

  INLA:::inla.sparse.det.bym(Q + (phi/(1-phi) + eps) * 

      Diagonal(n), rankdef = rdef,

    adjust.for.con.comp = adjust.for.con.comp) -

  (log.q1.det - (n-rankdef) * log(phi))

c(aa, bb)

if (aa >= bb) sqrt(aa-bb) else NA

[Håvard Rue]

thanks; very weird. I'll check

> -- 
> 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 post to this group, send email to r-inla-discussion-group@googlegr
> oups.com.
> Visit this group at https://groups.google.com/group/r-inla-discussion
> -group.
> For more options, visit https://groups.google.com/d/optout.

[INLA help]

Patrick,

can you try again now? I havn't updated the stable one. as far as I'm
able to check, the results now should be ok. the numerics comes into
play here, as two approaches, which are mathematically the same, gives
slightly different results due to almost singular matrices (one based
on eigenvalues, and one based on a determinant identity). but its almos
only for \phi close to zero and the err is on the good side. some of
the details are also like that we need to define what is correct, as
its not obivous.

thanks!!!!!!
H

> > send an email to r-inla-discussion-group+unsubscribe@googlegroups.c

> > om
> > .
> > To post to this group, send email to r-inla-discussion-group@google
> > gr
> > oups.com.
> > Visit this group at https://groups.google.com/group/r-inla-discussi
> > on
> > -group.
> > For more options, visit https://groups.google.com/d/optout.
>
> -- 
> Håvard Rue
> Department of Mathematical Sciences
> Norwegian University of Science and Technology
> N-7491 Trondheim, Norway
> Voice:  +47-7359-3533    URL  : http://www.math.ntnu.no/~hrue  
> Mobile: +47-9260-0021    Email: havar...@math.ntnu.no
>
> R-INLA: www.r-inla.org
>

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