ANOVA Table for Full Linear Model

I'm teaching Multiple Regression again this semester, so I'm once again reminded that the ANOVA output from R is not the same as from other statistical software. A quick google turned up this post from r-help which has essentially my exact question, so I'll quote from it here (data changed for better reproducibility):

I am trying to figure out how to get an ANOVA table that shows the sum of squares. degrees of freedom, etc, for the full model versus the error (aka residuals).

Here is an example of the kind of table I'd like to get:

Analysis of Variance
Source           DF          SS          MS         F        P
Regression       2        0.9715       0.485728    52.13    0.000
Error            29       0.2702       0.009318
Total            31       1.2417

This kind of table is prevalent throughout my statistics textbook, and can apparently be easily obtained in other statistical software tools. I'm not saying this as a gripe, but just as evidence that I'm not trying to do something obviously bizarre.

Here is an example of the only kind of ANOVA table for a single linear model that I've been able to get using R:

library(Stat2Data)
data("NFLStandings2016")
m1 = lm(WinPct ~ PointsFor + PointsAgainst, data = NFLStandings2016)
anova(m1)
#> Analysis of Variance Table
#> 
#> Response: WinPct
#>               Df  Sum Sq Mean Sq F value    Pr(>F)    
#> PointsFor      1 0.41262 0.41262  44.280 2.692e-07 ***
#> PointsAgainst  1 0.55884 0.55884  59.972 1.537e-08 ***
#> Residuals     29 0.27023 0.00932                      
#> ---
#> Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Created on 2019-10-11 by the reprex package (v0.2.1)

Is there a way to get an ANOVA table with the full linear regression model considered as a whole rather than broken down into each additional predictor variable? In other words, is there a way to get the former kind of table?

The response on r-help includes quotes like "the logical thing to do would be" and also "if you really want the column of MS, you have a little extra work to do." Is there really no command in R that produces this kind of summary table? My guess is it's hiding in some modeling package and I haven't stumbled across it, but I also wouldn't be surprised if there was some base print method I'm missing.

Hi Amelia,

Does the following code achieve what you were hoping for?

library(Stat2Data)
library(rms)
data("NFLStandings2016")

anova(rms::ols(WinPct ~ PointsFor + PointsAgainst, data = NFLStandings2016))
                Analysis of Variance          Response: WinPct 

 Factor        d.f. Partial SS MS          F     P     
 PointsFor      1   0.3895631  0.389563129 41.81 <.0001
 PointsAgainst  1   0.5588400  0.558839981 59.97 <.0001
 REGRESSION     2   0.9714563  0.485728140 52.13 <.0001
 ERROR         29   0.2702337  0.009318404           

Here is a modification of the base R anova call. I believe this should work for general lm/aov objects, but I haven't tested it thoroughly. Note that the reg_collapse argument will collapse to Source by default, but by setting it to FALSE you can get the original rows back, just with a Total SS row.

anova_alt = function (object, reg_collapse=TRUE,...) 
{
  if (length(list(object, ...)) > 1L) 
    return(anova.lmlist(object, ...))
  if (!inherits(object, "lm")) 
    warning("calling anova.lm(<fake-lm-object>) ...")
  w <- object$weights
  ssr <- sum(if (is.null(w)) object$residuals^2 else w * object$residuals^2)
  mss <- sum(if (is.null(w)) object$fitted.values^2 else w * 
               object$fitted.values^2)
  if (ssr < 1e-10 * mss) 
    warning("ANOVA F-tests on an essentially perfect fit are unreliable")
  dfr <- df.residual(object)
  p <- object$rank
  if (p > 0L) {
    p1 <- 1L:p
    comp <- object$effects[p1]
    asgn <- object$assign[stats:::qr.lm(object)$pivot][p1]
    nmeffects <- c("(Intercept)", attr(object$terms, "term.labels"))
    tlabels <- nmeffects[1 + unique(asgn)]
    ss <- c(vapply(split(comp^2, asgn), sum, 1), ssr)
    df <- c(lengths(split(asgn, asgn)), dfr)
    if(reg_collapse){
      if(attr(object$terms, "intercept")){
        collapse_p<-2:(length(ss)-1)
        ss<-c(ss[1],sum(ss[collapse_p]),ss[length(ss)])
        df<-c(df[1],sum(df[collapse_p]),df[length(df)])
        tlabels<-c(tlabels[1],"Source")
      } else{
        collapse_p<-1:(length(ss)-1)
        ss<-c(sum(ss[collapse_p]),ss[length(ss)])
        df<-c(df[1],sum(df[collapse_p]),df[length(df)])
        tlabels<-c("Regression")
      }
    }
  }else {
    ss <- ssr
    df <- dfr
    tlabels <- character()
    if(reg_collapse){
      collapse_p<-1:(length(ss)-1)
      ss<-c(sum(ss[collapse_p]),ss[length(ss)])
      df<-c(df[1],sum(df[collapse_p]),df[length(df)])
    }
  }
 
  ms <- ss/df
  f <- ms/(ssr/dfr)
  P <- pf(f, df, dfr, lower.tail = FALSE)
  table <- data.frame(df, ss, ms, f, P)
  table <- rbind(table, 
                 colSums(table))
  table$ms[nrow(table)]<-table$ss[nrow(table)]/table$df[nrow(table)]
  table[length(P):(length(P)+1), 4:5] <- NA
  dimnames(table) <- list(c(tlabels, "Error","Total"), 
                          c("Df","SS", "MS", "F", 
                            "P"))
  if (attr(object$terms, "intercept")){
    table <- table[-1, ]
    table$MS[nrow(table)]<-table$MS[nrow(table)]*(table$Df[nrow(table)])/(table$Df[nrow(table)]-1)
    table$Df[nrow(table)]<-table$Df[nrow(table)]-1
  }
  structure(table, heading = c("Analysis of Variance Table\n"), 
            class = c("anova", "data.frame"))
}

## Warpbreaks example
fm1 <- lm(breaks ~ wool*tension, data = warpbreaks)
anova_alt(fm1)
#> Analysis of Variance Table
#> 
#>        Df    SS     MS      F         P
#> Source  5  3488 697.54 5.8279 0.0002772
#> Error  48  5745 119.69                 
#> Total  53 52018 981.47

Created on 2019-10-11 by the reprex package (v0.3.0)