Binomial GLM in R: the same data, but two different models

Question

Consider a logistic regression on these data:

X1 X2 Y
1  0  0
1  0  1
0  1  0
0  1  0
0  1  0
0  1  1
1  1  1

R accepts three different representations of the data: one row per table entry, and two condensed representations (one with weights, one with successes and failures). In my mind, these three specifications should all be the same mathematically: the data is the same 7 observations, and they are presented to R in different formats.

data1 <- data.frame(x1=c(1,1,0,0,0,0,1), x2=c(0,0,1,1,1,1,1), y=c(0,1,0,0,0,1,1))
data2 <- data.frame(x1=c(0,1,0,1), x2=c(0,0,1,1), y=c(0,0.5,0.25,1), w=c(0,2,4,1))
data3x <- data.frame(x1=c(0,1,0,1), x2=c(0,0,1,1))
data3y <- cbind(c(0,1,1,1), c(0,1,3,0))

model1 <- glm(y~x1+x2, data=data1, family="binomial")
model2 <- glm(y~x1+x2, data=data2, family="binomial", weight=w)
model3 <- glm(data3y~data3x$x1+data3x$x2, family="binomial")

Models 2 and 3 are the same, which makes sense. But Model 1 is different from model 2 and 3 and I can't reason out why the same data should return different model statistics (coefficients, null and residual deviance) than the others. Models 2 and 3 just use a different representation of the same data.

This might be a red herring, but Model 1 has its coefficients shifted by 4 units compared to Model 2, which is exactly the difference in the number of (populated) rows/residual degrees of freedom between the two.

> model1

Call:  glm(formula = y ~ x1 + x2, family = "binomial", data = data1)

Coefficients:
(Intercept)           x1           x2  
     -19.66        19.66        18.57  

Degrees of Freedom: 6 Total (i.e. Null);  4 Residual
Null Deviance:      9.561 
Residual Deviance: 7.271    AIC: 13.27
> model2

Call:  glm(formula = y ~ x1 + x2, family = "binomial", data = data2, 
    weights = w)

Coefficients:
(Intercept)           x1           x2  
     -23.66        23.66        22.57  

Degrees of Freedom: 2 Total (i.e. Null);  0 Residual
Null Deviance:      2.289 
Residual Deviance: 3.167e-10    AIC: 9.112

Answer 1

The model is

$$\operatorname{E}Y = \frac{1}{1+ \exp[-(\beta_0 + \beta_1 x_1 + \beta_2 x_2)]}$$

& it's saturated, having as many parameters to estimate as the no. distinct covariate patterns. So the equations to solve are as follows:

For $x_1=1$, $x_2 =0$, $\operatorname{E}Y=\frac{1}{2}$

$\beta_0 + \beta_1 = 0$

For $x_1=0$, $x_2 =1$, $\operatorname{E}Y=\frac{1}{4}$

$\beta_0 + \beta_2 = -\log 3$

For $x_1=1$, $x_2=1$, $\operatorname{E}Y=1$

$\beta_0 + \beta_1 + \beta_2 = \infty$

There is quasi-complete separation (if $x_1+x_2>1$ then $E{Y}=1$), so maximum-likelihood estimates of the coefficients are unbounded. But any sufficiently large value $c$ can stand in for infinity, giving the solutions:

$\beta_0 = -(c + \log 3)$

$\beta_1 = c + \log3$

$\beta_2 = c$

I don't know why glm gives up trying to maximize the likelihood at different values for $c$ depending on the data structure, but it's of no practical consequence: predictions & differences in likelihoods will be almost the same.

Answer 2

Despite the converge fail that was illustrated in this example, it should be noted there are indeed some key differences in these applications. A weighted GLM has number of observations equal to the number of response levels, even when the weights are frequency weights. On the other hand, if you replicate the factor levels according to the frequency weights, the number of observations is equal to the sum of the weights (appropriately). Ultimately, they will converge to the same thing, but interesting behavior is observed when you inspect the properties of the one-step estimators:

set.seed(123)
x <- 0:2
y <- c(1,0,2)/2
w <- 1:3*10

## weighted and unweighted one step glms
summary(glm(y ~ x, family=binomial, weights=w, control=list(maxit = 1)))
summary(glm(y ~ x, family=binomial, data.frame('y'=rep.int(y, w), 'x'=rep.int(x,w)), control=list(maxit = 1)))

Give the following (different) results:

Call:
glm(formula = y ~ x, family = binomial, weights = w, control = list(maxit = 1))

Deviance Residuals: 
      1        2        3  
 0.8269  -7.0855   2.3210  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)  
(Intercept)  -0.5260     0.6210  -0.847   0.3970  
x             1.4456     0.7484   1.932   0.0534 .
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67.640  on 2  degrees of freedom
Residual deviance: 56.275  on 1  degrees of freedom
AIC: 63.079

Number of Fisher Scoring iterations: 1

Warning message:
glm.fit: algorithm did not converge 
> 

Call:
glm(formula = y ~ x, family = binomial, data = data.frame(y = rep.int(y, 
    w), x = rep.int(x, w)), control = list(maxit = 1))

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-1.1496  -1.1496   0.5946   0.5946   0.8376  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)    
(Intercept)  -1.7747     0.5502  -3.226  0.00126 ** 
x             1.7089     0.3700   4.618 3.87e-06 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 67.640  on 59  degrees of freedom
Residual deviance: 44.055  on 58  degrees of freedom
AIC: 44.171

Number of Fisher Scoring iterations: 1

Warning messages:
1: In eval(expr, envir, enclos) :
  non-integer #successes in a binomial glm!
2: glm.fit: algorithm did not converge 
> 

So to answer OP's question, the reason why these are irreconcilable results (despite the converge fail) is that the actual trace of the Fisher Scoring differs for weighted and unweighted analyses because in the weighted case, the Fisher information is based on the 3 observation weighted sample, in the unweighted case, the Fisher Information is based on the 60 observation unweighted information. The 3 observation weighted and 60 observation unweighted likelihoods only agree when Fisher scoring actually obtains a beta estimate giving a 0 sum-score solution.