# Posterior prediction in Turing

**URL:** <https://discourse.julialang.org/t/posterior-prediction-in-turing/97355>\
**Category:** Probabilistic Programming\
**Tags:** turing\
**Created:** [April 11, 2023, 5:11pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355 "2023-04-11T17:11:41Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![jfhawkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jfhawkin/32/48098_2.png) [@jfhawkin](https://discourse.julialang.org/u/jfhawkin)\
**Post date:** [April 11, 2023, 5:11pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355/1 "2023-04-11T17:11:41Z")

</div>

I’m looking to replicate the below code snippet from a Stan example. I’m having a hard time figuring out the syntax to replicate the Turing predict() example for my logistic regression model and translating the Stan code to Julia/Turing. The issue seems to be translating the blank dependent, y, variable from the linear regression example into an equivalent for my logistic case. If I simply use, `y = Vector{Union{Missing, Int64}}(undef, n)`, converting the `Float64` into `Int64` from the example to match the binary 0/1 for a logistic model, it gives me a `MethodError` because there isn’t a Bernoulli model for ````Union{Missing, Int64}‘’'. I tested the linear model example seems to work on my version of Turing.

```julia
# Posterior_epred returns the posterior estimates for the different subgroups stored in the
# poststrat_df dataframe. 
epred_mat <- posterior_epred(fit, newdata = poststrat_df, draws = 1000)
mrp_estimates_vector <- epred_mat %*% poststrat_df$n / sum(poststrat_df$n)
mrp_estimate <- c(mean = mean(mrp_estimates_vector), sd = sd(mrp_estimates_vector))
cat("MRP estimate mean, sd: ", round(mrp_estimate, 3))

```

I’m working on the `predict()` example in Inference.jl, with key code from the example:

```julia
m_test = linear_reg(xs_test, Vector{Union{Missing, Float64}}(undef, length(ys_test)), σ);
predictions = predict(m_test, chain_lin_reg)
ys_pred = vec(mean(Array(group(predictions, :y)); dims = 1));

```

My Julia/Turing code to reproduce the Stan snippet is:

```julia
y = Vector{Union{Missing, Int64}}(undef, n) # This doesn't work (MethodError) as y but is based on Turing example in documentation for Inference.jl
y = ones(n) # This works as y
post_model = varying_intercept(
        state_idx, eth_idx, male_eth_idx, educ_idx, educ_age_idx, educ_eth_idx, y, X, n, standardize = false
) # Model runs for training data set
epred_mat = predict(post_model, chains)

vec(mean(Array(group(epred_mat, :y)); dims = 1)) # Not the correct output

```

---

<div class="post-metadata">

**Author:** ![Christopher\_Fisher](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/christopher_fisher/32/26132_2.png) [@Christopher\_Fisher](https://discourse.julialang.org/u/Christopher_Fisher)\
**Post date:** [April 12, 2023, 4:00pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355/2 "2023-04-12T16:00:19Z")

</div>

Can you post a runnable example of your problem? That would help others provide a solution.

---

<div class="post-metadata">

**Author:** ![jfhawkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jfhawkin/32/48098_2.png) [@jfhawkin](https://discourse.julialang.org/u/jfhawkin)\
**Post date:** [April 13, 2023, 2:00pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355/3 "2023-04-13T14:00:45Z")

</div>

Putting together a runnable example, I saw that the Turing case study for logistic regression has a prediction function I could adapt. The example can be run after installing RDatasets.jl: [https://turing.ml/dev/tutorials/02-logistic-regression/](https://turing.ml/dev/tutorials/02-logistic-regression/)

This gets me most of the way to replicating my Stan example. The output from the Stan function `posterior_epred` is an N \times M matrix where N is the number of chains and M is the total observations in the prediction data. A simplified version of my function that excludes variables that aren’t germane to the problem, based on the logistic example, is as follows:

```julia
function prediction(X::Matrix, chains;predictors=size(X, 2))
            
    # Compute the mean parameter value across chains/iterations
    inter = DataFrame(mean(group(chains,:inter)))[!, :mean]
    # Select parameter values for betas
    β = DataFrame(mean(group(chains,:β)))[!,:mean]

    v = inter .+ X * β

     return logistic.(v)
end;

```

I call this with:

```julia
sub_chains = chains[1001:end]
post_model = prediction(X, sub_chains)

```

This give the correct return in terms of format, but it isn’t quite right. The issue is that `mean` (and its parent function `summarize()`) aggregates across chains and samples. I need the mean across chains only, so \beta and inter should both be 1000 \times 1 vectors and not single values. Then, `logistic(v)` would return a matrix of probabilities, or I would loop through samples to generate vectors of probabilities. My thought was there might be a way to adapt `mean()` to only aggregate across the chains and not samples, but I couldn’t find a way to do it in the source code. Any ideas (probably from the developers of that function)?

---

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [April 13, 2023, 2:59pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355/4 "2023-04-13T14:59:58Z")

</div>

I’m not exactly sure if the following code does what you are looking for, but I thought I’d provide it just in case. I usually convert the chain to a `DataFrame` because I find it easier to work with. So, using the logistic regression tutorial, I would do something like this:

```julia
m = logistic_regression(train, train_label, n, 1)
chain = sample(m, NUTS(), MCMCThreads(), 1000, 4)
chaindf = DataFrame(chain)

```

This gives me a `DataFrame` that looks like this:

```julia
4000×18 DataFrame
  Row │ iteration chain intercept student balance income lp n_steps is_accept acceptance_rate log_density hamiltonian_energy ⋯
      │ Int64 Int64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 Float64 ⋯
──────┼───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
    1 │ 501 1 -4.78973 0.440116 1.8584 0.486572 -56.5809 7.0 1.0 0.962241 -56.5809 56.9433 ⋯
    2 │ 502 1 -4.87119 0.811197 1.98672 0.385192 -58.3785 3.0 1.0 0.893347 -58.3785 58.8032  
    3 │ 503 1 -4.14025 -1.30256 1.83855 -0.161454 -56.879 23.0 1.0 0.987357 -56.879 59.4473  
    4 │ 504 1 -4.45475 -0.122842 1.97072 0.35409 -55.7659 7.0 1.0 0.956346 -55.7659 57.3812  
    5 │ 505 1 -4.453 -0.654141 1.87193 0.693527 -57.6155 3.0 1.0 0.733881 -57.6155 58.1016 ⋯
    6 │ 506 1 -4.01181 -0.810257 1.89205 0.419372 -57.1889 7.0 1.0 0.994691 -57.1889 57.7374  
  ⋮ │ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋮ ⋱
 3996 │ 1496 4 -4.01397 0.274545 1.27778 0.783362 -59.1164 3.0 1.0 1.0 -59.1164 60.3322  
 3997 │ 1497 4 -4.01407 -0.260389 1.92296 0.460915 -57.7695 7.0 1.0 0.698575 -57.7695 63.512   
 3998 │ 1498 4 -4.19074 -0.636925 1.59508 0.471109 -56.8324 7.0 1.0 0.996946 -56.8324 58.3185 ⋯
 3999 │ 1499 4 -4.66654 -0.347037 2.1733 0.467141 -56.4114 7.0 1.0 0.953491 -56.4114 58.7883  
 4000 │ 1500 4 -4.84313 -0.751735 2.0547 0.0107546 -57.1793 15.0 1.0 0.921095 -57.1793 57.7937  
                                                                                                                       6 columns and 3989 rows omitted

```

Then, I would proceed with something like this:

```julia
p = [
    logistic(row.intercept + row.student * test[i, 1] + row.balance * test[i, 2] + row.income * test[i, 3])
    for row in eachrow(chaindf), i in 1:size(test, 1)
]

4000×9500 Matrix{Float64}:
 0.465484 0.00344039 0.00508513 0.00137203 0.261555 0.000450104 … 0.00421789 0.00069823 0.00792776 0.314582 0.0921833 0.013099
 0.488537 0.0032657 0.00365576 0.000998941 0.298431 0.000396564 0.00399184 0.000555768 0.00674533 0.476042 0.151211 0.0193158
 0.440289 0.0101318 0.00349736 0.00157082 0.400781 0.001867 0.0114167 0.0015384 0.0104211 0.229439 0.058257 0.00721661        
 0.574727 0.00507123 0.00532339 0.00150388 0.384724 0.00063604 0.00616453 0.000864896 0.0100408 0.35092 0.0966909 0.0118373
 0.611662 0.00417768 0.00978534 0.0022429 0.335778 0.000487701 0.00526206 0.000925815 0.0124358 0.150412 0.038205 0.00529938        
 0.654961 0.00772059 0.00966926 0.00271075 0.452906 0.00101429 … 0.00940061 0.00147208 0.0163597 0.248997 0.0660443 0.00869838        
 0.544748 0.00612234 0.00897271 0.00257814 0.335795 0.000884618 0.00743218 0.00134973 0.0135882 0.326597 0.103618 0.0164337
 0.351807 0.00370288 0.00620632 0.00185623 0.18436 0.000599739 0.00447022 0.000950212 0.00846804 0.154772 0.0469492 0.00833413        
 ⋮ ⋮ ⋱ ⋮
 0.422075 0.00426458 0.0135892 0.00354346 0.187335 0.0006838 0.00527318 0.00140718 0.0133127 0.367753 0.15215 0.0358649
 0.509708 0.00680076 0.0312533 0.00795952 0.222816 0.00117337 0.00846254 0.00278392 0.0247064 0.145091 0.0568705 0.0149864
 0.419215 0.00726688 0.0258561 0.00754298 0.18959 0.0014455 … 0.0088408 0.00302761 0.0223092 0.145294 0.0598602 0.0168975
 0.679729 0.0074264 0.0100926 0.00269873 0.467505 0.000927589 0.0091118 0.00140119 0.016628 0.369742 0.109152 0.0145805
 0.465057 0.00682179 0.0107435 0.00336271 0.274343 0.00116335 0.00817516 0.00178612 0.014866 0.139324 0.0427343 0.00780718        
 0.659689 0.00358769 0.00450103 0.00105151 0.431728 0.000349499 0.00449446 0.000528723 0.00839054 0.325864 0.0756088 0.0073663
 0.431515 0.00421881 0.00197272 0.000694956 0.33773 0.00057689 0.00494311 0.000559974 0.00552086 0.254671 0.056975 0.00570127

```

Finally,

```julia
threshold = 0.07
predictions = [mean(c) ≥ threshold ? 1 : 0 for c in eachcol(p)]

9500-element Vector{Int64}:
 1
 0
 0
 0
 1
 0
 0
 0
 ⋮
 0
 0
 0
 0
 1
 1
 0

```

---

<div class="post-metadata">

**Author:** ![jfhawkin](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jfhawkin/32/48098_2.png) [@jfhawkin](https://discourse.julialang.org/u/jfhawkin)\
**Post date:** [April 13, 2023, 6:10pm UTC](https://discourse.julialang.org/t/posterior-prediction-in-turing/97355/5 "2023-04-13T18:10:09Z")

</div>

Thanks! Recognizing that conversion to a `DataFrame` flattens the Array into 2D was the key. Looping over i,j pairs is fine. This is exactly what I needed to replicate the Stan implementation.