Linear regression with STAN

Make a STAN model. Specify priors.

require(rstan)
data(cars)

modelcode <- "
data {
int<lower=0> N;
vector[N] dist;
vector[N] speed;
}
parameters {
real intercept;
real slope;
real<lower=0> sigma;
}
model {
intercept ~ normal(0,30);
slope ~ normal(1,2);
sigma ~ exponential(0.3);
dist ~ normal(intercept+slope*speed, sigma);
}
"

aModel<-stan_model(model_code = modelcode)

Priors predictive check of the model. Make some predictions based on prior distributions only.

set.seed(2021)
nsim <- 200
simprior <- data.frame(
  intercept = rnorm(nsim,0,30),
  slope = rnorm(nsim,1,2),
  sigma = rexp(nsim,0.3),
  speed = rep(20, nsim)
)
simprior <- within(simprior, 
  dist<-rnorm(nsim, intercept+slope*speed, sigma)
)
head(simprior)

par(mfrow=c(1,3))
with(simprior, {
  hist(intercept, freq=F)
  curve(dnorm(x,0,30), col="gray", add=T)
  hist(slope, freq=F)
  curve(dnorm(x,1,2), col="gray", add=T)
  hist(sigma, freq=F)
  curve(dexp(x,0.3), col="gray", add=T)
}
)

options(mc.cores=4)
fit<-sampling(aModel,data=list(N=dim(cars)[1],speed=cars$speed,dist=cars$dist),iter=1000,chains=4)
print(fit)

Inference for Stan model: 24e22daa4d287a2e1cc029098e963c65.
4 chains, each with iter=1000; warmup=500; thin=1; 
post-warmup draws per chain=500, total post-warmup draws=2000.

             mean se_mean   sd    2.5%     25%     50%     75%   97.5% n_eff Rhat
intercept  -15.20    0.22 6.26  -27.20  -19.33  -15.51  -11.23   -2.58   840    1
slope        3.78    0.01 0.39    3.00    3.54    3.79    4.04    4.49   842    1
sigma       15.04    0.05 1.44   12.56   14.00   14.98   15.96   18.12  1016    1
lp__      -165.00    0.05 1.25 -167.95 -165.50 -164.72 -164.14 -163.66   764    1

Samples were drawn using NUTS(diag_e) at Fri Feb 14 14:19:19 2020.
For each parameter, n_eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor on split chains (at 
convergence, Rhat=1).

Make posterior predictive check by drawing samples of model parameter values from fitted model. The actual steps of the STAN sampler (No U-Turn Sampler, NUTS) are used for sampling.

fit_ss <- extract(fit, pars=c('intercept','slope','sigma'), permuted=F)
isamp <- sample(1:dim(fit_ss)[1],nsim)
simpost <- data.frame(fit_ss[,1,][isamp,])
simpost$speed <- rep(20, nsim)
simpost <- within(simpost, 
  dist <- rnorm(nsim,intercept+slope*speed, sigma)
)
head(simpost)

Histograms of posterior distributions of model parameters versus priors (in gray).

#stan_plot(fit, ci_level=0.8, outer_level=0.95, show_density=T)
#stan_hist(fit, bins=30)
par(mfrow=c(1,3))
with(simpost, {
  hist(intercept, xlim=c(-100,100), freq=F)
  curve(dnorm(x,0,30), col="gray", add=T)
  hist(slope, xlim=c(-4,6), freq=F)
  curve(dnorm(x,1,2), col="gray", add=T)
  hist(sigma, xlim=c(0,20), freq=F)
  curve(dexp(x,0.3), col="gray", add=T)
  }
)

Are the estimated Intercept and Slope correlated?

cor(fit_ss[,1,1], fit_ss[,1,2])

[1] -0.9550463

Plot Intercept vs. Slope. Overlay the 95% ellipse from linear regression.

par(mfrow=c(1,2))
require(ellipse)
with(simpost, plot(x=intercept, y=slope))
lines(ellipse(lm(dist ~ speed, cars), n=50, level=0.95), col="gray")

fullsample<-extract(fit, pars=c('intercept','slope'))
inter_c<-cut(fullsample$intercept,50)
slope_c<-cut(fullsample$slope,50)
z<-table(inter_c,slope_c)
persp(x=z,col="lightgray",xlab="Intercept",ylab="Slope",zlab="Density")