R中由群组运行的数百个线性回归[关闭]

3 回答

• 1

  我认为最简单的答案仍然缺失 . 您可以使用嵌套和映射的组合 . 我'll show you how it works for linear regression. I think you'能够将相同的原理应用于 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)
  )
  

  如果我们想为每个组 Build 一个回归模型，研究IQ得分和时间点之间的关系，我们只需要五行代码 .

  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) 在 tibble 中为每个组创建 tibbles ，包含 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将汇总统计信息外包给负责该模型的特定软件包开发人员 .

  回复于 2024-04-29T15:02:06+08:00
• 2

  虽然AndS给出的代码可以工作，但每个组运行lm函数4次，这使得效率有点低 . 您可以使用以下内容 . 我试图将其分解为更简单的步骤：

  假设您的数据框（df）有三个变量：“Group”，“Dep”，“Indep”：

  #Getting the unique list of groups
  groups <- unique(df$Group)
  
  #Creating a model summary list to combine the model summary of each model
  model_summaries = list()
  
  #Running the models
  for(i in 1:length(groups)){
    model <- lm(Dep ~ Indep, df[df$Group == Groups[i], c("Dep", "Indep")])
    model_summaries[i] <- summary(model)
  }
  

  在每个模型摘要中，您有以下元素RSQ，系数（包含p值和截距）

  如果这有帮助，请告诉我 .

  回复于 2024-04-29T15:02:06+08:00
• 0

  如果我正确理解您的问题，您希望对多个组进行多次回归 . 以下是如何使用mtcars数据执行此操作的示例 .

  library(dplyr)
  mtcars %>% group_by(cyl) %>% 
      summarise_at(vars(disp:wt), funs(
          r.sqr = summary(lm(mpg~.))$r.squared,
          intercept = summary(lm(mpg~.))$coefficients[[1]],
          slope = summary(lm(mpg~.))$coefficients[[2]],
          p.value = summary(lm(mpg~.))$coefficients[[8]]
      ))
  

  这将针对每个变量运行每组的回归并提取您要求的信息 . 如果您的公式始终相同，则可以简化如下 .

  mtcars %>% group_by(cyl) %>% 
      summarise(
          r.sqr = summary(lm(mpg~wt))$r.squared,
          intercept = summary(lm(mpg~wt))$coefficients[[1]],
          slope = summary(lm(mpg~wt))$coefficients[[2]],
          p.value = summary(lm(mpg~wt))$coefficients[[8]]
      )
  

  这实际上是运行回归4次（每个感兴趣的值一次） . 如果您的实际数据需要太长时间，您可以尝试这样做：

  df <- mtcars %>% group_by(cyl) %>% summarise(model = list(summary(lm(mpg~wt))))
  

  它只是每组运行一次模型，然后提取出你想要的信息 . 问题是以这种方式提取值可能是一种痛苦

  df$model[[1]]$coefficients[[1]]
  [1] 39.5712
  

  回复于 2024-04-29T15:02:06+08:00