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.
fonte
2013-04-16 06:16:30
Può essere off-topic .. e buona [cross-validated] (http://stats.stackexchange.com/). –
Perché? Si tratta di come implementare una tecnica statistica nella lingua [r] (http://stackoverflow.com/tags/r/info) che ha circa 30.000 domande. Se preferisci, potrei rimuovere gli elementi statistici della domanda e concentrarmi solo sull'implementazione R? –
Dai un'occhiata a http://www.statmethods.net/stats/regression.html – NPE