r语言 - run lm()函数在数据框中的2个不同变量的子集上

Question

我想创建一个函数（最好使用apply family函数，该功能在数据框架中从两个单独的变量中运行lm函数。参考，我的数据是当前构造如下（这是一个虚拟数据框架，通常代表我正在使用的真实数据框架）：

District = c(rep("A", times = 25),
             rep("B", times = 25),
             rep("C", times = 25),
             rep("D", times = 25))
Year = c(1971:1975, 1971:1975, 1971:1975, 1971:1975,
         1971:1975, 1971:1975, 1971:1975, 1971:1975,
         1971:1975, 1971:1975, 1971:1975, 1971:1975,
         1971:1975, 1971:1975, 1971:1975, 1971:1975,
         1971:1975, 1971:1975, 1971:1975, 1971:1975)
Crop = c(rep("Wheat",  times = 5),
         rep("Maize",  times = 5),
         rep("Rice",   times = 5),
         rep("Barley", times = 5))
set.seed(100)
Yield = rnorm(100, 2, 0.5)
df <- data.frame(District, Year, Crop, Yield)

我需要生成一个线性模型（lm），该模型（CC_3）预测Yield作为每个Crop的Year的函数。因此，我需要每个地区的小麦模型，每个地区的大麦等等。我需要尽可能地自动化这一点，因为我的实际数据将导致大约7000个线性模型。

我发现詹姆斯·邦德（James Bond）的这个答案非常有帮助，这是一个出发点。它非常优雅，使用lapply，如果可能的话，我也想做的（7000个线性型号=非常慢）。但是，它只在数据集中子集中一个列，而我需要为每个模型子集两个变量。这是我当前（非工作）代码修改以在上述虚拟数据集上运行的代码：

df$Dist <- as.factor(df$Dist)
df$Crop <- as.factor(df$Crop)
for (i in 1:length(levels(df$Crop))) {
  x <- levels(df$Crop)[i]
  dat <- df[df$Crop == x, ]
  out <- lapply(levels(dat$Dist), function(z) {
            data.frame(District = z, 
                       Slope = lm(Yield ~ Year, data = dat[dat$Dist == z, ])$coef[2], 
                       Crop = x,
                       row.names=NULL)
  })
}
do.call(rbind ,out)

不幸的是，运行上述代码仅生成数据集（小麦）中第一个Crop级别的模型。请参阅下面的输出：

  District       Slope  Crop
1        A  0.03125866 Wheat
2        B -0.08108222 Wheat
3        C  0.17172314 Wheat
4        D -0.11278486 Wheat

能够在Crop和District变量上循环的任何帮助都将不胜感激。好像我很亲密，但似乎错过了有关for循环的根本性。

如果可以将2个参数传递给lapply函数并避免for循环，那将是惊人的。预先感谢！

Answer 1

使用 dplyr 的一个选项：

df_lm <- df %>%
  group_by(District,Crop) %>%
  do(mod = lm(Yield ~ Year,data = .))
df_coef <- df_lm %>%
  do(data.frame(
    District = .$District,
    Crop = .$Crop,
    var = names(coef(.$mod)),
    coef(summary(.$mod)))
    )
> df_coef
Source: local data frame [32 x 7]
Groups: <by row>
# A tibble: 32 × 7
   District   Crop         var      Estimate   Std..Error    t.value   Pr...t..
*    <fctr> <fctr>      <fctr>         <dbl>        <dbl>      <dbl>      <dbl>
1         A Barley (Intercept) -407.66953514 378.49788671 -1.0770722 0.36034462
2         A Barley        Year    0.20771336   0.19183872  1.0827499 0.35818046
3         A  Maize (Intercept)  159.81133118 212.90233600  0.7506321 0.50738515
4         A  Maize        Year   -0.08002266   0.10790790 -0.7415830 0.51211787
5         A   Rice (Intercept)  -68.01125454 117.60578244 -0.5782986 0.60361684
6         A   Rice        Year    0.03552764   0.05960758  0.5960255 0.59313364
7         A  Wheat (Intercept)  -59.61828825 134.67806297 -0.4426726 0.66972726
8         A  Wheat        Year    0.03125866   0.06826053  0.4579317 0.65918309
9         B Barley (Intercept) -319.99755207  57.14553545 -5.5996947 0.01125215
10        B Barley        Year    0.16332436   0.02896377  5.6389189 0.01103509
# ... with 22 more rows

要看的另一件事是 nlme 中的lmList函数。