R中的二项式GLM:相同的数据,但两个不同的模型

机器算法验证 r 物流
2022-03-11 00:57:11

考虑对这些数据进行逻辑回归:

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 接受三种不同的数据表示:每个表条目一行,以及两种压缩表示(一种带有权重,一种带有成功和失败)。在我看来,这三个规范在数学上应该都是相同的:数据是相同的 7 个观察值,它们以不同的格式呈现给 R。

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

模型 2 和模型 3 是相同的,这是有道理的。但是模型 1 与模型 2 和 3 不同,我无法解释为什么相同的数据应该返回与其他数据不同的模型统计数据(系数、空值和残差)。模型 2 和 3 只是使用相同数据的不同表示。

这可能是一个红鲱鱼,但与模型 2 相比,模型 1 的系数移动了 4 个单位,这正是两者之间(填充)行数/剩余自由度的差异。

> 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
2个回答

模型是

EY=11+exp[(β0+β1x1+β2x2)]

& 它是饱和的,有尽可能多的参数来估计。不同的协变量模式。所以要求解的方程如下:

为了x1=1,x2=0,EY=12

β0+β1=0

为了x1=0,x2=1,EY=14

β0+β2=log3

为了x1=1,x2=1,EY=1

β0+β1+β2=

存在准完全分离(如果x1+x2>1然后EY=1),因此系数的最大似然估计是无界的。但是任何足够大的值c可以代表无穷大,给出解决方案:

β0=(c+log3)

β1=c+log3

β2=c

我不知道为什么glm放弃尝试最大化不同值的可能性c取决于数据结构,但它没有实际影响:预测和可能性的差异几乎相同。

尽管此示例中说明了收敛失败,但应注意这些应用程序中确实存在一些关键差异。加权 GLM 的观测数等于响应水平数,即使权重是频率权重也是如此。另一方面,如果您根据频率权重复制因子水平,则观测数等于权重之和(适当地)。最终,它们将收敛到相同的东西,但是当您检查一步估计器的属性时会观察到有趣的行为:

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

给出以下(不同的)结果:

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 
> 

因此,要回答 OP 的问题,这些是不可调和的结果(尽管收敛失败)的原因是,Fisher Scoring 的实际轨迹对于加权和未加权分析是不同的,因为在加权情况下,Fisher 信息基于 3 个观察加权样本,在未加权的情况下,Fisher 信息是基于 60 个观察未加权信息。3 个观察加权和 60 个观察未加权似然仅在 Fisher 评分实际获得给出 0 总分解决方案的 beta 估计时才一致。