LGCP with non-separable covariance function?

[Seth Flaxman]

Hi all. Is it possible to fit a log-Gaussian Cox Process with a non-separable covariance function in INLA for a spatiotemporal point process? I'm agnostic as to the parametric form, I just need something in which space and time aren't separable.

Best,

Seth

[Havard Rue]

Working on it.... Later this year.

--
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-disc...@googlegroups.com.
Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Finn Lindgren]

On 26 Mar 2015, at 20:44, Seth Flaxman <fla...@gmail.com> wrote:

Hi all. Is it possible to fit a log-Gaussian Cox Process with a non-separable covariance function in INLA for a spatiotemporal point process? I'm agnostic as to the parametric form, I just need something in which space and time aren't separable.

In addition to Håvard's reply:

You probably have some idea about what _kind_ of non-separability that you need?

We are working on _some_ types, but there are many possibilities.  E.g. spatial drift is a lot harder than non-separability without spatial drift.

Finn

[Seth Flaxman]

I'm using one of the forms in Gneiting 2002 at the moment, because
these have the advantage of being relatively simple, with a single
parameter beta controlling "space/time interaction" which can be 0
meaning none--very little is known in my application domain
(criminology) on the question we're trying to answer, so the first
step is descriptive statistics, model comparison between different
models, etc. Thus, I simply want to be able to compare a separable and
nonseparable model for my point process, keeping everything else the
same (so I don't want to adopt a mechanistic model).

I don't think I need spatial drift. Heat equations would be very neat,
but probably too fancy. What other types of nonseparability is/isn't
hard or in the works?

Thanks for the responses and very nice work on INLA by the way,
Seth

[Seth Flaxman]

Quick followup: for the moment I'm happy investigating a separable model since that's the baseline anyway. I followed the tutorial here and I have a spatial model working (impressive how easy that was!) 

Is there a tutorial or do you otherwise have pointers on adding time? (This is possible, right?) I don't see an inla.mesh.3d, for example.

Thanks!

Seth

[Daniel Simpson]

In the long tutorial there's a chapter on this (I think). There is also the Cameletti paper.

If you've implemented the spatial model with inla.stack and inla.spde.make.A then it's not a lot more work to include the temporal component.

D

[Seth Flaxman]

Hello. I'm getting back to this now, trying to follow the SPDE tutorial (Chapter 15, pp 104-106). I'm confused about the code and explanation on the top of page 105:

"The data stack can be made considering the ideas for the purerly spatial model. In

addition, we have to consider the temporal volume at each temporal node. As we have

equally spaced time knots and setted the time range equals 1, it is jus the one over the

number of time knots."

What are y and expected and what's the logic behind the code for creating them?

m <- spde$n.spde

n <- nrow(burkitt)

y <- rep(0:1, c(k * m, n))

st.volume <- diag(kronecker(Diagonal(n=k), spde$param.inla$M0))

expected <- c(st.volume, rep(0, n))

stk <- inla.stack(data=list(y=y, expect=expected),

A=list(rBind(Diagonal(n=k*m), Ast), 1), effects=list(idx, list(a0=rep(1, k*m + n))))

i.e. why is y length k * m + n with n ones and why is expected padded with n zeros? (Since I've got a mesh, I thought that inference would be cheap because I'd only be considering the nodes of the mesh, but n is the number of observed points...)

Also, there's code for Voronoi polygons in Chapter 12 for the spatial only LGCP. Do I not need to do that for some reason for the spatiotemporal LGCP?

Sorry I'm a little all over the place with these questions. Happy to create a self-contained example if that'll help!

Thanks,

Seth

[Samir Bhatt]

Hi Seth,

Someone else will know more clearly what is going on, but i have been working on the similar models so ill give you my 2 cents

The discretised spde mesh in space or space time (through the Kronecker), in a log Gaussian cox model, are considered integration points in INLA. So the area/volume of these mesh nodes is proportional to the expected number of events or the intensity.

So when paramaterising your y's -  you need zeros for your mesh nodes and ones for you actual data

But for the expect/E part where you have the areas/volumes of your mesh, you have a vector with the volumes of the mesh and zeros for your data locations

Finally your observation matrix has the standard projection from mesh to data, and also a diagonal matrix of 1's for the integration points of your mesh.

 Finn or Dan correct me if i am massively off?

> email to r-inla-discussion-group+unsub...@googlegroups.com.

> To post to this group, send 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/d/optout.

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

[Finn Lindgren]

Yes, Samir's explanation looks right to me.

One minor issue:

The example code in the earlier email used ...$M0 to obtain the local integration area for each mesh node.  That doesn't always give the right values, especially for non-default spde models.  It's much safer to use

  diag(inla.mesh.fem(mesh)$c0)

which always gives the local areas.  (For a kronecker space time model, you may need to do the kronecker for c0 as well.)

Finn

To unsubscribe from this group and stop receiving emails from it, send an email to r-inla-discussion...@googlegroups.com.

[Seth Flaxman]

Thanks guys, this is starting to make sense. 

But what's the implication for the overall runtime then, if y and expected are length actual data + # integration points? I've got a big dataset (n > 100,000) and 5000 integration points, and my call to inla(...) has been going for more than a day (or is hung). Should I expect to be able to fit a problem of this size?

Is there some other way of setting this up that should be more efficient (or some other reason things are really slow that I'm not aware of?) 

Also, does anyone happen to have example code of a Kronecker space/time LGCP (i.e. discretised spde mesh)?

Thanks,

Seth

> email to r-inla-discussion...@googlegroups.com.
> To post to this group, send 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/d/optout.

--
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-disc...@googlegroups.com.
Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

--
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-disc...@googlegroups.com.
Visit this group at http://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Elias T Krainski]

Hi Seth,

We've done an aggregation approach for n>2,000,000. We aggregate the
data into counts in each Voronoi polygons built from the mesh locations.
Time discrete as well. Counting using the Voronoi polygons was chosen
because it made the fitting process easier because the locations and
integration points are the same.

Elias

[Elias T Krainski]

On 26/08/15 03:19, Seth Flaxman wrote:
> Also, does anyone happen to have example code of a Kronecker
> space/time LGCP (i.e. discretised spde mesh)?

I have updated the spde-tutorial with improvements in the spatio
temporal point process example. Thanks for motivating it and starting
this thread :)

[Seth Flaxman]

Aggregation sounds like the way to go for my case. 

Anyone know if there's an efficient way to do this aggregation using the Voronoi polygons? Something like quadratcount in spatstat? See attached, based on Elias's tutorial, on line 49 I want y to be the number of points in each space-time cell.

Thanks,

Seth

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

[Elias T Krainski]

I just added an example with an aggregation approach.
Elias.

[Seth Flaxman]

OK, getting closer I think. I found a small bug in your code, which isn't a problem for your synthetic data, but will be an issue whenever a grid cell is empty:

agg.dat <- as.data.frame(table(area, time))

I think this fixes it because then table() knows about missing values:

time=ordered(time,1:(length(t.breaks)-1))

area=ordered(area,1:length(tiles))

agg.dat <- as.data.frame(table(area,time))

I'd like to fit a model that has a latent field decomposed as f1(s) + f2(t) + f12(s,t). Does this code seem right to do that?

stk.data=inla.stack.data(stk)

stk.data$s1 = stk.data$s

stk.data$s.group1 = stk.data$s.group

formula <- y ~ 0 + b0 +  f(s,model=spde) +

  f(s.group1, model='ar1') +

  f(s1, model=spde, group=s.group, control.group=list(model='ar1'))

Thanks!

Seth

On Mon, Aug 31, 2015 at 8:28 PM, Elias T Krainski <eliask...@gmail.com> wrote:

I just added an example with an aggregation approach.
Elias.

[Elias T Krainski]

On 25/09/15 18:15, Seth Flaxman wrote:

OK, getting closer I think. I found a small bug in your code, which isn't a problem for your synthetic data, but will be an issue whenever a grid cell is empty:

agg.dat <- as.data.frame(table(area, time))

I think this fixes it because then table() knows about missing values:

time=ordered(time,1:(length(t.breaks)-1))

area=ordered(area,1:length(tiles))

agg.dat <- as.data.frame(table(area,time))

in this case you forced to always guarantee to have
  nrow(agg.dat)==(length(t.breaks)*length(tiles))
thanks for noticing it! I've updated with this correction.

I'd like to fit a model that has a latent field decomposed as f1(s) + f2(t) + f12(s,t). Does this code seem right to do that?

yes. However, main effects + interactions needs some care...

Elias

[Seth Flaxman]

Any suggestions for dealing with main effects+interaction? So far I've had things crash/take a very, very long time to run with small datasets. Has anyone had luck with a two stage approach? I.e. fit main effects first, then fix those hyperparameters (or constrain then with informative priors) and fit interactions...

To post to this group, send email to r-inla-disc...@googlegroups.com.