Yrurtj

TL;DR: Which $N$ is correct for BIC in logistic regression, the aggregated binomial or Bernoulli $N$?

Suppose I have a data set to which I'd like to apply logistic regression. For the sake of example, suppose there are $j=5$ groups with $m=100$ participants each, for a total $n=500$. The outcome is 0 or 1. For example, the following data set (R code):

library(dplyr)



set.seed(44)

d <- tibble(y = rbinom(500, 1, .5),

            x = factor(rep(LETTERS[1:5], each = 100))) 

There are two ways I can represent this: as is, above, treating every observation as a Bernoulli random variable, or aggregating observations within groups and treating each observation as Binomial. The number of rows in the data set will be 500 in the first instance, and 5 in the second.

I can construct the aggregated data set:

d %>% 

  group_by(x, y) %>% 

  summarise(n = n()) %>%

  spread(y, n) %>%

  rename(f = `0`, s = `1`) %>%

  mutate(n = s + f) -> d_agg

I can then fit the logistic regression using both data sets in R:

g_bern  <- glm(y ~ x,          data=d,     family=binomial)

g_binom <- glm(cbind(s,f) ~ x, data=d_agg, family=binomial)

and compute the AIC:

AIC(g_bern)   # [1] 693.8487

AIC(g_binom)  # [1] 35.21523

which of course differ by a constant

2*sum(lchoose(d_agg$n, d_agg$s))  # [1] 658.6335

as expected (see: Logistic Regression: Bernoulli vs. Binomial Response Variables).

However, the BICs differ by that constant AND a factor that depends on the "number of observations", and the number of observations differ in each:

BIC(g_bern)    # [1] 714.9217

BIC(g_binom)   # [1] 33.26242

nobs(g_bern)   # [1] 500

nobs(g_binom)  # [1] 5

Just to confirm, we can recalculate BIC for both:

-2*logLik(g_bern) + attr(logLik(g_bern),"df")*log(nobs(g_bern))

# 'log Lik.' 714.9217 (df=5)

-2*logLik(g_binom) + attr(logLik(g_binom),"df")*log(nobs(g_binom))

# 'log Lik.' 33.26242 (df=5)

and indeed the only place these two numbers differ is $N$.

This surprises me, since I would think that R would "know" which of the two to use to prevent ambiguity. It has the same information in both cases.

Which one is "right"? Or is BIC really this arbitrary?

The BIC (and the AIC) are relative measures for comparing models. However, it makes no sense to compare what is otherwise the same model between using an aggregated vs. a disaggregated response. Nor would it make sense to compare models that would otherwise be different (e.g., different regressors), but where one model uses an aggregated response and the other model uses a disaggregated version of the response. As long as the two models being compared both represent the response variable in the same format, everything will be fine. Note that the two formats are ultimately equivalent—they contain the same information and mostly just look different on the outside, see: Input format for response in binomial glm in R.

Interesting question! Coming at this from an applied setting, I think you need to remember that both BIC and AIC are measures of relative model fit.

In other words, these measures don't tell you much when you examine them for a single model, but can help you to select an appropriate model among a set of competing models. In particular:

1. If your goal is to find the 'best' among those competing models for prediction of the outcome variable, then select the model with the lowest AIC value;


2. If your goal is to find the 'best' among those competing models for understanding and describing the effects of the predictor variables included in the model on the outcome variable, then select the model with the lowest BIC value.

In defining your set of competing models, you would have to make sure the models follow the same conceptual framework. Thus, you would either compare several binomial logistic regression models or several binary logistic models, but not a mixture of both. (It is important to compare like with like, otherwise you won't know if a model won the competition based on its own merits or simply because you changed the model specification/fitting procedure.)

From this perspective, the only thing that matters is that R is consistent when computing the AIC and BIC across models of the same type (e.g., binomial logistic regression models).

Just to clarify: g_bern is a binary logistic regression model, whereas g_binom is a binomial logistic regression model. While they both model the probability of success in one trial, you wouldn't mix together variations of these models when defining your set of competing models (for the reasons explained above and also covered by @gung).