Получение значений p из оставляющего один в R
У меня есть данные 96 наблюдений (пациентов) и 1098 переменных (генов). Ответ является двоичным (Y и N), а предикторы являются числовыми. Я пытаюсь выполнить перекрестную проверку с пропуском, но меня интересует не стандартная ошибка, а p-значения для каждой переменной из каждой из 95 моделей логистической регрессии, созданных из LOOCV. Это мои попытки до сих пор:
#Data frame 96 observations 1098 variables
DF2
fit <- list()
for (i in 1:96){
df <- DF2[-i,]
fit[[i]] <- glm (response ~., data= df, family= "binomial")
}
model_pvalues <- data.frame(model = character(), p_value = numeric())
Эти выходные данные подходят как большой список с 16 элементами и списком из 30: коэффициенты $, $residuals, $fit.values ....
Попытка 1:
for (i in length(fit)){
model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]])))
}
Это выводит в "model_pvalues" 95 наблюдений (Перехват и 94 переменных) и 4 переменные: Estimate, Std. Ошибка, значение z, Pr(>|z|). Однако то, что я действительно пытаюсь получить, это p-значение для всех 1097 переменных, для 95 моделей, построенных путем перекрестной проверки.
Попытка 2:
for (i in length(fit)){
model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]]))[4])
}
Когда я запускаю это, я получаю одно число (не уверен откуда, принимая бета) для одной переменной.
Попытка 3:
for (i in 1:96){
df <- DF2[-i,]
fit[[i]] <- glm (response ~., data= df, family= "binomial")
model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]])))
}
Когда я запускаю это, я получаю набор данных из 1520 наблюдений 4 переменных: Estimate, Std. Ошибка, значение z, Pr(>|z|). Наблюдения начинаются с (Перехват), за которым следуют 82 переменных. После этого он повторяет этот шаблон с (Intercept1) и теми же 82 переменными вплоть до (Intercept15).
Поэтому моя конечная цель - создать 95 моделей через LOOCV и получить p-значения для всех 1097 переменных, используемых во всех моделях. Любая помощь будет очень высоко ценится!
Изменить: пример данных (реальные наблюдения DF 96 для 1098 переменных)
Response X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
P1 N 1 1 1 0 1 0 1 0 2 2
P2 N 2 1 1 0 2 2 1 2 2 2
P3 N 2 1 2 1 1 0 1 1 0 1
P4 Y 1 1 2 0 1 0 0 1 1 1
P5 N 2 2 1 1 1 0 0 0 1 1
P6 N 2 1 2 1 1 0 0 0 2 1
P7 Y 2 1 1 0 2 0 0 0 2 0
P8 Y 2 1 1 0 2 0 0 1 0 2
P9 N 1 1 1 0 2 0 0 0 1 0
P10 N 2 1 2 1 1 0 1 0 0 2
1 ответ
За n наблюдения (96 для ваших реальных данных, 10 в данных примера) и p переменных (1098 для ваших реальных данных, 10 в примере данных), код ниже должен извлечь p ряд за n матрица столбцов р-значений. Я чувствую себя обязанным предупредить вас о том, что n<<p случай (очень мало наблюдений относительно количества параметров), вероятно, будет иметь крайне плохие статистические свойства, и, возможно, даже будет невозможен, если вы не используете метод, такой как штрафная регрессия... это также, вероятно, причина, почему так много ваших параметров отсутствуют в оценках (т.е. вы получаете только 94 из возможных 1097 переменных) - тем более что ваши шаблоны выражений просты (только 0, 1 или 2), большое количество параметров коллинеарны и не могут оценивать совместно (вы бы видели много NAс вашей оригинальной модели подойдет тоже).
Получить пример данных:
DF2 <- read.table(row.names=1,header=TRUE,text="
Resp. X1 X2 X3 X4 X5 X6 X7 X8 X9 X10
P1 N 1 1 1 0 1 0 1 0 2 2
P2 N 2 1 1 0 2 2 1 2 2 2
P3 N 2 1 2 1 1 0 1 1 0 1
P4 Y 1 1 2 0 1 0 0 1 1 1
P5 N 2 2 1 1 1 0 0 0 1 1
P6 N 2 1 2 1 1 0 0 0 2 1
P7 Y 2 1 1 0 2 0 0 0 2 0
P8 Y 2 1 1 0 2 0 0 1 0 2
P9 N 1 1 1 0 2 0 0 0 1 0
P10 N 2 1 2 1 1 0 1 0 0 2")
Подходящие модели
n <- nrow(DF2)
fit <- vector(mode="list",n) ## best to pre-allocate objects
for (i in 1:n) {
df <- DF2[-i,]
fit[[i]] <- glm (Resp. ~., data= df, family= "binomial")
}
В этом случае мы должны быть немного осторожнее при извлечении p-значений, потому что из-за коллинеарности некоторые из них отсутствуют - R оставляет NA в векторе коэффициентов (coef()) для не оцененных параметров, но аналогичным образом не заполняет строки таблицы коэффициентов в сводке.
tmpf <- function(x) {
## extract coef vector - has NA values for collinear terms
## [-1] is to drop the intercept
r1 <- coef(x)[-1]
## fill in values from p-value vector; leave out intercept with -1,
r2 <- coef(summary(x))[-1,"Pr(>|z|)"]
r1[names(r2)] <- r2
return(r1)
}
pvals <- sapply(fit,tmpf)
Конечно, для игрушечного примера все значения p по существу равны 1 ...
## round(pvals,4)
## [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
## X1 0.9998 0.9998 0.9998 0.9998 0.9998 0.9998 0.9999 0.9998 0.9999 0.9998
## X2 0.9999 0.9999 0.9999 0.9999 NA 0.9999 0.9999 0.9999 0.9999 0.9999
## X3 0.9999 0.9999 0.9999 0.9999 0.9999 0.9998 0.9999 0.9999 0.9999 0.9999
## X4 0.9998 0.9998 0.9998 NA 0.9998 0.9998 0.9998 0.9998 0.9998 0.9998
## X5 NA 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 NA 1.0000
## X6 0.9999 NA 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999 0.9999
## X7 1.0000 1.0000 1.0000 1.0000 1.0000 NA 1.0000 1.0000 1.0000 1.0000
## X8 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000
## X9 1.0000 1.0000 NA 1.0000 1.0000 1.0000 NA NA 1.0000 NA
## X10 NA NA NA NA NA NA NA NA NA NA