7

我在 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？
4

3 回答 3

4

因此，接下来是对@NPR 从 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)

并获得生成的原始和交叉验证的 r-square：

  raw_rsq    cv_rsq 
0.7237907 0.7050297

警告：虽然raw_rsq显然是正确的，并且cv_rsq在我期望的范围内，但请注意，我还没有检查过该函数的确切crosval作用。因此，使用风险自负，如果有人有任何反馈，我们将非常欢迎。它也仅适用于具有截距和标准主效应符号的线性模型。

于 2013-04-16T06:16:30.670 回答
1

我为此编写了一个函数。它也适用于名义预测变量。它只适用于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
于 2016-04-13T08:00:00.597 回答
0

关于 stats.stackexchange（例如，链接 1和链接 2）的讨论认为应该使用均方误差 (MSE) 而不是R^2.

留一法交叉验证（k-folds cv 的特殊情况，其中 k=N）具有允许使用简单公式快速计算线性模型的 CV MSE 的属性。请参阅“R 中的统计学习简介”的第 5.1.2 节。以下代码应计算lm模型的 RMSE 值（使用同一部分的公式 5.2）：

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

您可以将其与“常规”RMSE 进行比较：

summary(fit)$sigma

或者 RMSE 从 5 倍或 10 倍交叉验证中获得，我想。

于 2017-11-09T21:10:20.803 回答