# Simple linear regression using Approximate Bayesian Computation (ABC)

To understand how ABC works, I like to use it in a simple linear regression model. I am using EasyABC package in R. My problem is that the results I get are not good. In particular, the posterior distributions are extremely wide. I think this is partly because I am using prior distributions with a very large variance. In fact, if I use different priors (e.g., Unif(0,20) --- something much tighter), results are better. But I am using those priors to mean something commonly used as non-informative priors in Bayesian method.

My question is that is there anything wrong in the way ABC is used in this example (the code is given below)? Why is the posterior so bad?

library(EasyABC)

## data generation
n <- 100
intercept <- 10
slope     <- 1
s         <- 2.5
x <- runif(n,10,20)
y <- rnorm(n,intercept+slope*x,s)

model <- lm(y~x)
intercept.data <- coef(model)[1]
slope.data <- coef(model)[2]
sd.data <- summary(model)$sigma sum_stat_obs <- c(intercept.data,slope.data,sd.data) sim.abc <- function(par){ samples <- rnorm(n, mean =par[1]+par[2]*x, sd = par[3]) model0 <- lm(samples~x) intercept.sim <- coef(model0)[1] slope.sim <- coef(model0)[2] sd.sim <- summary(model0)$sigma
return(c(intercept.sim,slope.sim,sd.sim))
}

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
prior=list(c("normal",0,1000),c("normal",0,1000),c("exponential",1)),
summary_stat_target=sum_stat_obs, n_rec = 10000)

par(mfrow=c(2,2))
plot(model.abc$param[5000:10000,1],type="l") hist(model.abc$param[5000:10000,1], main = "Posterior for intercept")
plot(model.abc$param[5000:10000,2],type="l") hist(model.abc$param[5000:10000,2], main = "Posterior for slope")


The main problem is the tolerance together with the terrible choice of priors (recall that ABC does not work as other MCMC methods since it is based on simulations from the priors). You need to reduce the tolerance_quantile in order to reduce the distance between simulated and observed summary statistics.

See for instance (better tolerance):

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
prior=list(c("normal",0,100),c("normal",0,100),c("exponential",1)),
summary_stat_target=sum_stat_obs, n_rec = 10000, tolerance_quantile=0.001)

apply(model.abc$param[5000:10000,],2,median) apply(model.abc$param[5000:10000,],2,mean)


and (better tolerance and priors):

model.abc <-ABC_mcmc(method="Marjoram", model=sim.abc,
prior=list(c("normal",10,5),c("normal",1,3),c("exponential",2)),
summary_stat_target=sum_stat_obs, n_rec = 10000, tolerance_quantile=0.0005)

apply(model.abc$param[5000:10000,],2,median) apply(model.abc$param[5000:10000,],2,mean)

par(mfrow=c(3,2))
plot(model.abc$param[5000:10000,1],type="l") hist(model.abc$param[5000:10000,1], main = "Posterior for intercept")
plot(model.abc$param[5000:10000,2],type="l") hist(model.abc$param[5000:10000,2], main = "Posterior for slope")
plot(model.abc$param[5000:10000,3],type="l") hist(model.abc$param[5000:10000,3], main = "Posterior for sigma")