Comparing a Bayesian model with a Classical model for linear regression

Question

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[1] + beta[2]*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[1]   9.937366e-03   0.09806290  0.002191658    0.002238168
    beta[2]   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.

bayesian jags r regression

Sergio · Accepted Answer

Try:
bayes_model="model {

for(i in 1:N){
dist[i] ~ dnorm(mu[i],tau)
mu[i] = beta[1] + beta[2]*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[1]   -16.88098  6.2008860 0.138586750     0.53168398
beta[2]     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.

Comparing a Bayesian model with a Classical model for linear regression

One Answer

Add your own answers!

Ask a Question