lm()：遍历多个线性模型，导出 F 统计量的 p 值

Question

我有一个大数据集，我需要运行一个线性模型比较组。我需要使用线性模型找到组比较的 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 值。帮助将不胜感激！谢谢！

Answer 1

我认为这至少可以让您入门。它使用 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))

Answer 2

(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 同意打印的摘要。

你的循环

我建议你在computation/update期间使用向量而不是数据框。

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) returns是一个向量，因此var[i] <- colnames(data)会触发错误。
i~group 看起来很奇怪。你的 data 里有 i 吗？

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

跟进 (2016-06-19)

Thank you. This has been extremely helpful. I don't have "i" in my data, but I was hoping that I could use that to represent the different column names, so that it goes through all of them. Is there a way to assign column names numbers so that this would work?

是的，但我需要知道每一列的内容。

Column 1 has a group number. The following columns have data for different factors I am looking at.

好的，所以我认为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().

相比，循环开销非常小

lm()：遍历多个线性模型，导出 F 统计量的 p 值

lm(): loop through multiple linear models exporting p-value of F-statistic

loops

regression

r

linear-regression

lm