Failing to get a simple spde model to work
541 views
Skip to first unread message
Bob O'Hara
Feb 5, 2014, 11:09:38 AM2/5/14
to r-inla-disc...@googlegroups.com
I've spent the last 2 days swaering at my computer whilst trying to get a simple spatial model to work. All I want is this model:
formula = Y ~ -1 + intercept + Covariate + f(ugh, model=spde)
which is even easier than the examples given in the tutorials, which is a problem because the extra complexity makes it more difficult to work out what is going on. Anywax, when I run the model with this:
inla.result = inla(formula, family="normal", data=inla.stack.data(Aagh.stack,spde=Spde))
I get the following error message:
Error in eval(expr, envir, enclos) : object 'intercept' not found
My full example code is below. Can anyone explain the error (and also why it has to be so complicated!)?
Thanks
Bob
(BTW, in the tutorials and code, calling an object spde is confusing, because model=spde and spde=spde refer to different things)
# Simulate some data...
N=100
Y=rnorm(N)
X.coord=runif(N)
Y.coord=runif(N)
Covariate=rnorm(N)
# create mesh
Mesh=inla.mesh.2d(cbind(X.coord, Y.coord), max.edge=0.2)
# plot(Mesh); points(X.coord, Y.coord, col=2)
# create SPDE object
Spde = inla.spde2.matern(Mesh)
# Create (sparse) matrix to represent spatial model
A.Aagh <- inla.spde.make.A(Mesh, loc=cbind(X.coord,Y.coord), index=1:N)
# Create an index (why? Can't I just pass an index from A.Aagh to hte effects list in inla.stack?)
Mesh.index=inla.spde.make.index(name="ugh", n.spde=Mesh$n, n.repl=1)
# stack the data
Aagh.stack=inla.stack(data=list(Y=Y), A=list(A.Aagh, 1),
effects=list(c(Mesh.index, list(intercept=1)), list(Covariate=Covariate)),
tag="AArg")
# (fail to) fit the model
formula = Y ~ -1 + intercept + Covariate + f(ugh, model=spde)
inla.result = inla(formula, family="normal", data=inla.stack.data(Aagh.stack,spde=Spde))
formula = Y ~ -1 + intercept + Covariate + f(ugh, model=spde)
which is even easier than the examples given in the tutorials, which is a problem because the extra complexity makes it more difficult to work out what is going on. Anywax, when I run the model with this:
inla.result = inla(formula, family="normal", data=inla.stack.data(Aagh.stack,spde=Spde))
I get the following error message:
Error in eval(expr, envir, enclos) : object 'intercept' not found
My full example code is below. Can anyone explain the error (and also why it has to be so complicated!)?
Thanks
Bob
(BTW, in the tutorials and code, calling an object spde is confusing, because model=spde and spde=spde refer to different things)
# Simulate some data...
N=100
Y=rnorm(N)
X.coord=runif(N)
Y.coord=runif(N)
Covariate=rnorm(N)
# create mesh
Mesh=inla.mesh.2d(cbind(X.coord, Y.coord), max.edge=0.2)
# plot(Mesh); points(X.coord, Y.coord, col=2)
# create SPDE object
Spde = inla.spde2.matern(Mesh)
# Create (sparse) matrix to represent spatial model
A.Aagh <- inla.spde.make.A(Mesh, loc=cbind(X.coord,Y.coord), index=1:N)
# Create an index (why? Can't I just pass an index from A.Aagh to hte effects list in inla.stack?)
Mesh.index=inla.spde.make.index(name="ugh", n.spde=Mesh$n, n.repl=1)
# stack the data
Aagh.stack=inla.stack(data=list(Y=Y), A=list(A.Aagh, 1),
effects=list(c(Mesh.index, list(intercept=1)), list(Covariate=Covariate)),
tag="AArg")
# (fail to) fit the model
formula = Y ~ -1 + intercept + Covariate + f(ugh, model=spde)
inla.result = inla(formula, family="normal", data=inla.stack.data(Aagh.stack,spde=Spde))
Elias T Krainski
Feb 5, 2014, 11:42:10 AM2/5/14
to r-inla-disc...@googlegroups.com
Hi Bob,
you need to set the control.predictor into the inla() call
control.predictor=list(A=inla.stack.A(Aagh.stack))
It is not easy to make a balance between simplification and to including
all possibilities. It maybe improved on forward versions. We wrote a
smaller document where the main purpose is the simplification. It is at
http://www.math.ntnu.no/inla/r-inla.org/tutorials/spde/inla-spde-howto.pdf
best,
Elias.
you need to set the control.predictor into the inla() call
control.predictor=list(A=inla.stack.A(Aagh.stack))
It is not easy to make a balance between simplification and to including
all possibilities. It maybe improved on forward versions. We wrote a
smaller document where the main purpose is the simplification. It is at
http://www.math.ntnu.no/inla/r-inla.org/tutorials/spde/inla-spde-howto.pdf
best,
Elias.
Finn Lindgren
Feb 5, 2014, 11:48:51 AM2/5/14
to Bob O'Hara, r-inla-disc...@googlegroups.com
Elias reply should solve the problem.
> (BTW, in the tutorials and code, calling an object spde is confusing,
> because model=spde and spde=spde refer to different things)
I don't see what you mean here; in both cases two different function
parameters (model and spde) are set to the same object, which is
called spde. Avoiding using names for variables that coincide with
names of function parameters would be an impossible task.
In the f() function, "model" specifies any of a large number of
models, usually as a string, but here as a special object.
For the functions dealing specifically with spde:s on the other hand,
it makes sense to call the function parameter spde.
Finn
> --
> 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 an email to
> r-inla-disc...@googlegroups.com.
> Visit this group at http://groups.google.com/group/r-inla-discussion-group.
> For more options, visit https://groups.google.com/groups/opt_out.
--
Finn Lindgren
email: finn.l...@gmail.com
> (BTW, in the tutorials and code, calling an object spde is confusing,
> because model=spde and spde=spde refer to different things)
parameters (model and spde) are set to the same object, which is
called spde. Avoiding using names for variables that coincide with
names of function parameters would be an impossible task.
In the f() function, "model" specifies any of a large number of
models, usually as a string, but here as a special object.
For the functions dealing specifically with spde:s on the other hand,
it makes sense to call the function parameter spde.
Finn
> 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 an email to
> r-inla-disc...@googlegroups.com.
> Visit this group at http://groups.google.com/group/r-inla-discussion-group.
> For more options, visit https://groups.google.com/groups/opt_out.
--
Finn Lindgren
email: finn.l...@gmail.com
Bob O'Hara
Feb 6, 2014, 9:56:58 AM2/6/14
to Elias T Krainski, r-inla-disc...@googlegroups.com
Wonderful, thanks! It's not clear to me why I need a control.predictor, but I'll work that out later.
I appreciate the balance is difficult to find, and I guess it's in different places for different people. The smaller document is very useful, because it's simpler - now I might find it easier to follow the other tutorials.--
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-group+unsub...@googlegroups.com.
To post to this group, send an email to r-inla-discussion-group@googlegroups.com.
Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/groups/opt_out.
--
Bob O'Hara
Biodiversity and Climate Research Centre
Senckenberganlage 25
D-60325 Frankfurt am Main,
Germany
Tel: +49 69 798 40226
Mobile: +49 1515 888 5440
WWW: http://www.bik-f.de/root/index.php?page_id=219
Blog: http://occamstypewriter.org/boboh/
Journal of Negative Results - EEB: www.jnr-eeb.org
Bob O'Hara
Feb 6, 2014, 9:57:12 AM2/6/14
to Finn Lindgren, r-inla-disc...@googlegroups.com
Ah, I've got it. Thanks!
I read model=spde as saying that the model was an spde, but I guess that should have been model="spde". If I've got it right, model=spde sets model to be the object called spde produced by the inla.stack.data() function, and this happens to be the same thing as the object already created called spde. Is this right? And this it's why if I have already created Spde, this works:
I read model=spde as saying that the model was an spde, but I guess that should have been model="spde". If I've got it right, model=spde sets model to be the object called spde produced by the inla.stack.data() function, and this happens to be the same thing as the object already created called spde. Is this right? And this it's why if I have already created Spde, this works:
formula=y ~ -1 + intercept + covar + f(field, model = Spde2, replicate = field.repl)
inla.result =inla(formula, data=inla.stack.data(stack,Spde2=Spde), family="normal", control.predictor=list(A=inla.stack.A(stack),compute=TRUE))
Finn Lindgren
Feb 6, 2014, 10:00:58 AM2/6/14
to Bob O'Hara, r-inla-disc...@googlegroups.com
On 06/02/14 14:57, Bob O'Hara wrote:
> Ah, I've got it. Thanks!
>
> I read model=spde as saying that the model was an spde, but I guess that
> should have been model="spde". If I've got it right, model=spde sets
> model to be the object called spde produced by the inla.stack.data()
> function, and this happens to be the same thing as the object already
> created called spde. Is this right? And this it's why if I have already
> created Spde, this works:
>
> formula=y ~ -1 + intercept + covar + f(field, model = Spde2, replicate
> = field.repl)
> inla.result =inla(formula, data=inla.stack.data(stack,Spde2=Spde),
> family="normal", control.predictor=list(A=inla.stack.A(stack),compute=TRUE))
Yes, precisely! For simpler models, one just specifies the name of the
> Ah, I've got it. Thanks!
>
> I read model=spde as saying that the model was an spde, but I guess that
> should have been model="spde". If I've got it right, model=spde sets
> model to be the object called spde produced by the inla.stack.data()
> function, and this happens to be the same thing as the object already
> created called spde. Is this right? And this it's why if I have already
> created Spde, this works:
>
> formula=y ~ -1 + intercept + covar + f(field, model = Spde2, replicate
> = field.repl)
> inla.result =inla(formula, data=inla.stack.data(stack,Spde2=Spde),
> family="normal", control.predictor=list(A=inla.stack.A(stack),compute=TRUE))
_type_ of model, and inla fills in the rest later, but since the spde
models need more preprocessing they are supplied as a whole object instead.
Finn
> <mailto:r-inla-disc...@googlegroups.com>.
> > Visit this group at
> http://groups.google.com/group/r-inla-discussion-group.
> > For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
> --
> Finn Lindgren
> email: finn.l...@gmail.com <mailto:finn.l...@gmail.com>
> http://groups.google.com/group/r-inla-discussion-group.
> > For more options, visit https://groups.google.com/groups/opt_out.
>
>
>
> --
> Finn Lindgren
>
>
>
>
> --
> Bob O'Hara
>
> Biodiversity and Climate Research Centre
> Senckenberganlage 25
> D-60325 Frankfurt am Main,
> Germany
>
> Tel: +49 69 798 40226
> Mobile: +49 1515 888 5440
> WWW: http://www.bik-f.de/root/index.php?page_id=219
> Blog: http://occamstypewriter.org/boboh/
> Journal of Negative Results - EEB: www.jnr-eeb.org <http://www.jnr-eeb.org>
>
>
>
> --
> Bob O'Hara
>
> Biodiversity and Climate Research Centre
> Senckenberganlage 25
> D-60325 Frankfurt am Main,
> Germany
>
> Tel: +49 69 798 40226
> Mobile: +49 1515 888 5440
> WWW: http://www.bik-f.de/root/index.php?page_id=219
> Blog: http://occamstypewriter.org/boboh/
Reply all
Reply to author
Forward
0 new messages