glm() 模型的交叉验证答案

【问题标题】：Cross validation for glm() modelsglm() 模型的交叉验证
【发布时间】：2014-02-18 06:06:51
【问题描述】：

我正在尝试对我之前在 R 中构建的一些 glm 模型进行 10 倍交叉验证。我对 boot 包中的 cv.glm() 函数有点困惑，尽管我已经阅读大量帮助文件。当我提供以下公式时：

library(boot)
cv.glm(data, glmfit, K=10)

这里的“数据”参数是指整个数据集还是仅指测试集？

到目前为止，我看到的示例提供了“数据”参数作为测试集，但这并没有真正的意义，例如为什么在同一个测试集上进行 10 折？它们都会给出完全相同的结果（我假设！）。

不幸的是?cv.glm 解释得模糊不清：

data：包含数据的矩阵或数据框。行应该是案例和列对应于变量，其中之一是回应

我的另一个问题是关于 $delta[1] 结果。这是 10 次试验的平均预测误差吗？如果我想得到每个折叠的错误怎么办？

这是我的脚本的样子：

##data partitioning
sub <- sample(nrow(data), floor(nrow(x) * 0.9))
training <- data[sub, ]
testing <- data[-sub, ]

##model building
model <- glm(formula = groupcol ~ var1 + var2 + var3,
        family = "binomial", data = training)

##cross-validation
cv.glm(testing, model, K=10)

【问题讨论】：

查看boot:::cv.glm的示例部分。您应该输入整个数据，模型和CV的折叠。
感谢您的回复@RomanLuštrik。听起来很棒。不过，我仍然想知道一些事情。此函数是否在交叉验证中使用所有提供的数据？假设我为cv.glm(data, glm, K=10) 提供了一个 1000 行的数据框，它是否对数据进行了 10 个分区，每个分区为 100 个并进行交叉验证？抱歉，我已经通过了 ?cv.glm 但我没有在那里找到。
如果您要进行 2 倍 CV，该函数将获取 50% 的数据并拟合模型。它将使用其他 50% 的数据来查看模型对数据的描述程度。或者，在留一法 CV 中，它会将模型拟合到除一个数据“点”之外的所有数据“点”，并查看单出“点”的效果如何。重复 N 次即可得到结果。
嗨@RomanLuštrik。你说如果我做一个 2-fold CV，函数会根据 50% 的数据拟合模型，并用另外 50% 作为测试集。如果函数这样做，那么为什么它需要一个参数“glmfit”，它是一个先前拟合的模型？
如果你有关于交叉验证的问题，我建议你在 crossvalidated.com 上打开一个线程。

标签： r partitioning prediction glm cross-validation

【解决方案1】：

我对使用各种包的 10 折交叉验证方法总是有点谨慎。我有自己的简单脚本，可以为任何机器学习包手动创建测试和训练分区：

#Randomly shuffle the data
yourData<-yourData[sample(nrow(yourData)),]

#Create 10 equally size folds
folds <- cut(seq(1,nrow(yourData)),breaks=10,labels=FALSE)

#Perform 10 fold cross validation
for(i in 1:10){
    #Segement your data by fold using the which() function 
    testIndexes <- which(folds==i,arr.ind=TRUE)
    testData <- yourData[testIndexes, ]
    trainData <- yourData[-testIndexes, ]
    #Use test and train data partitions however you desire...
}

【讨论】：

谢谢杰克·德鲁。出于测试目的，我将您上面的代码与 cv.glm 的结果进行了比较，结果是相同的。感谢您的帖子，我现在可以信任 cv.glm ;-)

【解决方案2】：

@Roman 在他的 cmets 中提供了一些答案，但是，您的问题的答案是通过使用 cv.glm 检查代码来提供的：

我相信这段代码将随机设置的数据拆分为 K-folds，如果 K 不整除 n，则根据需要安排舍入：

if ((K > n) || (K <= 1)) 
    stop("'K' outside allowable range")
K.o <- K
K <- round(K)
kvals <- unique(round(n/(1L:floor(n/2))))
temp <- abs(kvals - K)
if (!any(temp == 0)) 
    K <- kvals[temp == min(temp)][1L]
if (K != K.o) 
    warning(gettextf("'K' has been set to %f", K), domain = NA)
f <- ceiling(n/K)
s <- sample0(rep(1L:K, f), n)

这里的这一位表明增量值不是均方根误差。就像帮助文件中所说的那样The default is the average squared error function. 这是什么意思？我们可以通过检查函数声明来看到这一点：

function (data, glmfit, cost = function(y, yhat) mean((y - yhat)^2), 
    K = n)

这表明在每个折叠中，我们计算误差平方的平均值，其中误差通常指预测响应与实际响应之间的误差。

delta[1] 只是每个折叠的所有这些术语的总和的weighted average，请参阅cv.glm 代码中的我的内联 cmets：

for (i in seq_len(ms)) {
    j.out <- seq_len(n)[(s == i)]
    j.in <- seq_len(n)[(s != i)]
    Call$data <- data[j.in, , drop = FALSE]
    d.glm <- eval.parent(Call)
    p.alpha <- n.s[i]/n #create weighted average for later
    cost.i <- cost(glm.y[j.out], predict(d.glm, data[j.out, 
        , drop = FALSE], type = "response"))
    CV <- CV + p.alpha * cost.i # add weighted average error to running total
    cost.0 <- cost.0 - p.alpha * cost(glm.y, predict(d.glm, 
        data, type = "response"))
}

【讨论】：