CODEWITHSUNDEEP
rlinear-regressionlm
Joshua Rosenberg
Joshua Rosenberg

Reputation: 4236

Linear regression in R using lm: Different summary output in a function

Brief question regarding linear regression in R using the lm function. I noticed that the output is different when using the summary command as part of a function.

When I enter:

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

The follow is returned in the console:

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

But, when writing a function in order to be able to produce output for multiple models and to be able to produce results for multiple dependent variables, I enter:

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)))}

And I get the same output compared to the output for model 1, but with a lot of additional output for these three models (model2, model2.1, and model2.2). Specifically, the output contains the residuals for each case for each of the three models and information about every case with missing data. Advice would be much appreciated. Thanks.

Upvotes: 3

Views: 2793

Answers (2)

Molx
Molx

Reputation: 6931

An alternative is to use David Robinson's broom package, which will convert the fit output into a dataframe. For a summary of the fit (without coefficients), using 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

For more information on the fit, you should use tidy, but the number of rows will vary and the end data.frame should be created carefully.

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

I'm not an user of this for a very long time, so there might be better ways to join this dataframes together.

Upvotes: 2

MrFlick
MrFlick

Reputation: 206616

Note that lm() returns an object of class "lm" and summary() on that object produces a "summary.lm" object. There are custom print.lm() and print.summary.lm() objects. So what ever is printed to the console may be different than what's in the object itself.

When you manually concatenate (c()) two summary.lm objects, you create one aggregated list and lose the proper class. You probably want to return a list of objects instead

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

Upvotes: 2

Related Questions