Allowing data to remain unclassified in Bayesian model selection

Question

This is a follow up to a previous question. In my analysis, I am using a Bayesian model selection (using MCMC) to determine which model best fits several data sets. However, because I want to allow my model to also leave some datasets unclassified, I also included a complex model that is more general than either of the two models (see here for more details on this approach).

Unfortunately, for some model parameters, the unclassified model seems to fit much better than both of the actual models. I can replicate this effect in a simulation:

library(tidyverse)
library(ggplot2)

N <- 10

reps <- 1000

P1 <- function(k1,k2,N) {
  beta(N-k2+k1+1,N-k1+k2+1)
}

P2 <- function(k1,k2,N) {
  beta(k2+k1+1,2*N-k1-k2+1)
}

P3 <- function(k1,k2,N) {
  beta(k1+1,N-k1+1)*beta(k2+1,N-k2+1)
}


v <- seq(0.01,0.99,by=0.01)

df <- expand.grid(v = v, i = seq_len(reps), mdl=c(1,2))

df$K1 <- map_dbl(df$v,function(v) rbinom(1,N,v))
df$K2 <- map2_dbl(df$v,df$mdl, function(v,mdl) ifelse(mdl == 1, rbinom(1,N,1-v), rbinom(1,N,v)) )


df$P1 <- P1(df$K1,df$K2,N)
df$P2 <- P2(df$K1,df$K2,N)
df$P3 <- P3(df$K1,df$K2,N)

# roulette wheel selection based on probabilities
select.rw <- function(p) {
  iv <- p/sum(p)
  pv <- runif(1,0, 1 )
  for(i in seq_along(iv)) {
    pv <- pv - iv[i]
    if(pv <= 0) {
      return( i )
    }
  }
}

# We do a roulette wheel selection to simulate model selection using MCMC
# Also randomly selecting removes artifacts that appear when P1==P2
df$ret <- df[,c("P1","P2","P3")] %>% array_branch(margin = 1) %>% map_dbl(select.rw)

df$ACC <- (df$mdl == df$ret) * 1
df %>% group_by(mdl,v) %>% summarize(M=mean(ACC)) %>% ungroup() %>% ggplot(aes(x=v,y=M,color=factor(mdl))) + geom_line()

df %>% group_by(mdl,v) %>% summarize(M1 = mean(ret==1),M2 = mean(ret==2),M3 = mean(ret==3)) %>%
  ggplot(aes(x=v,linetype=factor(mdl))) + 
  geom_line(aes(y=M1),color="green") + geom_line(aes(y=M2),color="blue") + geom_line(aes(y=M3),color="red")

If I look at the rate of misclassification in relation to the parameter $v$, then I find that models are more likely to become misclassified if $v$ is close to $0.5$:

In this simulation, this is no problem, because in the end I still get a $50/50$ classification rate for both actual models after classification, which fits the initial distribution.

Unfortunately, during analysis, I found that the parameters are not distributed equally for both classes. I can simulate this by choosing the parameter from a different distribution for the two models:

library(tidyverse)
library(ggplot2)

N <- 10

reps <- 1000

P1 <- function(k1,k2,N) {
  beta(N-k2+k1+1,N-k1+k2+1)
}

P2 <- function(k1,k2,N) {
  beta(k2+k1+1,2*N-k1-k2+1)
}

P3 <- function(k1,k2,N) {
  beta(k1+1,N-k1+1)*beta(k2+1,N-k2+1)
}


v <- seq(0.01,0.99,by=0.01)

df <- expand.grid(i = seq_len(reps*100), mdl=c(1,2))

# Unequal probabilities in either group
df$v <- ifelse(df$mdl == 1, rbeta(nrow(df),2,2), rbeta(nrow(df),0.5,0.5))

df$K1 <- map_dbl(df$v,function(v) rbinom(1,N,v))
df$K2 <- map2_dbl(df$v,df$mdl, function(v,mdl) ifelse(mdl == 1, rbinom(1,N,1-v), rbinom(1,N,v)) )


df$P1 <- P1(df$K1,df$K2,N)
df$P2 <- P2(df$K1,df$K2,N)
df$P3 <- P3(df$K1,df$K2,N)

# roulette wheel selection based on probabilities
select.rw <- function(p) {
  iv <- p/sum(p)
  pv <- runif(1,0, 1 )
  for(i in seq_along(iv)) {
    pv <- pv - iv[i]
    if(pv <= 0) {
      return( i )
    }
  }
}

# We do a roulette wheel selection to simulate model selection using MCMC
# Also randomly selecting removes artifacts that appear when P1==P2
df$ret <- df[,c("P1","P2","P3")] %>% array_branch(margin = 1) %>% map_dbl(select.rw)

(df$ret %>% factor() %>% summary())/nrow(df)

However, in this case I find that the final rate differs from the initial 50/50 rate for the two models because model 1 has a much higher chance of remaining unclassified:

> (df$ret %>% factor() %>% summary())/nrow(df)
      1       2       3 
0.29076 0.41503 0.29421 

So my initial attempt of allowing the MCMC to keep some datasets as unclassified creates artifacts in the data.

Is there a way to either remove these artifacts (see also my previous question on this) or to handle possible unclassified datasets differently that does not create these artifacts?