从 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 个变量：Estimate、Std。误差，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 案例（相对于参数数量的观察很少）可能具有极差的统计特性，甚至可能是不可能的，除非您使用诸如惩罚回归之类的技术。 ..这也可能是您的许多参数从估计中丢失的原因（即您只从可能的 1097 个变量中得到 94 个） - 特别是因为您的表达式模式很简单（只有 0、1 或 2 )，大量的参数是共线的，不能联合估计（你应该在你的原始模型拟合中看到很多NAs）。

获取示例数据：

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 在系数向量 (coef()) 中留下了一个 NA 用于非估计参数，但不会类似地填写摘要中系数表的行。

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