आर में रैखिक मॉडल से क्रॉस-मान्य आर-वर्ग कैसे प्राप्त करें?

Question

मुझे लगता है मैं नमूना आर वर्ग से बाहर के एक अनुमान प्राप्त करने के लिए चाहते हैं आर

में

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) 

एक रेखीय मॉडल की है। मैं कुछ रूप के-गुना क्रॉस सत्यापन का उपयोग करने के बारे में सोच रहा था।

• आर में कौन सा कोड रैखिक मॉडल फिट लेता है और एक क्रॉस-मान्य आर-वर्ग देता है?
• या आर का उपयोग कर क्रॉस-मान्य आर-वर्ग प्राप्त करने के लिए कोई अन्य दृष्टिकोण है?

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

Figures.stackexchange पर चर्चा (उदाहरण के लिए, link 1, और link 2) तर्क देते हैं किके बजाय औसत-स्क्वायर त्रुटि (एमएसई) का उपयोग किया जाना चाहिए।

छोड़ें-एक-आउट क्रॉस-सत्यापन (के-फ़ोल्डर्स सीवी का विशेष मामला जहां के = एन) एक संपत्ति है जो एक सरल सूत्र का उपयोग करके रैखिक मॉडल के लिए सीवी एमएसई की त्वरित गणना की अनुमति देता है। "आर में सांख्यिकीय सीखने के परिचय" की धारा 5.1.2 देखें।

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

जो तुम "नियमित" RMSE के लिए तुलना कर सकते हैं:

summary(fit)$sigma 

या RMSE 5 से प्राप्त निम्नलिखित कोड lm मॉडल (एक ही अनुभाग से समीकरण 5.2 का प्रयोग करके) के लिए RMSE मूल्य की गणना करना चाहिए या 10 गुना क्रॉस-सत्यापन, मुझे लगता है।