# Am I wrong in assuming that the bayesian estimate has to converge to the true mean?

I wrote a small simulation in R wherein I simulate two slightly different groups that are normally distributed. I then simulated a series of ten identical experiments where 25 new participants are sampled for each group for each experiment.

I start with a simple bayesian model with a noniformative prior for the effect of group on Y. Then:

• I take the posterior mean and sd of the model
• construct a new prior from the posterior mean and sd
• feed this prior into the model for the next experiment

and so on and so on. basically, the posterior of experiment t-1 becomes the prior for experiment t.

My intuition would suggest that over time, the estimate should slowly converge towards the true mean of the population. However, this is not what happens. Am I missing something?

#########################################################################################
#########################################################################################

if (!require(tidyverse)) install.packages('tidyverse')
#if (!require(cowplot)) install.packages('cowplot')
if (!require(brms)) install.packages('brms')
library(tidyverse)
#library(cowplot)
library(brms)

#########################################################################################
# Simulating Data
#########################################################################################

# difference between means of two groups
difference <- 0.2

# variation within group
variation <- 1

# number of experiments
experiments <- 10

# samples per experiment
n <- 25

#------------------------------------------------------------------------------------
#simulation of new data for the two groups
#------------------------------------------------------------------------------------

mysample <- c(rnorm(n,0,variation),rnorm(n,difference,variation))
df <- data.frame(mysample,group=rep(c(1,0),each=n))

#------------------------------------------------------------------------------------
#variables reserved for the for loop
#------------------------------------------------------------------------------------

mofs <- rep(NA,experiments)
mypmean <- rep(NA,experiments)
mypsd <- rep(NA,experiments)

# dimensions of the plot
minva <- difference-(variation*3)
maxva <- difference+(variation*3)

#########################################################################################
# Kickstart the Model
#########################################################################################

# We first need to kickstart the model with some rather broad prior:
bmodel1 <- brm(formula=mysample~group,
prior=c(
set_prior("normal(0,50)",class="b",coef="group") #set an uninformative prior
),
data = df,
warmup = 200, iter = 1000, chains = 4,
cores = getOption("mc.cores", 3L))

#########################################################################################
# Extract the posterior summary
#########################################################################################

# extract posterior
param <- posterior_summary(bmodel1)

#create a new Prior from the posterior
newprior <- paste0("normal(",param[[2]],",",param[[3]],")")

#########################################################################################
# Run the specified number of experiments
#########################################################################################

for (i in c(1:experiments)) {

# two groups with mean difference as specified and a variation as specified. combined as a data frame
mysample <- c(rnorm(n,0,variation),rnorm(n,difference,variation))
df <- data.frame(mysample,group=rep(c(1,0),each=n))

# the mean difference between the groups of the sample "mofs" is saved just to see how it varied in each experiment
mofs[i] <- mean(df$mysample[which(df$group==1)]-
df$mysample[which(df$group==0)])

#------------------------------------------------------------------------------------
#The kickstarted model us updated with new data AND a new prior (the posterior of the experiment before)
#------------------------------------------------------------------------------------
tmp <- update(bmodel1,prior=
set_prior(newprior),newdata = df)

# save the estimated difference between both groups  and the sd
mypmean[i] <- param[[2]]
mypsd[i] <- param[[3]]

#------------------------------------------------------------------------------------
# create a new prior from the posterior
#------------------------------------------------------------------------------------
param <- posterior_summary(tmp)
newprior <- paste0("normal(",param[[2]],",",param[[3]],")")
}
#####################################################################################
# Plot the results
####################################################################################

df <- data.frame(estimatedDif=mypmean,
index=c(1:experiments),
sd=mypsd,
sampledDif=mofs)

x11()
ggplot(df,aes(index,estimatedDif))+
geom_point(aes(index,sampledDif),size=3,alpha=0.3,color="lightblue")+
geom_point(alpha=1)+
ylim(c(minva,maxva))+
geom_hline(aes(yintercept=difference*-1),color="green")+
ggtitle(paste0("Sampled mean (blue) and Estimated (black) for ",experiments," Experiments and n=",n))

• Could it be that you are only updating your prior based on a model that incorporates new data from only n = 25, and that you just haven't updated your prior enough? There must be an analytical solution here to compare to since it's only two normals (really one normal distribution, that of the differences) you're estimating. – Brash Equilibrium Jun 13 '18 at 16:31
• The updated prior is only based on n=25, that was kind of what I wanted to see: will a (bunch of) scientist(s) using bayes get a proper idea of the true difference if they have noisy data, small effects, and small n by just running an identical experiment over and over again that modifies their prior believe. – Braino Jun 13 '18 at 16:37
• Try the experiments with (a) a larger sample size per iteration, (b) more iterations, and see what happens. – Brash Equilibrium Jun 13 '18 at 16:46
• i just did. I ran 300 experiments with n=25. it did not converge. When I run 300 experiments with n=300 even then. It does not seem to converge – Braino Jun 13 '18 at 16:59

Okay. So there where three things that I did wrong:

1. I did not include a prior for the intercept and sigma
2. I added the wrong parameter (sigma) to the prior for the group instead if its SE
3. The "update" function needs to recompile the model to incoorporate the new priors, otherwise it will just proceed with the old uniformative prior.

Of course, letting brms recompile the model for each subsequent update takes a lot more time so that 300 simulations are not feasible for my very, very old laptop. So I just proceeded with ten experiments, an effect of 0.2 n=25 and a variation of 1:

The changes in the code are

    tmp <- update(bmodel1,prior=c(
set_prior(newpriorInt,class="Intercept"),
set_prior(newpriorG,class="b",coef="group"),
set_prior(newpriorSig,class="sigma")
),newdata = df)

#------------------------------------------------------------------------------------
# create a new prior from the posterior
#------------------------------------------------------------------------------------
param <- posterior_summary(tmp)

newpriorInt <- paste0("normal(",param[[1]],",",param[[1,2]],")")
newpriorG <- paste0("normal(",param[[2,1]],",",param[[2,2]],")")
newpriorSig <- paste0("normal(",param[[3,1]],",",param[[3,2]],")")

• Andrew Gelman has famously said and I paraphrase that when the model output makes zero sense it is probably a coding error. Then again he also says he’s routinely amazed by how robust some of his analysis results are to some pretty egregious coding errors. I’m sorry I didn’t look at your code closely enough to catch this, but I’m glad you did. – Brash Equilibrium Jun 21 '18 at 5:00