R中的子集回归[重复]

Posted 2023-03-12

【中文标题】R中的子集回归[重复]

【英文标题】：Regression by subset in R [duplicate]

【发布时间】：2013-08-04 05:38:23

【问题描述】：

我是 R 新手，正在尝试对单个文件中的多个数据子集（“案例”）运行线性回归。我有 50 种不同的案例，所以我不想运行 50 种不同的回归……很高兴将其自动化。我已经找到并尝试了ddply 方法，但出于某种原因，这会为每种情况返回相同的系数。我使用的代码如下：

ddply(MyData, "Case", function(x) coefficients(lm(Y~X1+X2+X3, MyData)))

我得到的结果再次是每个“案例”的相同系数。关于如何改进我的代码以便回归为每种情况运行一次并为每种情况提供唯一系数的任何想法？

【参考方案1】：

虽然你可以用 (d)plyr 走得更远，但我通过一个简单的 for 循环获得了最大的灵活性（因为我对 do() 不太有信心，因为我想要的不仅仅是系数的输出()，例如）

让我们从加载数据和加载一些包开始：

library(dplyr)
library(broom) # get lm output (coefficients) as a dataframe
library(reshape2) # Melt / DCAST for swapping between wide / long format
data(warpbreaks)  

接下来，创建一个我们将要迭代的组列表。

GROUPS <- unique(warpbreaks$tension) # Specify groups, one unique model per group.

在这种情况下，回答您问题的代码可能是这样的：

for (i in 1:length(GROUPS))
  CURRENT_GROUP <- GROUPS[i]
  df <- filter(warpbreaks, tension==CURRENT_GROUP) # subset the dataframe
  fit <- lm(breaks ~ wool, data = df) # Build a model
  coeff <- tidy(fit) # Get a pretty data frame of the coefficients & p values
  coeff <- coeff[,c(1,2,5)] # Extract P.Value & Estimate

  # Rename (intercept) to INT for pretty column names
  coeff[coeff$term=="(Intercept)", ]$term <- "INT"

  # Make it into wide format with reshape2 package.
  coeff <- coeff %>% melt(id.vars=c("term"))

  # Defactor the resulting data.frame
  coeff <- mutate(coeff,
                  variable=as.character(variable),
                  term=as.character(term))

  # Rename for prettier column names later
  coeff[coeff$variable=="estimate", ]$variable <- "Beta"
  coeff[coeff$variable=="p.value", ]$variable <- "P"

  coeff <- dcast(coeff, .~term+variable)[,-1]

  rsquared <- summary(fit)$r.squared

  # Create a df out of what we just did.
  row <- cbind(
    data.frame(
      group=CURRENT_GROUP,
      rsquared=rsquared),
    coeff
  )

  # If first iteration, create data.frame -- otherwise: rowbind
  if (i==1)
    RESULT_ROW = row
  
  else
    RESULT_ROW = rbind(RESULT_ROW, row)
   # End if.
 # End for loop

在编写内部for循环代码时，只需通过模拟循环来测试一些东西（首先运行i

RESULT_ROW  
##   group   rsquared INT_Beta        INT_P woolB_Beta    woolB_P
## 1     L 0.26107601 44.55556 9.010046e-08 -16.333333 0.03023435
## 2     M 0.07263228 24.00000 5.991177e-07   4.777778 0.27947884
## 3     H 0.12666292 24.55556 9.234773e-08  -5.777778 0.14719571

我并不是说这比上面描述的答案更好，我只是喜欢这样一个事实，即这让我可以灵活地从任何模型中提取任何东西（不仅仅是具有漂亮的 coefficients() 函数的 lm 模型）。

【讨论】：

对于具有许多组的大型数据帧，for 循环方法非常慢，因为 a) 您在每个组上运行单独的过滤器（尝试 library(hflights); for (g in unique(hflights$TailNum)) filter(hflights, TailNum == g): ~ 35 秒，它甚至什么都不做，只是过滤），b）一系列连续的rbinds 很慢，因为它每次都必须将数据帧的内容复制到一个新的。因此，这不是一个好习惯。

其次，没有必要。为什么不直接编写自己的函数并使用do？这可以使用func <- function(df) 来完成，然后是fit <- lm(... 行到row <- cbind( 之间的所有内容，然后是return(row) 。然后你可以做warpbreaks %>% group_by(tension) %>% do(func(.))

我绝对可以看出它根本没有效率。您能否提供一个示例，说明如何在没有 for 循环的情况下执行与我刚刚发布的内容相同的操作？非常感谢学习经验。

当然，请看我的第二条评论 :)

@David Robinson 您的建议如何与您建议包含在row <- cbind(data.frame( group=CURRENT_GROUP,rsquared=rsquared), coeff) 中的代码的group=CURRENT_GROUP 参数一起使用

【参考方案2】：

你也可以只使用基本函数：

# load data
data(warpbreaks)

# fit linear model to each subset
fits <- by(warpbreaks, warpbreaks[,"tension"],
           function(x) lm(breaks ~ wool, data = x))

# Combine coefficients from each model
do.call("rbind", lapply(fits, coef))

【讨论】：

谢谢，迈克尔。这也做得很好。现在我试图将每个子集的系数附加到原始表中。 （显然，这些在每个案例中都是重复的）有什么提示吗？

如果附加的意思是在原始数据框中添加一列，给出每个子集的系数，您可以使用基础 R 中的 merge() 函数或 plyr 包中的 join() 函数。

正确，迈克尔。我正在寻找使用模型系数向原始数据框添加一列。这两个功能都运行良好。我正在慢慢学习这门语言的来龙去脉！

【参考方案3】：

ddply 将 data.frames（来自分割输入 data.frame）传递给函数。你可能想要这个：

ddply(MyData, "Case", function(df) coefficients(lm(Y~X1+X2+X3, data=df)))

（未测试，因为您没有提供可重现的示例。）

您将每个组的整个输入 data.frame 传递给 lm..

【讨论】：

太棒了！它工作得很好。现在……刚刚发生了什么？！ df 在 function() 内部时是什么意思？我猜它只代表每个案例对应的数据？

我刚刚将函数参数命名为df 以突出它是一个data.frame。如果您愿意，可以将其称为 x。

我认为您应该在stats.stackexchange.com/questions/66390/… 的相同问题中感谢@Roland，而不仅仅是“找到解决方案”，并链接到此答案。

谢谢你让我诚实，Henrik。归功于交叉验证帖子中的 Roland。