Logistic Regression: Bernoulli vs. Binomial Response Variables

Question

I want to perform logistic regression with the following binomial response and with $X_1$ and $X_2$ as my predictors.

I can present the same data as Bernoulli responses in the following format.

The logistic regression outputs for these 2 data sets are mostly the same. The deviance residuals and AIC are different. (The difference between the null deviance and the residual deviance is the same in both cases - 0.228.)

The following are the regression outputs from R. The data sets are called binom.data and bern.data.

Here is the binomial output.

    Call:
    glm(formula = cbind(Successes, Trials - Successes) ~ X1 + X2, 
        family = binomial, data = binom.data)

    Deviance Residuals: 
    [1]  0  0  0

    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)  -2.9649    21.6072  -0.137    0.891
    X1Yes        -0.1897     2.5290  -0.075    0.940
    X2            0.3596     1.9094   0.188    0.851

    (Dispersion parameter for binomial family taken to be 1)

    Null deviance:  2.2846e-01  on 2  degrees of freedom
    Residual deviance: -4.9328e-32  on 0  degrees of freedom
    AIC: 11.473

    Number of Fisher Scoring iterations: 4

Here is the Bernoulli output.

    Call:
    glm(formula = Success ~ X1 + X2, family = binomial, 
               data = bern.data)

    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -1.6651  -1.3537   0.7585   0.9281   1.0108  

    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)  -2.9649    21.6072  -0.137    0.891
    X1Yes        -0.1897     2.5290  -0.075    0.940
    X2            0.3596     1.9094   0.188    0.851

    (Dispersion parameter for binomial family taken to be 1)

    Null deviance: 15.276  on 11  degrees of freedom
    Residual deviance: 15.048  on  9  degrees of freedom
    AIC: 21.048

    Number of Fisher Scoring iterations: 4

My questions:

1. I can see that the point estimates and standard errors between the 2 approaches are equivalent in this particular case. Is this equivalence true in general?

2. How can the answer for Question #1 be justified mathematically?

3. Why are the deviance residuals and AIC different?

Answer 1

1) Yes. You can aggregate/de-aggregate (?) binomial data from individuals with the same covariates. This comes from the fact that the sufficient statistic for a binomial model is the total number of events for each covariate vector; and the Bernoulli is just a special case of the binomial. Intuitively, each Bernoulli trial that makes up a binomial outcome is independent, so there shouldn't be a difference between counting these as a single outcome or as separate individual trials.

2) Say we have $n$ unique covariate vectors $x_1, x_2, \ldots, x_n$, each of which has a binomial outcome on $N_i$ trials, i.e. $$Y_i \sim \mathrm{Bin}(N_i, p_i)$$ You've specified a logistic regression model, so $$\mathrm{logit}(p_i) = \sum_{k=1}^K \beta_k x_{ik}$$ although we'll see later that this isn't important.

The log-likelihood for this model is $$\ell(\beta; Y) = \sum_{i=1}^n \log {N_i \choose Y_i} + Y_i \log(p_i) + (N_i - Y_i) \log(1-p_i)$$ and we maximise this with respect to $\beta$ (in the $p_i$ terms) to get our parameter estimates.

Now, consider that for each $i = 1, \ldots, n$, we split the binomial outcome into $N_i$ individual Bernoulli/binary outcomes, as you have done. Specifically, create $$Z_{i1}, \ldots, Z_{iY_i} = 1$$ $$Z_{i(Y_i+1)}, \ldots, Z_{iN_i} = 0$$ That is, the first $Y_i$ are 1s and the rest are 0s. This is exactly what you did - but you could equally have done the first $(N_i - Y_i)$ as 0s and the rest as 1s, or any other ordering, right?

Your second model says that $$Z_{ij} \sim \mathrm{Bernoulli}(p_i)$$ with the same regression model for $p_i$ as above. The log-likelihood for this model is $$ \ell(\beta; Z) = \sum_{i=1}^n \sum_{j=1}^{N_i} Z_{ij}\log(p_i) + (1-Z_{ij})\log(1-p_i) $$ and because of the way we defined our $Z_{ij}$s, this can be simplified to $$ \ell(\beta; Y) = \sum_{i=1}^n Y_i \log(p_i) + (N_i - Y_i)\log(1-p_i) $$ which should look pretty familiar.

To get the estimates in the second model, we maximise this with respect to $\beta$. The only difference between this and the first log-likelihood is the term $\log {N_i \choose Y_i}$, which is constant with respect to $\beta$, and so doesn't affect the maximisation and we'll get the same estimates.

3) Each observation has a deviance residual. In the binomial model, they are $$ D_i = 2\left[Y_i \log \left( \frac{Y_i/N_i}{\hat{p}_i} \right) + (N_i-Y_i) \log \left( \frac{1-Y_i/N_i}{1-\hat{p}_i} \right)\right] $$ where $\hat{p}_i$ is the estimated probability from your model. Note that your binomial model is saturated (0 residual degrees of freedom) and has perfect fit: $\hat{p}_i = Y_i/N_i$ for all observations, so $D_i = 0$ for all $i$.

In the Bernoulli model, $$ D_{ij} = 2\left[Z_{ij} \log \left( \frac{Z_{ij}}{\hat{p}_i} \right) + (1-Z_{ij}) \log \left(\frac{1-Z_{ij}}{1-\hat{p}_i} \right)\right] $$ Apart from the fact that you will now have $\sum_{i=1}^n N_i$ deviance residuals (instead of $n$ as with the binomial data), these will each be either $$D_{ij} = -2\log(\hat{p}_i)$$ or $$D_{ij} = -2\log(1-\hat{p}_i)$$ depending on whether $Z_{ij} = 1$ or $0$, and are obviously not the same as the above. Even if you sum these over $j$ to get a sum of deviance residuals for each $i$, you don't get the same: $$ D_i = \sum_{j=1}^{N_i} D_{ij} = 2\left[Y_i \log \left( \frac{1}{\hat{p}_i} \right) + (N_i-Y_i) \log \left( \frac{1}{1-\hat{p}_i} \right)\right] $$

The fact that the AIC is different (but the change in deviance is not) comes back to the constant term that was the difference between the log-likelihoods of the two models. When calculating the deviance, this is cancelled out because it is the same in all models based on the same data. The AIC is defined as $$AIC = 2K - 2\ell$$ and that combinatorial term is the difference between the $\ell$s:

$$AIC_{\mathrm{Bernoulli}} - AIC_{\mathrm{Binomial}} = 2\sum_{i=1}^n \log {N_i \choose Y_i} = 9.575$$

Answer 2

I just want make comments on the last paragraph, “The fact that the AIC is different (but the change in deviance is not) comes back to the constant term that was the difference between the log-likelihoods of the two models. When calculating the change in deviance, this is cancelled out because it is the same in all models based on the same data." Unfortunately, this is not correct for the change in deviance. The deviance does not include the constant term Ex (extra constant term in the log-likelihood for the binomial data). Therefore, the change in deviance does nothing to do with the constant term EX. The deviance compares a given model to the full model. The fact that the deviances are different from Bernoulli/binary and binomial modelling but change in deviance is not is due to the difference in the full model log-likelihood values. These values are cancelled out in calculating the deviance changes. Therefore, Bernoulli and binomial logistic regression models yield an identical deviance changes provided the predicted probabilities pij and pi are the same. In fact, that is true for the probit and other link functions.

Let lBm and lBf denote the log-likelihood values from fitting model m and full model f to Bernoulli data. The deviance is then

    DB=2(lBf - lBm)=-2(lBm – lBf).

Although the lBf is zero for the binary data, we have not simplified the DB and kept it as is. The deviance from the binomial modelling with the same covariates is

    Db=2(lbf+Ex – (lbm+Ex))=2(lbf – lbm) = -2(lbm – lbf)

where the lbf+Ex and lbm+Ex are the log-likelihood values by the full and m models fitted to the binomial data. The extra constant term (Ex) is disappeared from the right hand side of the Db. Now look at change in deviances from Model 1 to Model 2. From Bernoulli modelling, we have change in deviance of

    DBC=DB2-DB1=2(lBf – lBm2)-2(lBf – lBm1) =2(lBm1 – lBm2).

Similarly, change in deviance from binomial fitting is

    DbC=DB2-DB1=2(lbf – lbm2)-2(lbf – lbm1) =2(lbm1 – lbm2).

It is immediately follows that the deviance changes are free from the log-likelihood contributions from full models, lBf and lbf. Therefore, we will get the same change in deviance, DBC = DbC, if lBm1 = lbm1 and lBm2 = lbm2. We know that is the case here and that why we are getting the same deviance changes from Bernoulli and binomial modelling. The difference between lbf and lBf leads to the different deviances.

Answer 3

The AIC are different since the response data changes. The important thing is that differences (i.e. delta) in the AIC, within models with the same response, are the same for the binomial versus Bernoulli:

    binom.data <- data.frame(Successes = c(2, 3, 3), Trials = c(3, 4, 
                   5), X1 = c("Yes", "No", "Yes"), 
                   X2 = c(10.7, 11.3, 9.9)) 

    summary(bi_Full <- glm(formula = cbind(Successes, Trials - 
           Successes) ~ X1 + X2, family = binomial, 
           data = binom.data))
    
    Call:
    glm(formula = cbind(Successes, Trials - Successes) ~ X1 + X2, 
        family = binomial, data = binom.data)
    
    Deviance Residuals: 
    [1]  0  0  0
    
    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)  -2.9649    21.6072  -0.137    0.891
    X1Yes        -0.1897     2.5290  -0.075    0.940
    X2            0.3596     1.9094   0.188    0.851
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 2.2846e-01  on 2  degrees of freedom
    Residual deviance: 4.4409e-16  on 0  degrees of freedom
    AIC: 11.473
    
    Number of Fisher Scoring iterations: 4

    summary(bi_X1 <- glm(formula = cbind(Successes, Trials - 
           Successes) ~ X1,family = binomial, data = binom.data))
    
    Call:
    glm(formula = cbind(Successes, Trials - Successes) ~ X1, 
         family = binomial, 
        data = binom.data)
    
    Deviance Residuals: 
          1        2        3  
     0.1503   0.0000  -0.1150  
    
    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)   1.0986     1.1547   0.951    0.341
    X1Yes        -0.5878     1.3663  -0.430    0.667
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 0.22846  on 2  degrees of freedom
    Residual deviance: 0.03581  on 1  degrees of freedom
    AIC: 9.5087
    
    Number of Fisher Scoring iterations: 4
        
    AIC(bi_Full) - AIC(bi_X1)
    [1] 1.96419
        
    bern.data <- data.frame(Successes = c("Yes", "Yes", "No", "Yes", 
          "Yes", "Yes", "No", "Yes", "Yes", "Yes", "No", "No"), X1 = 
          c("Yes", "Yes", "Yes", "No", "No", "No", "No", "Yes", 
            "Yes", "Yes", "Yes", "Yes"), X2 = c(10.7, 10.7, 10.7, 
            11.3, 11.3, 11.3, 11.3, 9.9, 9.9, 9.9, 9.9, 9.9))
        
    summary(bern_Full <- glm(formula = bern.data$Success == 'Yes' ~ 
            X1 + X2, family = binomial, data = bern.data))
    
    Call:
    glm(formula = bern.data$Success == "Yes" ~ X1 + X2, 
        family = binomial, 
        data = bern.data)
    
    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -1.6651  -1.3537   0.7585   0.9281   1.0108  
    
    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)  -2.9649    21.6072  -0.137    0.891
    X1Yes        -0.1897     2.5290  -0.075    0.940
    X2            0.3596     1.9094   0.188    0.851
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 15.276  on 11  degrees of freedom
    Residual deviance: 15.048  on  9  degrees of freedom
    AIC: 21.048
    
    Number of Fisher Scoring iterations: 4
    
    summary(bern_X1 <- glm(formula = bern.data$Success == 'Yes' ~ X1, 
            family = binomial, data = bern.data))
    
    Call:
    glm(formula = bern.data$Success == "Yes" ~ X1, family = binomial, 
        data = bern.data)
    
    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -1.6651  -1.4006   0.7585   0.9695   0.9695  
    
    Coefficients:
                Estimate Std. Error z value Pr(>|z|)
    (Intercept)   1.0986     1.1547   0.951    0.341
    X1Yes        -0.5878     1.3663  -0.430    0.667
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 15.276  on 11  degrees of freedom
    Residual deviance: 15.084  on 10  degrees of freedom
    AIC: 19.084
    
    Number of Fisher Scoring iterations: 4
    
    AIC(bern_Full) - AIC(bern_X1)
    [1] 1.96419