dplyr 版本的数据框分组然后在每个组上创建回归模型

Question

谁能建议dplyr 回答以下问题？ Split data.frame by country, and create linear regression model on each subset

为了完整起见，链接中的问题和答案包含在下面。

问题

作为参考，这是乔希的问题：

我有一个来自世界银行的数据的 data.frame，看起来像这样；

  country date BirthRate     US.   
4   Aruba 2011    10.584 25354.8
5   Aruba 2010    10.804 24289.1
6   Aruba 2009    11.060 24639.9
7   Aruba 2008    11.346 27549.3
8   Aruba 2007    11.653 25921.3
9   Aruba 2006    11.977 24015.4

在这个数据框中总共有 70 个国家的子集，我想对其进行线性回归。如果我使用以下内容，我会为一个国家/地区获得不错的 lm；

andora = subset(high.sub, country == "Andorra")

andora.lm = lm(BirthRate~US., data = andora)

anova(andora.lm)
summary(andora.lm)

但是当我尝试在 for 循环中使用相同类型的代码时，我会收到一个错误，我将在代码下方打印该错误；

high.sub = subset(highInc, date > 1999 & date < 2012)
high.sub <- na.omit(high.sub)
highnames <- unique(high.sub$country)

for (i in highnames) {
  linmod <- lm(BirthRate~US., data = high.sub, subset = (country == "[i]"))  
}

#Error message:
Error in lm.fit(x, y, offset = offset, singular.ok = singular.ok, ...) : 
  0 (non-NA) cases

如果我可以让这个循环运行，我理想情况下希望将每个模型的系数甚至更好的 r 平方值附加到一个空的 data.frame。任何帮助将不胜感激。

回答

作为参考，这里是 jlhoward 的答案（包含 BondedDust 的评论），它使用了这个优秀问题中的 *apply 函数： R Grouping functions: sapply vs. lapply vs. apply. vs. tapply vs. by vs. aggregate

models <- sapply(unique(as.character(df$country)),
                 function(cntry)lm(BirthRate~US.,df,subset=(country==cntry)),
                 simplify=FALSE,USE.NAMES=TRUE)

# to summarize all the models
lapply(models,summary)
# to run anova on all the models
lapply(models,anova)

#This produces a named list of models, so you could extract the model for Aruba as:
models[["Aruba"]]

Answer 1

看看nlme包的lmList函数：

library(nlme)
lmList(BirthRate ~ US. | country, df)

这里，| country 用于为每个国家/地区创建回归。

Answer 2

目前尚无法从dplyr 返回列表。如果您只需要截距和斜率 @jazzurro 的答案就是方法，但是如果您需要整个模型，则需要执行类似

的操作

library(dplyr)
models <- df %>% group_by(country) %>% do(mod = lm(BirthRate ~ US., data = .))

如果你想对每个拟合模型进行方差分析，你可以使用rowwise

models %>% rowwise %>% do(anova(.$mod))

但结果再次被强制转换为数据帧，与lapply(models$mod, anova) 的操作不完全相同。

目前（即直到下一个版本的dplyr）如果您需要将整个结果存储在一个列表中，您可以使用plyr 中的dlply，例如plyr::dlply(df, "country", function(d) anova(lm(BirthRate ~ US., data = d)))，或者当然如果您不必绝对使用dplyr，您可以选择@SvenHohenstein 的答案，无论如何这看起来是一种更好的方法。

Answer 3

这是使用dplyr 的一种方式。我使用您的数据创建了Brazil。所以你对这两个国家有相同的结果。您会在最终的数据框中看到截距和斜率。

数据和代码

structure(list(country = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 
2L, 2L, 2L, 2L, 2L, 2L), .Label = c("Aruba", "Brazil"), class = "factor"), 
date = c(2011L, 2010L, 2009L, 2008L, 2007L, 2006L, 2011L, 
2010L, 2009L, 2008L, 2007L, 2006L), BirthRate = c(10.584, 
10.804, 11.06, 11.346, 11.653, 11.977, 10.584, 10.804, 11.06, 
11.346, 11.653, 11.977), US. = c(25354.8, 24289.1, 24639.9, 
27549.3, 25921.3, 24015.4, 25354.8, 24289.1, 24639.9, 27549.3, 
25921.3, 24015.4)), .Names = c("country", "date", "BirthRate", 
"US."), class = "data.frame", row.names = c("4", "5", "6", "7", 
"8", "9", "10", "11", "12", "13", "14", "15"))

#   country date BirthRate     US.
#4    Aruba 2011    10.584 25354.8
#5    Aruba 2010    10.804 24289.1
#6    Aruba 2009    11.060 24639.9
#7    Aruba 2008    11.346 27549.3
#8    Aruba 2007    11.653 25921.3
#9    Aruba 2006    11.977 24015.4
#10  Brazil 2011    10.584 25354.8
#11  Brazil 2010    10.804 24289.1
#12  Brazil 2009    11.060 24639.9
#13  Brazil 2008    11.346 27549.3
#14  Brazil 2007    11.653 25921.3
#15  Brazil 2006    11.977 24015.4

group_by(mydf, country) %>%
    do({model = lm(BirthRate ~ US., data = .);
       data.frame(int = coef(model)[1], slope = coef(model)[2])})

#  country      int        slope
#1   Aruba 11.02503 8.393295e-06
#2  Brazil 11.02503 8.393295e-06