1# File src/library/stats/R/dummy.coef.R 2# Part of the R package, https://www.R-project.org 3# 4# Copyright (C) 1998-2020 The R Core Team 5# Copyright (C) 1998 B. D. Ripley 6# 7# This program is free software; you can redistribute it and/or modify 8# it under the terms of the GNU General Public License as published by 9# the Free Software Foundation; either version 2 of the License, or 10# (at your option) any later version. 11# 12# This program is distributed in the hope that it will be useful, 13# but WITHOUT ANY WARRANTY; without even the implied warranty of 14# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 15# GNU General Public License for more details. 16# 17# A copy of the GNU General Public License is available at 18# https://www.R-project.org/Licenses/ 19 20dummy.coef <- function(object, ...) UseMethod("dummy.coef") 21 22dummy.coef.lm <- function(object, use.na=FALSE, ...) 23{ 24 xl <- object$xlevels 25 if(!length(xl)) # no factors in model 26 return(as.list(coef(object))) 27 Terms <- terms(object) 28 tl <- attr(Terms, "term.labels") 29 int <- attr(Terms, "intercept") 30 facs <- attr(Terms, "factors")[-1, , drop=FALSE] 31 Terms <- delete.response(Terms) 32 mf <- object$model %||% model.frame(object) 33 vars <- dimnames(facs)[[1]] # names 34 xtlv <- lapply(mf[,vars, drop=FALSE], levels) ## levels 35 nxl <- pmax(lengths(xtlv), 1L) ## (named) number of levels 36 lterms <- apply(facs, 2L, function(x) prod(nxl[x > 0])) 37 nl <- sum(lterms) 38 ## dummy: data frame of vars 39 args <- sapply(vars, function(i) 40 if (nxl[i] == 1) rep.int(1, nl) 41 else factor(rep.int(xtlv[[i]][1L], nl), levels = xtlv[[i]]), 42 simplify=FALSE) 43 ## dummy <- as.data.frame(args) # slightly more efficiently: 44 dummy <- do.call(data.frame, args); names(dummy) <- vars 45 pos <- 0L 46 rn <- rep.int(tl, lterms) 47 rnn <- character(nl) # all "" --- will be names of rows 48 for(j in tl) { 49 i <- vars[facs[, j] > 0] 50 ifac <- i[nxl[i] > 1] 51 lt.j <- lterms[[j]] 52 if(length(ifac) == 0L) { # quantitative factor 53 rnn[pos+1L] <- j 54 } else { 55 p.j <- pos + seq_len(lt.j) 56 if(length(ifac) == 1L) { # main effect 57 dummy[p.j, ifac] <- x.i <- xtlv[[ifac]] 58 rnn[p.j] <- as.character(x.i) 59 } else { # interaction 60 tmp <- expand.grid(xtlv[ifac], KEEP.OUT.ATTRS=FALSE) 61 dummy[p.j, ifac] <- tmp 62 rnn[p.j] <- apply(as.matrix(tmp), 1L, paste, collapse = ":") 63 } 64 } 65 pos <- pos + lt.j 66 } 67 attr(dummy,"terms") <- attr(mf,"terms") 68 lcontr <- object$contrasts 69 lci <- vapply(dummy, is.factor, NA) 70 lcontr <- lcontr[names(lci)[lci]] ## factors with 1 level have disappeared (?) 71 mm <- model.matrix(Terms, dummy, lcontr, xl) 72 if(anyNA(mm)) { 73 warning("some terms will have NAs due to the limits of the method") 74 mm[is.na(mm)] <- NA 75 } 76 coef <- object$coefficients 77 if(!use.na) coef[is.na(coef)] <- 0 78 asgn <- attr(mm,"assign") 79 res <- setNames(vector("list", length(tl)), tl) 80 if(isM <- is.matrix(coef)) { # isM is true for "mlm", multivariate lm (incl manova) 81 for(j in seq_along(tl)) { 82 keep <- which(asgn == j) 83 cf <- coef[keep, , drop=FALSE] 84 ij <- rn == tl[j] 85 cf <- 86 if (any(na <- is.na(cf))) { 87 if(ncol(cf) >= 2) 88 stop("multivariate case with missing coefficients is not yet implemented") 89 ## else: 1 column --> treat 'cf' as vector 90 rj <- t( mm[ij, keep[!na], drop=FALSE] %*% cf[!na]) 91 rj[apply(mm[ij, keep[ na], drop=FALSE] != 0, 1L, any)] <- NA 92 rj 93 } else 94 t(mm[ij, keep, drop = FALSE] %*% cf) 95 dimnames(cf) <- list(colnames(coef), rnn[ij]) 96 res[[j]] <- cf 97 } 98 } else { ## regular univariate lm case 99 for(j in seq_along(tl)) { 100 keep <- which(asgn == j) 101 cf <- coef[keep] 102 ij <- rn == tl[j] 103 res[[j]] <- 104 if (any(na <- is.na(cf))) { 105 rj <- setNames(drop(mm[ij, keep[!na], drop = FALSE] %*% 106 cf[!na]), rnn[ij]) 107 rj[apply(mm[ij, keep[na], drop=FALSE] != 0, 1L, any)] <- NA 108 rj 109 } else 110 setNames(drop(mm[ij, keep, drop = FALSE] %*% cf), rnn[ij]) 111 } 112 } 113 if(int > 0) 114 res <- c(list("(Intercept)" = if(isM) coef[int, ] else coef[int]), res) 115 structure(res, class = "dummy_coef", matrix = isM) 116} 117 118## NB: This is very much duplication from dummy.coef.lm -- keep in sync ! 119dummy.coef.aovlist <- function(object, use.na = FALSE, ...) 120{ 121 xl <- attr(object, "xlevels") 122 if(!length(xl)) # no factors in model 123 return(as.list(coef(object))) 124 Terms <- terms(object, specials="Error") 125 err <- attr(Terms,"specials")$Error - 1 126 tl <- attr(Terms, "term.labels")[-err] 127 int <- attr(Terms, "intercept") 128 facs <- attr(Terms, "factors")[-c(1,1+err), -err, drop=FALSE] 129 stopifnot(length(names(object)) == (N <- length(object))) 130 mf <- object$model %||% model.frame(object) 131 vars <- dimnames(facs)[[1]] # names 132 xtlv <- lapply(mf[,vars, drop=FALSE], levels) ## levels 133 nxl <- pmax(lengths(xtlv), 1L) ## (named) number of levels 134 lterms <- apply(facs, 2L, function(x) prod(nxl[x > 0])) 135 nl <- sum(lterms) 136 args <- setNames(vector("list", length(vars)), vars) 137 for(i in vars) 138 args[[i]] <- if(nxl[[i]] == 1) rep.int(1, nl) 139 else factor(rep.int(xl[[i]][1L], nl), levels = xl[[i]]) 140 ## dummy <- as.data.frame(args) # slightly more efficiently: 141 dummy <- do.call(data.frame, args); names(dummy) <- vars 142 pos <- 0L 143 rn <- rep.int(tl, lterms) 144 rnn <- character(nl) # all "" --- will be names of rows 145 for(j in tl) { 146 i <- vars[facs[, j] > 0] 147 ifac <- i[nxl[i] > 1] 148 lt.j <- lterms[[j]] 149 if(length(ifac) == 0L) { # quantitative factor 150 rnn[pos+1L] <- j 151 } else { 152 p.j <- pos + seq_len(lt.j) 153 if(length(ifac) == 1L) { # main effect 154 dummy[p.j, ifac] <- x.i <- xtlv[[ifac]] 155 rnn[p.j] <- as.character(x.i) 156 } else { # interaction 157 tmp <- expand.grid(xtlv[ifac], KEEP.OUT.ATTRS=FALSE) 158 dummy[p.j, ifac] <- tmp 159 rnn[p.j] <- apply(as.matrix(tmp), 1L, paste, collapse = ":") 160 } 161 } 162 pos <- pos + lt.j 163 } 164 form <- paste0("~", paste0(tl, collapse = " + "), if(!int) "- 1") 165 lcontr <- object$contrasts 166 lci <- vapply(dummy, is.factor, NA) 167 lcontr <- lcontr[names(lci)[lci]] ## factors with 1 level have disappeared 168 mm <- model.matrix(terms(formula(form)), dummy, lcontr, xl) 169 tl <- c("(Intercept)", tl) 170 res <- setNames(vector("list", N), names(object)) 171 allasgn <- attr(mm, "assign") 172 for(i in names(object)) { 173 coef <- object[[i]]$coefficients 174 if(!use.na) coef[is.na(coef)] <- 0 175 asgn <- object[[i]]$assign 176 uasgn <- unique(asgn) 177 tll <- tl[1L + uasgn] 178 mod <- setNames(vector("list", length(tll)), tll) 179### FIXME --- npk.aovE --- fails : "N" gets length 0 !!!! 180 if((isM <- is.matrix(coef))) { # "mlm", multivariate lm (incl manova) 181 for(j in uasgn) { 182 keep <- which(asgn == j) 183 cf <- coef[keep, , drop=FALSE] 184 ij <- rn == tl[j] 185 cf <- 186 if (any(na <- is.na(cf))) { 187 if(ncol(cf) >= 2) 188 stop("multivariate case with missing coefficients is not yet implemented") 189 if(j == 0) { 190 structure(cf[!na], names="(Intercept)") 191 } else { 192 ## else: 1 column --> treat 'cf' as vector 193 rj <- t( mm[ij, keep[!na], drop=FALSE] %*% cf[!na]) 194 rj[apply(mm[ij, keep[ na], drop=FALSE] != 0, 1L, any)] <- NA 195 rj 196 } 197 } else { # no NA's 198 if(j == 0) 199 structure(cf, names="(Intercept)") 200 else 201 t(mm[ij, keep, drop=FALSE] %*% cf) 202 } 203 dimnames(cf) <- list(colnames(coef), rnn[ij]) 204 mod[[tl[j+1L]]] <- cf 205 } 206 } else { ## regular univariate lm case 207 for(j in uasgn) { 208 keep <- which(asgn == j) 209 cf <- coef[keep] 210 mod[[tl[j+1L]]] <- 211 if(j == 0) { 212 structure(cf, names="(Intercept)") 213 } else { 214 ij <- rn == tl[j+1L] 215 if (any(na <- is.na(cf))) { 216 rj <- setNames(drop(mm[ij, keep[!na], drop = FALSE] %*% 217 cf[!na]), rnn[ij]) 218 rj[apply(mm[ij, keep[na], drop=FALSE] != 0, 1L, any)] <- NA 219 rj 220 } else 221 setNames(drop(mm[ij, allasgn == j, drop=FALSE] %*% cf), 222 rnn[ij]) 223 } 224 } 225 } 226 res[[i]] <- structure(mod, matrix = isM) 227 } ## for( i ) 228 229 structure(res, class = "dummy_coef_list") 230} 231 232print.dummy_coef <- function(x, ..., title) 233{ 234 terms <- names(x) 235 n <- length(x) 236 isM <- attr(x, "matrix") 237 nr.x <- if(isM) vapply(x, NROW, 1L) else lengths(x) 238 line <- 0L # 'lineEnd' 239 if(isM) { # "mlm" - multivariate case 240 ansnrow <- sum(1L + nr.x) 241 addcol <- max(nr.x) - 1L 242 nm <- addcol + if(isM) max(vapply(x, NCOL, 1L)) else 1L 243 ans <- matrix("", ansnrow , nm) 244 rn <- character(ansnrow) 245 for (j in seq_len(n)) { 246 this <- as.matrix(x[[j]]) 247 nr1 <- nrow(this) 248 nc1 <- ncol(this) 249 dn <- dimnames(this) 250 dimnames(this) <- 251 list(dn[[1]] %||% character(nr1), 252 dn[[2]] %||% character(nc1)) 253 if(nc1 > 1L) { 254 lin0 <- line + 2L 255 line <- line + nr1 + 1L 256 ans[lin0 - 1L, addcol + (1L:nc1)] <- colnames(this) 257 ans[lin0:line, addcol + (1L:nc1)] <- format(this, ...) 258 rn[lin0 - 1L] <- paste0(terms[j], ": ") 259 } else { 260 lin0 <- line + 1L 261 line <- line + nr1 262 ans[lin0:line, addcol + 1L] <- format(this, ...) 263 rn[lin0] <- paste0(terms[j], ": ") 264 } 265 if(addcol > 0) ans[lin0:line, addcol] <- rownames(this) 266 } 267 } else { ## regular univariate lm case 268 nm <- max(nr.x) 269 ans <- matrix("", 2L*n, nm) 270 rn <- character(2L*n) # "" 271 for (j in seq_len(n)) { 272 this <- x[[j]] 273 n1 <- length(this) 274 if(n1 > 1) { 275 line <- line + 2L 276 ans[line-1L, 1L:n1] <- names(this) 277 ans[line, 1L:n1] <- format(this, ...) 278 rn [line-1L] <- paste0(terms[j], ": ") 279 } else { 280 line <- line + 1L 281 ans[line, 1L:n1] <- format(this, ...) 282 rn[line] <- paste0(terms[j], ": ") 283 } 284 } 285 } 286 dimnames(ans) <- list(rn, character(nm)) 287 cat(if(missing(title)) "Full coefficients are" else title, "\n") 288 print(ans[1L:line, , drop=FALSE], quote=FALSE, right=TRUE) 289 invisible(x) 290} 291 292print.dummy_coef_list <- function(x, ...) 293{ 294 for(strata in names(x)) 295 print.dummy_coef(x[[strata]], ..., title=paste("\n Error:", strata)) 296 invisible(x) 297} 298