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-fold cross validation을 사용할 생각이었다.

• R의 어떤 코드가 선형 모델에 적합하고 교차 유효성이 검사 된 r-square를 반환합니까?
• 또는 R을 사용하여 교차 유효성이 검사 된 r-square를 얻는 다른 방법이 있습니까?

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) 

을 그리고 그 결과 원료 및 교차 검증 연구를 얻을 수 정사각형 :

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 

그리고 영업 샘플 :

> MOD_k_fold_r2(lmfit) 
#raw_r2 cv_r2 
# 0.724 0.718 

Answer 3

stats.stackexchange에 대한 토론 (예 : link 1 및 link 2)은대신 평균 제곱 오류 (MSE)를 사용해야한다고 주장합니다..

k- 폴드 cv의 특별한 경우 (k = N 인 경우)는 간단한 수식을 사용하여 선형 모델의 CV MSE를 빠르게 계산할 수있는 속성이 있습니다. "R의 통계 학습 소개"섹션 5.1.2를 참조하십시오.

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

는 "일반"RMSE에 비교할 수 :

summary(fit)$sigma 

또는 RMSE는 5에서 얻은 다음 코드는 (같은 섹션에서 식 5.2을 사용하여) lm 모델 RMSE 값을 계산한다 또는 10 배 교차 유효성 검사, 나는 생각한다.