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

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
• Data2 is incorrectly factored. The [1, 0, .5] response level receives a weight of 2 indicating 2 levels with y taking 0 and 1 as an average response. However, there are no [1,0,.5] response levels in the data you show. – AdamO Oct 31 '16 at 16:08
• These models are giving you the same result in the sense that they both exploded :) – AdamO Oct 31 '16 at 20:15
• @AdamO that's fair. I should have thought through these results more before posting. I was just so shocked that the estimates were so very different that I had to ask for help. – Sycorax Oct 31 '16 at 20:25

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.

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.