使用 lm 在 R 中进行线性回归：函数中的不同汇总输出答案

【问题标题】：Linear regression in R using lm: Different summary output in a function使用 lm 在 R 中进行线性回归：函数中的不同汇总输出
【发布时间】：2014-12-14 22:38:21
【问题描述】：

关于使用 lm 函数在 R 中进行线性回归的简短问题。我注意到在函数中使用 summary 命令时输出不同。

当我进入时：

model1 <- lm (PostVal_Ave ~ Int)
summary(model1)

控制台返回如下：

Call:
lm(formula = PostVal_Ave ~ Int)

Residuals:
    Min      1Q  Median      3Q     Max 
-3.9871 -0.8897  0.4853  1.0129  1.5129 

Coefficients:
        Estimate Std. Error t value Pr(>|t|)
(Intercept)   5.4871     0.1426  38.491   <2e-16
Int           0.2776     0.1988   1.396    0.164

(Intercept) ***
Int            
---
Signif. codes:  
0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.322 on 175 degrees of freedom
(35 observations deleted due to missingness)
Multiple R-squared:  0.01102,   Adjusted R-squared:  0.005366 
F-statistic: 1.949 on 1 and 175 DF,  p-value: 0.1644

但是，当编写一个函数以便能够为多个模型生成输出并能够为多个因变量生成结果时，我输入：

allModels <- function(x){
    model2 <- lm (x ~ Int)
    model2.1 <- lm (x ~ Int + cPreEff)
    model2.2 <- lm (x ~ Int + cPreEff + Gender + Grade)
    return(c(summary(model1), summary(model1.1), summary(model1.2)))}

与模型 1 的输出相比，我得到了相同的输出，但是这三个模型（模型 2、模型 2.1 和模型 2.2）有很多额外的输出。具体来说，输出包含三个模型中每一个的每个案例的残差，以及关于每个缺少数据的案例的信息。建议将不胜感激。谢谢。

【问题讨论】：

那么这里的问题到底是什么？请注意，lm() 返回一个“lm”类的对象，而该对象上的summary() 会生成一个“summary.lm”对象。有自定义的 print.lm() 和 print.summary.lm() 对象。因此，打印到控制台的内容可能与对象本身的内容不同。当您手动连接 (c()) 两个 summary.lm 对象时，整个列表将失去其类别。您可能希望返回对象列表 return(list(summary(model1), summary(model1.1), summary(model1.2)))
谢谢 - 返回一个对象列表就可以了。

标签： r linear-regression lm

【解决方案1】：

请注意，lm() 返回一个“lm”类的对象，而该对象上的summary() 会生成一个“summary.lm”对象。有自定义的 print.lm() 和 print.summary.lm() 对象。因此，打印到控制台的内容可能与对象本身的内容不同。

当您手动连接 (c()) 两个 summary.lm 对象时，您会创建一个聚合列表并丢失正确的类。您可能希望返回一个对象列表

return(list(summary(model1), summary(model1.1), summary(model1.2)))

【讨论】：

【解决方案2】：

另一种方法是使用 David Robinson 的 broom package，它将拟合输出转换为数据框。对于拟合的总结（没有系数），使用glance：

library(broom)

fitIris <- function() {
  fit1 <- glance(lm(Sepal.Length ~ Species, data = iris))
  fit2 <- glance(lm(Petal.Length ~ Species, data = iris))
  fit3 <- glance(lm(Sepal.Width ~ Species, data = iris))
  out <- rbind(fit1, fit2, fit3)
  cbind(Fit = c("fit1", "fit2", "fit3"), out)
}

fitIris()

#    Fit r.squared adj.r.squared     sigma  statistic      p.value df    logLik      AIC      BIC deviance df.residual
# 1 fit1 0.6187057     0.6135181 0.5147894  119.26450 1.669669e-31  3 -111.7260 231.4520 243.4945  38.9562         147
# 2 fit2 0.9413717     0.9405741 0.4303345 1180.16118 2.856777e-91  3  -84.8467 177.6934 189.7359  27.2226         147
# 3 fit3 0.4007828     0.3926302 0.3396877   49.16004 4.492017e-17  3  -49.3663 106.7326 118.7751  16.9620         147

有关拟合的更多信息，您应该使用tidy，但行数会有所不同，应谨慎创建结束data.frame。

fitIris2 <- function() {
  fit1 <- tidy(lm(Sepal.Length ~ Species, data = iris))
  fit2 <- tidy(lm(Petal.Length ~ Species, data = iris))
  fit3 <- tidy(lm(Sepal.Width  ~ Species, data = iris))
  out <- rbind(fit1, fit2, fit3)
  cbind(Fit = rep(c("fit1", "fit2", "fit3"), each=3), out)
}

fitIris2()

#    Fit              term estimate  std.error statistic       p.value
# 1 fit1       (Intercept)    5.006 0.07280222 68.761639 1.134286e-113
# 2 fit1 Speciesversicolor    0.930 0.10295789  9.032819  8.770194e-16
# 3 fit1  Speciesvirginica    1.582 0.10295789 15.365506  2.214821e-32
# 4 fit2       (Intercept)    1.462 0.06085848 24.022945  9.303052e-53
# 5 fit2 Speciesversicolor    2.798 0.08606689 32.509597  5.254587e-69
# 6 fit2  Speciesvirginica    4.090 0.08606689 47.521176  4.106139e-91
# 7 fit3       (Intercept)    3.428 0.04803910 71.358540 5.707614e-116
# 8 fit3 Speciesversicolor   -0.658 0.06793755 -9.685366  1.832489e-17
# 9 fit3  Speciesvirginica   -0.454 0.06793755 -6.682608  4.538957e-10

我很长一段时间都不是这个用户，所以可能有更好的方法可以将这些数据框连接在一起。

【讨论】：