K-кратная или удерживающая перекрестная проверка для регрессии гребня с использованием R

Я работаю над перекрестной проверкой прогноза моих данных с 200 субъектами и 1000 переменных. Меня интересует регрессия гребня, поскольку число переменных (которые я хочу использовать) больше, чем количество выборок. Поэтому я хочу использовать оценки усадки. Ниже приведены примеры данных:

 #random population of 200 subjects with 1000 variables 
    M <- matrix(rep(0,200*100),200,1000)
    for (i in 1:200) {
    set.seed(i)
      M[i,] <- ifelse(runif(1000)<0.5,-1,1)
    }
    rownames(M) <- 1:200

    #random yvars 
    set.seed(1234)
    u <- rnorm(1000)
    g <- as.vector(crossprod(t(M),u))
    h2 <- 0.5 
    set.seed(234)
    y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

    myd <- data.frame(y=y, M)
myd[1:10,1:10]

y X1 X2 X3 X4 X5 X6 X7 X8 X9
1   -7.443403 -1 -1  1  1 -1  1  1  1  1
2  -63.731438 -1  1  1 -1  1  1 -1  1 -1
3  -48.705165 -1  1 -1 -1  1  1 -1 -1  1
4   15.883502  1 -1 -1 -1  1 -1  1  1  1
5   19.087484 -1  1  1 -1 -1  1  1  1  1
6   44.066119  1  1 -1 -1  1  1  1  1  1
7  -26.871182  1 -1 -1 -1 -1  1 -1  1 -1
8  -63.120595 -1 -1  1  1 -1  1 -1  1  1
9   48.330940 -1 -1 -1 -1 -1 -1 -1 -1  1
10 -18.433047  1 -1 -1  1 -1 -1 -1 -1  1

Я хотел бы сделать следующее для перекрестной проверки -

(1) разделить данные на две части - использовать первую половину в качестве тренировки, а вторую половину - в качестве теста

(2) Кросс-валидация в K-кратном порядке (скажем, 10-кратное или предложение о любом другом подходящем для моего случая случае приветствуется)

Я могу просто сэмплировать данные на две части (получить и проверить) и использовать их:

# using holdout (50% of the data) cross validation 
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

 myd_train <- myd[training.id,]
 myd_test  <- myd[test.id,]   

Я использую lm.ridgeиз MASSпакета R

library(MASS)
out.ridge=lm.ridge(y~., data=myd_train, lambda=seq(0, 100,0.001))
plot(out.ridge)
select(out.ridge)

lam=0.001
abline(v=lam)

out.ridge1 =lm.ridge(y~., data=myd_train, lambda=lam)
hist(out.ridge1$coef)
    out.ridge1$ym
hist(out.ridge1$xm)

У меня два вопроса -

(1) Как я могу предсказать набор тестов и рассчитать точность (как соотношение прогнозируемого и фактического значений)?

(2) Как я могу выполнить K-кратную проверку? скажете в 10 раз?

Вы можете использовать caret пакет (виньетки , бумага ) для этого типа вещей, который может обернуть несколько моделей машинного обучения или вы можете использовать свои собственные индивидуальные модели . Поскольку вас интересует регрессия гребня, здесь представлены только пользовательские коды для регрессии гребня, вы можете более точно адаптироваться к вашей ситуации.

Для простого разделения данных:

set.seed(107)
# stratified random split of the data
inTrain <- createDataPartition(y = myd$y, p = .5,list = FALSE)
training <- myd[ inTrain,]
testing <- myd[-inTrain,]

Для проверки K-fold и других типов резюме, включая загрузку по умолчанию

ridgeFit1 <- train(y ~ ., data = training,method = 'ridge', 
preProc = c("center", "scale"), metric = "ROC")
plot(ridgeFit1)

Вот обсуждение того, как использовать trainфункцию. Обратите внимание, что метод ridge зависит от elasticnetфункций пакета (и его зависимости от того lars, должен или должен быть установлен). Если не установлен в системе, вас спросят, хотите ли вы это сделать.

используемый тип повторной выборки. По умолчанию используется простая начальная загрузка. Для изменения метода повторной выборки используется функция trainControl.

Метод option контролирует тип пересэмплирования и по умолчанию «boot». Другой метод, "repeatcv", используется для указания повторной перекрестной проверки в K-кратном порядке (и аргумент повторяет контролирует количество повторений). K контролируется аргументом числа и по умолчанию 10.

 ctrl <- trainControl(method = "repeatedcv", repeats = 5)

 ridgeFit <- train(y ~ ., data = training,method = 'ridge',
preProc = c("center", "scale"),trControl = ctrl, metric = "ROC")

plot(ridgefit)

Для прогнозов:

plsClasses <- predict(ridgeFit, newdata = testing)

Это продолжение предложения Фрэнка в комментариях. Доктор Харрель, пожалуйста, исправьте, если я ошибаюсь (оцените исправления).

Ваши данные:

#random population of 200 subjects with 1000 variables 
    M <- matrix(rep(0,200*100),200,1000)
    for (i in 1:200) {
    set.seed(i)
      M[i,] <- ifelse(runif(1000)<0.5,-1,1)
    }
    rownames(M) <- 1:200

    #random yvars 
    set.seed(1234)
    u <- rnorm(1000)
    g <- as.vector(crossprod(t(M),u))
    h2 <- 0.5 
    set.seed(234)
    y <- g + rnorm(200,mean=0,sd=sqrt((1-h2)/h2*var(g)))

    myd <- data.frame(y=y, M)

Установите rmsпакет и загрузите его.

require(rms)

ols функция используется для оценки линейной модели с использованием обычных наименьших квадратов, где можно указать срок штрафа.

Как предлагается ниже в комментариях, я добавил petraceфункцию. Эта функция отслеживает AIC и BIC против казни.

# using holdout (50% of the data) cross validation 
training.id <- sample(1:nrow(myd), round(nrow(myd)/2,0), replace = FALSE)
test.id <- setdiff(1:nrow(myd), training.id)

 myd_train <- myd[training.id,]
 myd_test  <- myd[test.id,] 

frm <- as.formula(paste("y~",paste(names(myd_train)[2:100],collapse="+")))

Важное примечание Я не смог использовать все 1000 переменных, поскольку программа жалуется, если число переменных превышает 100. Также y~.не сработало обозначение формулы типа. Итак, см. Выше способ сделать то же самое создание объекта формулыfrm

f <- ols(frm, data = myd_train, method="qr", x=TRUE, y=TRUE)


p <- pentrace(f, seq(.2,1,by=.05))

Error in array(x, c(length(x), 1L), if (!is.null(names(x))) list(names(x),  : 
'data' must be of a vector type, was 'NULL'

 plot(p)

«Для обычного непенализованного подбора из lrm или ols и для вектора или списка штрафов подходит серия логистических или линейных моделей, использующих штрафную оценку максимального правдоподобия, и сохраняет эффективные степени свободы, информационный критерий Акаике (AIC), Шварц Байесовский Information Criterion (BIC) и исправленный AIC Хурвича и Цая (AIC_c). При желании pentrace может использовать функцию nlminb для определения оптимального штрафного коэффициента или комбинации факторов, штрафующих различные типы терминов в модели ». из rmsпакета руководства.

calibrateФункция предназначена для калибровки модели с передискретизацией и использует начальную загрузку или перекрестную проверку, чтобы получить исправленные смещения (с поправкой на переоснащение) оценки прогнозируемых и наблюдаемых значений на основе подстановки прогнозов в интервалы. validateФункция не передискретизации проверки регрессионной модели, с или без обратного переменного удаления понижающего. B = количество повторений. Для метода = "перекрестная проверка" - это число групп пропущенных наблюдений.

cal <- calibrate(f, method = "cross validation", B=20)  
plot(cal)

Вы можете использовать Predictфункцию для вычисления прогнозируемых значений и доверительных интервалов. Я не уверен, что это работает в тестовой ситуации.

Пакет R glmnet( виньетка ) имеет функцию-обертку, которая делает именно то, что вы хотите, называется cv.glmnet( док ). Я только вчера использовал его, он работает как сон.