glm 和 LogitModelFit 之间的区别

Question

我对 R 中的 glm 函数有疑问。

具体来说，我不确定如何包含标称变量。

运行使用 glm 函数后，我在 R 中得到的结果如下：

> df

   x1 x2 y
1  a  2  0
2  b  4  1
3  a  4  0
4  b  2  1
5  a  4  1
6  b  2  0

> str(df)
'data.frame':   6 obs. of  3 variables:
 $ x1: Factor w/ 2 levels "a","b": 1 2 1 2 1 2
 $ x2: num  2 4 4 2 4 2
 $ y: Factor w/ 2 levels "0","1": 1 2 1 2 2 1

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

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -39.132  15208.471  -0.003    0.998
x1b            19.566   7604.236   0.003    0.998
x2              9.783   3802.118   0.003    0.998

但是，当我运行 LogitModelFit 函数 Wolfram Mathematica 时，我得到了不同的参数。

Wolfram 中的代码如下：

data = {{a, 2, 0}, {b, 4, 1}, {a, 4, 0}, {b, 2, 1}, {a, 4, 1}, {b, 2, 0}};

model = LogitModelFit[data, {x, y}, {x, y}, NominalVariables -> x]

model["BestFitParameters"]

这些是我估计的参数：

{-18.5661, -18.5661, 9.28303}

model // Normal

1/(1 + E^(18.5661 - 9.28303 y + 18.5661 DiscreteIndicator[x, a, {a, b}]))

那么，这里有什么不同？为什么结果差异如此之大？

我在 R 或 Wolfram 中做错了什么吗？

Answer 1

好像在你的 LogitModelFit

1/(1 + E^(18.5661 - 9.28303 y + 18.5661 DiscreteIndicator[x, a, {a, b}]))

DiscreteIndicator指离散变量匹配条件x1 == 'a',

而在您的 glm 拟合结果中有一个离散变量 x1b 匹配条件 x1 == 'b':

> str(df)
'data.frame':   6 obs. of  3 variables:
 $ x1: Factor w/ 2 levels "a","b": 1 2 1 2 1 2
 $ x2: num  2 4 4 2 4 2
 $ y: Factor w/ 2 levels "0","1": 1 2 1 2 2 1

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

Coefficients:
             Estimate Std. Error z value Pr(>|z|)
(Intercept)   -39.132  15208.471  -0.003    0.998
x1b            19.566   7604.236   0.003    0.998
x2              9.783   3802.118   0.003    0.998

所以差异似乎是由于 LogitModelFit 和 glm 排除一个依赖类别的方式不同。 LogitModelFit 排除从属类别 x=='a' 而 glm 排除它的补语 x=='b'.

Answer 2

您实际上有 4 个组，您正试图为其估算 3 个参数：

library(dplyr)
df %>% group_by(x1, x2) %>% summarise(n = n(), y = mean(y))

从巨大的标准误差可以看出，参数估计值不稳定。 wolfram 的标准误差也应该非常大（如果给定的话）。

其次，wolfram，似乎使用了不同的参考组，对于x1:

> df$x1 <- relevel(df$x1, "b")
> m <- glm(y ~ x1 + x2, family = binomial(), data = df, control = list(maxit = 100))
> summary(m)

Call:
glm(formula = y ~ x1 + x2, family = binomial(), data = df, control = list(maxit = 100))

Deviance Residuals: 
       1         2         3         4         5         6  
-0.00008   0.00008  -1.17741   1.17741   1.17741  -1.17741  

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  -19.566   7604.236  -0.003    0.998
x1a          -19.566   7604.236  -0.003    0.998
x2             9.783   3802.118   0.003    0.998

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 8.3178  on 5  degrees of freedom
Residual deviance: 5.5452  on 3  degrees of freedom
AIC: 11.545

Number of Fisher Scoring iterations: 18

这更接近于wolfram的结果（这其实和你找到的是同一个模型，我只是选择了另一个参考组）。

两种模型（glm 和 wolfram）的预测实际上是相等的。实际上，前两个参数非常小（最好的模型是 -Inf）且第三个参数等于前两个参数的一半 (9.783*2 = 19.566) 的任何模型都会给出几乎相同的结果。

因子二源于 x2 取值 2 和 4，两者相差 2。

glm 和 LogitModelFit 之间的区别

Difference between glm and LogitModelFit

parameters

wolfram-mathematica

r

logistic-regression

coefficients