[½»Á÷] ÇóÖú ÈçºÎʵÏÖleava one out

cv.lm <- function(obj, loo = TRUE) {
#
# The Leave-One-Out (LOO) and/or Leave-Group-Out (LGO) Cross-Validation in R for (Multiple) Linear Regression.
#
# Input:
#   obj: the model of MLR
#   loo: logic, if FALSE, do LGO CV
#
# Note:
#   the following two varible should be assigned according to the number of observation in the data set.
#   m: the number of loop when doing validation.
#   n: the number of removed observations.
#
# Output:
#   q.squared: cross-validation relation coefficient.
#   SDEP: Standard Deviation of Error of Predictions
#   newsq: variance in Y explained only for LOO CV
#
# Usage:
# loo <- cv.lm(obj)
# lgo <- cv.lm(obj, loo = FALSE)
#
# Copyright (C) 2005
#
        data <- data.frame(y = obj$model[,1], x = obj$model[, -1])
        N <- nrow(data)
        if (loo == TRUE) {
                ytest <- data[, 1]
                newrsq <- numeric(N)
                ypred <- numeric(N)
                for(i in 1:N){
                        newtrain <- data[-i, ]
                        xtest <- data[i, -1]
                        newfm <- lm(y ~., data = newtrain)
                        newrsq[i] <- summary(newfm)$r.squared
                        ypred[i] <- predict(newfm, xtest)
                }
        }
        else {
                m <- 40
                n <- round(N / 10)
                newrsq <- numeric(m)
                ytest <- numeric(m*n)
                ypred <- numeric(m*n)
                for (i in 1:m) {
                        j <- sample(N, n, replace = FALSE)
                        newtrain <- data[-j, ]
                        xtest <- data[j, -1]
                        ytest[((i-1) * n + 1):(i * n)] <- data[j, 1]
                        newfm <- lm(y ~., data = newtrain)
                        ypred[((i-1) * n + 1):(i * n)] <- predict(newfm, xtest)
                }
        }
        q.squared <- 1 - sum((ytest - ypred)^2) / sum((ytest - mean(ytest))^2)
        if (loo == TRUE) {
                SDEP <- sqrt(sum((ytest - ypred)^2) / N)
        }
        else {
                SDEP <- sqrt(sum((ytest - ypred)^2) / (m * n))
        }
        if (loo == TRUE) {
                return(list(q.squared = q.squared, SDEP = SDEP, newrsq = newrsq))
        }
        else {
                return(list(q.squared = q.squared, SDEP = SDEP))
        }
}