Làm cách nào để có được hình chữ nhật r được xác thực chéo từ mô hình tuyến tính trong R?

Question

Tôi có một mô hình tuyến tính trong 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) 

Tôi muốn để có được một ước lượng ra khỏi mẫu r-vuông. Tôi đã nghĩ đến việc sử dụng một số hình thức xác nhận chéo k-fold.

• Mã nào trong R có mô hình tuyến tính phù hợp và trả về quảng cáo r được xác thực chéo?
• Hoặc có cách tiếp cận nào khác để nhận được hình chữ nhật được xác thực chéo bằng R không?

Answer 1

Vì vậy, điều sau đây là một chút thích ứng với the example that @NPR linked to from statsmethods. Về cơ bản tôi đã điều chỉnh ví dụ để biến nó thành một hàm.

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

Vì vậy, bằng cách sử dụng dữ liệu từ trước

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

Chúng ta có thể phù hợp với một mô hình tuyến tính và gọi hàm kiểm chứng chéo:

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

Và có được kết quả r liệu và cross-xác nhận -square:

raw_rsq cv_rsq 
0.7237907 0.7050297 

Nhớ lại: Mặc dù raw_rsq rõ ràng là chính xác và cv_rsq nằm trong công viên bóng mà tôi mong đợi, lưu ý rằng tôi chưa kiểm tra chính xác chức năng của crosval. Vì vậy, sử dụng có nguy cơ của riêng bạn và nếu bất cứ ai có bất kỳ thông tin phản hồi, nó sẽ được chào đón nhất. Nó cũng chỉ được thiết kế cho các mô hình tuyến tính với một đánh chặn và tiêu chuẩn hiệu ứng chính ký hiệu.

Answer 2

Tôi đã viết một chức năng để thực hiện việc này. Nó cũng hoạt động cho những người dự đoán danh nghĩa. Nó chỉ hoạt động cho lm đối tượng (tôi nghĩ), nhưng có thể dễ dàng được mở rộng để glm, vv

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

kiểm tra nó:

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

Và trên mẫu OP:

> MOD_k_fold_r2(lmfit) 
#raw_r2 cv_r2 
# 0.724 0.718 

Answer 3

Thảo luận về số liệu thống kê.stackexchange (ví dụ: link 1 và link 2) cho rằng lỗi bình phương bình phương (MSE) nên được sử dụng thay vì R^2.

Bỏ qua xác thực chéo một lần (trường hợp đặc biệt k-folds cv trong đó k = N) có thuộc tính cho phép tính toán nhanh MSE CV cho các mô hình tuyến tính bằng cách sử dụng công thức đơn giản. Xem phần 5.1.2 của "Giới thiệu về học tập thống kê trong R". Các mã sau đây cần tính toán giá trị RMSE cho lm mô hình (sử dụng Equation 5.2 từ cùng một phần):

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

nào bạn có thể so sánh với RMSE "thường xuyên":

summary(fit)$sigma 

hoặc RMSE thu được từ 5 hoặc xác thực chéo 10 lần, tôi cho là vậy.