Р - к-кратная кросс-проверка для линейной регрессии со стандартной погрешностью оценки

Question

Я хотел бы выполнить к-кратная кросс-проверка в R для модели линейной регрессии и протестировать правило один стандартная ошибка:

https://stats.stackexchange.com/questions/17904/one-standard-error-rule-for-variable-selection

Таким образом, мне нужна функция, которая возвращает мне оценку кросс-проверки ошибки предсказания и стандартная ошибка этой оценки (или, по меньшей мере, СКО для каждого раза, так что я могу вычислить стандартную ошибку себя). Во многих пакетах есть функции, которые вычисляют ошибку перекрестной проверки (например, cv.glm в пакете boost), но обычно они возвращают только оценку CV ошибки предсказания, а не ее стандартную ошибку, или MSE для каждой складки.

Я пробовал использовать пакет DAAG, функция которого CVlm должна давать более богатый результат, чем cv.glm. Однако я не могу заставить его работать! Вот мой код:

a=c(0.0056, 0.0088, 0.0148, 0.0247, 0.0392, 0.0556, 0.0632, 0.0686, 0.0786, 0.0855, 0.0937) 
b=c(6.0813, 9.5011, 15.5194, 23.9409, 32.8492, 40.8399, 43.8760, 45.5270, 46.7668, 46.1587, 43.4524) 
dataset=data.frame(x=a,y=b) 
CV.list=CVlm(df=dataset,form.lm = formula(y ~ poly(x,2)), m=5) 

я получаю едва информативную ошибку

Error in xy.coords(x, y, xlabel, ylabel, log) : 
'x' and 'y' lengths differ 

, который не имеет большого смысла для меня. x и y - это ту же длину (11), так что, очевидно, функция жалуется на некоторые другие переменные, создаваемые внутри.

Я с радостью принимаю решения с другими пакетами (например, caret). Кроме того, было бы здорово, если бы я мог указать несколько повторений для кросс-валидации k-fold. Благодаря!

Answer 1

CVlm не нравится poly(x,2) в вашей формуле. Вы можете легко обойти, что при добавлении результата poly(x,2) первый в вашем DataTable и вызова CVlm на этих новых переменных:

dataset2 <- cbind(dataset,poly(dataset$x,2)) 
names(dataset2)[3:4] <- c("p1","p2") 
CV.list=CVlm(df=dataset2,form.lm = formula(y ~ p1+p2)) 

И как вы заинтересованы значениями напечатанных, которые, к сожалению, не сохранены в любом месте вы можете использовать что-то вроде :

# captures the printed output 
printOut <- capture.output(CV.list=CVlm(df=dataset2,form.lm = formula(y ~ p1+p2))) 

# function to parse the output 
# to be adapted if necessary for your needs 
GetValues <- function(itemName,printOut){ 
    line <- printOut[grep(itemName,printOut)] 
    items <- unlist(strsplit(line,"[=]| +")) 
    itemsMat <- matrix(items,ncol=2,byrow=TRUE) 
    vectVals <- as.numeric(itemsMat[grep(itemName,itemsMat[,1]),2]) 
    return(vectVals) 
} 

# get the Mean square values as a vector 
MS <- GetValues("Mean square",printOut) 

Answer 2

Среднее значение MSE хранится как атрибут объекта модели. attributes(CV.list)$ms дает вам то, что вы ищете.