R 中的多项式逻辑回归:nnet 包中的多项式结果与 mlogit 包中的 mlogit 不同?
multinomial logistic regression in R: multinom in nnet package result different from mlogit in mlogit package?
两个 R 函数,multinom
(包 nnet
)和 mlogit
(包 mlogit
)均可用于多项逻辑回归。但是为什么这个例子 returns 系数 p 值的不同结果?
#prepare data
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
mydata$gre[1:10] = rnorm(10,mean=80000)
#multinom
:
test = multinom(admit ~ gre + gpa + rank, data = mydata)
z <- summary(test)$coefficients/summary(test)$standard.errors
# For simplicity, use z-test to approximate t test.
pv <- (1 - pnorm(abs(z)))*2
pv
# (Intercept) gre gpa rank2 rank3 rank4
# 0.00000000 0.04640089 0.00000000 0.00000000 0.00000000 0.00000000
#mlogit
:
mldata = mlogit.data(mydata,choice = 'admit', shape = "wide")
mlogit.model1 <- mlogit(admit ~ 1 | gre + gpa + rank, data = mldata)
summary(mlogit.model1)
# Coefficients :
# Estimate Std. Error t-value Pr(>|t|)
# 1:(intercept) -3.5826e+00 1.1135e+00 -3.2175 0.0012930 **
# 1:gre 1.7353e-05 8.7528e-06 1.9825 0.0474225 *
# 1:gpa 1.0727e+00 3.1371e-01 3.4195 0.0006274 ***
# 1:rank2 -6.7122e-01 3.1574e-01 -2.1258 0.0335180 *
# 1:rank3 -1.4014e+00 3.4435e-01 -4.0697 4.707e-05 ***
# 1:rank4 -1.6066e+00 4.1749e-01 -3.8482 0.0001190 ***
为什么 multinorm
和 mlogit
的 p 值如此不同?我想这是因为我使用 mydata$gre[1:10] = rnorm(10,mean=80000)
添加了异常值。如果异常值是不可避免的问题(例如在基因组学、代谢组学等方面),我应该使用哪个 R 函数?
这里的区别是 Wald $z$ 检验(你在 pv
中计算的)和似然比检验(summary(mlogit.model)
返回的东西)之间的区别。Wald 检验在计算上是更简单,但通常具有不太理想的属性(例如,其 CI 不是 scale-invariant)。您可以阅读有关这两个程序的更多信息 here.
要对 nnet
模型系数执行 LR 测试,您可以加载 car
和 lmtest
包并调用 Anova(test)
(尽管您必须这样做单个 df 测试需要多做一些工作)。
作为替代方案,您可以使用 broom
,它为 multinom
class 模型输出整齐的格式。
library(broom)
tidy(test)
它会 return data.frame
具有 z 统计量和 p 值。
查看 tidy
文档以获取更多信息。
P.S.: 因为我无法从您发布的 link 中获取数据,所以我无法复制结果
两个 R 函数,multinom
(包 nnet
)和 mlogit
(包 mlogit
)均可用于多项逻辑回归。但是为什么这个例子 returns 系数 p 值的不同结果?
#prepare data
mydata <- read.csv("http://www.ats.ucla.edu/stat/data/binary.csv")
mydata$rank <- factor(mydata$rank)
mydata$gre[1:10] = rnorm(10,mean=80000)
#multinom
:
test = multinom(admit ~ gre + gpa + rank, data = mydata)
z <- summary(test)$coefficients/summary(test)$standard.errors
# For simplicity, use z-test to approximate t test.
pv <- (1 - pnorm(abs(z)))*2
pv
# (Intercept) gre gpa rank2 rank3 rank4
# 0.00000000 0.04640089 0.00000000 0.00000000 0.00000000 0.00000000
#mlogit
:
mldata = mlogit.data(mydata,choice = 'admit', shape = "wide")
mlogit.model1 <- mlogit(admit ~ 1 | gre + gpa + rank, data = mldata)
summary(mlogit.model1)
# Coefficients :
# Estimate Std. Error t-value Pr(>|t|)
# 1:(intercept) -3.5826e+00 1.1135e+00 -3.2175 0.0012930 **
# 1:gre 1.7353e-05 8.7528e-06 1.9825 0.0474225 *
# 1:gpa 1.0727e+00 3.1371e-01 3.4195 0.0006274 ***
# 1:rank2 -6.7122e-01 3.1574e-01 -2.1258 0.0335180 *
# 1:rank3 -1.4014e+00 3.4435e-01 -4.0697 4.707e-05 ***
# 1:rank4 -1.6066e+00 4.1749e-01 -3.8482 0.0001190 ***
为什么 multinorm
和 mlogit
的 p 值如此不同?我想这是因为我使用 mydata$gre[1:10] = rnorm(10,mean=80000)
添加了异常值。如果异常值是不可避免的问题(例如在基因组学、代谢组学等方面),我应该使用哪个 R 函数?
这里的区别是 Wald $z$ 检验(你在 pv
中计算的)和似然比检验(summary(mlogit.model)
返回的东西)之间的区别。Wald 检验在计算上是更简单,但通常具有不太理想的属性(例如,其 CI 不是 scale-invariant)。您可以阅读有关这两个程序的更多信息 here.
要对 nnet
模型系数执行 LR 测试,您可以加载 car
和 lmtest
包并调用 Anova(test)
(尽管您必须这样做单个 df 测试需要多做一些工作)。
作为替代方案,您可以使用 broom
,它为 multinom
class 模型输出整齐的格式。
library(broom)
tidy(test)
它会 return data.frame
具有 z 统计量和 p 值。
查看 tidy
文档以获取更多信息。
P.S.: 因为我无法从您发布的 link 中获取数据,所以我无法复制结果