Intepretation of crossvalidation result - cv.glm()

Question

My logistic model has been suspicious due to enormous coefficients, so I tried to do a crossvalidation, and also do a crossvalidation of simplified model, to confirm the fact that the original model is overspecified, as James suggested. However, I don't know how to interpret the result (this is the model from the linked question):

> 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

I guess the delta is the mean fitting error, but how to interpret it? Is it a good or bad fit? BTW, the algorithm did not converge, maybe due to the enormous coefficients (?)

I tried a simplified model:

> 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

and it's crossvalidation:

> 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

Now it converged. But the delta seems more or less the same, despite of the fact that this model looks much more sane! I'm confused by the crossvalidation now... please give me a hint on how interpret it.

Answer 1

I started digging through the code for the boot package and found the function cv.glm() at https://github.com/cran/boot/blob/5b1e0fea4d1ab1716f2226d673e981d669495b75/R/bootfuns.q#L825, as well as going through Introduction to Statistical Learning by James et al. I haven't gotten to the $K$-fold CV section yet, but here's my understanding...

The first component of delta is the average mean-squared error that you obtain from doing $K$-fold CV.

The second component of delta is the average mean-squared error that you obtain from doing $K$-fold CV, but with a bias correction. How this is achieved is, initially, the residual sum of squares (RSS) is computed based on the GLM predicted values and the actual response values for the entire data set. As you're going through the $K$ folds, you generate a training model, and then you compute the RSS between the entire data set of $y$-values (not just the training set) and the predicted values from the training model. These resulting RSS values are then subtracted from the initial RSS. After you're done going through your $K$ folds, you will have subtracted $K$ values from the initial RSS. This is the second component of delta.

I'm hoping this is right, as this is how I'm interpreting the code.

Here is the code snippet, for your reference. Thankfully, it appears that this code is mostly self-contained.

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)
}

Answer 2

Similar to what mambo said, the delta values are useful to compare this model with alternative models. You might, for example, plot the delta values of this vs. comparable models to see which produce the lowest MSE (delta). The first value of delta is the standard k-fold estimate and the second is bias corrected.

Answer 3

This might be useful for the understanding of prediction errors (delta):

R crossvalidation cv.glm: prediction error and confidence interval

This answer from AdamO was particularly helpful:

"Prediction errors are different from standard errors in two critical ways.

1. Prediction errors provide intervals for predicted values, i.e. values which could be observed in the outcome controlling for some or all of the variation (through conditioning) in the predictors. Standard errors provide intervals for estimated statistics, e.g. parameters which are never truly observed. Continuously valued parameters such as log odds ratios in a logistic regression model can create "prediction intervals" for binary outcomes in the form of a confusion matrix (this is natural for Bayesians).

2. Prediction errors do not vanish in large n whereas confidence intervals do. This is because no amount of sampling will reduce the variability inherent in a single observation drawn from the data generating mechanism. Prediction errors do decrease in large however, since the precision of the estimated predictive model improves. Confidence intervals do vanish in large as a result of the central limit theorem (usu.). This is because sampling the universe repeatedly would yield the exact same thing with 0 variation."

Answer 4

I found this online, which helps explain what delta is: http://home.strw.leidenuniv.nl/~jarle/IAC/Tasks/IAC-lecture4-homework.pdf

It seems to me that that comparative values of delta between models are of importance rather than the absolute values.