使用 r 中 LogisticDx 包中的 gof() 函数时出错

Error using gof() function from LogisticDx package in r

您好,我在尝试使用 LogisticDx 包中的 gof() 函数时不断收到以下错误消息。

因子错误(G, 标签 = dx1[ 格式(max(P), 数字 = 3), by = G]$V1) : 无效'labels';长度 6 应该是 1 或 5

我不知道是什么原因导致了这个错误。您可以在下面找到导致错误的代码以及有效的代码。 如果您需要更多信息,请告诉我

H

R 版本和加载的库

R 版本 3.2.3 (2015-12-10) 平台:x86_64-pc-linux-gnu(64 位) 运行下:Ubuntu15.04

附加基础包: [1] stats graphics grDevices utils 数据集方法基础

其他附包: [1] ROCR_1.0-7 gplots_2.17.0 LogisticDx_0.2 xtable_1.8-2 pander_0.5.2
[6] plyr_1.8.3 Amelia_1.7.4 mice_2.25 Rcpp_0.11.6 knitr_1.11

产生错误的代码

PD <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
        1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
        0, 0, 0)
E <- c(4, 4, 3, 1, 0, 5, 1, 5, 5, 5, 1, 5, 4, 2, 1, 1, 1, 1, 5,
       0, 5, 5, 5, 5, 0, 0, 4, 5, 2, 4, 5, 4, 5, 3, 0, 5, 3, 3,
       5, 2, 4, 0, 0, 5, 1, 0, 3, 2, 5, 1, 2, 4, 0, 2, 5, 5, 4,
       3, 5, 2, 0, 2, 0, 1, 0, 5, 2, 0, 3, 0, 0, 0, 1, 0, 0, 2,
       3, 0, 0, 1, 0, 0, 0, 0, 0, 2, 0, 1, 3, 0, 2, 0, 3, 0, 0,
       0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 2, 2, 0, 0, 2, 0, 0,
       0, 0, 1, 0, 0, 0, 1, 1, 0, 0, 0, 2, 1, 0, 0)

test.data <- data.frame(PD,E)
factor(test.data$PD)

g.error <- gof(glm(PD ~ E,family=binomial,data=test.data),plot=FALSE)

有效的代码

data(uis)
g.works <- gof(glm(RACE ~ NDRGTX,family=binomial,data=uis),plot=FALSE)

这里的问题与正在执行的 Hosmer–Lemeshow 测试中的组数有关。 gof 的默认参数似乎是 g=10,它似乎不适用于 129 个观察值(如果您再添加一行它确实有效)。 如果将组数设置为 9,它应该可以工作。

gof(glm(PD ~ E,family=binomial,data=test.data), plot=FALSE, g=9)

也许其他人可以回答为什么会这样?

问题出现在源码第397行 https://github.com/dardisco/LogisticDx/blob/master/R/gof.R

### Hosmer-Lemeshow
    ## sort by probability
    dx1 <- dx1[order(P), ]
    ## base group size
    gs1 <- dx1[, sum(n) %/% g] # Line 393
    g1 <- rep(gs1, g)
    ## remainer
    mod1 <- dx1[, sum(n) %% g]
    g1[seq(1, g-1, by=g/mod1)] <- gs1 + 1 # Line 397

这里seq(1, g-1, by=g/mod1)的长度不一定等于第393行的除法余数使得sum(g1)不等于n(ex n= 129,g=10)。这会导致第 399 行出现错误

dx1[, "G" := factor(G,
              labels=dx1[, format(max(P), digits=3), by=G]$V1)]

因为当 n=129 和 g=10 时,最后一个 dx1[129,G] 是 NA。解决方案可能是

g1[seq(1, mod1, by = 1)] <- gs1+1