dummy.coef() and print.dummy_coef() modified to deal with multivariate lm fits

dummy.coef_mult.R

# modified against R-Devel build 2016-01-22 r69971
#  File src/library/stats/R/dummy.coef.R
#  Part of the R package, https://www.R-project.org
#
#  Copyright (C) 1998 B. D. Ripley
#  Copyright (C) 1998-2015 The R Core Team
#
#  This program is free software; you can redistribute it and/or modify
#  it under the terms of the GNU General Public License as published by
#  the Free Software Foundation; either version 2 of the License, or
#  (at your option) any later version.
#
#  This program is distributed in the hope that it will be useful,
#  but WITHOUT ANY WARRANTY; without even the implied warranty of
#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#  GNU General Public License for more details.
#
#  A copy of the GNU General Public License is available at
#  https://www.R-project.org/Licenses/

dummy.coef <- function(object, ...) UseMethod("dummy.coef")

dummy.coef.lm <- function(object, use.na=FALSE, ...)
{
    xl <- object$xlevels
    if(!length(xl)) # no factors in model
	return(as.list(coef(object)))
    Terms <- terms(object)
    tl <- attr(Terms, "term.labels")
    int <- attr(Terms, "intercept")
    facs <- attr(Terms, "factors")[-1, , drop=FALSE]
    Terms <- delete.response(Terms)
    vars <- all.vars(Terms) # e.g. drops I(.), ...
    nxl <- setNames(rep.int(1, length(vars)), vars)
    tmp <- lengths(xl)
    nxl[names(tmp)] <- tmp
    lterms <- apply(facs, 2L, function(x) prod(nxl[x > 0]))
    nl <- sum(lterms)
    args <- setNames(vector("list", length(vars)), vars)
    for(i in vars)
	args[[i]] <- if(nxl[[i]] == 1) rep.int(1, nl)
	else factor(rep.int(xl[[i]][1L], nl), levels = xl[[i]])
    dummy <- do.call("data.frame", args)
    pos <- 0
    rn <- rep.int(tl, lterms)
    rnn <- rep.int("", nl)
    for(j in tl) {
	i <- vars[facs[, j] > 0]
	ifac <- i[nxl[i] > 1]
	if(length(ifac) == 0L) {        # quantitative factor
	    rnn[pos+1] <- j
	} else if(length(ifac) == 1L) {	# main effect
	    dummy[ pos+1L:lterms[j], ifac ] <- xl[[ifac]]
	    rnn[ pos+1L:lterms[j] ] <- as.character(xl[[ifac]])
	} else {			# interaction
	    tmp <- expand.grid(xl[ifac])
	    dummy[ pos+1L:lterms[j], ifac ] <- tmp
	    rnn[ pos+1L:lterms[j] ] <-
		apply(as.matrix(tmp), 1L, function(x) paste(x, collapse=":"))
	}
	pos <- pos + lterms[j]
    }
    ## some terms like poly(x,1) will give problems here, so allow
    ## NaNs and set to NA afterwards.
    mf <- model.frame(Terms, dummy, na.action=function(x)x, xlev=xl)
    mm <- model.matrix(Terms, mf, object$contrasts, xl)
    if(anyNA(mm)) {
        warning("some terms will have NAs due to the limits of the method")
        mm[is.na(mm)] <- NA
    }
    coef <- as.matrix(object$coefficients)
    if(!use.na) coef[is.na(coef)] <- 0
    asgn <- attr(mm,"assign")
    res <- setNames(vector("list", length(tl)), tl)
    for(j in seq_along(tl)) {
	keep <- asgn == j
	ij <- rn == tl[j]
      dummyCoeffs <- t(mm[ij, keep, drop=FALSE] %*% coef[keep, ])
      dimnames(dummyCoeffs) <- list(colnames(coef), rnn[ij])
      res[[j]] <- dummyCoeffs
    }
    if(int > 0) {
	res <- c(list("(Intercept)" = coef[int]), res)
    }
    class(res) <- "dummy_coef"
    res
}

print.dummy_coef <- function(x, ..., title)
{
    terms <- names(x)
    n <- length(x)
    xnrow <- vapply(x, function(y) nrow(as.matrix(y)), 1L)
    ansnrow <- sum(1L + xnrow)
    addcol <- max(xnrow)-1
    nm <- addcol + max(vapply(x, function(y) ncol(as.matrix(y)), 1L))
    ans <- matrix("", ansnrow , nm)
    rn <- rep.int("", ansnrow)
    lineBeg <- 0
    lineEnd <- 0
    for (j in seq_len(n)) {
	this <- as.matrix(x[[j]])
      dnthis <- dimnames(this)
      rnthis <- if(is.null(dnthis[[1]])) rep.int("", nrow(this)) else dnthis[[1]]
      cnthis <- if(is.null(dnthis[[2]])) rep.int("", ncol(this)) else dnthis[[2]]
      dimnames(this) <- list(rnthis, cnthis)
      nr1 <- nrow(this)
	nc1 <- ncol(this)
	if(nc1 > 1L) {
	  lineBeg <- lineEnd + 2
        lineEnd <- lineBeg + nr1 - 1
	  ans[lineBeg-1, addcol + (1L:nc1)] <- colnames(this)
	  ans[lineBeg:lineEnd, addcol + (1L:nc1)] <- format(this, ...)
	  rn[lineBeg-1] <- paste0(terms[j], ":   ")
	} else {
	  lineBeg <- lineEnd + 1
        lineEnd <- lineBeg + nr1 - 1
	  ans[lineBeg:lineEnd, addcol + 1L] <- format(this, ...)
	  rn[lineBeg] <- paste0(terms[j], ":   ")
	}
      if(addcol > 0) ans[lineBeg:lineEnd, addcol] <- rownames(this)
    }
    rownames(ans) <- rn
    colnames(ans) <- rep.int("", nm)
    cat(if(missing(title)) "Full coefficients are" else title, "\n")
    print(ans[1L:lineEnd, , drop=FALSE], quote=FALSE, right=TRUE)
    invisible(x)
}

Line #83 is not quite correct, missing a , but that was easy to fix. The bigger issue is that the other patch I had deals with NA coefficients nicely, and that's even a bit more complicated in the multivariate case. For now, I will omit the combined presence of
NA coefficients and multivariate response. No need, to post anything further... before I commit my version to the R source repository (soonish).