# Simple PPC for logistic regression

**URL:** <https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749>\
**Category:** Modeling\
**Created:** [August 30, 2020, 7:31pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749 "2020-08-30T19:31:56Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![dkaplan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dkaplan/32/13835_2.png) [@dkaplan](https://discourse.mc-stan.org/u/dkaplan)\
**Post date:** [August 30, 2020, 7:31pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/1 "2020-08-30T19:31:56Z")

</div>

Greetings all,

A while back I was asking about the simplest way to obtain a posterior predictive check in logistic regression. The ppc\_error\_binned function in Bayesplot doesn’t seem to work well as I think it requires proportions and not binary 1/0s. Any suggestions would be most appreciated.

Thanks,

David

My modelString is

modelString ="  
data {  
int \<lower=0\> n;  
vector[n] age;  
vector[n] male;  
int \<lower=0,upper=1\> chd[n];  
}  
parameters {  
real alpha;  
real beta1;  
real beta2;  
}

transformed parameters {  
real \<lower=0\> oddsbeta1;  
real \<lower=0\> oddsbeta2;

oddsbeta1 = exp(beta1);  
oddsbeta2 = exp(beta2);

}

model {  
for (i in 1:n) {  
chd[i] ~ bernoulli\_logit(alpha + beta1_male[i] + beta2_age[i]);  
}

// Priors  
alpha ~ normal(0, .01);  
beta1 ~ normal(2, .01);  
beta2 ~ normal(.5, .01);  
}

// To obtain the loo for logistic regression

generated quantities {  
vector[n] log\_lik;  
int chd\_rep[n];  
for (i in 1:n) {  
chd\_rep[i] = bernoulli\_logit\_rng(alpha + beta1_male[i] + beta2_age[i]);  
log\_lik[i] = bernoulli\_logit\_lpmf(chd[i] | alpha + beta1_male[i] + beta2_age[i]);  
}  
}  
"

```

```

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [August 31, 2020, 5:17am UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/2 "2020-08-31T05:17:19Z")

</div>

Hi David,

Do the responses in this thread cover what you’re after?

> [@Posterior Prediction from Logit Regression](https://discourse.mc-stan.org/t/posterior-prediction-from-logit-regression/12217):
>
> 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 predictiv…](https://discourse.mc-stan.org/t/posterior-predictive-check-for-binomial-regression/7203)

---

<div class="post-metadata">

**Author:** ![dkaplan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dkaplan/32/13835_2.png) [@dkaplan](https://discourse.mc-stan.org/u/dkaplan)\
**Post date:** [September 1, 2020, 1:06am UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/3 "2020-09-01T01:06:53Z")

</div>

Hi Andrew,

The code I sent you got a bit messed up in the copying. In any case, the links that were in the previous post about PPCs for logistic regression were not helpful. Here is the code

modelString ="  
data {  
int \<lower=0\> n;  
vector[n] age;  
vector[n] male;  
int \<lower=0,upper=1\> chd[n];  
}  
parameters {  
real alpha;  
real beta1;  
real beta2;  
}

transformed parameters {  
real \<lower=0\> oddsbeta1;  
real \<lower=0\> oddsbeta2;

oddsbeta1 = exp(beta1);  
oddsbeta2 = exp(beta2);

}

model {  
for (i in 1:n) {  
chd[i] ~ bernoulli\_logit(alpha + beta1_age[i] + beta2_male[i]);  
}

// Priors  
alpha ~ normal(0, .01);  
beta1 ~ normal(2, .01);  
beta2 ~ normal(.5, .01);  
}

generated quantities {  
real chd\_rep[n];  
for (i in 1:n) {  
chd\_rep[i] = bernoulli\_logit\_rng(alpha + beta1_age[i] + beta2_male[i]);  
}

}  
"

Everything runs just fine, but when I run the command

ppc\_stat(y = chd, yrep = chd\_rep, stat = “mean”)

I get the following message

Error: object ‘chd\_rep’ not found

Thanks for looking at this.

David

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 1, 2020, 2:07am UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/4 "2020-09-01T02:07:45Z")

</div>

Hi David,

The error:

```
Error: object ‘chd_rep’ not found

```

Implies that you haven’t extracted the `chd_rep` estimates from the stanfit object before calling `ppc_stat`.

So given a stan model fitted via:

```
fitt_mod = stan(model, data)

```

You need to extract `chd_rep` into a matrix before passing it to `ppc_stat`, you can do this via:

```
chd_rep = extract(fitt_mod,"chd_rep")$chd_rep
ppc_stat(y=chd,yrep=chd_rep)

```

---

<div class="post-metadata">

**Author:** ![dkaplan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dkaplan/32/13835_2.png) [@dkaplan](https://discourse.mc-stan.org/u/dkaplan)\
**Post date:** [September 1, 2020, 1:29pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/5 "2020-09-01T13:29:49Z")

</div>

Hi Andrew,

That did the trick, but I’m still left with the problem of finding the correct way to plot chd against chd\_rep. Every plot I try in ppc\_ gives very strange figures. chd is coded 1/0, but the model should give the logit of chd and the logit of chd\_rep, so a nice simple ppc density or histogram should come out. For example, when I run the following

chd\_rep = extract(myfitCHD,“chd\_rep”)$chd\_rep  
ppc\_dens\_overlay(y = chd, yrep = chd\_rep)

I get the following plot which seems strange to me.

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/8/87eebb38b394a1c1c9f89aeee68114bccf9d1ad2.png)

This is what I meant in my first email about not finding a good way to display ppc’s for logistic regression.

Thanks so much for your help and patience.

David

---

<div class="post-metadata">

**Author:** ![andrjohns](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/andrjohns/32/15297_2.png) [@andrjohns](https://discourse.mc-stan.org/u/andrjohns)\
**Post date:** [September 1, 2020, 2:03pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/6 "2020-09-01T14:03:16Z")

</div>

No worries! This is an area where @jonah would have some good advice. Jonah is there a recommended way to set up the ppc plotting for binary data from a binomal model?

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [September 1, 2020, 6:40pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/7 "2020-09-01T18:40:50Z")

</div>

Yeah if you have a binary outcome it can be tricky to think of how to plot it. It’s easier to make plots on the logit scale if you have a binomial outcome with number of trials \> 1. But for binary data there are some other options.

Using an rstanarm model as an example:

```
library(bayesplot)
library(rstanarm)

# switch is 0/1
fit <- stan_glm(switch ~ I(dist/100) + log(arsenic), data = wells, family = binomial)
y <- wells$switch
yrep <- posterior_predict(fit) # draws from the posterior predictive distribution 

```

One thing you can look at is Pr(y = 1) by comparing the proportion of 1s in the data vs the proportions of 1s in the posterior predictive distribution:

```
ppc_stat(y, yrep, stat = "mean") 

```

You can look at the distribution of other functions too using `ppc_stat`. The `stat` argument can be any function that takes in a vector and returns a scalar so, e.g. `stat = "sd"` would work too.

If you want to look at the distribution of 0s and 1s itself you can use `ppc_bars()`:

```
ppc_bars(y, yrep)

```

There are also the `ppc_stat_grouped()` and `ppc_bars_grouped()` functions that let you make each of these plots for separate levels of a grouping variable, if that exists in your model.

---

<div class="post-metadata">

**Author:** ![dkaplan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dkaplan/32/13835_2.png) [@dkaplan](https://discourse.mc-stan.org/u/dkaplan)\
**Post date:** [September 1, 2020, 8:14pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/8 "2020-09-01T20:14:56Z")

</div>

Hi Jonah,

That really helps. So, when I see a plot like this

 ![image](https://canada1.discourse-cdn.com/flex030/uploads/mc_stan/original/2X/0/0ca4d43d97fb94cdbb273b41c1c4ced1e2ea16e7.png)

it is telling me that while have about 43% of the sample with heart disease, the model predicts 100% of the sample to have chd. Fairly high number of false positives.

Thanks

David

---

<div class="post-metadata">

**Author:** ![jonah](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/jonah/32/200_2.png) [@jonah](https://discourse.mc-stan.org/u/jonah)\
**Post date:** [September 1, 2020, 9:15pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/9 "2020-09-01T21:15:13Z")

</div>

Glad that’s helpful. Yeah it looks like the sample is somewhat balanced but the model is predicting all 1s and no 0s. Presumably the value of

```stan
inv_logit(alpha + beta1*age[i] + beta2*male[i])

```

is coming out to 1 or nearly 1 every time it’s calculated. One thing to consider is that this prior with mean 2 and SD 0.01

> [@dkaplan](#):
>
> beta1 ~ normal(2, .01);

is super informative, implying that you pretty much know for sure that `beta1` is between 1.97 and 2.03. If that’s true from the prior information you have then that’s fine, but 2 is pretty big on the logit scale so this is probably contributing to the probability being so close to 1.

---

<div class="post-metadata">

**Author:** ![dkaplan](https://yyz2.discourse-cdn.com/flex030/user_avatar/discourse.mc-stan.org/dkaplan/32/13835_2.png) [@dkaplan](https://discourse.mc-stan.org/u/dkaplan)\
**Post date:** [September 1, 2020, 9:35pm UTC](https://discourse.mc-stan.org/t/simple-ppc-for-logistic-regression/17749/10 "2020-09-01T21:35:22Z")

</div>

And that’s because those are not the write priors. :-/ Everything looks more sensible now. Thanks.

David