# Posterior Prediction from Logit Regression

**URL:** <https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217>\
**Category:** General\
**Created:** [December 7, 2019, 2:19pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217 "2019-12-07T14:19:50Z")\
**Posts on this page:** 9\
**Page:** 1

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [December 7, 2019, 2:19pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/1 "2019-12-07T14:19:50Z")

</div>

I can certainly sympathise with the apparent desperation in this post ([https://stats.stackexchange.com/questions/281235/posterior-predictive-checking-for-bernoulli-distributed-data](https://stats.stackexchange.com/questions/281235/posterior-predictive-checking-for-bernoulli-distributed-data)). It seems like when you ask how to do posterior predictive checks for a logistic regression you get ‘look at this vignette’ or snippets of code that get you close to where you need to be by showing you how to generate the draws from the posterior (see [PPC for simple logistic regression in Stan](https://discourse.mc-stan.org/t/ppc-for-simple-logistic-regression-in-stan/6274) and [Posterior predictive check for binomial regression](https://discourse.mc-stan.org/t/posterior-predictive-check-for-binomial-regression/7203)) but only take you that far. I fully understand that the wonderful people at the Stan project administer this website in their spare time and thus don’t have time to work up a full answer. So I have no right to expect details, but I sure would like some.

**Can someone help me by actually showing me how to graph a PPC with some data and code I provide?**

In this dataset there are three groups of caffeine-deprived coffee drinkers. One group was given water, one was given decaf, and one was given caffeinated coffee. Each participant indicated whether their caffeine withdrawal was alleviated after they drank the beverage they were given (response = “feelBetter” = 1) or not (response = “noEffect”).

Here is the data

```
library(rstan)
library(dplyr)
df <- data.frame(id = 1:62, 
                 group = factor(c(rep("water", 20),
                                  rep("decaf", 22),
                                  rep("caffeine", 20)), 
                                levels = c("water", "decaf", "caffeine")),
                 effect = factor(c(rep("noEffect",16), rep("feelBetter",4),
                                   rep("noEffect",16), rep("feelBetter",6),
                                   rep("noEffect",4), rep("feelBetter",16)),
                                 levels = c("noEffect", "feelBetter")))

```

So, taking group into account, how many people expected that their beverage would make them feel better?

```
df %>% group_by(group, effect) %>% 
       summarise(count = n()) %>% 
       mutate(perc = round(100*count/sum(count),2))

# output

# group effect count perc
# 1 water noEffect 16 80  
# 2 water feelBetter 4 20  
# 3 decaf noEffect 16 72.7
# 4 decaf feelBetter 6 27.3
# 5 caffeine noEffect 4 20  
# 6 caffeine feelBetter 16 80 

```

I ran the model as a logistic regression with levels-means coding of group factors (each group getting their own intercept parameter). Here is the model code.

```
dList <- list(N = nrow(df),
              nGroup = nlevels(df$group),
              gIndex = as.integer(df$group),
              effect = as.integer(df$effect) - 1) 

ms <- "
  data{
    int<lower=0> N;
    int<lower=0> nGroup;
    int<lower=1> gIndex[N];
    int<lower=0> effect[N];
  }

  parameters{
    vector [nGroup] a;
  }

  model{
    a ~ normal(0,10);
    for (i in 1:N) {
    effect[i] ~ bernoulli_logit(a[gIndex[i]]);
    }
  }
 
  generated quantities{
    real y_draw[N];
    for (i in 1:N) {
      y_draw[i] = bernoulli_logit_rng(a[gIndex[i]]);
    }
  }
"

postExp <- stan(model_code = ms,
                data = dList,
                warmup = 1e3,
                iter = 4e3,
                cores = 1,
                chains = 1, 
                seed = 2334)

```

The summary…

```
print(postExp, pars = "a", probs = c(0.025, 0.975))

# output
# mean se_mean sd 2.5% 97.5% n_eff Rhat
# a[1] -1.50 0.01 0.58 -2.72 -0.46 2172 1
# a[2] -1.04 0.01 0.48 -2.02 -0.12 2415 1
# a[3] 1.46 0.01 0.58 0.44 2.68 1892 1

```

…all looks about right. I know this because when I convert these group log odds to probabilities the model predicted proportions of people saying their beverage made them feel better matches the raw probabilities from the descriptive statistics above pretty closely

```
logOdds <- as.data.frame(summary(postExp, pars = "a", probs = c(0.025, 0.975))$summary)$mean
logitToProb <- function(lo) exp(lo)/(1+exp(lo))
groupProbs <- sapply(1:3, function (i) round(logitToProb(logOdds[i])*100,1))
names(groupProbs) = c("water", "decaf", "caffeine")
groupProbs

# output
# water decaf caffeine 
# 18.2 26.2 81.1 

```

So the model seems to be working. But I need PPCs. So how to do them?

In order to generate the draws from the joint posterior I extracted the draws from the generated quantities block.

```
post_draw <- data.frame(rstan::extract(postExp, pars = "y_draw"))

```

There are 3000 draws from the joint posterior for each observed data point, which in this case is each participant. That is 3000 0s and 1s, hopefully in the correct group proportions for each participant within each group

```
head(post_draw[,1:6])

# output
 # y_draw.1 y_draw.2 y_draw.3 y_draw.4 y_draw.5 y_draw.6
 # 1 0 0 0 0 0
 # 0 0 0 0 1 0
 # 0 0 0 1 0 0
 # 1 0 1 0 0 1
 # 0 1 1 0 0 0
 # 0 0 0 0 0 0

```

So great, I have what I need, but **what do I do with this matrix of draws from the posterior?**

I can take means and get upper and lower HPDIs for each row of the matrix and then bind that to the actual data, but how to graph it?

I tried the `ppc_error_binned` function from the bayesplot package but to be honest I didn’t really know what I was doing, since the vignette’s all seem to do with questions that don’t really suit my problem (e.g. poisson regression or logit regression where each value of the outcome is a count of successes within several trials rather than a single trial [0,1] value for each participant).

Again an example PPC figure based on my own data with the code and an explanation would be very much appreciated. Again I have no right to expect this, just a hope.

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 7, 2019, 2:35pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/2 "2019-12-07T14:35:42Z")

</div>

> [@llewmills](#):
>
> I can certainly sympathise with the apparent desperation in this post ([https://stats.stackexchange.com/questions/281235/posterior-predictive-checking-for-bernoulli-distributed-data](https://stats.stackexchange.com/questions/281235/posterior-predictive-checking-for-bernoulli-distributed-data)).

The answer to their question is “yes”—you want to look at the y\_rep[i], and with binary data, all you get are counts.

> [@llewmills](#):
>
> the wonderful people at the Stan project administer this website in their spare time

Thanks for the vote of confidence, but we don’t quite deserve that much sympathy. Many of us, including me, are paid to work full-time on Stan (well, not quite full-time for me—I have to spend a fraction of my time writing grant proposals to keep the team employed at Columbia Uni).

I’d suggest using parallel naming for PPC simulations. In this case, use `effect` and `effect_sim` or `y` and `y_sim`. Each run of the Stan program will give you `y_sim` values. What I would do is plot `sum(y_sim)` as a histogram over multiple runs, then place a vertical line at `sum(y)`. The posterior p-value will be the proportion of `sum(y_sim)` greater than (or equal to—this is discretes, so it matters) `sum(y)`. That shouldn’t be too close to 0 or 1.

> [@llewmills](#):
>
> Can someone help me by actually showing me how to graph a PPC

You need an R expert for help with R. @jonah will know how you can do it.

---

<div class="post-metadata">

**Author:** ![FJCC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fjcc/32/6653_2.png) [@FJCC](https://discourse.mc-stan.org/u/FJCC)\
**Post date:** [December 7, 2019, 6:06pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/3 "2019-12-07T18:06:27Z")

</div>

Try appending these lines to your code.

```
library(bayesplot)
postDF <- as.data.frame(postExp)

ppc_stat_grouped(y = as.integer(df$effect) - 1, yrep = as.matrix(postDF[4:65]), group = df$group, stat = mean)

```

The postDF[4:65] picks out all ot the y\_draw[n] columns from the data frame of the posterior, so _yrep_ is a 3000 x 62 matrix. I should have used your post\_draw matrix in that place! The _group_ parameter labels each of those columns with its corresponding treatment.

 ![PostPredChk](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/3/34dc8995c4ce4fb35ccc636a1533068e3992318c.jpeg)

---

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [December 7, 2019, 7:11pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/4 "2019-12-07T19:11:15Z")

</div>

Thank you @FJCC!. Can you explain what I’m looking at here?

---

<div class="post-metadata">

**Author:** ![FJCC](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/fjcc/32/6653_2.png) [@FJCC](https://discourse.mc-stan.org/u/FJCC)\
**Post date:** [December 7, 2019, 8:24pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/5 "2019-12-07T20:24:49Z")

</div>

The histograms show the distributions of the fraction of simulated data that resulted in feelBetter in each of the 3000 posterior draws. In your post\_draw matrix, the first 20 columns are simulations using the _a_ value for water (_a[1]_), the next 22 come from the _a_ value for decaf, and the remaining columns are caffeine results. If we take the mean value in one row of the first 20 columns, that gives the fraction of feelBetter results produced by simulating data with the _a[1]_ from that particular draw. Doing that for all 3000 rows provides a distribution of the fraction of feelBetter outcomes given all of the _a[1]_ values explored in the posterior, the inherent spread of a bernoulli process, and the fact that you observe 20 individuals.  
The graph shows that your observed results, the black vertical lines, are near the most likely results that the model predicts but that quite a range of results is reasonably probable.  
Below is a manual calculation of the same result starting with the postDF data frame from my code. I hope it is clearer than my attempt to explain the calculation. Note that the plot order is different than in the plot from bayesplot.

```
Ydraw <- postDF[4:65]
GroupMeans <- data.frame(water = rowMeans(Ydraw[, 1:20]), 
                         decaf= rowMeans(Ydraw[, 21:42]), 
                         caffeine = rowMeans(Ydraw[, 43:62]))
GroupMeansTall <- tidyr::gather(GroupMeans, key = "group", value = "Value")

ggplot(GroupMeansTall, aes(x = Value)) + geom_histogram(fill = "skyblue", color = "darkblue") +
  facet_wrap(~group)

```

 ![PostPredChk_manual](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/b/b006088107ea005f7e914fb84663ff08eeb4113e.jpeg)

---

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [December 7, 2019, 9:00pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/6 "2019-12-07T21:00:57Z")

</div>

Thank you so much @FJCC! You have made my day.

So the function takes the row means for all the columns of 3000 simulated draws assigned to observations from each group (the mean of all the 0s and 1s being the same as the proportion of 1s) and then plots each as a histogram? That makes sense!

Thank you again

---

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [December 7, 2019, 9:07pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/7 "2019-12-07T21:07:11Z")

</div>

Thank you @Bob Carpenter, both for your answer to my question and for shining a light on the recondite inner workings of the Stan project. I wish I could give the solution to both yours and FJCC’s answers.

Most grant applications involve justifying how the project will contribute to the greater good. I think that should be easy when it comes to stan and what you guys are doing at Columbia, it is such amazing work. Of course I am assuming the presence on the grant review panel of assessors people can take the broader view and who can see how better inferential statistics practices can benefit all of mankind.

---

<div class="post-metadata">

**Author:** ![llewmills](https://avatars.discourse-cdn.com/v4/letter/l/4da419/32.png) [@llewmills](https://discourse.mc-stan.org/u/llewmills)\
**Post date:** [December 8, 2019, 6:53am UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/8 "2019-12-08T06:53:55Z")

</div>

@Bob Carpenter FJCC suggested taking the group means of [0,1] draws from the posterior, but in your answer you suggested summing these. Is summing a better technique for posterior prediction than averaging in this case?

---

<div class="post-metadata">

**Author:** ![Bob\_Carpenter](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/bob_carpenter/32/9230_2.png) [@Bob\_Carpenter](https://discourse.mc-stan.org/u/Bob_Carpenter)\
**Post date:** [December 14, 2019, 7:16pm UTC](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217/9 "2019-12-14T19:16:33Z")

</div>

> [@llewmills](#):
>
> Is summing a better technique for posterior prediction than averaging in this case?

All averaging does is divide by the number of items, so the only effect is on the axis labels. The posterior p-values, for example, don’t change.