Come ottenere un r-quadrato convalidato incrociato dal modello lineare in R?

Question

Ho un modello lineare in 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) 

desidero ottenere una stima dei fuori di campione r-quadrato. Stavo pensando di usare qualche forma di convalida incrociata k-fold.

• Che codice in R prende un modello di interpolazione lineare e restituisce un cross-validato r-quadrato?
• Oppure c'è qualche altro approccio per ottenere un r-quadro convalidato incrociato usando R?

Answer 1

Quindi ciò che segue è un leggero adattamento a the example that @NPR linked to from statsmethods. Essenzialmente ho adattato l'esempio per renderlo una funzione.

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

Quindi, utilizzando i dati da prima

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

possiamo inserire un modello lineare e chiamare la funzione di validazione croce:

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

E ottenere la risultante r crudo e cross-validati -square:

raw_rsq cv_rsq 
0.7237907 0.7050297 

Avvertenza: Mentre raw_rsq è chiaramente corretto e cv_rsq è nel parco palla che mi aspetto, si noti che non ho ancora esaminato esattamente cosa fa la funzione crosval. Quindi usa a tuo rischio e se qualcuno ha dei feedback, sarebbe il benvenuto. Inoltre è progettato solo per modelli lineari con notazione di intercettazione e effetti principali standard.

Answer 2

Ho scritto una funzione per farlo. Funziona anche per i predittori nominali. Funziona solo per lm oggetti (credo), ma potrebbe essere facilmente ampliato per glm ecc

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

provarla:

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

E sul campione OP:

> MOD_k_fold_r2(lmfit) 
#raw_r2 cv_r2 
# 0.724 0.718 

Answer 3

La discussione su stats.stackexchange (ad esempio, link 1 e link 2) sostiene che l'errore medio quadrato (MSE) deve essere utilizzato anziché R^2.

Convalida incrociata di leave-one-out (caso speciale di k-folds cv dove k = N) ha una proprietà che consente il calcolo rapido del CV MSE per i modelli lineari utilizzando una formula semplice. Vedere la sezione 5.1.2 di "Introduzione all'apprendimento statistico in R". Il seguente codice dovrebbe calcolare il valore RMSE per lm modelli (usando l'equazione 5.2 da stessa sezione):

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

cui si potrebbe paragonare alla RMSE "regolare":

summary(fit)$sigma 

o RMSE ottenuti da 5- o 10 volte la convalida incrociata, suppongo.