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.