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

Question

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

因子错误（G，标签 = dx1[，格式（max（P），数字 = 3），by = G]$V1）： 无效的“标签”；长度 6 应该是 1 或 5

我无法弄清楚导致此错误的原因。您可以在下面找到导致错误的代码以及有效的代码。 如果您需要更多信息，请告诉我

H

R 版本和加载的库

R 版本 3.2.3 (2015-12-10) 平台：x86_64-pc-linux-gnu（64位） 运行于：Ubuntu 15.04

附加的基础包： [1] stats graphics grDevices utils datasets methods base

其他附加包： [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)

Answer 1

问题出现在源码第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 (例如 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

Answer 2

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

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

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