交叉验证结果的解释 - cv.glm()

机器算法验证 r 回归 物流 交叉验证
2022-04-22 06:34:37

由于系数很大,我的逻辑模型一直很可疑,所以我尝试进行交叉验证,并对简化模型进行交叉验证,以确认原始模型被过度指定的事实,正如James 建议的那样。但是,我不知道如何解释结果(这是来自链接问题的模型):

> summary(m5)

Call:
glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj + resid_usili2 + 
    rok:obdobi + rok:kraj + obdobi:kraj + kraj:resid_usili2 + 
    rok:obdobi:kraj, family = "quasibinomial")
[... see https://stats.stackexchange.com/q/48739/5509 for complete summary output ]

> cv.glm(na.omit(data.frame(orel, resid_usili2)), m5, K = 10)
$call
cv.glm(data = na.omit(data.frame(orel, resid_usili2)), glmfit = m5, 
    K = 10)

$K
[1] 10

$delta
[1] 0.2415355 0.2151626

$seed
  [1]         403         271  1234892862 -1124595763  -489713400  1566924080   147612843
  [8]  1879282918  -694084381  1171051622  2063023839 -1307030905  -477709428  1248673977
 [15]  -746898494   420363755  -890078828   460552896  -758793089  -913500073  -882355605
[....]
Warning message:
glm.fit: algorithm did not converge

我猜 delta 是平均拟合误差,但如何解释呢?它是好还是坏?顺便说一句,算法没有收敛,可能是因为系数很大(?)

我尝试了一个简化模型:

> summary(m)

Call:
glm(formula = cbind(ml, ad) ~ rok + obdobi + kraj, family = "quasibinomial")

Deviance Residuals: 
    Min       1Q   Median       3Q      Max  
-2.7335  -1.2324  -0.1666   1.0866   3.1788  

Coefficients:
              Estimate Std. Error t value Pr(>|t|)    
(Intercept) -107.60761   48.06535  -2.239 0.025335 *  
rok            0.05381    0.02393   2.249 0.024683 *  
obdobinehn    -0.26962    0.10372  -2.599 0.009441 ** 
krajJHC        0.68869    0.27617   2.494 0.012761 *  
krajJHM       -0.26607    0.28647  -0.929 0.353169    
krajLBK       -1.11305    0.55165  -2.018 0.043828 *  
krajMSK       -0.61390    0.37252  -1.648 0.099593 .  
krajOLK       -0.49704    0.32935  -1.509 0.131501    
krajPAK       -1.18444    0.35090  -3.375 0.000758 ***
krajPLK       -1.28668    0.44238  -2.909 0.003691 ** 
krajSTC        0.01872    0.27806   0.067 0.946322    
krajULKV      -0.41950    0.61647  -0.680 0.496315    
krajVYS       -1.17290    0.39733  -2.952 0.003213 ** 
krajZLK       -0.38170    0.36487  -1.046 0.295698    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for quasibinomial family taken to be 1.304775)

    Null deviance: 2396.8  on 1343  degrees of freedom
Residual deviance: 2198.6  on 1330  degrees of freedom
AIC: NA

Number of Fisher Scoring iterations: 4

这是交叉验证:

> cv.glm(orel, m, K = 10)
$call
cv.glm(data = orel, glmfit = m, K = 10)

$K
[1] 10

$delta
[1] 0.2156313 0.2154078

$seed
  [1]         403         526   300751243  -244464717  1066448079  1971573706 -1154513152
  [8]   634841816 -1521293072 -1040655077   505710009  -323431793 -1218609191  1060964279
 [15]  1349082996   -32847357 -1387496845   821178952  -971482876  1295018851  1380491861

现在它收敛了。但是 delta 看起来或多或少是一样的,尽管这个模型看起来更加理智!我现在对交叉验证感到困惑......请给我一个关于如何解释它的提示。

4个回答

与 mambo 所说的类似,delta 值对于将此模型与替代模型进行比较很有用。例如,您可以绘制此模型与可比较模型的 delta 值,以查看产生最低 MSE (delta) 的模型。delta 的第一个值是标准的 k 倍估计值,第二个是偏差校正值。

我开始挖掘包的代码,并https://github.com/cran/boot/blob/5b1e0fea4d1ab1716f2226d673e981d669495b75/R/bootfuns.q#L825boot找到了该函数,并通过 James et 的Introduction to Statistical Learning人。我还没到cv.glm()K-折叠简历部分,但这是我的理解......

的第一个组成部分delta是您从中获得的平均均方误差K-折叠简历。

的第二个组成部分delta是您从中获得的平均均方误差K-fold CV,但带有偏差校正。这是如何实现的,最初,残差平方和 (RSS) 是根据 GLM 预测值和整个数据集的实际响应值计算的。当你正在经历K折叠,生成一个训练模型,然后计算整个数据集之间的 RSSy-values(不仅仅是训练集)和训练模型的预测值。然后从初始 RSS 中减去这些得到的 RSS 值。在你完成后K折叠,你会减去K来自初始 RSS 的值。这是 的第二个组成部分delta

我希望这是正确的,因为这就是我解释代码的方式。

这是代码片段,供您参考。值得庆幸的是,这段代码似乎大部分是独立的。

sample0 <- function(x, ...) x[sample.int(length(x), ...)]
cv.glm <- function(data, glmfit, cost=function(y,yhat) mean((y-yhat)^2),
           K=n)
{
# cross-validation estimate of error for glm prediction with K groups.
# cost is a function of two arguments: the observed values and the
# the predicted values.
    call <- match.call()
    if (!exists(".Random.seed", envir=.GlobalEnv, inherits = FALSE)) runif(1)
    seed <- get(".Random.seed", envir=.GlobalEnv, inherits = FALSE)
    n <- nrow(data)
    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)
    n.s <- table(s)
#    glm.f <- formula(glmfit)
    glm.y <- glmfit$y
    cost.0 <- cost(glm.y, fitted(glmfit))
    ms <- max(s)
    CV <- 0
    Call <- glmfit$call
    for(i in seq_len(ms)) {
        j.out <- seq_len(n)[(s == i)]
        j.in <- seq_len(n)[(s != i)]
        ## we want data from here but formula from the parent.
        Call$data <- data[j.in, , drop=FALSE]
        d.glm <- eval.parent(Call)
        p.alpha <- n.s[i]/n
        cost.i <- cost(glm.y[j.out],
                       predict(d.glm, data[j.out, , drop=FALSE],
                               type = "response"))
        CV <- CV + p.alpha * cost.i
        cost.0 <- cost.0 - p.alpha *
            cost(glm.y, predict(d.glm, data, type = "response"))
    }
    list(call = call, K = K,
         delta = as.numeric(c(CV, CV + cost.0)),  # drop any names
         seed = seed)
}

这可能有助于理解预测误差(delta):

R交叉验证cv.glm:预测误差和置信区间

AdamO 的回答特别有帮助:

“预测误差在两个关键方面不同于标准误差。

  1. 预测误差为预测值提供区间,即在控制预测变量中的部分或全部变化(通过调节)的结果中可以观察到的值。标准误差为估计的统计数据提供区间,例如从未真正观察到的参数。逻辑回归模型中的对数优势比等连续值参数可以以混淆矩阵的形式为二元结果创建“预测区间”(这对于贝叶斯来说是很自然的)。

  2. 预测误差不会在大 n 中消失,而置信区间会。这是因为采样量不会减少从数据生成机制中提取的单个观测值所固有的可变性。然而,由于估计的预测模型的精度提高了,因此预测误差确实会大大降低。由于中心极限定理(usu.),置信区间确实消失了。这是因为重复对宇宙进行采样会产生完全相同的结果,但变异为 0。”

我在网上找到了这个,这有助于解释 delta 是什么:http: //home.strw.leideuniv.nl/~jarle/IAC/Tasks/IAC-lecture4-homework.pdf

在我看来,模型之间 delta 的比较值比绝对值更重要。