# Comparing a Bayesian model with a Classical model for linear regression

I have fitted following Bayesian linear regression model using rjags package in R, with the help of cars data set. I used some weakly informative priors for parameters.

require(rjags)
dim(cars)
N=length(cars$speed) bayes_model="model { for(i in 1:N){ dist[i] ~ dnorm(mu[i],tau) mu[i] = beta + beta*speed[i] } for (l in 1:2) { beta[l] ~dnorm(0, 100) } tau ~ dgamma(.001,.001) sigma_tau = 1/tau }" model2 <- jags.model(textConnection(bayes_model), data = list(dist=cars$$dist,N=N,speed=cars$$speed), n.chains=2) params <- c('beta','sigma_tau') samps.1 <- coda.samples(model2, params, n.iter = 2000) burn.in=1000 summary.model.1=summary(window(samps.1, start = burn.in)) Stat.model.1=as.data.frame(summary.model.1$statistics)


This is how the results looks like .

 > Stat.model.1
Mean           SD     Naive SE Time-series SE
beta   9.937366e-03   0.09806290  0.002191658    0.002238168
beta   1.650041e-01   0.09903592  0.002213404    0.002330977
sigma_tau 2.341437e+03 522.81381343 11.684631408   11.700676273


When I fit a classical linear regression model , following results are obtained .

  summary(lm(dist~speed ,data=cars))

Call:
lm(formula = dist ~ speed, data = cars)

Residuals:
Min      1Q  Median      3Q     Max
-29.069  -9.525  -2.272   9.215  43.201

Coefficients:
Estimate Std. Error t value Pr(>|t|)
(Intercept) -17.5791     6.7584  -2.601   0.0123 *
speed         3.9324     0.4155   9.464 1.49e-12 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 15.38 on 48 degrees of freedom
Multiple R-squared:  0.6511,    Adjusted R-squared:  0.6438
F-statistic: 89.57 on 1 and 48 DF,  p-value: 1.49e-12


It can be observed that the results based on Bayesian and classical method is not similar. What may be the reason for this ? Is there any problem with my prior distributions ?

Are there any diagnostic plots that should look into ? Also please feel free to let me know if I have missed any important steps in this analysis.

I am relatively new to Bayesian and I am working with different sort of examples to learn how to assign correct prior distributions.

• I do not know R but I will interpret what I can! Looks like your coefficients were shrunk. Try a uniform prior, those results should be very similar and if they are then your weakly informative priors are not weakly informative and need more work. Aug 27 '20 at 15:58

Try:

bayes_model="model {

for(i in 1:N){
dist[i] ~ dnorm(mu[i],tau)
mu[i] = beta + beta*speed[i]
}

for (l in 1:2) { beta[l] ~dnorm(0, 0.001) }
tau ~ dgamma(.01,.01)
sigma_tau = 1/tau
}"


You'll get something like:

> Stat.model.1
Mean         SD    Naive SE Time-series SE
beta   -16.88098  6.2008860 0.138586750     0.53168398
beta     3.89631  0.3812761 0.008521332     0.03262728
sigma_tau 244.69856 53.4120568 1.193733180     1.25568735


Why? Because in JAGS the parameters of the normal distribution are mean $$\mu$$ and precision $$\tau$$, and $$\tau=\sigma^{-2}$$, so

beta[l] ~ dnorm(0, 100)
is a very strong, not a weakly informative prior, it forces beta[l] to be as close to zero as possible.

• Thank you for your answer. I understood it now. So after getting these results , how can we so sure that the results are accurate ? I mean how to validate the results from the Bayesian model ? Like do I need to look at some plots ? As I told in the question , I am kind of new to Bayesian analysis. So learning these steps will be a great help for me. Aug 27 '20 at 17:56
• It's a long story :) A clear, easy, and entertaining introduction could be Richard McElreath, Statistical Rethinking. Aug 27 '20 at 18:02