从 R 中的留一法中获取 p 值

Question

我有一个包含 96 个观察值（患者）和 1098 个变量（基因）的数据框。响应是二进制的（Y 和 N），预测变量是数字。我正在尝试执行留一法交叉验证，但我的兴趣不是标准错误，而是从 LOOCV 创建的 95 个逻辑回归模型中每个变量的每个变量的 p 值。这些是我迄今为止的尝试：

#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 个元素的大列表：$coefficients、$residuals、$fitted.values...

尝试 1：

for (i in length(fit)){ 
  model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]])))
}

此输出为 "model_pvalues" 95 个观察值（截距和 94 个变量）和 4 个变量：估计值、标准值。错误，z 值，Pr(>|z|)。然而，我真正想要得到的是所有 1097 个变量的 p 值，对于通过留一法交叉验证构建的 95 个模型。

尝试 2：

for (i in length(fit)){ 
  model_pvalues <- rbind(model_pvalues, coef(summary(fit[[i]]))[4])
}

当我运行这个时，我为一个变量得到一个数字（不确定从哪里来，假设是 beta）。

尝试 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]])))
}

当我运行这样做时，我得到了一个包含 4 个变量的 1520 个观察值的数据框：估计值、标准值。错误，z 值，Pr(>|z|)。观察以 (Intercept) 开头，后跟 82 个变量。之后，它使用 (Intercept1) 和相同的 82 个变量重复此模式，直到 (Intercept15)。

所以我的最终目标是通过 LOOCV 创建 95 个模型，并获得所有模型中使用的所有 1097 个变量的 p 值。非常感谢任何帮助！

编辑：示例数据（1098 个变量的真实 DF 96 观察）

  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

Answer 1

对于 n 个观测值（96 个用于您的真实数据，10 个在示例数据中）和 p 个变量（1098 个用于您的真实数据，10 个在示例数据中），下面的代码应该通过 p 值的 n 列矩阵提取 p 行。我觉得有必要警告你，尝试拟合 n<<p 案例（相对于参数数量的观察很少）可能具有极差的统计特性，甚至可能是不可能的，除非你使用像 penalized 这样的技术回归......这也可能是为什么你的许多参数在估计中丢失的原因（即你只从可能的 1097 个变量中得到 94 个） - 特别是因为你的表达模式很简单（只有 0, 1, or 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()) 对于非估计参数，但不会类似地填写摘要中系数 table 的行。

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

从 R 中的留一法中获取 p 值

Getting p-values from leave-one-out in R

r

bioinformatics

cross-validation