# How to specify the unit information prior

Continuing on from this question and this question re BIC and its approximation to the Bayes factor with a unit information prior (Kass & Wasserman, 1995), I'm trying to quantify this relationship as a stepping stone into Bayesian stats. So far, my calculation of the BIC approximation of the Bayes factor (based upon my impression of Wagenmakers 2007) is linearly related to my Bayes factor that is calculated from my interpretation of the unit information prior using the INLA package in R. Good start! However, my BIC Bayes factor is ~ 3 times smaller than the Bayes factor calculated with INLA and I'm not sure why.

The prior I've used in the "inla" function is N(0, 1/(variance * n)) and this seems to me the likely place where I'm out. I'm not sure how I got the multiply by n in the formula, but it appears to work... roughly. Kass and Wasserman have N(0, variance / n) which when converted to precision would be N(0, n / variance), but this gives me a less good relationship.

Help based on other Bayesian packages is also welcome.

EDIT

*Deleted code, see below answer instead*


EDIT

So I'm pretty sure I've figured out the one sample case. I would still appreciate help for the two sample case and the regression case (which I'll start working on now).

I think I have figured it out for the one sample case after updating my knowledge with Doing Bayesian Data Analysis and the rjags code contained within. The model specification I've used below gives me a slope of ~ 1 and intercept of ~ 0 when the BIC Bayes factor approximation is regressed against the Bayes factor from the rjags model.

model {
# Likelihood:
for( i in 1 : N ) {
y[i] ~ dnorm( mu , tau ) # tau is precision, not SD
}
# Prior:
tau ~ dnorm( yprec , (yprec^2)/2 )
mu ~ dnorm( amu , atau )
#Hyperprior
atau <- tauModel[ modelIndex ]
amu  <- muModel[ modelIndex ]
tauModel[1] <- 10000
tauModel[2] <- yprec
muModel[1] <- 0
muModel[2] <- ymean

#Hyperhyperprior
modelIndex ~ dcat( modelProb[] )
modelProb[1] <- 0.5
modelProb[2] <- 0.5
}


And the input data:

y = rnorm( n=200 , mean=1 , sd=1 )
yprec = 1/var(y)
ymean = mean(y)


And for linear regression:

model {
for( i in 1 : Ndata ) {
y[i] ~ dnorm( mu[i] , tau )
mu[i] <- beta0 + beta1 * x[i]
}
beta0 ~ dnorm( b0mu , b0tau )
beta1 ~ dnorm( b1mu , b1tau )
tau ~ dgamma( yrprec , (yrprec^2)/2 )

#Hyperprior
b0tau <- b0tauModel[ modelIndex ]
b0mu  <- b0muModel[ modelIndex ]
b1tau <- b1tauModel[ modelIndex ]
b1mu  <- b1muModel[ modelIndex ]

b0muModel[1] <- zymean
b0tauModel[1] <- yrprec
b1muModel[1] <- 0
b1tauModel[1] <- 10000

b0muModel[2] <- xbeta0
b0tauModel[2] <- yrprec
b1muModel[2] <- xbeta1
b1tauModel[2] <- yrprec

#Hyperhyperprior
modelIndex ~ dcat( modelProb[] )
modelProb[1] <- 0.5
modelProb[2] <- 0.5
}


Regression input data:

x <- rnorm(50, sd=30) + 1:50  ;  y <- rnorm(50, sd=30) + 1:50
nSubj <- length(x)
xM = mean( x ) ; xSD = sd( x ); yM = mean( y ) ; ySD = sd( y )
zx = ( x - xM ) / xSD         ; zy = ( y - yM ) / ySD
#
lm1 <- lm(zy ~ zx); slm1 <- summary(lm1)
yrprec <- 1/var(resid(lm1))
zymean <- mean(zy)
xbeta0 <- coef(slm1)[1,1]      ; xbeta1 <- coef(slm1)[2,1]


Framework for these models is from Kruschke

• Could you please explain the code rather then providing code-only answer? – Tim Nov 21 '17 at 9:14