1
$\begingroup$

I am relatively new to statistical learning and need some advice regarding the use of nested cross-validation for model selection and performance evaluation for binary logistic regression with LASSO in R. My understanding is that the outer loop will go through the model parameters (alpha and lambda). For each alpha, I have picked lambda.min using cv,glmnet LOOCV, and then performed 10-cross-validation on the selected parameters to evaluate the model performance (sensitivity and specificity at predefined threshold for example). Here is the full description.

    # create a vector that allocates each observation to one of k-folds
    # and create a matrix in which we store the results
    k = 10 
    A = 10 # alpha
    set.seed(1234)
    folds = sample(1:k, nrow(x), replace=TRUE)
    cv.sen = matrix(NA, k , A) # k x A matrix
    cv.spe = matrix(NA, k , A) # k x A matrix


    # Nested cross validation
    # in the jth fold, the elements of folds that equal j are in the test set, and the remainder are in the training set
    # we make our prediction for each alpha and lambda.min, compute test errors on the appropriate test set,
    # store them in the appropriate slot in the matrix cv.sen etc. This gives kxA matrix of which the (i,j) element
    # corresponds to the test set sensitivity for the ith cross-validation fold for the best model at ath alpha model.
    for (a in 1:A){
      alpha = a/10

      for (j in 1:k){
        cvfit = cv.glmnet(x[folds!=j,], y=y[folds!=j], family = "binomial", alpha=alpha, nfolds=length(y[folds!=j]))
        best.fit = glmnet(x[folds!=j,], y=y[folds!=j], family = "binomial", alpha=alpha, lambda=cvfit$lambda.min)

        fitted.results <- predict(best.fit, newx = x[folds==j,],  type = "response")
        fitted.results_threshold <- ifelse(fitted.results > threshold, 1, 0)
        table <- table(pred = fitted.results_threshold, true = y[folds==j])

        sen = round(table[4]/(table[3]+table[4]),4)
        spe = round(table[1]/(table[1]+table[2]),4)
        cv.sen[j,a] = sen
        cv.spe[j,a] = spe

      }
    }

    mean.cv.sen = apply(cv.sen, 2, mean, na.rm=TRUE)
    mean.cv.spe = apply(cv.spe, 2, mean, na.rm=TRUE)

enter image description here

The mean sensitivity and specificity plots can then be used to select the final model, for instance Alpha = 0.6 seems appropriate. The final model was then applied to a hidden test set and the generalisation error was confirmed to be very low.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.