如何在R中循环/重复执行线性回归

16

我已经找到了如何在R中使用4个变量创建表格,用于多元线性回归。每次回归的因变量(Lung)都从一个包含22,000列的csv表格的一列中获取。其中一个自变量(Blood)则从类似表格的相应列中获取。

每列代表特定基因的水平,这就是为什么有这么多列的原因。此外还有两个附加变量(年龄性别)。当我输入线性回归方程时,我使用lm(Lung[,1] ~ Blood[,1] + Age + Gender),它适用于一个基因。

我想找到一种输入此方程并让R计算所有其余LungBlood列的方法,并希望将系数输出到一个表格中。

任何帮助都将不胜感激!

5个回答

26
你想要运行22,000个线性回归并提取系数?从编程角度来看,这很简单。
set.seed(1)

# number of columns in the Lung and Blood data.frames. 22,000 for you?
n <- 5 

# dummy data
obs <- 50 # observations
Lung <- data.frame(matrix(rnorm(obs*n), ncol=n))
Blood <- data.frame(matrix(rnorm(obs*n), ncol=n))
Age <- sample(20:80, obs)
Gender  <- factor(rbinom(obs, 1, .5))

# run n regressions
my_lms <- lapply(1:n, function(x) lm(Lung[,x] ~ Blood[,x] + Age + Gender))

# extract just coefficients
sapply(my_lms, coef)

# if you need more info, get full summary call. now you can get whatever, like:
summaries <- lapply(my_lms, summary)
# ...coefficents with p values:
lapply(summaries, function(x) x$coefficients[, c(1,4)])
# ...or r-squared values
sapply(summaries, function(x) c(r_sq = x$r.squared, 
                                adj_r_sq = x$adj.r.squared))

模型被存储在一个列表中,第三个模型(使用DV Lung[, 3]和IVs Blood[,3] + Age + Gender)存储在my_lms[[3]]中等等。您可以对该列表使用apply函数进行摘要,从中提取您想要的数字。


1
这太棒了!我能问一下你是如何编写年龄和性别方程的吗?我之前做法是制作一个包含年龄/性别的表格,并设置AgeGender <- read.csv("AgeGender.csv", header=TRUE)。然后我设置Age <- AgeGender[,1]和Gender <- AgeGender[,2]。性别也被编码为1表示男性,0表示女性,那么我应该将你的性别函数中的0.5更改为0吗? - user4438232
尝试获取R平方值时,我尝试使用summary(lm(Lung[,x] ~ Blood[,x] + Age + Gender))$r.squared,但它一直显示“维度数量不正确”。有什么建议吗? - user4438232
@JHall1020 我知道已经过了很长时间,但你有找到计算回归系数和R平方的p值的方法吗? - Guu
嗨,我正在浏览这个页面,因为我遇到了与OP类似的问题。为什么你需要在这里设置种子?随机数是什么? - BonnieKlein
只需使虚拟数据部分可重现(即 rnorm() 调用)。 - arvi1000
显示剩余4条评论

3
这个问题似乎是关于如何在循环内部修改公式并调用回归函数的。以下是您可以使用钻石数据集实现的方法:
attach(ggplot2::diamonds)
strCols = names(ggplot2::diamonds)

formula <- list(); model <- list()
for (i in 1:1) {
  formula[[i]] = paste0(strCols[7], " ~ ", strCols[7+i])
  model[[i]] = glm(formula[[i]]) 

  #then you can plot or do anything else with the result ...
  png(filename = sprintf("diamonds_price=glm(%s).png", strCols[7+i]))
  par(mfrow = c(2, 2))      
  plot(model[[i]])
  dev.off()
  }

2

无论是否合理,为了使循环至少在某种程度上起作用,您需要:

y<- c(1,5,6,2,5,10) # response 
x1<- c(2,12,8,1,16,17) # predictor 
x2<- c(2,14,5,1,17,17) 
predictorlist<- list("x1","x2") 
for (i in predictorlist){ 
  model <- lm(paste("y ~", i[[1]]), data=df) 
  print(summary(model)) 
} 

粘贴功能将解决该问题。

变量(y / x1 / x2)需要作为数据框存储,以使示例正常工作,我认为:df <- data.frame( y = c(1,5,6,2,5,10), x1 = c(2,12,8,1,16,17), x2 = c(2,14,5,1,17,17) ) - Tfsnuff

2

一个tidyverse的补充 - 使用map()

另一种方法 - 使用purrr包中的map2()

library(purrr)

xs <- anscombe[,1:3] # Select variables of interest
ys <- anscombe[,5:7]

map2_df(ys, xs,
        function(i,j){
          m <- lm(i ~j + x4 , data = anscombe)
          coef(m)
        })

输出是一个数据框(tibble),包含所有系数:
  `(Intercept)`     j      x4
1          4.33 0.451 -0.0987
2          6.42 0.373 -0.253 
3          2.30 0.526  0.0518

如果有更多的变量在改变,可以使用pmap()函数来实现。

0

以下方法适用于多元模型,其中有多个结果和预测变量。我将使用一些示例数据来说明这个想法以及它是如何工作的。

df <- data.frame(y1=sample(1:5, size=50, replace=TRUE),
                 y2=sample(1:5, size=50, replace=TRUE),
                 x1=sample(1:5, size=50, replace=TRUE),
                 x2=sample(1:5, size=50, replace=TRUE),
                 x3=sample(1:5, size=50, replace=TRUE),
                 x4=sample(1L:2L, size=50, replace=TRUE))
df

该函数需要为dv指定一个命名参数,但使用省略号表示您可以拥有任意数量的预测变量。在函数内部,您使用deparse()substitute()处理预测变量,并将它们与dv一起传递给reformulate()函数。我在函数内部包含了dv=dv,以便可以看到模型输出与哪个dv相关联。
# custom lm function
lm_func <- function(dv, ...){
  x = sapply(substitute(...()), deparse)
  f = reformulate(termlabels=x, response=dv)
  model = eval(lm(f, data=df))
  list(dv=dv, model_summary=summary(model))
}

在下一步中,我们从目标数据框中选择dvs并对其进行命名。
# select the dvs and set names
dvs <- names(df)[1:2]
dvs <- purrr::set_names(dvs)

最后,对dvs运行一个循环并存储结果。
# run a for loop and save the output for each loop
lm_out = list()
for (i in 1:length(dvs)){
  lm_out[[i]] = (lm_func(dvs[i], x1, x2))
  }
lm_out

注意:在lm_func中可以做更多的事情;例如,提取模型摘要的哪些部分。

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接