Extend a spatial code to be appropriate

[Silvia Shawky]

Dear INLA experts,

I am running a spatial and space-time coregionalization model that are explained in chapter (8) in Advanced Spatial Modelling with SPDE using R and INLA. 

In the spatial coregionalization model I have used the code below to predict the mean response on the spatial domain and it works well. (I found the code in chapter 8 in Spatial and spatiotemporal bayesian models book. I need help to extend the code below to obtain the mean response on spatial domain and across time for a spatiotemporal coregionalization model.

#Plot the prediction data
s1 <- inla.posterior.sample(n=1, result = output.apr.NP)
names(s1[[1]])
grep("intercept", rownames(s1[[1]]$latent), fixed = TRUE)
ids <- lapply(c('intercept', 'u.field', 'x.field'),
              function(x) grep(x, rownames(s1[[1]]$latent),fixed = TRUE))
pred.y.f <- function(s)
  (s$latent[ids[[1]],1]+ s$latent[ids[[2]],1] +s$latent[ids[[3]],1])
s1000 <- inla.posterior.sample(1000, output.apr.NP)
prd.y <- sapply(s1000, pred.y.f)
prd <- list(y.mean = inla.mesh.project(projgrid, field = rowMeans(prd.y)))
prd$y.sd <- inla.mesh.project(projgrid, field = apply(prd.y,1,sd))
for (j in 1:2) prd[[j]][!xy.in] <- NA
par(mfcol=c(2,2), mar=c(1,1,1,1))
for (i in 1:2) {
  image.plot(list(x=projgrid$x, y=projgrid$y, z=prd[[i]]),
             asp=1, axes=FALSE, nlevel=101, col=viridis(10),
             legend.mar=7,legend.args=list(text=names(prd)[i], line=0))
  

Thanks in Advance 

 Sylvia Alber

[Elias T. Krainski]

Hi,

you can use the following code example (at the end it plots for one of the outcomes). But, notice that it E(y_0|y_observed), not actually in terms of p_y(y_0|y_observed). You can add the noise term afterworks.

### ADDED FROM HERE
### (with added control.compute=list(config=TRUE) in the previous inla(...) call)

### 1. samples from the joint posterior
system.time(s1 <- inla.posterior.sample(1, result))
system.time(ss <- inla.posterior.sample(1e3, result, add.names=FALSE))

### 2. identify the components of each linear predictor
tail(xn <- rownames(s1[[1]]$latent))

ii1 <- list(grep('intercept1', xn),
            grep('s1:', xn))
str(ii1)

ii2 <- list(grep('intercept2', xn),
            grep('s2:', xn),
            grep('s12:', xn))
str(ii2)

ii3 <- list(grep('intercept3', xn),
            grep('s3:', xn),
            grep('s13:', xn),
            grep('s23:', xn))
str(ii3)

### 1. define a (spatial) grid projector
nxy <- c(99, 101)
proj1grid <- inla.mesh.projector(mesh, dims=nxy)
np <- nrow(proj1grid$proj$A)

### 2. build the spacetime projector matrix
stloc0 <- kronecker(
    matrix(1, k, 1), ### repeat coords each time
    proj1grid$lattice$loc)
A0st <- inla.spde.make.A(mesh, stloc0, n.group = k,
  group = rep(1:k, each = np))

### 3. evaluate the sum of each linear predictor term
pred1s <- sapply(ss, function(x)
    drop(x$latent[ii1[[1]],1] + A0st%*%x$latent[ii1[[2]], ]))
pred2s <- sapply(ss, function(x)
    drop(x$latent[ii2[[1]],1] + A0st%*%x$latent[ii2[[2]], ] +
         A0st%*%x$latent[ii2[[3]], ]))
pred3s <- sapply(ss, function(x)
    drop(x$latent[ii3[[1]],1] + A0st%*%x$latent[ii3[[2]], ] +
         A0st%*%x$latent[ii3[[3]], ] +
         A0st%*%x$latent[ii3[[4]], ]))

### 4. compute summaries
pred1.m <- matrix(rowMeans(pred1s), np)
pred1.sd <- matrix(apply(pred1s, 1, sd), np)

### 5. visualize
par(mfrow=c(2,2), mar=c(1,1,0.1,0.1), mgp=c(2,0.5,0))
for(j in 1:4)
    image.plot(proj1grid$x, proj1grid$y, matrix(pred1.m[,j], nxy[1]))

par(mfrow=c(2,2), mar=c(1,1,0.1,0.1), mgp=c(2,0.5,0))
for(j in 1:4)
    image.plot(proj1grid$x, proj1grid$y, matrix(pred1.sd[,1], nxy[1]))

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/CAHABYd1hHEWqL8U%2BNP6KyjD6jOTsDSOFiAVtCjoyCn%2BBwWtm5A%40mail.gmail.com.

[Silvia Shawky]

Thank you, I will try it. 

[Silvia Shawky]

Dear Team

I have tried the code above and it gives me the error below.

pred1s <- sapply(ss, function(x)
    drop(x$latent[ii1[[1]],1] + A0st%*%x$latent[ii1[[2]], ]))

Error: cannot allocate vector of size 8.9 GB

pred2s <- sapply(ss, function(x)
    drop(x$latent[ii2[[1]],1] + A0st%*%x$latent[ii2[[2]], ] +
         A0st%*%x$latent[ii2[[3]], ]))

Error: cannot allocate vector of size 8.9 GB

Any suggestion to solve this problem?

[Elias T. Krainski]

It happened because your computer didn't have enough memory. You can try removing large (unused) objects or subsetting the latent field previously and then doing the matrix multiplication.

To view this discussion on the web, visit https://groups.google.com/d/msgid/r-inla-discussion-group/6ed4642b-1643-4b24-a421-173ae603c178n%40googlegroups.com.