r - R 中按组运行的数百个线性回归

Question

已关闭。此问题需要 debugging details 。目前不接受答案。

---

编辑问题以包含 desired behavior, a specific problem or error, and the shortest code necessary to reproduce the problem 。这将有助于其他人回答问题。

已关闭 5 年前。

Improve this question

我有一个包含 3,000 多行和 10 多个变量的表。我正在尝试使用一个变量作为预测变量，使用另一个变量作为 300 个不同组的响应来运行线性回归。我需要每个回归的斜率、p 值和 r 平方。单独进行每个回归并记录汇总变量将花费数小时甚至数天的时间。

我使用了以下包来获取每个组的截距和斜率，但我不知道如何获取每个组相应的 p 值和 r 平方:

library(lme4)
groupreg<-lmList(logpop ~ avgp | id, data=data)
groupreg

我实现了下面的列表示例，其中“Adams #”是 id 值。 NA 的存在是因为并非所有组都有多个点来绘制和比较:

Coefficients:
                (Intercept)          avgp
Adams 6           4.0073332            NA
Adams 7           6.5177389 -7.342443e+00
Adams 8           4.7449321            NA
Adams 9                  NA            NA

但是，该表不包含任何显着性统计数据。我仍然需要 p 值和 r 平方统计量。如果有一个代码可以一次性完成所有组值的所有操作，或者有一个代码可以只提取剩余的值，那将会很有帮助。

还有没有办法对所有组的斜率输出求幂？我的结果是对数转换的。

谢谢大家!!

Answer 1

我认为最简单的答案仍然缺失。您可以结合使用嵌套和映射。我将向您展示它如何用于线性回归。我认为您可以将相同的原理应用于 lme4 的模型包。

让我们创建一个玩具数据集，我们在其中测量了三个不同群体在两个不同时间点的 IQ 分数。

library(tidyverse)
library(broom)

df <- tibble(
  id = seq_len(90),
  IQ = rnorm(90, 100, 15),
  group = rep(c("A", "B", "C"), each = 30),
  time = rep(c("T1", "T2"), 45)
)

如果我们想为每个组建立一个回归模型，研究智商分数和时间点之间的关系，我们只需要五行代码。

df %>% 
  nest(-group) %>% 
  mutate(fit = map(data, ~ lm(IQ ~ time, data = .)),
         results = map(fit, glance)) %>% 
  unnest(results) %>% 
  select(group, r.squared, p.value) 

将会返回

 # A tibble: 3 x 3
  group r.squared p.value
  <chr>     <dbl>   <dbl>
1 A       0.0141    0.532
2 B       0.0681    0.164
3 C       0.00432   0.730

哪里nest(-group)创建tibbles在您的tibble内对于每个组，包含相应的变量 id , IQ和time 。然后添加一个新列 fit与 mutate()在其中为每个组应用回归模型以及包含结果的新列，我们unnest()不久后访问值 glance()正确返回。最后一步我们select()兴趣的三个值。

要获取坡度，您需要调用 tidy()此外。也许可以以某种方式缩短代码，但一种解决方案是

df %>% 
  nest(-group) %>% 
  mutate(fit = map(data, ~ lm(IQ ~ time, data = .)),
         results1 = map(fit, glance),
         results2 = map(fit, tidy)) %>% 
  unnest(results1) %>% 
  unnest(results2) %>% 
  select(group, term, estimate, r.squared, p.value) %>% 
  mutate(estimate = exp(estimate)) 

要对斜率求幂，您只需添加另一个 mutate()陈述。终于返回了

# A tibble: 6 x 5
  group term        estimate r.squared p.value
  <chr> <chr>          <dbl>     <dbl>   <dbl>
1 A     (Intercept) 3.34e+46   0.0141    0.532
2 A     timeT2      3.31e- 2   0.0141    0.532
3 B     (Intercept) 1.17e+47   0.0681    0.164
4 B     timeT2      1.34e- 3   0.0681    0.164
5 C     (Intercept) 8.68e+43   0.00432   0.730
6 C     timeT2      1.25e- 1   0.00432   0.730

请注意，估计值已取幂。如果没有求幂，您可以使用 base R 仔细检查斜率和 p 值。调用

summary(lm(IQ ~ time, data = filter(df, group == "A")))

如果您使用更复杂的模型 ( lme4 )，有一个名为 lmerTest 的包它为 lme4 提供包装函数它返回 p 值(至少对于我已经使用过的混合模型)。

关于使用 glance() 的警告对于 lme4模型应该被讲出来，因为 broom 的维护者包，会尝试 new concept他们将汇总统计数据外包给负责模型的特定包开发人员。