机器算法验证 - 在 R（nnet 包）中获取“multinom”的 p 值 - 吾爱随笔录

在 R（nnet 包）中获取“multinom”的 p 值

机器算法验证 r 回归 p 值多项分布

2022-01-20 07:49:43

如何使用package in的multinom函数获取 p 值？nnetR

我有一个数据集，其中包含“病理评分”（缺席、轻度、重度）作为结果变量，以及两个主要影响：年龄（两个因素：二十 / 三十天）和治疗组（四个因素：没有 ATB 感染；感染 + ATB1；感染 + ATB2；感染 + ATB3)。

首先，我尝试拟合一个序数回归模型，考虑到我的因变量（序数）的特征，这似乎更合适。然而，赔率比例的假设被严重违反（图形上），这促使我改用多项模型，使用nnet包。

首先，我选择了需要用作基线类别的结果级别：

Data$Path <- relevel(Data$Path, ref = "Absent")

然后，我需要为自变量设置基线类别：

Data$Age <- relevel(Data$Age, ref = "Twenty")
Data$Treat <- relevel(Data$Treat, ref="infected without ATB")

该模型：

test <- multinom(Path ~ Treat + Age, data = Data) 
# weights:  18 (10 variable) 
initial value 128.537638 
iter 10 value 80.623608 
final  value 80.619911 
converged

输出：

Coefficients:
         (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate   -2.238106   -1.1738540      -1.709608       -1.599301        2.684677
Severe     -1.544361   -0.8696531      -2.991307       -1.506709        1.810771

Std. Errors:
         (Intercept)    infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate   0.7880046    0.8430368       0.7731359       0.7718480        0.8150993
Severe     0.6110903    0.7574311       1.1486203       0.7504781        0.6607360

Residual Deviance: 161.2398
AIC: 181.2398

有一段时间，我找不到一种方法来获取模型的值和使用. 昨天我遇到了一篇文章，其中作者提出了一个关于估计系数值的类似问题（如何在 R 中建立和估计多项式 logit 模型？）。在那里，一位博主建议从结果中值非常容易，首先获取值如下： $p$ nnet:multinom $p$ $p$ summarymultinom $t$

pt(abs(summary1$coefficients / summary1$standard.errors), df=nrow(Data)-10, lower=FALSE) 

         (Intercept)   infected+ATB1   infected+ATB2   infected+ATB3    AgeThirty
Moderate 0.002670340   0.08325396      0.014506395     0.02025858       0.0006587898
Severe   0.006433581   0.12665278      0.005216581     0.02352202       0.0035612114

值，至少有 2 个因子缺失分布用于真正的统计量通常是错误的；对于聚合数据，它可能是一个非常严重的错误。” 根据 Brian Ripley 的说法，“使用 Wald 检验进行拟合也是一个错误，因为它们与二项式拟合存在相同（可能很严重）的问题。使用轮廓似然置信区间（软件包确实为此提供了软件），或者如果你必须测试，似然比测试（同上）。” $p$ $t$ $z$ multinom

我只需要能够得出可靠的值。 $p$

4个回答

怎么用

z <- summary(test)$coefficients/summary(test)$standard.errors
# 2-tailed Wald z tests to test significance of coefficients
p <- (1 - pnorm(abs(z), 0, 1)) * 2
p

基本上，这将基于相对于其标准误差的估计系数，并且将使用 az test 来测试基于双尾检验的零显着差异。因子二纠正了上面提到的问题 Peter Dalgaard（你需要它是因为你想要一个双尾检验，而不是一个尾检验），它使用 z 检验而不是 t 检验来解决另一个问题你提到的问题。

您也可以使用获得相同的结果（Wald z-tests）

library(AER)
coeftest(test)

似然比检验通常被认为比 Wald z 检验更准确（后者使用正态近似，LR 检验不使用），这些可以使用

library(afex)
set_sum_contrasts() # use sum coding, necessary to make type III LR tests valid
library(car)
Anova(test,type="III")

如果您想进行成对的 Tukey 事后测试，那么可以使用我在另一篇文章中解释lsmeans的包来获得这些测试！

您也可能对似然比检验 p 值感兴趣，如下所示：

http://thestatsgeek.com/2014/02/08/wald-vs-likelihood-ratio-test/

你可以像这样提取（对不起，它是一个自定义函数：D）

likehoodmultinom_p <- function(model_lmm) 
{

  i <- 1

  variables <-c("No funciona")
  valores <- c("No funciona") 


  for (var in model_lmm$coefnames[-1]) { # Qutiamos el -1 de coefnames para no obener un NA

  variables[i] =paste(var)
  valores[i]= lrtest(model_lmm, var)[[5]][2]
  i=i+1
   ## Contributed to stack at:  https://stackoverflow.com/questions/23018238/assesing-the-goodness-of-fit-for-the-multinomial-logit-in-r-with-the-nnet-packag/60835647#60835647
  }
  return (data.frame(variables,valores))
}

L_iris= likehoodmultinom_p(iris_fit)

在我的函数中，您获得了一个带有因子的 df，因此您可能需要稍微更改它们以提取它们。我还没有纠正我原来的功能：

L_iris= likehoodmultinom_p(iris_fit)

L_iris$valores = as.character(L_iris$valores) # Pass them as chr
L_iris$valores = as.numeric(L_iris$valores) # And as numeric.

然后，您可以轻松访问它们。我通常也根据 p 值对它们进行排序。

正如 OP（他引用 B Ripley 的话）已经说过的那样，wald 测试对于多项模型来说并不是很好，我们真的应该使用似然比测试。下面我展示了一种通过包中的函数获取它的简单方法MASS，使用nnet::multinom. 使用的主力功能是MASS::dropterm：

> library(nnet)
> example(birthwt)
(bwt.mn <- multinom(low  ~ . , bwt) )


brthwt> bwt <- with(birthwt, {
brthwt+ race <- factor(race, labels = c("white", "black", "other"))
brthwt+ ptd <- factor(ptl > 0)
.
.
.
Call:
multinom(formula = low ~ ., data = bwt)

Coefficients:
(Intercept)         age         lwt   raceblack   raceother   smokeTRUE 
 0.82320102 -0.03723828 -0.01565359  1.19240391  0.74065606  0.75550487 
    ptdTRUE      htTRUE      uiTRUE        ftv1       ftv2+ 
 1.34375901  1.91320116  0.68020207 -0.43638470  0.17900392 

Residual Deviance: 195.4755 
AIC: 217.4755 
> confint(bwt.mn)   
                  2.5 %       97.5 %
(Intercept) -1.61649875  3.262900795
age         -0.11309745  0.038620883
lwt         -0.02953168 -0.001775495
raceblack    0.14190092  2.242906899
raceother   -0.16438896  1.645701076
smokeTRUE   -0.07755089  1.588560638
ptdTRUE      0.40173272  2.285785295
htTRUE       0.50053490  3.325867418
uiTRUE      -0.22990670  1.590310835
ftv1        -1.37601313  0.503243725
ftv2+       -0.71550657  1.073514417
> MASS::dropterm(bwt.mn, trace=FALSE, test="Chisq") 
# weights:  11 (10 variable)
initial  value 131.004817 
iter  10 value 98.297550
.
.
.
Single term deletions

Model:
low ~ age + lwt + race + smoke + ptd + ht + ui + ftv
       Df    AIC    LRT  Pr(Chi)   
<none>    217.48                   
age     1 216.42 0.9419 0.331796   
lwt     1 220.95 5.4739 0.019302 * 
race    2 219.23 5.7513 0.056380 . 
smoke   1 218.67 3.1982 0.073717 . 
ptd     1 223.58 8.1085 0.004406 **
ht      1 222.93 7.4584 0.006314 **
ui      1 217.59 2.1100 0.146342   
ftv     2 214.83 1.3582 0.507077   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

还有许多其他方法……但是将函数 fromMASS与函数 from 一起使用nnet似乎是谨慎的。

考虑 p 值的一种方法是进行可能性测试，即您的拟合比一些更简单的更少项（或可能没有项，恒定拟合）的拟合更好。下面是一些代码。

# Multinomial fit
fit <- nnet::multinom(cyl ~ mpg + hp, data=datasets::mtcars)

# Multinomial fit with one or more terms dropped
base_fit <- nnet::multinom(cyl ~ 1, data=datasets::mtcars)
base_fit2 <- nnet::multinom(cyl ~ mpg, data=datasets::mtcars)

# p-value that the fit is better than the base_fit
result <- lmtest::lrtest(fit, base_fit)
p_val1 <- result$`Pr(>Chisq)`[[2]]

# p-value that the fit is better than the base_fit2
result <- lmtest::lrtest(fit, base_fit2)
p_val2 <- result$`Pr(>Chisq)`[[2]]

结果：

> p_val1
[1] 6.250054e-14
> p_val2
[1] 0.0003148036

显然，下降hp会使拟合更差（p = 0.0003），而同时mpg下降会使拟合更hp差（p 接近于零）。

其它你可能感兴趣的问题

上一篇维数大于样本数时的 PCA 下一篇使用R的时间序列的STL趋势