can't figure out how to plot spatial random fields

[cp...@terpmail.umd.edu]

Hi!

I have a Bernoulli space-time model and am trying to plot the spatial random fields. However, I keep getting stuck.

Here is my code:

#I have 3 years of data

NYear <- length(unique(datafile.final2$YearSeq))

#create the mesh

mesh2 <- inla.mesh.create.helper(points=cbind(datafile.final2$x, datafile.final2$y),max.edge=c(20000,40000), cutoff=4000)

coords <- as.matrix(datafile.final2[,1:2])

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

A <- inla.spde.make.A(mesh2, loc=coords,group=datafile.final2$YearSeq,n.group = NYear)

time.index <- inla.spde.make.index("w", n.spde=spde$n.spde,n.group = NYear)

N <- nrow(datafile.final2)

X <- data.frame(Intercept=rep(1,N), datafile.final2$Phos.std,datafile.final2$salinity.std,datafile.final2$total.area.std, datafile.final2$distance.to.edge.std,datafile.final2$distance.between.beds.std,datafile.final2$bath.std)

colnames(X) <- c("Intercept","Phos.std","salinity.std","total.area.std", "distance.to.edge.std","distance.between.beds.std", "bath.std")

stack3 <- inla.stack(data=list(y=datafile.final2$pres, link=1,Ntrials=1),A=list(A,1),effects=list(w=time.index, X=as.data.frame(X)), tag="est")

hyper <- list(theta = list(prior="pc.prec", param=c(1,0.01)))

formula3 <- y ~ -1 + Intercept + Phos.std+salinity.std+total.area.std + distance.to.edge.std+distance.between.beds.std + bath.std+f(w, model=spde,group=w.group,control.group = list(model="ar1",scale.model = T),hyper=hyper)

data = inla.stack.data(stack3)

r3 <- inla(formula3, control.inla=list(reordering="metis"),family="binomial",Ntrials=1, data=data, control.compute=list(dic=T,cpo=T,config=T),control.predictor=list(A=inla.stack.A(stack3), compute=T, link=1, quantiles=c(0.025, 0.5, 0.975)),control.fixed=list(prec.intercept = 0.1))

#my attempt at prediction of the random field

w.pm <- r3$summary.random$coordinates$mean

projgrid <- inla.mesh.projector(mesh2, xlim=range(datafile.final2$x),ylim=range(datafile.final2$y))

#for loop for 3 years

xmean <- list()

for (j in 1:3)

  xmean[[j]] <- inla.mesh.project(

    projgrid, r3$summary.random$w$mean[time.index$w.group == j])

I'm just not sure where to go from here to plot the spatial fields!

I have no trouble getting the overall predictions by year; I would like to generate a plot something like this (see attached; which brings me to another question-is there a way to convert this plot from what I assume are log odds to probabilities?):

#overall predictions

idx <- inla.stack.data(stack3, tag = "pred")

preds <- r3$summary.linear.predictor

datafile.final2$pred <- preds$mean[1:nrow(datafile.final2)]

lattice::levelplot(pred ~ x + y|YearSeq, data = datafile.final2)

Thank you very much!!

-Carrie

[Elias T. Krainski]

Hi,

You can consider the examples from Chapter 7 of the SPDE book. The links for the online version and the code is available at https://www.r-inla.org/learnmore/books

Elias

--
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 view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/6412a3a9-9cc2-4873-b951-7f9d88630a7bn%40googlegroups.com.

[Carrie Perkins]

Thanks for your message, Elias!

That was actually the code I tried first. I get something really weird (attached) that looks nothing like my map. Here is my code, based on Chapter 7:

stepsize <- 10000 * 1 / 111
x.range <- diff(range(datafile.final2$x))
y.range <- diff(range(datafile.final2$y))
nxy <- round(c(x.range, y.range) / stepsize)

projgrid <- inla.mesh.projector(mesh2, xlim=range(datafile.final2$x),ylim=range(datafile.final2$y), dims = nxy)

xmean <- list()
for (j in 1:3)
  xmean[[j]] <- inla.mesh.project(
    projgrid, r3$summary.random$w$mean[time.index$w.group == j])

library(splancs)
xy.in <- inout(projgrid$lattice$loc,
               cbind(datafile.final2$x, datafile.final2$y))

par(mfrow = c(4,3), mar = c(0, 0, 0, 0))

for (j in 1:3) {

  xmean[[j]][!xy.in] <- NA
  book.plot.field(list(x = projgrid$x, y = projgrid$y, z = xmean[[j]]),
                  zlim = round(range(unlist(xmean), na.rm = TRUE), 1) )
}

What do you think I'm doing wrong? 

[Elias T. Krainski]

seems that you considered the locations coordinates instead of a polygon around. Can you check how it goes without the following line

  xmean[[j]][!xy.in] <- NA

[Carrie Perkins]

Yes, all the locations are coordinates, each representing a 100X100m raster pixel. When I try removing that line I get the attached image (Attempt 2). 

I'm not really sure where to go from here; I ended up trying a bunch of other things as well, such as creating a boundary raster from my coordinates and trying to use it the way they use the PRborder polygon in chapter 7, but no luck (Attempt 3).

[Carrie Perkins]

So far this is the best I can do (used a different approach):

gproj <- inla.mesh.projector(mesh2,  dims = c(300, 300))

xmean <- list()
for (j in 1:3)
  xmean[[j]] <- inla.mesh.project(

    gproj, r3$summary.random$w$mean[time.index$w.group == j])


grid.arrange(levelplot(xmean[[1]], scales=list(draw=F), xlab='', ylab='', main='mean',col.regions = heat.colors(16)),
             levelplot(xmean[[2]], scal=list(draw=F), xla='', yla='', main='mean' ,col.regions = heat.colors(16)),
             levelplot(xmean[[3]], scal=list(draw=F), xla='', yla='', main='mean' ,col.regions = heat.colors(16)))

Let me know if you think of something else I can try to get the other way to work! Still would love to limit the predictions to the map if possible.

[Finn Lindgren]

Hi,

There are tools for masking predictions to polygons and plotting with ggplot in the inlabru package. The examples we have are based on having used inlabru itself to do the predictions, but it also works with more generic objects. I’ll get back to you with an example.

Finn

On 30 Jul 2021, at 06:52, Carrie Perkins <cp...@terpmail.umd.edu> wrote:



To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/CAPtr_T6NUoi0N0g9szuQ2jhcTQvfV_4YvMdH_eWf%3D0Jf_SHp_w%40mail.gmail.com.

<attempt5_jul292021.png>

[Finn Lindgren]

Hi,

Here's an example for using ggplot, adapted from https://inlabru-org.github.io/inlabru/articles/web/2d_lgcp_plotsampling.html

The masking of the region requires you have the boundary as a SpatialPolygons object.

Otherwise it will show the whole mesh.

First a basic example, then one where I've tried to use your variables.

# Load inlabru

library(inlabru)

library(ggplot2)

# Generic example, with colors taken from the mesh coordinates (just to have something to plot):

mesh<-inla.mesh.2d(cbind(0,0),offset = 1, max.edge=0.1)

# Alt 1: Let inlabru::gg handle the conversion from mesh to pixels:

ggplot()+gg(mesh, color=mesh$loc[,1])

# Alt 2: Do your own conversion to pixels:

pix<-pixels(mesh)

proj <- inla.mesh.projector(mesh, coordinates(pix))

pix[["A"]] <- inla.mesh.project(proj, mesh$loc[,1])

pix[["B"]] <- inla.mesh.project(proj, mesh$loc[,2])
ggplot()+gg(pix, aes(fill=A))
ggplot()+gg(pix, aes(fill=B))

# Both side-by side with correct aspect ratio and common legend:

library(patchwork)

(ggplot()+gg(mesh, color=mesh$loc[,1]) |

 ggplot()+gg(mesh, color=mesh$loc[,2])) +

  coord_fixed() +

  plot_layout(guides = "collect") &
  theme(legend.position = "right") &
  scale_fill_continuous()

# Alt 1:

ggplot() +

  gg(mesh2, color = r3$summary.random$w$mean[time.index$w.group == j],

     nx = 300, ny =300, mask = boundary)

# Alt 2:

pix <- pixels(mesh2, mask = boundary)

gproj <- inla.mesh.projector(mesh2, coordinates(pix))

for (j in 1:3) {

  pix[[paste0("mean", j)]] <- inla.mesh.project(

    gproj, r3$summary.random$w$mean[time.index$w.group == j])

}

ggplot() + gg(pix, aes(fill = mean1))

# All three side-by side with patchwork and common legend:

library(patchwork)

(ggplot() + gg(pix, aes(fill = mean1)) |

  ggplot() + gg(pix, aes(fill = mean2)) |

  ggplot() + gg(pix, aes(fill = mean3))) +

  coord_fixed() +

  plot_layout(guides = "collect") &
  theme(legend.position = "bottom") &
  scale_fill_continuous()

Finn

--

Finn Lindgren
email: finn.l...@gmail.com

[Carrie Perkins]

Thank you so much, Finn!!