Dependence of Bayes factor on prior for nuisance parameters

lm_ranef_H0.stan

data {
  int<lower=0> N;
  vector[N] y;
  real<lower=0> intercept_prior_sd;
  real<lower=0> sigma_prior;
}

parameters {
  real intercept;
  real<lower=0> sigma_raw;
}

transformed parameters {
  real sigma;
  if(sigma_prior == 0) {
    sigma = sqrt(sigma_raw);
  } else {
    sigma = sigma_raw;
  }
}

model {
  if(intercept_prior_sd != 0) {
    target +=  normal_lpdf(intercept | 0, intercept_prior_sd);
  }
  if(sigma_prior != 0) {
    target += normal_lpdf(sigma | 0, sigma_prior) -normal_lccdf(0 | 0, 1); //norm. constant for the sigma prior;;
  } else {
    target += -1 * log(sigma_raw);
  }

  target += normal_lpdf(y | intercept, sigma);
}

lm_ranef_H1.stan

data {
  int<lower=0> N;
  vector[N] y;
  int<lower=1> N_groups;
  array[N] int<lower=1,upper=N_groups> groups;
  real<lower=0> intercept_prior_sd;
  real<lower=0> sigma_prior;
}

parameters {
  real intercept;
  real<lower=0> sigma_raw;
  vector[N_groups] group_z;
  real<lower=0> group_sd;
}

transformed parameters {
  vector[N_groups] group_r = group_z * group_sd;
  real sigma;
  if(sigma_prior == 0) {
    sigma = sqrt(sigma_raw);
  } else {
    sigma = sigma_raw;
  }
}

model {
  if(intercept_prior_sd != 0) {
    target +=  normal_lpdf(intercept | 0, intercept_prior_sd);
  }
  if(sigma_prior != 0) {
    target += normal_lpdf(sigma | 0, sigma_prior) -normal_lccdf(0 | 0, 1); //norm. constant for the sigma prior;;
  } else {
    target += -1 * log(sigma_raw);
  }
  target += std_normal_lpdf(group_z);
  target += normal_lpdf(group_sd | 0, 1)
            -normal_lccdf(0 | 0, 1); //norm. constant for the sigma prior;
  target += normal_lpdf(y | intercept + group_r[groups], sigma);
}

nuisance_parameters_BF.R

library(bridgesampling)
library(rstan)

# H0: intercept only
m_H0 <- stan_model("stan/lm_ranef_H0.stan")
# H1: intercept + single random intercept
m_H1 <- stan_model("stan/lm_ranef_H1.stan")

# The seed gives some variability, and not all seeds exhibit the problem
# but the problem is quite reliable
seed <- 14611187
set.seed(seed)

# Simulate dataset assuming random intercept present
N <- 40
N_groups <- 20

groups <- rep(1:N_groups, length.out = N)

intercept <- 1
group_sd <- 0.1
group_r <- rnorm(N_groups, sd = group_sd)

# Reducing sigma seems to result in stronger changes in BF, but not heavily tested
sigma <- 0.1

mu <- intercept + group_r[groups]
y <- rnorm(N, mu, sigma)

# Fit with default improper prior on intercept and residual sd
stan_data <- list(N = N,
                  N_groups = N_groups,
                  groups = groups,
                  y = y,
                  # Setting zeroes here make the prior improper - flat for intercept and 1/sigma^2 for residual sigma
                  intercept_prior_sd = 0,
                  sigma_prior = 0)

fit_H0 <- sampling(m_H0, data = stan_data,
                   iter = 15500, warmup = 500, chains = 4, seed = 1, refresh = 0)

fit_H1 <- sampling(m_H1, data = stan_data,
                   iter = 15500, warmup = 500, chains = 4, seed = 1, refresh = 0, control = list(adapt_delta = 0.95))

bridge_H0 <- bridge_sampler(fit_H0)
bridge_H1 <- bridge_sampler(fit_H1)

# Fit with weakly informative N(0, 1) prior on both intercept and residual sd
stan_data_prior <- stan_data
stan_data_prior$intercept_prior_sd <- 1
stan_data_prior$sigma_prior <- 1
fit_H0_prior <- sampling(m_H0, data = stan_data_prior,
                         iter = 15500, warmup = 500, chains = 4, seed = 1, refresh = 0)

fit_H1_prior <- sampling(m_H1, data = stan_data_prior,
                         iter = 15500, warmup = 500, chains = 4, seed = 1, refresh = 0, control = list(adapt_delta = 0.95))

bridge_H0_prior <- bridge_sampler(fit_H0_prior)
bridge_H1_prior <- bridge_sampler(fit_H1_prior)

# Check no big errors
error_measures(bridge_H0)$percentage
error_measures(bridge_H1)$percentage


error_measures(bridge_H0_prior)$percentage
error_measures(bridge_H1_prior)$percentage


# Compare Bayes factors under different priors for nuisance parameters
bf(bridge_H1, bridge_H0) # 4.23839
bf(bridge_H1_prior, bridge_H0_prior) # 2.83979