【问题标题】:Linear Regression and storing results in data frame [duplicate]线性回归并将结果存储在数据框中[重复]
【发布时间】:2015-03-17 18:30:10
【问题描述】:

我正在对数据框中的一些变量进行线性回归。我希望能够通过分类变量对线性回归进行子集化,为每个分类变量运行线性回归,然后将 t-stats 存储在数据框中。如果可能的话,我想在没有循环的情况下这样做。

这是我正在尝试做的一个示例:

  a<-  c("a","a","a","a","a",
         "b","b","b","b","b",
         "c","c","c","c","c")     
  b<-  c(0.1,0.2,0.3,0.2,0.3,
         0.1,0.2,0.3,0.2,0.3,
         0.1,0.2,0.3,0.2,0.3)
  c<-  c(0.2,0.1,0.3,0.2,0.4,
         0.2,0.5,0.2,0.1,0.2,
         0.4,0.2,0.4,0.6,0.8)
      cbind(a,b,c)

我可以从运行以下线性回归开始并非常轻松地提取 t 统计量:

  summary(lm(b~c))$coefficients[2,3]

但是,我希望能够在 a 列是 a、b 或 c 时运行回归。然后我想将 t-stats 存储在一个如下所示的表中:

variable t-stat
a        0.9
b        2.4
c        1.1

希望这是有道理的。如果您有任何建议,请告诉我!

【问题讨论】:

标签: r linear-regression lm


【解决方案1】:

这是使用broom 包中的dplyrtidy() 的解决方案。 tidy() 将各种统计模型输出(例如lmglmanova 等)转换成整齐的数据框。

library(broom)
library(dplyr)

data <- data_frame(a, b, c)

data %>% 
  group_by(a) %>% 
  do(tidy(lm(b ~ c, data = .))) %>% 
  select(variable = a, t_stat = statistic) %>% 
  slice(2)

#   variable     t_stat
# 1        a  1.6124515
# 2        b -0.1369306
# 3        c  0.8000000  

或者同时提取截距和斜率项的 t 统计量:

data %>% 
  group_by(a) %>% 
  do(tidy(lm(b ~ c, data = .))) %>% 
  select(variable = a, term, t_stat = statistic)

#   variable        term     t_stat
# 1        a (Intercept)  1.2366939
# 2        a           c  1.6124515
# 3        b (Intercept)  2.6325081
# 4        b           c -0.1369306
# 5        c (Intercept)  1.4572335
# 6        c           c  0.8000000

【讨论】:

  • 选择命令对我不起作用(变量、术语和 t_stat 未使用)。
  • @Ben:代码对我来说仍然是完全可重现的,并且一切都按照上面的描述工作。
【解决方案2】:

这是对plyr 包和ddply() 的投票。

plyrFunc <- function(x){
  mod <- lm(b~c, data = x)
  return(summary(mod)$coefficients[2,3])
  }

tStats <- ddply(dF, .(a), plyrFunc)
tStats
  a         V1
1 a  1.6124515
2 b -0.1369306
3 c  0.6852483

【讨论】:

  • 谢谢,这非常简单明了。
  • 没问题。 plyr 包非常方便。我最近才知道它,但它非常有用。
【解决方案3】:

您可以使用nlme 包中的lmList 函数将lm 应用于数据子集:

# the data
df <- data.frame(a, b, c)

library(nlme)
res <- lmList(b ~ c | a, df, pool = FALSE)
coef(summary(res))

输出:

, , (Intercept)

   Estimate Std. Error  t value   Pr(>|t|)
a 0.1000000 0.08086075 1.236694 0.30418942
b 0.2304348 0.08753431 2.632508 0.07815663
c 0.1461538 0.10029542 1.457233 0.24110393

, , c

     Estimate Std. Error    t value  Pr(>|t|)
a  0.50000000  0.3100868  1.6124515 0.2052590
b -0.04347826  0.3175203 -0.1369306 0.8997586
c  0.15384615  0.1923077  0.8000000 0.4821990

如果你只想要 t 值,你可以使用这个命令:

coef(summary(res))[, "t value", -1]
#          a          b          c 
#  1.6124515 -0.1369306  0.8000000  

【讨论】:

    【解决方案4】:

    使用split 对数据进行子集化并通过lapply 进行循环

    dat <- data.frame(b,c)
    dat_split <- split(x = dat, f = a)
    res <- sapply(dat_split, function(x){
      summary(lm(b~c, data = x))$coefficients[2,3]
    })
    

    根据您的需要重塑结果:

    data.frame(variable = names(res), "t-stat" = res) 
    
      variable     t.stat
    a        a  1.6124515
    b        b -0.1369306
    c        c  0.8000000
    

    【讨论】:

      【解决方案5】:

      你可以这样做:

      a<-  c("a","a","a","a","a",
             "b","b","b","b","b",
             "c","c","c","c","c")     
      b<-  c(0.1,0.2,0.3,0.2,0.3,
             0.1,0.2,0.3,0.2,0.3,
             0.1,0.2,0.3,0.2,0.3)
      c<-  c(0.2,0.1,0.3,0.2,0.4,
             0.2,0.5,0.2,0.1,0.2,
             0.4,0.2,0.4,0.6,0.8)
      df <- data.frame(a,b,c)
      
      
      t.stats <- t(data.frame(lapply(c('a','b','c'), 
                   function(x) summary(lm(b~c,data=df[df$a==x,]))$coefficients[2,3])))
      colnames(t.stats) <- 't-stat'
      rownames(t.stats) <- c('a','b','c')
      

      输出:

      > t.stats
            t-stat
      a  1.6124515
      b -0.1369306
      c  0.8000000
      

      除非我弄错了,否则您在输出中给出的值不是正确的。

      或者:

      t.stats <- data.frame(t.stats)
      t.stats$variable <- rownames(t.stats)
      
      > t.stats[,c(2,1)]
        variable     t.stat
      a        a  1.6124515
      b        b -0.1369306
      c        c  0.8000000
      

      如果你想要一个 data.frame 和一个单独的列。

      【讨论】:

        猜你喜欢
        • 2020-03-15
        • 1970-01-01
        • 2019-02-07
        • 2019-10-23
        • 2019-09-18
        • 1970-01-01
        • 1970-01-01
        • 1970-01-01
        相关资源
        最近更新 更多