Simple multi-parameter bayesian linear regression returning un-expected results!

Specific Domains Statistics
, ,
6

A few things I see:

  • rstanarm uses mean(y) as the mean of the intercept prior for a normal likelihood with identity link function, but you are using 0.
  • Stan and Distributions use different parameterizations of the exponential distribution. Use Exponential(sy).

In general, the model could be written more succinctly by writing x as Matrix(df[:, 2:8]):

@model function beauty_regression(
    y,
    x,
    μ_x = vec(mean(x; dims=1)),
    σ_x = vec(std(x; dims=1)),
    x_c = x .- μ_x',
    μ_y = mean(y),
    σ_y = std(y),
    βprior = arraydist(Normal.(0, 2.5*σ_y ./ σ_x))
)
    ϵ ~ Exponential(σ_y)
    α ~ Normal(μ_y, 2.5*σ_y)
    β ~ βprior
    μ = muladd(x_c, β, α)
    y .~ Normal.(μ, ϵ)
end

This implementation moves all of the computations that will not change from one run to another out of the body of the model. muladd fuses the matrix-vector multiplication and addition by a scalar. Another advantage to writing it this way is that the same model can be used for different numbers of predictors.

In case you’re not aware, there’s also TuringGLM, which is the rstanarm/brms of Turing.

EDIT: added centering

4 Likes