Can't predict the Linear predictor by using INLA

[Jin Hyung Lee]

In the spatial statistics setting for 100 locations, I aim to compare predictions of the linear predictor XB+W obtained from MCMC using Metropolis-Hastings and INLA. Here, X represents the covariates, B is the vector of regression coefficients (beta), and W denotes the spatial random effects. I have successfully implemented XB+W, which serves as the linear predictor for both count and binary data under this Bayesian Hierarchical Model.

However, when I applied this code to a Bayesian Hierarchical model for Gaussian data without W in the linear predictor, it caused issues. I also used summary.linear.predictor, but it did not perform as well.

1. How can I predict XB using INLA in my case?
2. My RMSPE from INLA is significantly worse than from MCMC. How can I obtain a more accurate RMSPE?

Thank you very much for your time and consideration. 

# INLA_Gaussian
getwd()
rm(list=ls())
load("G_phi1_Gen500_1.RData")
load("G_phi1_Gen500_1_MCMC.RData")
################################################################################
## (1) Creating Mesh using Coordinates
mesh <- inla.mesh.2d(X_trainingdata, max.edge=c(0.1,0.2))
plot(mesh)
dim(mesh$loc)
points(x=mesh$loc[,1], y=mesh$loc[,2],col="blue",pch=16,cex=1)
points(x=X_trainingdata[,1], y=X_trainingdata[,2],col="red",pch=16,cex=1)
################################################################################
## (2) SPDE object
spde = inla.spde2.matern(mesh=mesh, alpha=2)
################################################################################
## (3) the projector matrix (sparse matrix which is eta_star)
AMat <- inla.spde.make.A(mesh, loc=X_trainingdata)  
dim(AMat)
AMatCV <- inla.spde.make.A(mesh, loc=X_testdata)  
dim(AMatCV)
################################################################################
## Before (4) We might need to make an index function but not here.
## (4) Stack all the Relevant pieces together
## (4)-1 For the estimation
stk.e <- inla.stack(tag='est', ## tag
                    data=list(y=Z_trainingdata), ## response
                    A=list(AMat, 1), ## two projector matrix
                    effects=list(## two elements:
                      s=1:spde$n.spde, ## RF index
                      data.frame( x1=X_trainingdata[,1],x2=X_trainingdata[,2])))

## (4)-2 Prediction fitting
stk.pred <- inla.stack(tag='pred', A=list(AMatCV, 1),
                       data=list(y=NA), ## response as NA
                       effects=list(s=1:spde$n.spde,
                                    data.frame(x1=X_testdata[,1],x2=X_testdata[,2])))

################################################################################
#(5) Finally Join all the stack objects p.222
stk.full <- inla.stack(stk.e, stk.pred)

################################################################################
#(6) We need to define the formula
formula <- y ~ -1  + x1 + x2 +## fixed part
  f(s, model=spde) ## RF term

################################################################################
#(7) call INLA function
link = c(rep(NA, dim(X_trainingdata)[1]),rep(1, dim(X_testdata)[1]))

print("Begin INLA")
res<-inla(formula,
          data=inla.stack.data(stk.full, spde=spde),
          family="gaussian",  #binomial, poisson
          control.predictor=list(A=inla.stack.A(stk.full), compute=TRUE, link=link),
          control.compute = list(return.marginals.predictor=TRUE),
          keep=FALSE, verbose=TRUE)
summary(res) # Summary
print("End INLA")
print(res$cpu.used)
length(res$marginals.linear.predictor)

################################################################################
#(8) call INLA function
field_est_mean =
  res$summary.linear.predictor[inla.stack.index(stk.full,"est")$data, "mean"]
field_est_sd =
  res$summary.linear.predictor[inla.stack.index(stk.full,"est")$data, "sd"]
# Cross Validation
field_pred_mean =
  res$summary.linear.predictor[inla.stack.index(stk.full,"pred")$data, "mean"]
field_pred_sd =
  res$summary.linear.predictor[inla.stack.index(stk.full,"pred")$data, "sd"]
obs_RMSPE_INLA<-sqrt(mean((Z_testdata-field_pred_mean)^2)) #CVRMSPE to be used as gold standard
comptTime_INLA<-res$cpu.used # Computation Time

################################################################################
# (9) Compare with MCMC and INLA by plot
#MCMC
plot(density(Z_prediction_MCMC[20,]), main = "Density for Linear Predictor Prediction of test_Eta_20", col="red", xlim=c(-3,3))
# INLA
xSeq<-seq(-3,3,length.out=1000)
ySeq<-dnorm(x = xSeq , mean = field_pred_mean[20] , sd = field_pred_sd[20])
lines(x=xSeq , y = ySeq , typ="l" , col="purple", lty=1)
legend("topright", lty=c(1,1,1), legend = c("MCMC", "INLA"),
       col=c("red", "purple"))
################################################################################
# (10) Compare RMSPE
obs_RMSPE_INLA #RMSPE of INLA
Z_RMSPE_MCMC #RMSPE of MCMC