# Use mean squared error (MSE) for comparing model fits of Bayesian models

I want to use mean squared error (MSE) to assess/copmare the model fit of the Bayesian models. The formula for MSE is

$$MSE=\frac{1}{n}\sum^n_i{(y_i-\hat{y}_i)^2}$$

I'm not sure how MSE is used for Bayesian models, although I see researchers use MSE to compare Bayesian models in their papers. Suppose I have the following linear model

$$y_i=\beta*X_i+\epsilon_i, \epsilon \sim N(0,\sigma^2)$$

I use the following Bayesian code to estimate the model,

###simulate dataset
set.seed(66)
#true value
beta=c(2,3,4) #true beta
sigmasq=c(4) #true sigmasq

n=200
nvar=length(beta)
X=cbind(rep(1,n),runif(n),runif(n))
y=X%*%beta+rnorm(n,sd=sqrt(sigmasq))

##Bayesian conjugate linear model-------
#priors
betabar=rep(0,nvar)
A=diag(nvar)

nu=max(4,0.01*n)
ssq=1

####Mcmc set up
R=20000 #number of iterations
Data=list(y=y,X=X)
Prior=list(betabar=betabar,A=A,nu=nu,ssq=ssq)
Mcmc=list(R=R)

conjugate_linear<-function(Data,Prior,Mcmc){
y=Data$$y x=Data$$X
beta0=Prior$$betabar sigma0=Prior$$A
s0=Prior$$ssq v0=Prior$$nu

nvar=length(beta0)

iter=Mcmc$R result.beta=matrix(0,iter,ncol=nvar) result.sigma2=rep(0,iter) for(i in 1:iter){ RA=chol(sigma0) W=rbind(x,RA) z=c(y,as.vector(RA%*%beta0)) IR=backsolve(chol(crossprod(W)),diag(nvar)) beta_title=crossprod(t(IR))%*%crossprod(W,z) res=z-W%*%beta_title s=crossprod(res) ##draw sigma2 sigma2=(s+s0*v0)/rchisq(1,df=n+v0) ##draw beta|sigma2 beta=beta_title+as.vector(sqrt(sigma2))*IR%*%rnorm(nvar) result.beta[i,]<-beta result.sigma2[i]<-sigma2 } return(list(betadraw=result.beta,sigmasqdraw=result.sigma2)) } estimates<-conjugate_linear(Data,Prior,Mcmc) #estimate the model #estimation results ##results for beta summary(estimates$$betadraw) beta_mean=apply(estimates$$betadraw,2,mean) ##results for sigmasq summary(estimates$$sigmasqdraw) sigma_sq_mean=mean(estimates$$sigmasqdraw)  Here, estimates$betadraw is the MCMC draws for $$\beta$$ and beta_mean ($$\bar{\beta}$$) is the mean of $$\beta$$ draws. Similarly, estimates$sigmasqdraw is the MCMC draws for $$\sigma^2$$ and sigma_sq_mean($$\bar{\sigma}^2$$) is the mean of $$\sigma^2$$ draws. My core question is what quantities should I use to generate the predicted y (i.e., $$\hat{y_i}$$)? I have the following sub-questions and options. 1. Should the error term be ignored as in the frequentist way when generate predicted y? That is, $$\hat{y}_i=\bar{\beta}*X_i$$ #ignore error y_pred1=X%*%beta_mean MSE1=mean((y-y_pred1)^2) print(MSE1) # [1] 4.048716 #simulate error y_pred2=X%*%beta_mean+rnorm(n,sd=sqrt(sigma_sq_mean)) MSE2=mean((y-y_pred2)^2) print(MSE2) # [1] 8.853691  1. Unlike the above options which use the mean of MCMC draws, is it better to use individual MCMC draws and generate simulated datasets, then generate a vector of MSEs? #simualte multiple datasets using MCMC draws ##ignore errors MSE3s=sapply(1:R,function(i){ beta_i=estimates$betadraw[i,]
y_pred3=X%*%beta_i
MSE3=mean((y-y_pred3)^2)
MSE3
})
MSE3=mean(MSE3s)
# [1] 4.108187

##simulate errors
MSE4s=sapply(1:R,function(i){
beta_i=estimates$$betadraw[i,] sigmasq_i=estimates$$sigmasqdraw[i]

y_pred4=X%*%beta_i+rnorm(n,sd=sqrt(sigmasq_i))
MSE4=mean((y-y_pred4)^2)
MSE4
})
MSE4=mean(MSE4s)
print(MSE4)
# [1] 8.267055


1. MSE3 and MSE4 use the mean of MSEs from each simulated dataset as the final MSE measure. Alternatively, I can supply the mean of y across simulated datasets into the MSE formula. Is it better?
y_pred5s=sapply(1:R,function(i){
y_pred5=X%*%beta_i
y_pred5
})
y_pred5<-apply(y_pred5s,1,mean)
MSE5=mean((y-y_pred5)^2)
print(MSE5)
# [1] 4.048716

y_pred6s=sapply(1:R,function(i){
beta_i=estimates$$betadraw[i,] sigmasq_i=estimates$$sigmasqdraw[i]

y_pred6=X%*%beta_i+rnorm(n,sd=sqrt(sigmasq_i))
y_pred6
})
y_pred6<-apply(y_pred6s,1,mean)
MSE6=mean((y-y_pred6)^2)
print(MSE6)
# [1] 4.052652


Which MSE above is the most appropriate one (most widely used in practice)? Or they are all acceptable, as long as different Bayesian models use the same measure?