Rの線形モデルからクロスバリデーションされたr-squareを得るには？

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倍のクロスバリデーションを使用することを考えていました。

• Rのどのコードが線形モデルに適合し、クロスバリデーションされたr角形を返しますか？
• または、Rを使用してクロスバリデーションされたr-squareを取得するための他の方法がありますか？

Answer 1

次に、the example that @NPR linked to from statsmethodsにわずかに適合しています。基本的に私はこの例を関数にするように修正しました。

だから、私たちは、線形モデルをフィットし、クロスバリデーション機能を呼び出すことができます前に、

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

からのデータを使用して

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

：

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

し、得られた生と交差検定Rを得ます正方形：

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）がではなく使用されるべきであると主張している。

離散クロスバリデーション（k-fold cvのk = Nの特別な場合）は、単純な公式を使用して線形モデルのCV MSEの迅速な計算を可能にする特性を有する。 「Rによる統計学習の入門」の5.1.2節を参照してください。

summary(fit)$sigma 

または5から取得したRMSE：あなたは「通常の」RMSEとの比較ができ

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

：次のコードは、lmモデル（同じセクションから、式5.2を使用して）のためのRMSE値を計算しなければなりませんまたは10倍のクロスバリデーションが可能です。