贝叶斯预测，下标越界答案

【问题标题】：Bayes predict, subscript out of bounds贝叶斯预测，下标越界
【发布时间】：2014-06-16 15:57:50
【问题描述】：

在使用bayesglm 时，我遇到了predict 函数的一些问题。我读过一些帖子，说当样本外数据的级别多于样本数据时，可能会出现此问题，但我将相同的数据用于拟合和预测函数。预测适用于常规 glm，但不适用于bayesglm。示例：

control <- y ~ x1 + x2

# this works fine:
glmObject <- glm(control, myData, family = binomial())
predicted1 <- predict.glm(glmObject , myData, type = "response")

# this gives an error: 
bayesglmObject <- bayesglm(control, myData, family = binomial())
predicted2 <- predict.bayesglm(bayesglmObject , myData, type = "response") 
Error in X[, piv, drop = FALSE] : subscript out of bounds

# Edit... I just discovered this works. 
# Should I be concerned about using these results?
# Not sure why is fails when I specify the dataset
predicted3 <- predict(bayesglmObject, type = "response")

无法弄清楚如何使用bayesglm 对象进行预测。有任何想法吗？谢谢！

【问题讨论】：

标签： r prediction glm predict bayesglm

【解决方案1】：

其中一个原因可能与bayesglm 命令中参数“drop.unused.levels”的默认设置有关。默认情况下，此参数设置为 TRUE。因此，如果有未使用的关卡，它会在模型构建过程中被丢弃。但是，预测函数仍然使用因子变量中存在未使用水平的原始数据。这会导致用于模型构建的数据和用于预测的数据之间存在级别差异（即使它是相同的数据名 - 在您的情况下，myData）。我在下面给出了一个例子：

    n <- 100
    x1 <- rnorm (n)
    x2 <- as.factor(sample(c(1,2,3),n,replace = TRUE))

    # Replacing 3 with 2 makes the level = 3 as unused
    x2[x2==3] <- 2

    y <- as.factor(sample(c(1,2),n,replace = TRUE))

    myData <- data.frame(x1 = x1, x2 = x2, y = y)
    control <- y ~ x1 + x2

    # this works fine:
    glmObject <- glm(control, myData, family = binomial())
    predicted1 <- predict.glm(glmObject , myData, type = "response")

    # this gives an error - this uses default drop.unused.levels = TRUE
    bayesglmObject <- bayesglm(control, myData, family = binomial())
    predicted2 <- predict.bayesglm(bayesglmObject , myData, type = "response") 

    Error in X[, piv, drop = FALSE] : subscript out of bounds

    # this works fine - value of drop.unused.levels is set to FALSE
    bayesglmObject <- bayesglm(control, myData, family = binomial(),drop.unused.levels   = FALSE)
    predicted2 <- predict.bayesglm(bayesglmObject , myData, type = "response")

我认为更好的方法是使用 droplevels 预先从数据框中删除未使用的级别，并将其用于模型构建和预测。

【讨论】：