Edit: I have solved this. I was incorrectly simulating values for the dependent variable and not using a good summary statistic. I have changed my algorithm to:
- Simulate y value using simulated B0, B1, B2. Random normal around simulated sigma.
- Calculate summary statistic as sum of squared differences between simulated y and y from observed data at the same x
I am practicing ABC using a regression model. x is estimated as $\beta_0$ + $\beta_1$*predictor + $\beta_2$*predictor$^2$.
Using lm() in R I am getting this:
data.lm <- lm(dependent ~ predictor + I(predictor^2), data = data) Call: lm(formula = dependent ~ predictor + I(predictor^2), data = data) Residuals: Min 1Q Median 3Q Max -2.8411 -0.9694 0.0017 1.0181 3.3900 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) -7.5555495 1.4241091 -5.305 1.69e-05 *** Predictor 1.2716937 0.0757321 16.792 3.99e-15 *** I(Predictor ^2) -0.0145014 0.0008719 -16.633 4.97e-15 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 Residual standard error: 1.663 on 25 degrees of freedom Multiple R-squared: 0.9188, Adjusted R-squared: 0.9123 F-statistic: 141.5 on 2 and 25 DF, p-value: 2.338e-14
I am implementing the algorithm as follows:
- Set priors for $\beta_0$, $\beta_1$, $\beta_2$ each as $\mu$ = 0, $\sigma$ = 10
- Simulate values for each beta using rnorm(mean, sd)
Set a linear model using the simulated y values and x values from the data
n <- length(data_predictor) y <- rnorm(n,b0_r[i] + b1_r[i]*x + b0_r[i]*x^2, sigma_r[i]) sim.lm <- lm(y~data_predictor + I(data_predictor^2))
Calculate the sum of squared distance from the simulated model's error to the error of the model using the data
- If distance[i] < epsilon (currently 5), accept and put each beta into accepted values arrays
The problem with this is I am getting inaccurate simulated models. The predicted y values are either extremely negative or positive, which gives a high error in the simulated linear model. I know that ABC is heavily dependent on the priors, so do my priors for $\beta_0$, $\beta_1$, $\beta_2$ need to be close to the coefficients from the linear model using the data?
Is there a correct way to implement ABC on a simple linear model?