以下は、@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-2 乗を取得します。
raw_rsq cv_rsq
0.7237907 0.7050297
警告:raw_rsq
は明らかに正しく、私が期待する大まかな範囲内にありますが、関数がcv_rsq
何をするのかをまだ正確に調べていないことに注意してください。crosval
したがって、自己責任で使用し、フィードバックがあれば大歓迎です。また、切片と標準の主効果表記を使用する線形モデル用にのみ設計されています。