Mapping spatial effect with SPDE

[Mark Myer]

Hello, 

I am conducting a spatiotemporal analysis of binomial data observed on an island using the SPDE approach, and I am running into a strange error when attempting to map the spatial effect across my study area. I have cobbled together this approach from several guides and tutorials online, so apologies if there are any glaring omissions. 

I've made my mesh from a polygon and a set of points, using unionSpatialPolygons and inla.sp2segment to turn the polygon into an inla.mesh.segment object for use in inla.mesh.create.helper. 

loc<-cbind(DataWithDaymet$Long,DataWithDaymet$Lat) #Stores the lat/longs of each observation site

shape<-readOGR(dsn="LocationOfLayer", layer=("SuffolkPolygonDecDeg"))

ob <- unionSpatialPolygons(shape, rep(1, nrow(shape))) #This simplifies and unifies the boundary for faster processing

bound<-inla.sp2segment(ob) #Converts the polygon boundary to a inla.mesh.segment object, which is a series of line segments

(mesh<-inla.mesh.create.helper(points=cbind(DataWithDaymet$Long,DataWithDaymet$Lat),points.domain=bound$loc,offset=c(0.01,0.15),max.edge=c(0.1,0.2)))$n

Here is what the shapefile looks like with the observation locations printed on it:

plot(mesh,asp=1)

lines(bound,lwd=2)

points(loc,cex=.75,pch=20,col=2)

Here is the mesh I am using, with the outline of the island my study is on superimposed in blue, and my observation points superimposed in red for reference:

However, when I run the analysis using the following code:

A <-inla.spde.make.A(mesh, loc)

spde<-inla.spde2.matern(mesh,alpha=2)

formula<-Result ~ 0 + Intercept + AREIS_SepticCount1km + NDVI + Developed..Low.Intensity + Emergent.Herbaceous.Wetlands + WkMeanPrcp + WkMeanTmax + f(spatial, model=spde) #Defining the formula and variables

vars<-DataWithDaymet[,c(11,12,15,27,49,50)] #Takes out from the full dataset only the variables of interest mentioned in the formula

vars<-scale(vars)

stk<-inla.stack(data=list(Result=DataWithDaymet$Result),A=list(A,1),effect=list(list(spatial=1:spde$n.spde),data.frame(Intercept=1,vars)))  #Creating the stack for easy specification of INLA

r<-inla(formula, family="binomial", Ntrials=10596, data=inla.stack.data(stk), control.family=list(link="logit"), control.predictor=list(A=inla.stack.A(stk)),control.compute=list(dic=TRUE,cpo=TRUE)) #Running the INLA and checking marginals

summary(r)

and then use the following code to try and plot the spatial effect mean and SD onto the projector grid.

#Plotting the response on a grid

library(lattice)

projection<-inla.mesh.projector(mesh, dims=c(500,500)) #projects the mesh onto a square grid

table(xy.in <- inout(projection$lattice$loc,cbind(bound$loc[,1], bound$loc[,2]))) #Removes the grid points that are outside the boundary

#These are for the spatial field random effect, not the response

prd.m<-inla.mesh.project(projection,r$summary.ran$spatial$mean) #Projects the spatial effect mean onto the grid

prd.s<-inla.mesh.project(projection,r$summary.ran$spatial$s) #Projects the spatial effect SD onto the grid

prd.m[!xy.in] <- prd.s[!xy.in] <- NA #sets the projections outside the boundary to "NA"

library(gridExtra)

levelplot(prd.m, col.regions=terrain.colors(16))

levelplot(prd.s, col.regions=terrain.colors(16))

The resulting graph looks like this. My outline is chopped-up and distorted, despite using the boundary shapefile "bound" shown above which, when plotted, shows a normal shape. What might be causing this?

Thank you very much, 

---Mark

[Finn Lindgren]

Hi,

I'm not sure if this is the only issue, but the following line is not guaranteed

to do what you want, and in many cases it will definitely not do what you want:

table(xy.in <- inout(projection$lattice$loc,cbind(bound$loc[,1], bound$loc[,2])))

1) The inla.mesh.segment objects do not necessarily store the vertices of the curve in order.
2) bound$idx contains the information about which order to go through the vertices

3) even if the point were in order in loc, the code above would fail for curves with more

     than one connected component, which looks like exactly what you have.

If you found code on r-inla.org using the above, please let me know where so I can

 locate and modify it! Code elsewhere doing that would also be good to know about 

even if I can't directly correct it.

To do this correctly you need to use your original SpatialPolygons object to do the 

checking.  I believe rgeos::gIntersection can give you the wanted subset of points.

There is an alternate but also correct approach that only needs the inla package,

 but there is no ready-made function for it yet: the idea is to build a simple mesh with

the boundary with no extension:

simplemesh=inla.mesh.2d(boundary=bounds, max.edge=1e9)

ok = inla.mesh.projector(simplemesh, locations)$proj$ok

locations[ok,,drop=FALSE]

Please let me know if that approach works or not. If it does, I'll add a wrapper for it.

In the coming days& weeks, there will also be more and more Spatial object support

directly built into the mesh functions (the development version already lets you use

SpatialPolygons instead of inla.mesh.segment objects in most places).

Finn L, meshes-in-inla-developer

--
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 https://groups.google.com/group/r-inla-discussion-group.
For more options, visit https://groups.google.com/d/optout.

[Finn Lindgren]

Further to my previous email:

This is a common use case, and filtering of this type for inla.mesh.projector is already

used internally in the latest development code, so I'll add a "mask" parameter to that

interface to make it even easier!

Finn

On 10 Aug 2016, at 18:58, Mark Myer <thereal...@gmail.com> wrote:

[Mark Myer]

Finn, 

Thank you very much for getting back to me so quickly. I'm trying to use your solution:

simplemesh=inla.mesh.2d(boundary=bound, max.edge=1e9)

ok = inla.mesh.projector(simplemesh, locations )$proj$ok

locations[ok,,drop=FALSE]

but I am not sure if I understand what I need to do. 

The simple mesh is created fine:

but running the next line, ok is simply a logical vector of all "TRUE". I assumed that locations refers to the coordinates of my observation points loc. Was this supposed to refer to something else? I also don't understand the next line, locations[ok,,drop=FALSE] . Should I store this and use it to set the locations outside of my inla.mesh.project grids to NA?

Apologies for the lack of understanding, I am new to R. 

---Mark

[Finn Lindgren]

I should have been more precise, sorry!

In your code, "locations" should be your original

"projection$lattice$loc", and "ok" should

be used like you use the "inout" result in your code.

My third line was just an example of such use.

Finn

[Mark Myer]

It works well! Thank you so much for your help. 

As for where I had gotten the original code where I used the inout() function, I found it on Page 46 of the SPDE Tutorial:

Again, thanks for your help and all the work you've done on the R-INLA project!

---Mark

[Finn Lindgren]

> On 10 Aug 2016, at 21:29, Mark Myer <thereal...@gmail.com> wrote:
> It works well! Thank you so much for your help.

Ok, good, and thanks! I'll add something similar to the package when I get the chance.

> As for where I had gotten the original code where I used the inout() function, I found it on Page 46 of the SPDE Tutorial:

Ah, yes, I found that code too, but the PBborder in that code wasn't constructed
in the same way, and probably most importantly that border is a single closed loop,
which your boundary is not.

Finn