循环实现 Leave-One-Out 观察并运行 glm，一次一个变量答案

【问题标题】：Loop to implement Leave-One-Out observation and run glm, one variable at a time循环实现 Leave-One-Out 观察并运行 glm，一次一个变量
【发布时间】：2016-10-09 17:20:23
【问题描述】：

我有一个包含 96 个观察值和 1106 个变量的数据框。

我想对观察结果进行逻辑回归，方法是留下一个，一次一个。（因此，对于第一组观察，总共有 95 个观察被删除，第二组观察总共有 95 个，第二个观察被删除，依此类推，因此有 95 组观察，每组都有一个观察被忽略了。）
此外，我想一次只对一个变量运行每组这些观察结果。在对一个变量进行 95 次观察后，我想提取 p 值（省略截距 p 值）。
我已经能够手动完成所有这些，一次一个。但是这样做 96 次是非常乏味的，我确信必须有一种方法可以通过一个或多个循环来自动执行此操作。

这是我如何手动执行 10 次观察的演示。

## Create 10 data frames by removing one observation from each ##
di.1 <- mainDF [-1,]
di.2 <- mainDF [-2,]
di.3 <- mainDF [-3,]
di.4 <- mainDF [-4,]
di.5 <- mainDF [-5,]
di.6 <- mainDF [-6,]
di.7 <- mainDF [-7,]
di.8 <- mainDF [-8,]
di.9 <- mainDF [-9,]
di.10 <- mainDF [-10,]

## Create data frames to put each p-value result in ## 
dt.1 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.2 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.3 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.4 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.5 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.6 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.7 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.8 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.9 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)
dt.10 <- data.frame(intercept=numeric(), gene=numeric(), stringsAsFactors=FALSE)

## Run logistic regression on each data frame with one one obs. left out ##
## GLM run on one variable at a time##
## Extract p-values and put in separate dfs ##

for (i in 2:1106)
{
  formulas <- glm(response ~ di.1[,i], data=di.1, family= "binomial")
  dt.1[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.2[,i], data=di.2, family= "binomial")
  dt.2[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.3[,i], data=di.3, family= "binomial")
  dt.3[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.4[,i], data=di.4, family= "binomial")
  dt.4[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.5[,i], data=di.5, family= "binomial")
  dt.5[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.6[,i], data=di.6, family= "binomial")
  dt.6[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.7[,i], data=di.7, family= "binomial")
  dt.7[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.8[,i], data=di.8, family= "binomial")
  dt.8[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.9[,i], data=di.9, family= "binomial")
  dt.9[i,] <- coef(summary(formulas))[,4]
}
for (i in 2:1106)
{
  formulas <- glm(response ~ di.10[,i], data=di.10, family= "binomial")
  dt.10[i,] <- coef(summary(formulas))[,4]
}

## Remove intercept p-values ##
dt.1<- dt.1[-c(1)]
dt.2<- dt.2[-c(1)]
dt.3<- dt.3[-c(1)]
dt.4<- dt.4[-c(1)]
dt.5<- dt.5[-c(1)]
dt.6<- dt.6[-c(1)]
dt.7<- dt.7[-c(1)]
dt.8<- dt.8[-c(1)]
dt.9<- dt.9[-c(1)]
dt.10<- dt.10[-c(1)]

## Export data frames, then manually copy and paste them into one CSV ##
write.csv(dt.1, file = "MyData.csv")
write.csv(dt.2, file = "MyData2.csv")
write.csv(dt.3, file = "MyData3.csv")
write.csv(dt.4, file = "MyData4.csv")
write.csv(dt.5, file = "MyData5.csv")
write.csv(dt.6, file = "MyData6.csv")
write.csv(dt.7, file = "MyData7.csv")
write.csv(dt.8, file = "MyData8.csv")
write.csv(dt.9, file = "MyData9.csv")
write.csv(dt.10, file = "MyData10.csv")

我想知道我怎样才能完成所有这些工作，而不必一次检查每个观察。

这是我正在使用的数据块：

  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

非常感谢您的宝贵时间！

【问题讨论】：

每个变量都是一个基因，我正在计算每个基因在所有观察集中被发现显着的次数。这是我用来确定哪些基因对我的研究具有重要意义的不同方法之一。
是的，不要那样做。使用正则化方法，例如弹性网。你可以使用 glmnet 包。
也标记为cross-validation，因为 LOO 是 CV 的一种形式。

标签： r loops regression cross-validation glm

【解决方案1】：

正如我之前在评论中所说，我不会使用 glm 和 summary.glm，因为这对于您的任务来说太慢了，因为您将适合 96 * 1106 GLM。我将使用glm.fit，并自己计算回归系数的 p 值。下面的函数f 就是这样做的。它需要一个一维向量x 作为协变量（不允许NA）和另一个一维向量y 作为响应（不允许NA）。由于做了Logistic回归，所以要求y是两个水平（或0-1二进制值）的因子。

f <- function (x, y) {
  ## call `glm.fit`
  fit <- glm.fit(cbind(1,x), y, family = binomial())
  ## estimated regression coefficients
  beta <- unname(fit$coefficients)
  ## since there are only two coefficients, I don't bother using `chol2inv`
  ## then extract square root of diagonals for standard errors
  se <- sqrt(diag(chol2inv(fit$qr$qr, size = fit$qr$rank)))
  ## deal with possible rank-deficient case
  if (length(se) < 2L) se <- c(se, NA_real_)
  ## z-score
  z <- beta / se
  ## p-value (0.05 significance level)
  2 * pnorm(-abs(z))
  }

如果您不相信它的正确性，我们可以对此功能进行测试。以您的示例数据框dat 为例，我们做Response ~ X1：

dat <- 
structure(list(Response = structure(c(1L, 1L, 1L, 2L, 1L, 1L, 
2L, 2L, 1L, 1L), .Label = c("N", "Y"), class = "factor"), X1 = c(1L, 
2L, 2L, 1L, 2L, 2L, 2L, 2L, 1L, 2L), X2 = c(1L, 1L, 1L, 1L, 2L, 
1L, 1L, 1L, 1L, 1L), X3 = c(1L, 1L, 2L, 2L, 1L, 2L, 1L, 1L, 1L, 
2L), X4 = c(0L, 0L, 1L, 0L, 1L, 1L, 0L, 0L, 0L, 1L), X5 = c(1L, 
2L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 1L), X6 = c(0L, 2L, 0L, 0L, 0L, 
0L, 0L, 0L, 0L, 0L), X7 = c(1L, 1L, 1L, 0L, 0L, 0L, 0L, 0L, 0L, 
1L), X8 = c(0L, 2L, 1L, 1L, 0L, 0L, 0L, 1L, 0L, 0L), X9 = c(2L, 
2L, 0L, 1L, 1L, 2L, 2L, 0L, 1L, 0L), X10 = c(2L, 2L, 1L, 1L, 
1L, 1L, 0L, 2L, 0L, 2L)), .Names = c("Response", "X1", "X2", 
"X3", "X4", "X5", "X6", "X7", "X8", "X9", "X10"), row.names = c("P1", 
"P2", "P3", "P4", "P5", "P6", "P7", "P8", "P9", "P10"), class = "data.frame")

## code response into factor
dat[[1]] <- factor(dat[[1]])

## call `f`
f(dat[[2]], dat[[1]])
# [1] 0.8559137 0.8804148

## call `glm` + `summary.glm`
coef(summary(glm(Response ~ X1, data = dat, family = binomial())))
#              Estimate Std. Error    z value  Pr(>|z|)
#(Intercept) -0.4700036   2.588435 -0.1815783 0.8559137
#X1          -0.2231436   1.483239 -0.1504434 0.8804148

所以 p 值匹配！

我们现在需要另一个函数g 来组织您计划执行的双嵌套循环。外层循环控制“leave-one-out”，而内层循环由lapply 安排以循环遍历数据框列。在外循环的每次迭代结束时，将生成的 p 值数据框写入“.csv”文件。

g <- function (dat) {
  ## convert response to factor (if it is not readily is)
  y <- as.factor(dat[[1]])
  ## leave-one-out
  for (i in 1:nrow(dat)) {
    ## covariates data frame
    covariates <- dat[-i, -1]
    ## response vector
    response <- y[-i]
    ## call `f` to get a data frame of p-values
    result <- as.data.frame(lapply(covariates, f, y = response))
    ## write data frame to file
    write.csv(result, file = paste0(i,".csv"), row.names = FALSE)
    }
  }

当我运行g(dat) 时，我按预期得到了十个“.csv”文件。

跟进：

我仍在掌握如何在 R 中执行循环，因此看到这一点非常有帮助。在将其应用于我的数据时，我是否会将我想使用的数据框的名称放在dat 中？我需要在glm.fit函数部分指定数据框吗？

没有。 glm.fit（以及lm.fit）没有公式接口。直接矩阵计算只需要没有缺失值的数值矩阵来获得估计。这正是它比glm 快的原因。它不接受和消化数据帧。你可以阅读?glm.fit 看看它需要什么参数。

您的数据框dat 不必有列名。如上所述，我们在任何地方都没有使用过公式界面。函数g 假设dat 的第一列是响应，而所有其他列都是自变量。此外，g 不检查缺失值/NA，因此您应确保dat 没有不完整的情况。这些只是g 和f 的要求。

在dat 中包含列名的唯一好处是这些列名将作为标题写入导出的“.csv”文件中，这可能会增加可读性。

【讨论】：