About the "besag2" model

[Olivier Nuñez]

Dear all,

I'm trying to fit a shared component model.

According to the INLA documentation, the "besag2" model defines a model with length N = 2*graph (like the "bym" one).

This point is quite clear for me.

But, I do not understand the definition of  the region index "i" given in the example: 

i=1:N;formula = y ~ f(i, model="besag2", graph=g)

Actually, the graph "g" contain only n=N/2 different regions!

If the previous definition is correct, what would be the meaning of the following model?

i=rep(1:n,2);formula = y ~ f(i, model="besag2", graph=g) 

Thanks.

[INLA help]

Hi

the besag2, defines a model for x = c(a*z, z/a), where z is a besag
model. it means that the graph is defined on z, with length n say, while
the length of x is 2*N. so indexs 1:n will extract a*z, while indexs
1:n+n will extract z/a.

Best
H

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

[Olivier Nuñez]

Ok.

Many thanks.

Best regards.

[Olivier Nuñez]

Dear Håvard,

Another question: In my experience, whether the random effects have zero mean or not, seems critical when estimating the weight parameter "a" in the besag2 model with INLA.

Actually, in the example given in the documentation, the constant 1 is added to the structured effects obtained from the bym fit:

## simulate two new datasets. 'a' is the weighting between the log.rel.risk:

a = 2

xx = x[1:n]+1 # why the constant 1 is added??

x = c(a*xx, xx/a) 

Then, the besag2 model is fitted without constant:

formula = y ~ f(i, model="besag2", graph=g) -1

This way to model the random effects is justified by an identifiability issue?

Can you please give more details about this point.

Thanks. Best regards.

[INLA help]

Hi,

the data is much more informative about the 'a' parameter when there is
a difference in the log.rel.risk. I think that was the motivation for
the '+1'.

H

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

[Olivier Nuñez]

It will be because it's easiest to estimate "a" from the mean of the vector z=c(a*xx, xx/a) than from its variance components? But, when the mean of z is near zero the INLA estimate of "a" obtained from a besag2 model seems to be quite biased. Following the example without the "+1" and choosing a=10, I obtain:

> set.seed(1000)

> data(Oral)

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

> ## use data Oral to estimate a spatial field in order to simulate a

> ## 'realistic' dataset.

> formula = Y ~ f(region, model="bym", graph=g)

> result = inla(formula, data = Oral, family = "poisson", E = E)

> x = result$summary.random$region$mean

> n = length(x)/2

> ## simulate two new datasets. 'a' is the weighting between the

> ## log.rel.risk:

> a = 10

> xx = x[1:n]

> z = c(a*xx, xx/a)

> summary(z)

     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 

-5.316000 -0.165200 -0.003621  0.013430  0.038230  5.156000 

> 

> E = c(Oral$E, Oral$E)

> N = 2*n

> y = rpois(N, lambda = E*exp(z))

> 

> ## model='besag2' defines a model with length N = 2*graph->n, the

> ## first half is weighted with 'a' the other half is weighted with

> ## 1/a. here there is no unstructed terms.

> i = 1:N

> formula = y ~ f(i, model="besag2", graph=g) -1

> r = inla(formula, family = "poisson", data = data.frame(E,y,i), E=E)

> summary(r)

Call:

c("inla(formula = formula, family = \"poisson\", data = data.frame(E, ",  "    y, i), E = E)")

Time used:

 Pre-processing    Running inla Post-processing           Total 

         0.4056          1.3104          0.0936          1.8096 

The model has no fixed effects

Random effects:

Name      Model

 i   Besags ICAR model for joint models 

Model hyperparameters:

                mean   sd     0.025quant 0.5quant 0.975quant mode  

Precision for i 4.6888 0.3733 3.9960     4.6747   5.4624     4.6475

a for i         3.7095 0.0752 3.5647     3.7085   3.8598     3.7062

Expected number of effective parameters(std dev): 444.47(3.161)

Number of equivalent replicates : 2.448 

Marginal Likelihood:  -10364.36 

  

[INLA help]

On Mon, 2015-04-13 at 07:16 -0700, Olivier Nuñez wrote:
>
> It will be because it's easiest to estimate "a" from the mean of the
> vector z=c(a*xx, xx/a) than from its variance components? But, when
> the mean of z is near zero the INLA estimate of "a" obtained from a
> besag2 model seems to be quite biased. Following the example without
> the "+1" and choosing a=10, I obtain:

the concept of 'bias' is not relevant for us, we just want to get as
close as we can to the Bayes-solution....


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

[Olivier Nuñez]

You're right, the concept of bias is not appropriate in this context.

I was just surprised that the credible interval was so far from the value of "a" used in the simulation.

Thank you very much for this clarification.

[Olivier Nuñez]

Hi again, 

actually, the shared component model I would like to estimate was the following:

y1 = rpois(n, E*exp(z) )

y2 = rpois(n, E*exp(beta*x + a*z) )

where "x" is a vector of covariates,

"z" follows a Besag model with graph=G (length(G)=n)

and "a" is a weight-parameter.

I used the following inla-model for its estimation, which seems to work fine in the simulations.

Y=c(y1,y2)

E=c(E,E)

X = c(rep(NA,n),x)

i = c(1:n,rep(NA,n))

j= c(rep(NA,n),1:n)

formule = Y ~ -1 + X +

f(i,model="besag",graph=G) +

f(j,copy="i",fixed=FALSE,range=c(0,1))

fit<-inla(formule,family=c("poisson"),data=data.frame(Y,X,E,i,j),E=E)

It is assumed that the weight-parameter "a" belongs to the interval (0,1).

I used the "range" option to specify this prior.

Is this the best way to express this assumption?

Thanks. Best regards.

[INLA help]

you're doing everything correct ;-)

there is an alternative model, 'besag2', which code this as

a*z
and
1/a*z

in case you find that more natural.

[Olivier Nuñez]

It was my first approach! 

But, for some reason the Besag2 model did not provide good results in my simulations, when "z" has a mean near zero.

I had to do something wrong..... maybe with the prior.

Anyway, the above model works fine and I'm glad to know that it is correct.

Thanks for your feedback. Best regards.