R/tenFoldCV.glm.R

In PELVIS: Probabilistic Sex Estimate using Logistic Regression, Based on VISual Traits of the Human Os Coxae

tenFoldCV.glm <- function(dat.glm, mod, conf_level = 0.95) {
### dat.glm: dataframe on which to perform a 10-fold CV (with no missing values)
### mod: GLM to evaluate
### conf_level: 0.95 by default, confidence level needed to produce a sex estimate

### Prepare constants and data:
    N <- nrow(dat.glm) # sample size in the dataset
    ## Reorder the rows to introduce some randomness:
    dat.glm <- dat.glm[sample(1:N, size = N, replace = FALSE) , ]

######################################################
### 1. Prepare and initialize a dataframe of results #
######################################################
    MatRes <- matrix(nrow = N, ncol = 3)
    colnames(MatRes) <- c("True_sex", "Proba_10CV", "Predict_sex_10CV")
    MatRes <- data.frame(MatRes)
    MatRes[ , 1] <- as.character(dat.glm$Sex)

#############################
### 2. Perform a 10-fold CV #
#############################
    for (k in 1:10) { # for each of the ten folds ("splits", "folds"),
        pioch <- (1:N <= k*N/10) & (1:N > (k-1)*N/10) # draw at random 10% of individuals as a test sample
        app <- dat.glm[!pioch, ] # and train a GLM on the remaining 90%.
        val <- dat.glm[pioch, ]
        modTemp <- glm(mod$call$formula, family = binomial, data = app)
        MatRes[pioch, 2] <- predict(modTemp, newdata = val, type = "response")
    }
    MatRes[ , 3] <- finalSE(MatRes[ , 2], conf_level)
    
#######################
### 3. Model accuracy #
#######################
    indet_rate <- nrow(MatRes[MatRes$Predict_sex_10CV == "I", ]) / nrow(MatRes)
    if (indet_rate < 1) { 
        success_rate <- nrow(MatRes[(MatRes$Predict_sex_10CV=="F" & MatRes$True_sex=="F") | (MatRes$Predict_sex_10CV=="M" & MatRes$True_sex=="M"), ]) / nrow(MatRes[MatRes$Predict_sex_10CV!="I", ])
    } else { # if there is only indeterminate individuals (might happen if one single trait is used as predictor)
        success_rate <- NA
    }
    
###########################
### 4. Return the results #
###########################
    return(list(ClassifResults = MatRes, IndetRate = indet_rate, SuccessRate = success_rate))
}