# Translating a TuringGLM model to a Turing model

**URL:** <https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818>\
**Category:** Probabilistic Programming\
**Created:** [December 18, 2022, 7:48pm UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818 "2022-12-18T19:48:51Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![PeetoomHeida](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/peetoomheida/32/23865_2.png) [@PeetoomHeida](https://discourse.julialang.org/u/PeetoomHeida)\
**Post date:** [December 18, 2022, 7:48pm UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818/1 "2022-12-18T19:48:51Z")

</div>

I am trying to implement a hierarchical model, and am running into issues getting the model to return the parameters I used to generate a test data set. When I implement the model using TuringGLM I get the parameters I expect, but I think I must be doing something wrong when I try to set up the model by hand, rather than relying on a model formula. My TuringGLM formula is `response ~ treatment * species + (1|genotype)`. Response is a measured value (rather than a count), `treatment` consists of three categories, `species` consists of four categories, and `genotype` consists of three categories. Technically, each genotype is specific to the species, so could be nested, but I think the model works out the same either way.

The model I have tried to set up is:

```julia
@model function randomintercept_regression(y,X,idx; n_gr=length(unique(idx)), predictors=size(X, 2))
  #priors, depends on data that has been rescaled
  α ~ Normal(0, 10) # population-level intercept
  β ~ filldist(Normal(0,10), predictors)
  σ ~ Exponential(10) # residual SD
  #prior for variance of random intercepts
  #usually requires thoughtful specification
  τ ~ truncated(Cauchy(0, 10); lower=0) #group level SDs
  αⱼ ~ filldist(Normal(0, τ), n_gr) # group level intercepts
  #likelihood
  ŷ = α .+ X * β .+ αⱼ[idx] 
  y ~ MvNormal(ŷ, σ^2 * I)
end

```

where `idx` is the list of genotypes coded as integers, `X` is a matrix of the predictor variables for each individual in the study (coded as 1 and 0, with each column representing a different level of either `treatment`, `species`, or the interaction between the two) and `y` represents the measured response of the individual. The fact that my TuringGLM model returns the correct parameter estimates suggests to me that it isn’t an issue with the structure of my data, rather how the model is set up. I have all the code associated with this project [here](https://github.com/PeetoomHeida/cereal_induction), including code for data generation, should you want to do a really deep dive. Much of the code is adapted from [this resource](https://storopoli.github.io/Bayesian-Julia/pages/11_multilevel_models/).

Thanks for your help!

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [December 19, 2022, 9:38am UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818/2 "2022-12-19T09:38:53Z")

</div>

I don’t think you are miss-specifying your `Turing.jl` model.  
`TuringGLM.jl` uses different priors as default priors, which takes into account the underlying data, that might be the culprit here.

Can you post the estimates `TuringGLM.jl` versus the `Turing.jl` model you’ve posted?

One other interesting idea would be to normalize all you data `X` and response `y` using `ZScoreTransform` (take a look here: [Bayesian-Statistics/02-linear\_regression-kidiq.jl at main · storopoli/Bayesian-Statistics · GitHub](https://github.com/storopoli/Bayesian-Statistics/blob/main/turing/02-linear_regression-kidiq.jl#L16-L21)) and proceed with `Normal(0, 2.5)` priors for `β` and `α`, and `Exponential(1)` for `σ`.

---

<div class="post-metadata">

**Author:** ![PeetoomHeida](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/peetoomheida/32/23865_2.png) [@PeetoomHeida](https://discourse.julialang.org/u/PeetoomHeida)\
**Post date:** [December 22, 2022, 7:02pm UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818/3 "2022-12-22T19:02:39Z")

</div>

Sure thing. Here are the estimates for the `Turing.jl` model:

```julia
22×8 DataFrame
 Row │ parameters mean std naive_se mcse ess rhat ess_per_sec 
     │ Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64     
─────┼────────────────────────────────────────────────────────────────────────────────────────────────────
   1 │ α -3.45131 5.76711 0.0911859 0.252182 230.766 1.02397 1.63194
   2 │ treat1_Control 0.813562 4.82649 0.0763135 0.211874 255.586 1.01681 1.80746
   3 │ treat2_Triticale -1.7054 4.33897 0.0686052 0.193385 365.632 1.02321 2.58569
   4 │ treat2_Wheat -1.96195 4.34805 0.0687488 0.197807 318.722 1.027 2.25395
   5 │ treat2_Barley 2.43193 4.34733 0.0687373 0.197601 320.813 1.02696 2.26874
   6 │ treat2_Oats -2.23604 4.35545 0.0688658 0.19943 308.04 1.02776 2.17841
   7 │ treat1_MeJA -3.31079 4.82681 0.0763186 0.211433 256.962 1.01672 1.81719
   8 │ treat1_Insect -3.10834 4.82561 0.0762996 0.211556 257.02 1.01669 1.8176
   9 │ σ 1.07758 0.0523837 0.00082826 0.00407205 20.7791 1.28952 0.146946
  10 │ τ 0.152454 0.139013 0.002198 0.0155327 11.6907 1.71629 0.0826748
  11 │ αⱼ[1] -0.000636618 0.124812 0.00197345 0.00809698 52.625 1.08083 0.372155
  12 │ αⱼ[2] -0.0461119 0.107025 0.00169221 0.00487228 83.2762 1.05705 0.588916
  13 │ αⱼ[3] 0.0772738 0.152267 0.00240756 0.0137612 16.1005 1.41414 0.11386
  14 │ αⱼ[4] -0.0550869 0.149021 0.00235623 0.0141615 14.2971 1.51724 0.101107
  15 │ αⱼ[5] -0.0663971 0.170312 0.00269288 0.016902 13.7304 1.55641 0.0970988
  16 │ αⱼ[6] 0.0555591 0.125597 0.00198586 0.00879491 30.5549 1.17043 0.216079
  17 │ αⱼ[7] 0.0147474 0.11526 0.00182243 0.0077317 53.4461 1.08289 0.377962
  18 │ αⱼ[8] -0.0221508 0.161494 0.00255344 0.0157733 14.5065 1.49708 0.102588
  19 │ αⱼ[9] -0.00909111 0.099925 0.00157995 0.00275826 1996.17 1.00353 14.1166
  20 │ αⱼ[10] -0.0792182 0.156167 0.00246922 0.0152744 14.0818 1.49238 0.0995839
  21 │ αⱼ[11] 0.00108857 0.197207 0.00311812 0.0201408 15.2257 1.44926 0.107674
  22 │ αⱼ[12] -0.0331863 0.1292 0.00204284 0.0103636 19.8629 1.29653 0.140467

```

And here are the estimates for the `TuringGLM.jl` model:

```julia
20×8 DataFrame
 Row │ parameters mean std naive_se mcse ess rhat ess_per_sec 
     │ Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64     
─────┼──────────────────────────────────────────────────────────────────────────────────────────────
   1 │ α -0.30134 0.185182 0.00292799 0.00596053 1051.16 1.00139 32.5437
   2 │ β[1] -4.62083 0.216361 0.00342096 0.00692687 1030.89 1.00269 31.9162
   3 │ β[2] -4.0301 0.226085 0.00357472 0.00873237 682.446 1.00343 21.1284
   4 │ β[3] -4.32156 0.223349 0.00353146 0.00798252 881.987 1.00341 27.3061
   5 │ β[4] -3.88406 0.162535 0.00256991 0.00304312 3126.17 0.99978 96.7853
   6 │ β[5] -4.08865 0.164911 0.00260747 0.00377165 2972.37 1.00138 92.0239
   7 │ σ 1.05227 0.0476444 0.000753325 0.000765722 4253.41 1.0003 131.684
   8 │ τ 0.102047 0.0878335 0.00138877 0.00380113 529.434 1.0037 16.3911
   9 │ zⱼ[1] -0.0284821 0.919094 0.0145322 0.0130367 4590.23 0.999867 142.112
  10 │ zⱼ[2] -0.283288 0.981594 0.0155204 0.0151085 3973.92 0.999478 123.031
  11 │ zⱼ[3] 0.187311 0.937938 0.0148301 0.0128819 4123.02 1.00007 127.648
  12 │ zⱼ[4] -0.117119 0.958306 0.0151521 0.0140993 4466.12 0.999459 138.27
  13 │ zⱼ[5] -0.222428 0.966174 0.0152765 0.0150198 3850.42 1.00036 119.208
  14 │ zⱼ[6] 0.22366 0.952422 0.0150591 0.0138048 4170.4 0.999588 129.115
  15 │ zⱼ[7] 0.121774 0.956345 0.0151211 0.0188772 3589.15 0.999184 111.119
  16 │ zⱼ[8] 0.13658 0.942066 0.0148954 0.0166025 3821.37 1.0002 118.309
  17 │ zⱼ[9] -0.0184732 0.939687 0.0148578 0.0169186 3493.39 0.999854 108.154
  18 │ zⱼ[10] -0.236547 0.923681 0.0146047 0.0138926 4463.18 0.999422 138.179
  19 │ zⱼ[11] 0.151902 0.969382 0.0153273 0.0145629 4167.17 1.00028 129.014
  20 │ zⱼ[12] -0.0126014 0.933396 0.0147583 0.0146358 4166.88 0.999517 129.006

```

Based on my sleuthing it looks like the TuringGLM model is incorporating the `Barley` and `Control` as part of the intercept, and then doing comparative estimates. But I back-calculated them at the math was pretty close.

For reference the parameters used in the data generation were:

```julia
19×2 DataFrame
 Row │ parameter coeff     
     │ String15 Float64   
─────┼───────────────────────
   1 │ Control 0.0
   2 │ Insect -4.10827
   3 │ MeJA -4.17167
   4 │ Barley 0.0
   5 │ Oats -4.97833
   6 │ Triticale -4.11681
   7 │ Wheat -4.51245
   8 │ Barley1 0.634788
   9 │ Barley2 0.731347
  10 │ Barley3 0.620631
  11 │ Oats1 0.262506
  12 │ Oats2 0.480827
  13 │ Oats3 0.681654
  14 │ Triticale1 0.956313
  15 │ Triticale2 0.206814
  16 │ Triticale3 0.603935
  17 │ Wheat1 0.786783
  18 │ Wheat2 0.605044
  19 │ Wheat3 0.166197

```

How would normalization work with predictors that are purely 1 or 0? I looked into normalization, but wasnt sure if it would make a difference if the variables were dummy coded.

---

<div class="post-metadata">

**Author:** ![PeetoomHeida](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/peetoomheida/32/23865_2.png) [@PeetoomHeida](https://discourse.julialang.org/u/PeetoomHeida)\
**Post date:** [December 23, 2022, 8:22pm UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818/4 "2022-12-23T20:22:16Z")

</div>

I tried again, this time adjusting the X matrix to exclude `Barley` and `Control` so they would be incorporated into the intercept. The new parameter estimates are:

```julia
20×8 DataFrame
 Row │ parameters mean std naive_se mcse ess rhat ess_per_sec 
     │ Symbol Float64 Float64 Float64 Float64 Float64 Float64 Float64     
─────┼───────────────────────────────────────────────────────────────────────────────────────────────
   1 │ α -0.218272 0.176584 0.00279204 0.00410713 1599.31 1.00187 18.9299
   2 │ β[1] -4.10688 0.221082 0.00349562 0.00534045 1838.74 1.00322 21.7639
   3 │ β[2] -4.40813 0.216661 0.00342572 0.00576301 1662.02 1.00416 19.6721
   4 │ β[3] -4.69131 0.21747 0.0034385 0.00515242 1814.77 1.00265 21.4801
   5 │ β[4] -4.13625 0.166991 0.00264036 0.00626846 322.938 1.01617 3.82238
   6 │ β[5] -3.92708 0.161498 0.00255351 0.0046016 1530.07 1.00321 18.1104
   7 │ σ 1.05389 0.0513702 0.000812234 0.0023493 120.992 1.02733 1.4321
   8 │ τ 0.111642 0.0884269 0.00139815 0.00555255 241.995 1.01163 2.86432
   9 │ αⱼ[1] -0.000398632 0.115827 0.00183138 0.00189341 3743.89 1.0005 44.3137
  10 │ αⱼ[2] -0.0353279 0.126083 0.00199355 0.00288435 1894.11 1.00021 22.4192
  11 │ αⱼ[3] 0.0359535 0.124245 0.00196448 0.00310668 1929.02 0.999885 22.8324
  12 │ αⱼ[4] -0.00971979 0.114741 0.00181422 0.00177179 2993.38 0.99945 35.4304
  13 │ αⱼ[5] -0.0272379 0.126527 0.00200056 0.00251181 2555.22 1.00055 30.2444
  14 │ αⱼ[6] 0.0423152 0.126711 0.00200348 0.00344255 1478.26 1.00011 17.4971
  15 │ αⱼ[7] 0.00848105 0.114706 0.00181365 0.0021953 3374.93 0.999696 39.9467
  16 │ αⱼ[8] 0.0131566 0.116819 0.00184707 0.00195295 2963.07 0.999978 35.0717
  17 │ αⱼ[9] -0.0131799 0.116765 0.00184621 0.00276529 2158.81 1.00227 25.5522
  18 │ αⱼ[10] -0.0341412 0.119848 0.00189497 0.00311182 1546.84 1.00051 18.3088
  19 │ αⱼ[11] 0.0250325 0.116745 0.0018459 0.00275691 2733.87 1.00108 32.3589
  20 │ αⱼ[12] 0.000353223 0.116472 0.00184158 0.00219533 3462.41 0.999878 40.9821

```

Which are much closer to the TuringGLM model, but have some slight differences.

---

<div class="post-metadata">

**Author:** ![Storopoli](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/storopoli/32/209278_2.png) [@Storopoli](https://discourse.julialang.org/u/Storopoli)\
**Post date:** [December 25, 2022, 7:31pm UTC](https://discourse.julialang.org/t/translating-a-turingglm-model-to-a-turing-model/91818/5 "2022-12-25T19:31:59Z")

</div>

The number of parameters are different. I think that TuringGLM is coding the categorical species using K-1 categories, try to use the same categorical coding for the plain Turing data.