group option

878 views
Skip to first unread message

Olivier Nuñez

unread,
Feb 20, 2014, 6:27:12 AM2/20/14
to r-inla-disc...@googlegroups.com
Dear all,

I'm trying to fit a spatio-temporal model to count data and have some doubts about the corresponding INLA syntax.
The mathematical model is the following

For county i and period t, it is assumed that 
y_it ~ Poisson(lamda_it),

and log(lamda_it) = mu + u_i + v_t + w_it , 

where the latent variables are respectively

u_i ~ BYM 

v_t ~ RW1 (first order random walk)

and  for each i, 

w_it ~ RW1 with precision tau_i .

The corresponding code I try was 

formula = counts~1+f(i,model="bym", graph.file="mygraph")+ f(t,model="rw1") +
f(i,model="bym",group=t, control.group=list(model="rw1"))

fit = inla(formula,family="poisson",data=df)
         
I guess this syntax is not correct, because I get a lot of Inf value in the model predictions.

Thanks, Olivier.

INLA help

unread,
Feb 20, 2014, 7:51:40 AM2/20/14
to Olivier Nuñez, r-inla-disc...@googlegroups.com
Hi,

seems like you have unidentifiability issues in the model... also you
cannot have two equal arguments to f(), so do

i2 = i


formula = counts~1+f(i,model="bym", graph.file="mygraph")+
f(t,model="rw1") +
f(i2,model="bym",group=t, control.group=list(model="rw1"))

you should have got a error about this, right? and an error saying that
you need to define a graph for model 'bym' in the last entry?

also, what you do is not what you write,


w_it ~ RW1 with precision tau_i .

is this for each 'i' a RW1 in time 't'?

H


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

Olivier Nuñez

unread,
Feb 20, 2014, 9:14:15 AM2/20/14
to r-inla-disc...@googlegroups.com, Olivier Nuñez, he...@r-inla.org
Dear Håvard,

you're right.While simplifying notations, I forgot some terms in the code.
A more accurate version of my model is actually:

formula = counts ~ 1 +  f(i,model="bym", graph.file="mygraph") +  f(t,model="rw1") +
f(i2,model="bym", graph.file="mygraph",group=t, control.group=list(model="rw1"))

What I try to do, is to fit a RW1 process for each county taking into account the spatial dependence. It is likely that my interpretation of the "group" feature is not correct. Is the following syntax more accurate?

formula = counts~1 + f(i,model="bym", graph.file="mygraph") + f(t,model="rw1") +
f(t2,model="rw1",group=i, control.group=list(model="besag",graph.file="mygraph"))



Best regards.

INLA help

unread,
Feb 21, 2014, 1:01:37 AM2/21/14
to Olivier Nuñez, r-inla-disc...@googlegroups.com
On Thu, 2014-02-20 at 06:14 -0800, Olivier Nuñez wrote:
> Dear Håvard,
>
>
> you're right.While simplifying notations, I forgot some terms in the
> code.
> A more accurate version of my model is actually:
>
>
> formula = counts ~ 1 + f(i,model="bym", graph.file="mygraph") +
> f(t,model="rw1") +
> f(i2,model="bym", graph.file="mygraph",group=t,
> control.group=list(model="rw1"))
>
>
> What I try to do, is to fit a RW1 process for each county taking into
> account the spatial dependence. It is likely that my interpretation of
> the "group" feature is not correct. Is the following syntax more
> accurate?
>
>
> formula = counts~1 + f(i,model="bym", graph.file="mygraph") +
> f(t,model="rw1") +
> f(t2,model="rw1",group=i,
> control.group=list(model="besag",graph.file="mygraph"))

I attach some slides about the group option. both your options hsould be
the same, as 'a' grouped with 'b' is the same as 'b' grouped with 'a'.

the problem is that you will have a confounding



formula = counts ~ 1 + f(i,model="bym", graph.file="mygraph") +
f(t,model="rw1") +
f(i2,model="bym", graph.file="mygraph",group=t,
control.group=list(model="rw1"))

between f(i,...) and f(i2,...), so I suggest you take f(i,...) away. if
you could add linear constraints so that i2 sums to zero over all time
points for each index in i2, this would go away, but that is not a good
idea really.


Also, I think I would go away from the 'bym' model as this is a joint
model between the structured and unstructured term, and although you can
'group' with with 'rw1', it is questionable to connect the intrinsic
model and the iid model component in the same way in time.
you might just want to replace 'bym' with 'besag' ???

In any case, I would suggest to run a small simulation study where you
simulate from the model and estimate it back, to verify that all is
correct.

Best
group-models.pdf

Olivier Nuñez

unread,
Feb 21, 2014, 4:21:41 AM2/21/14
to r-inla-disc...@googlegroups.com, Olivier Nuñez, he...@r-inla.org
Dear Håvard,

I will simplify my model in the way you suggest.
Thank you very much for your comments and the slides, they are helpful.
Best regards. Olivier 
Reply all
Reply to author
Forward
0 new messages