# Posterior Predictive CARBayesST

I'm trying to use the CARBayesST package and I need to do Spatio-temporal predictions. In the vignette of the package on page 27 says " If there had been saying m missing values, then the Y component of the list would have contained m columns, with each one containing posterior predictive samples for one of the missing observations."

But I don't understand well how to obtain the posterior predictive values of Y, let's say I want to predict the value of Y for the next 3 periods for each zone. How should I do it?

This is the reproducible code (found in the vignette):

library(CARBayesST)
#################################################
#### Run the model on simulated data on a lattice
#################################################
#### set up the regular lattice
x.easting <- 1:10
x.northing <- 1:10
Grid <- expand.grid(x.easting, x.northing)
K <- nrow(Grid)
N <- 10
N.all <- N * K
#### set up spatial neighbourhood matrix W
distance <- as.matrix(dist(Grid))
W <-array(0, c(K,K))
W[distance==1] <-1
#### Simulate the elements in the linear predictor and the data
gamma <- rnorm(n=N.all, mean=0, sd=0.001)
x <- rnorm(n=N.all, mean=0, sd=1)
beta <- 0.1
Q.W <- 0.99 * (diag(apply(W, 2, sum)) - W) + 0.01 * diag(rep(1,K))
Q.W.inv <- solve(Q.W)
phi.temp <- mvrnorm(n=1, mu=rep(0,K), Sigma=(0.1 * Q.W.inv))
phi <- phi.temp
for(i in 2:N)
{
phi.temp2 <- mvrnorm(n=1, mu=(0.8 * phi.temp), Sigma=(0.1 * Q.W.inv))
phi.temp <- phi.temp2
phi <- c(phi, phi.temp)
}
LP <- 3 + x * beta + phi
mean <- exp(LP)
Y <- rpois(n=N.all, lambda=mean)
#### Run the model
model <- ST.CARar(formula=Y~x, family="poisson", W=W, burnin=10,
n.sample=50)

• following this post here, I believe you need to set some of the values of the outcome column to NA. Once these are set to NA, the parameter model$samples$y should have estimates in the plausible range of values. You could then carry out a standard 'score' evaluation comparing the predicted Y to the actual Y.
– user286537
May 27, 2020 at 15:24

Following my comment above, here is a reproducible answer:

            library(CARBayesST)
#################################################
#### Run the model on simulated data on a lattice
#################################################
#### set up the regular lattice
x.easting <- 1:10
x.northing <- 1:10
Grid <- expand.grid(x.easting, x.northing)
K <- nrow(Grid)
N <- 10
N.all <- N * K
#### set up spatial neighbourhood matrix W
distance <- as.matrix(dist(Grid))
W <-array(0, c(K,K))
W[distance==1] <-1
#### Simulate the elements in the linear predictor and the data
gamma <- rnorm(n=N.all, mean=0, sd=0.001)
x <- rnorm(n=N.all, mean=0, sd=1)
beta <- 0.1
Q.W <- 0.99 * (diag(apply(W, 2, sum)) - W) + 0.01 * diag(rep(1,K))
Q.W.inv <- solve(Q.W)
phi.temp <- mvrnorm(n=1, mu=rep(0,K), Sigma=(0.1 * Q.W.inv))
phi <- phi.temp
for(i in 2:N)
{
phi.temp2 <- mvrnorm(n=1, mu=(0.8 * phi.temp), Sigma=(0.1 * Q.W.inv))
phi.temp <- phi.temp2
phi <- c(phi, phi.temp)
}
LP <- 3 + x * beta + phi
mean <- exp(LP)
Y <- rpois(n=N.all, lambda=mean)

# Set last 51 observations to NA
Y2 <- Y
Y2[950:1000] <- NA
Y2

#### Run the model
# NOTES: We changed Y to Y2; We increasted n.sample to 5000
model <- ST.CARar(formula=Y2~x, family="poisson", W=W, burnin=10,
n.sample=5000)

# Access the predictions for the NA values
model$$samples$$Y
# Get mean (or median) of the column for the average posterior prediction
colMeans(model$$samples$$Y)


Note how this prints out values for those 51 values set to NA. We can calculate a naive accuracy as:

caret::RMSE(colMeans(model$samples$Y), Y[950:1000])

Note I took the colMeans, although you may want to take the median.

You can also use other accuracy metrics to your liking. Note that increasing the n.sample parameter greatly increased my accuracy (and most examples you will see will have 20,000+).

This is my shot at the problem - I am working with a similar model set. Let me know what you think! Hope it is helpful.