Random effect in posterior predictive distribution

200 views
Skip to first unread message

Richard

unread,
May 9, 2024, 5:09:33 AM5/9/24
to R-inla discussion group
Dear INLA community,

I am trying to use the function of INLA inla.posterior.sample and inla.posterior.sample.eval as I am interested to simulate the answers of the model.
I am having trouble adding the spde random effects, I have already read the vignette that INLA has for this topics and other answers in this group but it still unclear for me how to do it. Here is an example with artificial data replicating what I have been trying to do

Asume the mesh and the other processes are already completed

formula <- y ~ -1 + iIntercept + x1 + x2 + f(s, model = spde) 
#s is this inla.spde.make.index("s", spde$n.spde)

res <- inla(formula, family = "lognormal",
            data = inla.stack.data(stk.full),
            control.predictor = list(A = inla.stack.A(stk.full), compute = TRUE),
            control.compute = list(cpo = TRUE, dic = TRUE, waic = TRUE,
                                   return.marginals.predictor = TRUE,
                                   config = TRUE),
            control.fixed=list(mean = 0,prec = 0.01,
                               mean.intercept = 0,prec.intercept = 0.01),
            control.family = list(hyper=list(theta1 = list(
              prior = "normal",
              param = c(0.8, 0.3)))))


samples <- inla.posterior.sample(100, res)

#here is where I am starting to have problems as I am not sure how to introduce the #random effect term I defined previously in fun2 

fun2 <- function() return (b0 + ind*grid_huelva$Industria + mar*grid_huelva$Marsh + agric*grid_huelva$Agric + refi*grid_huelva$Refineria + urban*grid_huelva$Urbano + phospho*grid_huelva$Phospho  )

inla.posterior.sample.eval(fun2, samples)

Does anyone have any idea how can I do this? To have the posterior predictive distribution take on account f(s, model = spde)

Elias T. Krainski

unread,
May 10, 2024, 9:01:29 AM5/10/24
to R-inla discussion group
You need to have something like

Ap <- inla.spde.make.A(
    mesh = mesh,
    loc = cbind(pred.data$x, pred.data$y)
)

fun2 <- function()
    return(
        iIntercept +
        x1 * pred.data$x1 +
        Ap %*% s
     )

PS.: You can consider the inlabru package that facilitates working with SPDE models in INLA.

best regards,
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/9e7e8be6-9f13-47a6-9f0f-2ec82c8bb4fdn%40googlegroups.com.

Richard

unread,
May 10, 2024, 4:57:24 PM5/10/24
to R-inla discussion group
Thank you it worked perfect. Just a question about this:

Ap %*% s 

If I understand correctly this is what compute the spatial random effect at each location where one is making predictions but I also see that inside inla.posterior.sample object result there are hyperparameters for the Matern covariance function and I have as many of those parameters as the numbers of samples I selected so my question is

The spatial random effect Ap %*% uses the hyperparameters of the first sample and make predictions, that is one sample of the posterior predictive, then it uses the hyperparameters of the second sample and make predictions, that is another sample which had a covariance function with different parameters than the first one and so on, is that correct? Or are the hyperparametrs of the covariance function the same for all the samples in the simulation of the predictions?

Elias T. Krainski

unread,
May 12, 2024, 2:46:36 AM5/12/24
to R-inla discussion group
If you want to work with samples of hyperparameters, it is better to use inla.hyperpar.sample() instead.

The inla.posterior.sample() sample hyperparameters accounting for the int.strategy you have set to fit the model. Given these samples, then sample the latent variables. Therefore, with int.strategy = "ccd" or int.strategy = "grid" you will have samples with hyperparameters "integrated out". The result of this function is a list and the i-th element of such list contains the sample for the hyperparamters (h_i) and the latent variables sampled give h_i.

best regards,
Elias

Elias T. Krainski

unread,
May 12, 2024, 2:46:44 AM5/12/24
to R-inla discussion group
*given h_i.

Richard

unread,
May 12, 2024, 4:37:29 AM5/12/24
to R-inla discussion group
Thank you Elias all your answer were really helpful.
Reply all
Reply to author
Forward
0 new messages