About inla.posterior.sample and AR1 models

[Julia Quero Pérez]

Hello. 

I've got a spatiotemporal model with R-INLA and wanted to do some simulations of the prediction in space. I'm obtaining the parameter estimates with inla.posterior.sample and then build the linear predictor summing as needed.

The thing is, I've got an AR1 part of my model, and I'm not sure if I should remove the part regarding the first day or the first year. 

Here's the formula of my model:

formula.8.1<- AT ~  -1 + Intercept + HGHT +
  CoastDist*(cos(yearday/365)+sin(yearday/365)) +
  poly(year,degree=1,raw=TRUE):(CoastDist+HGHT) + poly(year,degree=1,raw=TRUE) +
  f(spatial,model=spde,
    group=spatial.group,control.group=list(model="iid")) +
  f(yearday,model="ar1",replicate = year)

Here's the code for the simulation of spatial predictions (my doubts are regarding the "Obtain linear predictor" paragraph, I have not removed any of the AR1 parameters)
# CCovGrid prediction covariables
# mmesh and ccoordsGrid needed to create projection matrix with tag="pred"
simulacion=function(CCovGrid=CovGrid, nnP=nP,nnT=nT, nnL=nL,
                    ccoordsGrid=coordsGrid,mmesh=mesh,groupidx=grp.idx,ngroup=ngr,
                    nSim=500,MSel=MSelT,mseed=NULL){
 
  #Defining projection matrix for defining the spatial random effects
  #Create Ap.s: matrix with nP*nT*nL rows and G=363 columns
  Ap.s<-inla.spde.make.A(mmesh, loc=ccoordsGrid,
                         group=groupidx,n.group=ngroup)
 
  # Take sample
  if (is.null(mseed)==F)
  {
    inla.seed = as.integer(runif(1)*.Machine$integer.max)
    set.seed(mseed)
    samples<-inla.posterior.sample(n=nSim, result=MSel, add.names=FALSE,
                                   seed = inla.seed, num.threads="1:1")
  }else
  {
    samples<-inla.posterior.sample(n=nSim, result=MSel, add.names=FALSE)
  }
 
  # Coger índices para obtener la matriz
  #  samples is a list with nSim elements, one for each simulation
  # each list has 3 elelments:names(samples[[1]]): "hyperpar" "latent"   "logdens"
  xnames <- rownames(samples[[1]]$latent) ### collect the names of each z_{tl}(s)
  # tags <- attr(samples, ".contents")$tag # obtén el nombre de las covariables y efectos
  idx <- lapply(c("Intercept","Cov","spatial","yearday"), function(nam)
    which(substr(xnames, 1, nchar(nam))==nam))
  idx[[2]]<-seq(idx[[1]]+1,length(samples[[1]]$latent),by=1) # covariables
  #idx is a list with 4 elements each of  different length: 1,9,23595=363*65, 9945=65*153
  # that give the indexes where those elements are in samples[[1]]$latent
 
  # Obtain parameter sample by coefficients
  mat.samples1 <- sapply(samples, function(spl)
    c(Intercept=spl$latent[idx[[1]]], Cov=spl$latent[idx[[2]]],
      spatial=spl$latent[idx[[3]]]))
  # yearday part
  mat.samplesaux <- sapply(samples, function(spl) yearday=spl$latent[idx[[4]]])
  rm(xnames,idx)
 
  # Obtain summands
  mu.samp<-cbind(Intercept=1, Cov=CCovGrid,s=Ap.s)%*%mat.samples1     # en matriz (ver abajo)
  #dim(mu.samp) #  nnL*nnT*nnP x nSim  92*61*844  x 500
  lmu.samp<-lapply(seq_len(nSim), function(i,mm) mm[,i], mm=mu.samp)  # en lista por columnas (cada simulación solo)
  rm(mu.samp)
  # yearday part
  lsamp2<-lapply(seq_len(nSim), function(i,mm) mm[,i], mm=mat.samplesaux)
 
 
  # Obtain linear predictor
  for (j in c(1:nSim))
  {
    for (i in seq(1,nnP*nnL*nnT, by=(nnL*nnT)))
    {
      lmu.samp[[j]][i:(i+nnL*nnT-1)]<-plogis(lmu.samp[[j]][i:(i+nnL*nnT-1)]+lsamp2[[j]])
    }
  }
 
  rm(mat.samplesaux, mat.samples1, lsamp2)
 
  # Return bernoulli response
  return(lapply(c(1:nSim),FUN=function(i) rbinom(n=nnL*nnT*nnP,size=1,prob=lmu.samp[[i]])))
}


I'm not providing files for reproducing the code as they're very heavy, and my doubt is more on the conceptual side.

Thanks in advance. Best regards,

Julia

[Helpdesk (Haavard Rue)]

the easy way out, is to add the predictions you need into the original data-
frame with NA response, and then the predictions will be part of the original
fit. would that be possible?

On Fri, 2025-08-01 at 03:18 -0700, 'Julia Quero Pérez' via R-inla discussion

> --
> 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, visit
> https://groups.google.com/d/msgid/r-inla-discussion-group/d90acd77-0f40-40a1-8942-2b0482a49aaan%40googlegroups.com
> .

--
Håvard Rue
he...@r-inla.org

[Julia Quero Pérez]

Hello,

The thing is that I don't only want one prediction, I want to do simulations, so I'll be "predicting"/needing the predictors for those points like 500 times. 
That's why my code is that way. I think what you're suggesting is something like adding a prediction stack on top of the estimation stack, right?

> https://urldefense.com/v3/__https://groups.google.com/d/msgid/r-inla-discussion-group/d90acd77-0f40-40a1-8942-2b0482a49aaan*40googlegroups.com__;JQ!!D9dNQwwGXtA!Ux2oyrrJwxpPYXR76x04nFmEm_8X4L1Es9XFRdzGPcEf5LPyWpYi8tYeJYp9Zf9973mlJmtZOwLXmH5e$