The problem is from chapter 10 question 7 of “Regression and Other Stories”, and the general idea is to use data which includes a continuous parameter “beauty” and then several binary predictors. I have solved the problem using stan_glm in R, so I know that out of the 7 predictors only three have sizable effects, beauty (continuous), gender(1,0), and nonnative English speaker (1,0). It is easy to recover these results using julia’s `GLM.jl package, and I show this below:
df = DataFrame(CSV.File("beauty.csv")) sanity_check = lm(@formula(eval~beauty+female+age+ minority+nonenglish+lower+course_id),df)
These are the results I would expect. However, when I create a Turing model for Bayesian linear regression I get that only beauty is significant.
My code for the Turing model is shown below. I use a small helper function (also shown) to standardize the predictors, HOWEVER, in this version of my turing model I only use the helper function calculate the standard deviations and means in one go. I leave the data un-transformed but center the priors because I am trying to emulate what stan_glm does here: Prior Distributions for rstanarm Models • rstanarm
function standardize(X) μ_X = mean(X) σ_X = std(X) X_std = (X .- μ_X)/σ_X return X_std, μ_X, σ_X end @model function beauty_regression(y,x1,x2,x3,x4,x5,x6,x7) x1_std,μ_x1,σ_x1 = standardize(x1) x2_std,μ_x2,σ_x2 = standardize(x2) x3_std,μ_x3,σ_x3 = standardize(x3) x4_std,μ_x4,σ_x4 = standardize(x4) x5_std,μ_x5,σ_x5 = standardize(x5) x6_std,μ_x6,σ_x6 = standardize(x6) x7_std,μ_x7,σ_x7 = standardize(x7) σ_y = std(y) ϵ ~ Exponential(1/σ_y) my = 0 β0 ~ Normal(my,2.5*σ_y) β1 ~ Normal(0,2.5*σ_y/σ_x1) β2 ~ Normal(0,2.5*σ_y/σ_x2) β3 ~ Normal(0,2.5*σ_y/σ_x3) β4 ~ Normal(0,2.5*σ_y/σ_x4) β5 ~ Normal(0,2.5*σ_y/σ_x5) β6 ~ Normal(0,2.5*σ_y/σ_x6) β7 ~ Normal(0,2.5*σ_y/σ_x7) μ = β0 .+ β1*(x1 .- μ_x1) + β2*(x2 .- μ_x2) + β3*(x3 .- μ_x3) + β4*(x4 .- μ_x4) + β5*(x5 .- μ_x5) + β6*(x6 .- μ_x6) + β7*(x7 .- μ_x7) y ~ MvNormal(μ,ϵ)
I run my simulation with the following code
model = beauty_regression(df[:,1],df[:,2],df[:,3],df[:,4],df[:,5],df[:,6],df[:,7],df[:,8]) chain = sample(model,NUTS(0.65),3_000);
And these are the results:
Same order but different names (sorry) We see that beauty (b1) is the only significant one…My guess is that I am doing something wrong specifying priors because I have a combination of continuous and binary, however, the stan_glm methodologies makes no such mention of doing anything different. Why am I getting pretty different results with my model? I appreciate any insights
EDIT1: Accidentally submitted the version where I DID standardize the data, should be the non-standardized version now
EDIT2: Fixed some embarrassing simple mistakes, but problem still persists. The code posted above has been updated to remove said mistakes.