lm()：循环通过多个线性模型导出 F 统计量的 p 值答案

【问题标题】：lm(): loop through multiple linear models exporting p-value of F-statisticlm()：循环通过多个线性模型导出 F 统计量的 p 值
【发布时间】：2016-10-19 20:04:43
【问题描述】：

我有一个大型数据集，我需要运行一个线性模型来比较组。我需要使用线性模型找到组比较的 p 值。有四组（所以我需要 1~2、1~3、1~4、2~3、2~4、3~4），并且有 130 列需要比较这些组的数据。任何帮助将不胜感激！

我有这个，这正是我所需要的。

fit<-lm(variable~group, data=data)
summary(fit)

但是，对于所有组和列，我要进行近 800 次比较，因此我想避免手动进行。我尝试编写一个 for 循环，但它不起作用。

k<-data.frame()
for (i in 1:130){
 [i,1]<-colnames(data)
 fit<- lm(i~group, data=data)
 [i,2] <- fit$p.value
}

但这给了我各种不同的错误。我真的只需要 p 值。帮助将不胜感激！谢谢！

【问题讨论】：

你能添加一个你的数据是什么样子的样本吗？即结构。这样可以更轻松地为您提供帮助。
数据有一列带有组号（例如，每一行的第一列有0、1、2或3）；这些是要比较的组。然后每个后续列都有一个我对值感兴趣的特定因素。我想单独比较每列中的组。我使用 read.csv 导入数据，然后使用 head(data)、string(data)。这有帮助吗？

标签： r loops regression linear-regression lm

【解决方案1】：

(2016-06-18) 您的问题目前无法完全回答。下面，我将指出几个问题。

如何正确获取p值

我假设您想要模型的 F 统计量的 p 值，作为拟合优度的指示。假设你的拟合模型是fit，我们应该这样做：

fstatistic <- summary(fit)$fstatistic
p_value <- unname(1 - pf(fstatistic[1], fstatistic[2], fstatistic[3]))

我以内置数据集trees为例进行演示。

fit <- lm(Height ~ Girth, trees)
## truncated output of summary(fit)
# > summary(fit)
# Residual standard error: 5.538 on 29 degrees of freedom
# Multiple R-squared:  0.2697,  Adjusted R-squared:  0.2445 
F-statistic: 10.71 on 1 and 29 DF,  p-value: 0.002758

fstatistic <- summary(fit)$fstatistic
p_value <- unname(1 - pf(fstatistic[1], fstatistic[2], fstatistic[3]))
## > p_value
# [1] 0.002757815

所以，p_value 同意打印的摘要。

你的循环

我建议您在计算/更新期间使用向量而不是数据框。

variable <- character(130)
p.value <- numeric(130)

您可以通过以下方式将最后的结果组合到一个数据框：

k <- data.frame(var = variable, p.value = p.value)

为什么？因为这是内存效率！现在，经过这些修正，我们得出：

variable <- character(130)
p.value <- numeric(130)
for (i in 1:130) {
  variable[i] <- colnames(data)
  fit <- lm(i~group, data=data)
  fstatistic <- summary(fit)$fstatistic
  p_value <- unname(1 - pf(fstatistic[1], fstatistic[2], fstatistic[3]))
  p.value[i] <- p_value
  }
k <- data.frame(var = variable, p.value = p.value)

其他问题

我仍然认为上面的代码行不通。因为我不确定以下是否正确：

  variable[i] <- colnames(data)
  fit <- lm(i~group, data=data)

在循环过程中，data 没有改变，所以colnames(data) 返回一个向量，因此var[i] <- colnames(data) 会触发错误。
i~group 看起来很奇怪。你的data 中有i 吗？

我无法帮助您解决这些问题。我不知道你的data 长什么样。但是，如果您可以放入数据的子集，那就没问题了。

跟进（2016-06-19）

谢谢。这非常有帮助。我的数据中没有“i”，但我希望我可以用它来表示不同的列名，以便它遍历所有列名。有没有办法分配列名编号，这样就可以了？

是的，但我需要知道你对每一列都有什么。

第 1 列有一个组号。以下列包含我正在查看的不同因素的数据。

好的，所以我想ncol(data) = 131，其中第一列是group，剩下的 130 列是你要测试的。那么这应该工作：

variable <- colnames(data)[-1]
p.value <- numeric(130)
for (i in 1:130) {
  fit <- lm(paste(variable[i], "group", sep = "~"), data=data)
  fstatistic <- summary(fit)$fstatistic
  p_value <- unname(1 - pf(fstatistic[1], fstatistic[2], fstatistic[3]))
  p.value[i] <- p_value
  }
k <- data.frame(var = variable, p.value = p.value)

可以使用sapply() 代替上面的for 循环。但我认为没有性能差异，因为与lm() 和summary() 相比，循环开销非常小。

【讨论】：

谢谢。这非常有帮助。我的数据中没有“i”，但我希望我可以用它来表示不同的列名，以便它遍历所有列名。有没有办法分配列名编号，这样就可以了？不幸的是，我无法上传我的数据，因为它在我的工作计算机上，并且我不允许共享它。我可以尝试在下一条评论中更好地解释它的外观：
第 1 列有一个组号（例如，0、1、2 或 3，取决于参与者分类）。以下列包含我正在查看的不同因素的数据（例如，特定位置的脑灰质衰减率）。然后我想浏览每一列，基本上做一个比较四组的方差分析。我认为循环是最好的方法，因为它是数百列中的相同测试。每行都是针对参与者的，包括他们的组分类和该特定区域的衰减率。这有帮助吗？
非常感谢！！这正是我要找的！

【解决方案2】：

我认为这至少可以帮助您入门。它使用 dplyr 和 broom 包。基本思想是将您想要的所有公式定义为字符，然后使用lapply() 将它们运行到lm()。

library(dplyr)
library(broom)

# Generate a vector of wanted formulas
forms <- c("mpg ~ cyl", "mpg ~ wt")

# Function to apply formula
lmit <- function(form){
  tidy(lm(as.formula(form), mtcars)) %>% 
    mutate(formula = form)
}

# Apply it and bind into a dataframe
results <- bind_rows(lapply(forms, lmit))

【讨论】：