如何從R中的線性模型獲得交叉驗證的r平方？

Question

我有一個線性模型中R.

set.seed(1234) 
x <- rnorm(100) 
z <- rnorm(100) 
y <- rnorm(100, x+z) 
mydata <- data.frame(x,y,z) 

fit <- lm(y ~ x + z, mydata) 

我想獲得了樣品的r-正方形的的估計。我正在考慮使用某種形式的k-fold交叉驗證。

• R中的哪些代碼需要線性模型擬合併返回經過交叉驗證的r平方？
• 或者還有其他的方法來獲得使用R的交叉驗證的R - 平方？

Answer 1

所以接下來是對the example that @NPR linked to from statsmethods的輕微調整。本質上我調整了這個例子來使它成爲一個函數。

library(bootstrap) 

k_fold_rsq <- function(lmfit, ngroup=10) { 
    # assumes library(bootstrap) 
    # adapted from http://www.statmethods.net/stats/regression.html 
    mydata <- lmfit$model 
    outcome <- names(lmfit$model)[1] 
    predictors <- names(lmfit$model)[-1] 

    theta.fit <- function(x,y){lsfit(x,y)} 
    theta.predict <- function(fit,x){cbind(1,x)%*%fit$coef} 
    X <- as.matrix(mydata[predictors]) 
    y <- as.matrix(mydata[outcome]) 

    results <- crossval(X,y,theta.fit,theta.predict,ngroup=ngroup) 
    raw_rsq <- cor(y, lmfit$fitted.values)**2 # raw R2 
    cv_rsq <- cor(y,results$cv.fit)**2 # cross-validated R2 

    c(raw_rsq=raw_rsq, cv_rsq=cv_rsq) 
} 

因此，使用從數據之前

# sample data 
set.seed(1234) 
x <- rnorm(100) 
z <- rnorm(100) 
y <- rnorm(100, x+z) 
mydata <- data.frame(x,y,z) 

我們可以適應線性模型，並調用交叉驗證功能：

# fit and call function 
lmfit <- lm(y ~ x + z, mydata) 
k_fold_rsq(lmfit, ngroup=30) 

並獲得所產生的原材料和交叉驗證[R -square：

raw_rsq cv_rsq 
0.7237907 0.7050297 

注意事項：雖然raw_rsq顯然是正確的，cv_rsq是我期待的球場，請注意，我還沒有仔細研究crosval函數的功能。因此，請自擔風險，如果有人有任何反饋意見，我們將非常歡迎。它也僅適用於帶截取和標準主效應符號的線性模型。

Answer 2

我寫了一個這樣做的函數。它也適用於名義預測。它僅適用於lm對象（我認爲），但可以很容易地擴展到glm等

# from 
# http://stackoverflow.com/a/16030020/3980197 
# via http://www.statmethods.net/stats/regression.html 

#' Calculate k fold cross validated r2 
#' 
#' Using k fold cross-validation, estimate the true r2 in a new sample. This is better than using adjusted r2 values. 
#' @param lmfit (an lm fit) An lm fit object. 
#' @param folds (whole number scalar) The number of folds to use (default 10). 
#' @export 
#' @examples 
#' fit = lm("Petal.Length ~ Sepal.Length", data = iris) 
#' MOD_k_fold_r2(fit) 
MOD_k_fold_r2 = function(lmfit, folds = 10, runs = 100, seed = 1) { 
    library(magrittr) 

    #get data 
    data = lmfit$model 

    #seed 
    if (!is.na(seed)) set.seed(seed) 

    v_runs = sapply(1:runs, FUN = function(run) { 
    #Randomly shuffle the data 
    data2 = data[sample(nrow(data)), ] 

    #Create n equally size folds 
    folds_idx <- cut(seq(1, nrow(data2)), breaks = folds, labels = FALSE) 

    #Perform n fold cross validation 
    sapply(1:folds, function(i) { 
     #Segement your data by fold using the which() function 

     test_idx = which(folds_idx==i, arr.ind=TRUE) 
     test_data = data2[test_idx, ] 
     train_data = data2[-test_idx, ] 

     #weights 
     if ("(weights)" %in% data) { 
     wtds = train_data[["(weights)"]] 
     } else { 
     train_data$.weights = rep(1, nrow(train_data)) 
     } 

     #fit 
     fit = lm(formula = lmfit$call$formula, data = train_data, weights = .weights) 

     #predict 
     preds = predict(fit, newdata = test_data) 

     #correlate to get r2 
     cor(preds, test_data[[1]], use = "p")^2 
    }) %>% 
     mean() 
    }) 

    #return 
    c("raw_r2" = summary(lmfit)$r.squared, "cv_r2" = mean(v_runs)) 
} 

測試它：

fit = lm("Petal.Length ~ Species", data = iris) 
MOD_k_fold_r2(fit) 
#> raw_r2  cv_r2 
#> 0.9413717 0.9398156 

而且在OP樣本：

> MOD_k_fold_r2(lmfit) 
#raw_r2 cv_r2 
# 0.724 0.718 

Answer 3

討論stats.stackexchange（例如，link 1和link 2）認爲應該使用均方差（MSE）而不是R^2。

一次性交叉驗證（k-fold cv的特殊情況，其中k = N）具有一個屬性，可以使用簡單公式爲線性模型快速計算CV MSE。參見「R中的統計學習入門」第5.1.2節。下面的代碼應該計算RMSE值lm模型（用公式5.2相同的部分）：

sqrt(sum((residuals(fit)/(1-hatvalues(fit)))^2)/length(fit$residuals)) 

，你可以比較 「常規」 RMSE：

summary(fit)$sigma 

或RMSE從5-獲得或者10倍交叉驗證，我想。