Binomial GLM in R: dieselben Daten, aber zwei verschiedene Modelle

Betrachten Sie eine logistische Regression dieser Daten:

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 akzeptiert drei verschiedene Darstellungen der Daten: eine Zeile pro Tabelleneintrag und zwei komprimierte Darstellungen (eine mit Gewichten, eine mit Erfolgen und Misserfolgen). Meiner Meinung nach sollten diese drei Spezifikationen mathematisch alle gleich sein: Die Daten sind die gleichen 7 Beobachtungen, und sie werden R in verschiedenen Formaten präsentiert.

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

Die Modelle 2 und 3 sind gleich, was Sinn macht. Aber Modell 1 unterscheidet sich von Modell 2 und 3, und ich kann nicht herausfinden, warum dieselben Daten andere Modellstatistiken (Koeffizienten, Null und Restabweichung) als die anderen zurückgeben sollten. Die Modelle 2 und 3 verwenden lediglich eine andere Darstellung derselben Daten.

Dies mag ein roter Hering sein, aber bei Modell 1 sind die Koeffizienten im Vergleich zu Modell 2 um 4 Einheiten verschoben. Dies ist genau der Unterschied in der Anzahl der (besiedelten) Reihen / verbleibenden Freiheitsgrade zwischen den beiden.

> 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

Das Modell ist

& es ist gesättigt und hat so viele Parameter zu schätzen wie die Nr. unterschiedliche kovariate Muster. Die zu lösenden Gleichungen lauten also wie folgt:

Für ist ,$x_1=1$$x_2 =0$$\operatorname{E}Y=\frac{1}{2}$

$\beta_0 + \beta_1 = 0$

Für , ,$x_1=0$$x_2 =1$$\operatorname{E}Y=\frac{1}{4}$

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

Für ist ,$x_1=1$$x_2=1$$\operatorname{E}Y=1$

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

Es gibt eine quasi vollständige Trennung (wenn dann ist ), so dass Schätzungen der Koeffizienten mit maximaler Wahrscheinlichkeit unbegrenzt sind. Aber jeder ausreichend große Wert kann für unendlich stehen und die Lösungen ergeben:$x_1+x_2>1$$E{Y}=1$$c$

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

$\beta_1 = c + \log3$

$\beta_2 = c$

Ich weiß nicht, warum glmes aufgegeben wird, die Wahrscheinlichkeit bei unterschiedlichen Werten für Abhängigkeit von der Datenstruktur zu maximieren , aber es hat keine praktische Konsequenz: Vorhersagen und Unterschiede in den Wahrscheinlichkeiten werden fast gleich sein.$c$

Trotz des in diesem Beispiel dargestellten Konvergenzfehlers sollte beachtet werden, dass es in diesen Anwendungen tatsächlich einige wesentliche Unterschiede gibt. Ein gewichtetes GLM hat eine Anzahl von Beobachtungen, die der Anzahl von Antwortniveaus entspricht, selbst wenn die Gewichte Frequenzgewichte sind. Wenn Sie dagegen die Faktorstufen gemäß den Frequenzgewichten replizieren, entspricht die Anzahl der Beobachtungen (entsprechend) der Summe der Gewichte. Letztendlich konvergieren sie zur gleichen Sache, aber ein interessantes Verhalten wird beobachtet, wenn Sie die Eigenschaften der Ein-Schritt-Schätzer untersuchen:

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

Geben Sie die folgenden (unterschiedlichen) Ergebnisse an:

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 
> 

Um die Frage von OP zu beantworten, ist der Grund, warum dies unvereinbare Ergebnisse sind (trotz des Konvergenzfehlers), dass die tatsächliche Spur der Fisher-Bewertung für gewichtete und ungewichtete Analysen unterschiedlich ist, da im gewichteten Fall die Fisher-Informationen auf den 3 gewichteten Beobachtungen basieren In dem ungewichteten Fall basieren die Fisher-Informationen auf den 60 ungewichteten Beobachtungsinformationen. Die 3 beobachtungsgewichteten und 60 beobachtungsungewichteten Wahrscheinlichkeiten stimmen nur überein, wenn die Fisher-Bewertung tatsächlich eine Beta-Schätzung erhält, die eine 0-Summen-Score-Lösung ergibt.