posterior predictive check for binomial glm with matrix response

[richardjtelford]

I've noticed that the posterior predictive check for a binomial glm is different depending on whether the response variable is a vector of proportions or a matrix of successes and failures.

Minimal reprex

tot <- rep(10, 100)
suc <- rbinom(100, prob = 0.9, size = tot)

df <- data.frame(tot, suc)

df$prop <- suc/tot

mod1 <- glm(cbind(suc, tot - suc) ~ 1 ,
            family = binomial,
            data = df)

#performance::check_model(mod1)
performance::check_posterior_predictions(mod1)

mod2 <- glm(prop ~ 1 ,
            family = binomial,
            data = df, 
            weights = tot)


performance::check_posterior_predictions(mod2)

I was, perhaps mistakenly, expecting these plots to look the same.

In performance::pp_check.glm, the result is calculated as

sapply(matrix_sim, function(i) i[, 1]/i[, 2], simplify = TRUE)

This is generating lots of Inf where the number of failures is zero.

Should this code not be calculating the proportion, with code such as

sapply(matrix_sim, function(i) i[, 1]/rowSums(i), simplify = TRUE)

[richardjtelford]

Thanks for having a look at this. I appear to have given a partial fix. Looking again line 288

  out$y <- response[, 1] / response[, 2]

would also need to become

  out$y <- response[, 1] / rowSums(response)

so that it was calculated in the same way

[richardjtelford]

Thank you. That looks how I would expect it to.

Would be nice to improve the x-axis label, but not sure what would work and be easy

[strengejacke]

After fixing a bug in insight, this is how it would look like with the current implementation, and your suggested fix.

set.seed(1)
tot <- rep(10, 100)
suc <- rbinom(100, prob = 0.9, size = tot)
df <- data.frame(tot, suc)
df$prop <- suc / tot

mod1 <- glm(cbind(suc, tot - suc) ~ 1,
  family = binomial,
  data = df
)

mod2 <- glm(prop ~ 1,
  family = binomial,
  data = df,
  weights = tot
)

mod3 <- glm(cbind(suc, tot) ~ 1,
  family = binomial,
  data = df
)

mod4 <- glm(am ~ 1,
  family = binomial,
  data = mtcars
)

Mod1

Curent (mod1)

New (mod1)

Mod2

Curent (mod2)

New (mod2)

Mod3

Curent (mod3)

New (mod3)

Mod4

Curent (mod4)

New (mod4)