GLM binomial em R: os mesmos dados, mas dois modelos diferentes

Considere uma regressão logística nesses dados:

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 aceita três representações diferentes dos dados: uma linha por entrada da tabela e duas representações condensadas (uma com pesos, outra com sucessos e falhas). Na minha opinião, essas três especificações devem ser todas matematicamente iguais: os dados são as mesmas 7 observações e são apresentadas ao R em diferentes formatos.

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")

Os modelos 2 e 3 são os mesmos, o que faz sentido. Mas o Modelo 1 é diferente dos modelos 2 e 3 e não consigo entender por que os mesmos dados devem retornar estatísticas diferentes do modelo (coeficientes, desvio nulo e residual) que os outros. Os modelos 2 e 3 usam apenas uma representação diferente dos mesmos dados.

Pode ser um arenque vermelho, mas o Modelo 1 tem seus coeficientes alterados em 4 unidades em comparação com o Modelo 2, que é exatamente a diferença no número de linhas (preenchidas) / graus residuais de liberdade entre as duas.

> 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

O modelo é

e está saturado, com tantos parâmetros para estimar quanto não. padrões covariáveis ​​distintos. Portanto, as equações a serem resolvidas são as seguintes:

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

$\beta_0 + \beta_1 = 0$

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

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

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

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

Existe uma separação quase completa (se $x_1+x_2>1$ então $E{Y}=1$), para que as estimativas de probabilidade máxima dos coeficientes sejam ilimitadas. Mas qualquer valor suficientemente grande$c$ pode substituir o infinito, fornecendo as soluções:

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

$\beta_1 = c + \log3$

$\beta_2 = c$

Não sei por que glmdesiste de tentar maximizar a probabilidade de diferentes valores para$c$ dependendo da estrutura dos dados, mas não tem conseqüências práticas: previsões e diferenças de probabilidade serão quase as mesmas.

Apesar da falha de convergência ilustrada neste exemplo, deve-se notar que existem algumas diferenças importantes nessas aplicações. Um GLM ponderado possui um número de observações igual ao número de níveis de resposta, mesmo quando os pesos são pesos de frequência. Por outro lado, se você replicar os níveis de fator de acordo com os pesos de frequência, o número de observações será igual à soma dos pesos (apropriadamente). Por fim, eles convergirão para a mesma coisa, mas um comportamento interessante é observado quando você inspeciona as propriedades dos estimadores de uma etapa:

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)))

Dê os seguintes resultados (diferentes):

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 
> 

Portanto, para responder à pergunta do OP, a razão pela qual esses resultados são inconciliáveis ​​(apesar da falha de convergência) é que o traço real da pontuação de Fisher difere nas análises ponderadas e não ponderadas, porque, no caso ponderado, as informações de Fisher são baseadas nas três observações ponderadas amostra, no caso não ponderado, a Informação Fisher é baseada nas 60 informações não ponderadas da observação. As três probabilidades de observação ponderada e 60 de observação não ponderada só concordam quando a pontuação de Fisher realmente obtém uma estimativa beta, fornecendo uma solução de soma-0.