How to specify the unit information prior

Question

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).

Answer 1

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