1
$\begingroup$

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){
  beta_i=estimates$betadraw[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?

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.