bym with pc.prec

[Tracey Wang]

Hi All,

I am currently using INLA to run a bym model and very interested in using pc priors. 

I read the examples in the published paper (Simpson et al 2015), but was not able the run the sample code bym.R on my own computer. (I was using the most updated stable version of R-INLA). 

However, my question is that, instead of using 'bym2' with prec and phi, can I specify the priors for prec.unstruct and prec.spatial in 'bym' model using pc.prec? Indeed it seems working as I have run the model. But I am a bit unsure if it makes sense since no one seems to mention about it. 

(It would make more sense if I could run both and compare the results, but I was not able to run 'bym2' for some unknown reason.)

Thanks,

Tracey

[INLA help]

this one discuss bym2 better:

@Article{art585,
  author =       {A. Riebler and S. H. S{\o}rbye and D. Simpson and H.
Rue},
  title =        {An intuitive {B}ayesian spatial model for disease
                  mapping that accounts for scaling},
  journal =      {Statistical Methods in Medical Research},
  keywords =     {year2016,my.own,my.cv},
  year =         2016,
  volume =       {25},
  number =       {4},
  pages =        {1145-1165}
}




I dont know what you mean with 'could not run' ?? Are you using the
testing version???


Best
H


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

[Xiaotian Wang]

Hi Håvard,

I am using the stable version. Do I have to use the testing version in order to run "bym2"?

The error I got is as following:

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

In addition: Warning messages:

1: In inla.model.properties.generic(inla.trim.family(model), (mm[names(mm) ==  :

  Model 'bym2' in section 'latent' is marked as 'experimental'; changes may appear at any time.

  Use this model with extra care!!! Further warnings are disabled.

2: running command '"C:/Program Files/R/R-3.3.3/library/INLA/bin/windows/64bit/inla.exe"  -b -s -v  "C:/Users/xiawang/AppData/Local/Temp/2/RtmpEvmBVo/file1ffc63f67c24/Model.ini"' had status 255 

                                        

I am using a binomial model with two covariates:

formula.bym2 <- SAC ~ Cycle + Age + f(RegionID, model = "bym2", graph = ON_adjmat, scale.model = TRUE, constr = TRUE,

                                                                        hyper=list(phi = list(prior = "pc", param = c(0.2, 0.95)), 

                                                                                       prec = list(prior = "pc.prec", param = c(0.5, 0.025))))

INLAResult.bym2 <- inla(formula.bym2, family = "binomial", data = cchs_work, control.family = list(link="logit"),

                                    control.compute = list(cpo = T, dic = T), control.inla = list(strategy = "simplified.laplace"), 

                                    control.predictor = list(compute = T), verbose = TRUE)      

Thanks,

Tracey

[INLA help]

On Thu, 2017-04-27 at 15:55 +0000, Xiaotian Wang wrote:
> Hi Håvard,
>
> I am using the stable version. Do I have to use the testing version
> in order to run "bym2"?
>
> The error I got is as following:
>
> Error in inla.inlaprogram.has.crashed() : 

>   The inla-program exited with an error. Unless you interupted it\
>


if you can send me code/data so I can rerun it here, it would be very
helpful. thanks.

[Xiaotian Wang]

Hi Håvard, 

Thank you very much for your offer. I am not sure if I could send you the data since it contains confidential information. However, I have tried to run "bym2" with another data set and it seems working! 

Another question here is I tried to use INLAResult.bym2$all.hyper$random to check my priors and the red part as following seems confusing to me since I am supposed to have "pc" and "pc.prec" priors?

Thanks,

Tracey

> str(INLAResult.bym2$all.hyper$random)

List of 1

 $ :List of 3

  ..$ hyperid    : chr "RegionID"

  ..$ hyper      :List of 2

  .. ..$ theta1:List of 9

  .. .. ..$ hyperid   : atomic [1:1] 10001

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ name      : atomic [1:1] log unstructured precision

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ short.name: atomic [1:1] prec.unstruct

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ prior     : atomic [1:1] loggamma

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ param     : atomic [1:2] 0.35 0.025

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ initial   : atomic [1:1] 4

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ fixed     : atomic [1:1] FALSE

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ to.theta  :function (x)  

  .. .. .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 1 2 6 1 6 1 2

  .. .. .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002ef... <truncated>

  .. .. ..$ from.theta:function (x)  

  .. .. .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 1 2 6 1 6 1 2

  .. .. .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002d0... <truncated>

  .. ..$ theta2:List of 9

  .. .. ..$ hyperid   : atomic [1:1] 10002

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ name      : atomic [1:1] log spatial precision

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ short.name: atomic [1:1] prec.spatial

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ prior     : atomic [1:1] loggamma

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ param     : atomic [1:2] 0.35 0.025

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ initial   : atomic [1:1] 4

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ fixed     : atomic [1:1] FALSE

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ to.theta  :function (x)  

  .. .. .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 1 2 6 1 6 1 2

  .. .. .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000002e1... <truncated>

  .. .. ..$ from.theta:function (x)  

  .. .. .. ..- attr(*, "srcref")=Class 'srcref'  atomic [1:8] 1 1 2 6 1 6 1 2

  .. .. .. .. .. ..- attr(*, "srcfile")=Classes 'srcfilecopy', 'srcfile' <environment: 0x000000004c2... <truncated>

  ..$ group.hyper:List of 1

  .. ..$ theta:List of 9

  .. .. ..$ hyperid   : atomic [1:1] 40001

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ name      : atomic [1:1] logit correlation

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ short.name: atomic [1:1] rho

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ initial   : atomic [1:1] 1

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ fixed     : atomic [1:1] FALSE

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ prior     : atomic [1:1] normal

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ param     : atomic [1:2] 0 0.2

  .. .. .. ..- attr(*, "inla.read.only")= logi FALSE

  .. .. ..$ to.theta  :function (x, REPLACE.ME.ngroup)  

  .. .. .. ..- attr(*, "inla.read.only")= logi TRUE

  .. .. ..$ from.theta:function (x, REPLACE.ME.ngroup)  

  .. .. .. ..- attr(*, "inla.read.only")= logi TRUE

[INLA help]

if I run this:

library(MASS)
data(Germany)
g = inla.read.graph(system.file("demodata/germany.graph",
package="INLA"))
source(system.file("demodata/Bym-map.R", package="INLA"))
Germany = cbind(Germany,region.struct=Germany$region)
Q = -inla.graph2matrix(g)
diag(Q) = g$nnbs
Q = Q * exp(mean(log(diag(ginv(as.matrix(Q))))))

tau = 1
rho = 0.3
n = dim(Q)[1]
u = inla.qsample(n = 1,  Q = Q + 1e-8 * Diagonal(n),
        constr=list(A = matrix(1, 1, n), e=0))
v = rnorm(n)
x = 1/sqrt(tau) * (sqrt(1-rho) * v + sqrt(rho) * u)
intercept = log(Germany$E)
y = rpois(n, lambda = exp(intercept + x))

formula = y ~ 1 + f(region,
        model = "bym2", graph = g, constr=TRUE)
r = inla(formula, family="poisson",
        data=data.frame(y, region=Germany$region), keep=T)

par(mfrow=c(1, 2))
germany.map(r$summary.random$region$mean[1:n]-x)
plot(r$summary.random$region$mean[1:n],x)
abline(a=0, b=1)






I get



> str(r$all.hyper$random[[1]]$hyper$theta1$prior)
 atomic [1:1] pc.prec

 - attr(*, "inla.read.only")= logi FALSE

> str(r$all.hyper$random[[1]]$hyper$theta2$prior)
 chr "table: -12.6923076923077 -12.6898382145907 -12.6873687368737
-12.6848992591567 -12.6824297814397 -12.6799603037227 -12.67749082"|
__truncated__
> 



don't you get this?

H

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

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

[Xiaotian Wang]

Hi Håvard,

The same issue as I mentioned first time came to me again when I was using the code you provided - inla crashed when running "bym2":

R version 3.3.3 (2017-03-06) -- "Another Canoe"

Copyright (C) 2017 The R Foundation for Statistical Computing

Platform: x86_64-w64-mingw32/x64 (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.

You are welcome to redistribute it under certain conditions.

Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.

Type 'contributors()' for more information and

'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or

'help.start()' for an HTML browser interface to help.

Type 'q()' to quit R.

> library(MASS)

> library(INLA)

Loading required package: sp

Loading required package: Matrix

This is INLA 0.0-1485844051, dated 2017-01-31 (09:14:12+0300).

See www.r-inla.org/contact-us for how to get help.

> data(Germany)

> g = inla.read.graph(system.file("demodata/germany.graph",

+ package="INLA"))

> source(system.file("demodata/Bym-map.R", package="INLA"))

Read 3921 items

> Germany = cbind(Germany,region.struct=Germany$region)

> Q = -inla.graph2matrix(g)

> diag(Q) = g$nnbs

> Q = Q * exp(mean(log(diag(ginv(as.matrix(Q))))))

> tau = 1

> rho = 0.3

> n = dim(Q)[1]

> u = inla.qsample(n = 1,  Q = Q + 1e-8 * Diagonal(n),

+ constr=list(A = matrix(1, 1, n), e=0))

> v = rnorm(n)

> x = 1/sqrt(tau) * (sqrt(1-rho) * v + sqrt(rho) * u)

> intercept = log(Germany$E)

> y = rpois(n, lambda = exp(intercept + x))

> formula = y ~ 1 + f(region,

+ model = "bym2", graph = g, constr=TRUE)

> r = inla(formula, family="poisson",

+ data=data.frame(y, region=Germany$region), keep=T)

Model and results are stored in working directory [inla.model]

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

In addition: Warning messages:

1: In inla.model.properties.generic(inla.trim.family(model), (mm[names(mm) ==  :

  Model 'bym2' in section 'latent' is marked as 'experimental'; changes may appear at any time.

  Use this model with extra care!!! Further warnings are disabled.

2: running command '"C:/Users/xiawang/Documents/R/win-library/3.3/INLA/bin/windows/64bit/inla.exe"  -b -s -v  "inla.model/Model.ini"' had status 5 

Best,

Tracey