1# File src/library/stats/R/proj.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 20proj <- function(object, ...) UseMethod("proj") 21 22proj.default <- function(object, onedf = TRUE, ...) 23{ 24 if(!inherits(object$qr, "qr")) 25 stop("argument does not include a 'qr' component") 26 if(is.null(object$effects)) 27 stop("argument does not include an 'effects' component") 28 RB <- c(object$effects[seq(object$rank)], 29 rep.int(0, nrow(object$qr$qr) - object$rank)) 30 prj <- as.matrix(qr.Q(object$qr, Dvec = RB)) 31 DN <- dimnames(object$qr$qr) 32 dimnames(prj) <- list(DN[[1L]], DN[[2L]][seq(ncol(prj))]) 33 prj 34} 35 36proj.lm <- function(object, onedf = FALSE, unweighted.scale = FALSE, ...) 37{ 38 if(inherits(object, "mlm")) 39 stop("'proj' is not implemented for multiple responses") 40 rank <- object$rank 41 if(rank > 0) { 42 prj <- proj.default(object, onedf = TRUE)[, 1L:rank, drop = FALSE] 43 if(onedf) { 44 df <- rep.int(1, rank) 45 result <- prj 46 } else { 47 asgn <- object$assign[object$qr$pivot[1L:object$rank]] 48 uasgn <- unique(asgn) 49 nmeffect <- c("(Intercept)", 50 attr(object$terms, "term.labels"))[1 + uasgn] 51 nterms <- length(uasgn) 52 df <- vector("numeric", nterms) 53 result <- matrix(0, length(object$residuals), nterms) 54 dimnames(result) <- list(rownames(object$fitted.values), nmeffect) 55 for(i in seq_along(uasgn)) { 56 select <- (asgn == uasgn[i]) 57 df[i] <- sum(select) 58 result[, i] <- prj[, select, drop = FALSE] %*% rep.int(1, df[i]) 59 } 60 } 61 } else { 62 result <- NULL 63 df <- NULL 64 } 65 if(!is.null(wt <- object$weights) && unweighted.scale) 66 result <- result/sqrt(wt) 67 use.wt <- !is.null(wt) && !unweighted.scale 68 if(object$df.residual > 0) { 69 res <- if(use.wt) object$residuals * sqrt(wt) else object$residuals 70 if(!is.matrix(result)) { 71 result <- matrix(res, length(res), 1L, 72 dimnames = list(names(res), "Residuals")) 73 } else { 74 dn <- dimnames(result) 75 d <- dim(result) 76 result <- setNames(c(result, res), NULL) 77 dim(result) <- d + c(0, 1) 78 dimnames(result) <- list(names(res), c(dn[[2L]], "Residuals")) 79 } 80 df <- c(df, object$df.residual) 81 } 82 names(df) <- colnames(result) 83 attr(result, "df") <- df 84 attr(result, "formula") <- object$call$formula 85 attr(result, "onedf") <- onedf 86 if(!is.null(wt)) attr(result, "unweighted.scale") <- unweighted.scale 87 result 88} 89 90proj.aov <- function(object, onedf = FALSE, unweighted.scale = FALSE, ...) 91{ 92 if(inherits(object, "maov")) 93 stop("'proj' is not implemented for multiple responses") 94 factors.aov <- function(pnames, tfactor) 95 { 96 if(!is.na(int <- match("(Intercept)", pnames))) 97 pnames <- pnames[ - int] 98 tnames <- setNames(lapply(colnames(tfactor), function(x, mat) 99 rownames(mat)[mat[, x] > 0], tfactor), 100 colnames(tfactor)) 101 if(!is.na(match("Residuals", pnames))) { 102 enames <- c(rownames(tfactor) 103 [as.logical(tfactor %*% rep.int(1, ncol(tfactor)))], 104 "Within") 105 tnames <- append(tnames, list(Residuals = enames)) 106 } 107 result <- tnames[match(pnames, names(tnames))] 108 if(!is.na(int)) result <- c("(Intercept)" = "(Intercept)", result) 109 ## should reorder result, but probably OK 110 result 111 } 112 projections <- NextMethod("proj") 113 t.factor <- attr(terms(object), "factors") 114 attr(projections, "factors") <- 115 factors.aov(colnames(projections), t.factor) 116 attr(projections, "call") <- object$call 117 attr(projections, "t.factor") <- t.factor 118 class(projections) <- "aovproj" 119 projections 120} 121 122 123proj.aovlist <- function(object, onedf = FALSE, unweighted.scale = FALSE, ...) 124{ 125 attr.xdim <- function(x) 126 { 127 ## all attributes except names, dim and dimnames 128 atrf <- attributes(x) 129 atrf[is.na(match(names(atrf), c("names", "dim", "dimnames")))] 130 } 131 "attr.assign<-" <- function(x, value) 132 { 133 ## assign to x all attributes in attr.x 134 ## attributes(x)[names(value)] <- value not allowed in R 135 for(nm in names(value)) attr(x, nm) <- value[nm] 136 x 137 } 138 factors.aovlist <- function(pnames, tfactor, 139 strata = FALSE, efactor = FALSE) 140 { 141 if(!is.na(int <- match("(Intercept)", pnames))) pnames <- pnames[-int] 142 tnames <- apply(tfactor, 2L, function(x, nms) 143 nms[as.logical(x)], rownames(tfactor)) 144 if(!missing(efactor)) { 145 enames <- 146 if(!is.na(err <- match(strata, colnames(efactor)))) 147 rownames(efactor)[as.logical(efactor[, err])] 148 else if(strata == "Within") 149 c(rownames(efactor)[as.logical(efactor %*% rep.int(1, ncol(efactor)))], 150 "Within") 151 ## else NULL 152 if(!is.null(enames)) 153 tnames <- append(tnames, list(Residuals = enames)) 154 } 155 result <- tnames[match(pnames, names(tnames))] 156 if(!is.na(int)) 157 result <- c("(Intercept)" = "(Intercept)", result) 158 ##should reorder result, but probably OK 159 result 160 } 161 if(unweighted.scale && is.null(attr(object, "weights"))) 162 unweighted.scale <- FALSE 163 err.qr <- attr(object, "error.qr") 164 Terms <- terms(object, "Error") 165 t.factor <- attr(Terms, "factors") 166 i <- attr(Terms, "specials")$Error 167 t <- attr(Terms, "variables")[[1 + i]] 168 error <- Terms 169 error[[3L]] <- t[[2L]] 170 e.factor <- attr(terms(formula(error)), "factors") 171 n <- nrow(err.qr$qr) 172 n.object <- length(object) 173 result <- setNames(vector("list", n.object), names(object)) 174 D1 <- seq_len(NROW(err.qr$qr)) 175 if(unweighted.scale) wt <- attr(object, "weights") 176 for(i in names(object)) { 177 prj <- proj.lm(object[[i]], onedf = onedf) 178 if(unweighted.scale) prj <- prj/sqrt(wt) 179 result.i <- matrix(0, n, ncol(prj), dimnames = list(D1, colnames(prj))) 180 select <- rownames(object[[i]]$qr$qr) %||% 181 rownames(object[[i]]$residuals) 182 result.i[select, ] <- prj 183 result[[i]] <- as.matrix(qr.qy(err.qr, result.i)) 184 attr.assign(result[[i]]) <- attr.xdim(prj) 185 D2i <- colnames(prj) 186 dimnames(result[[i]]) <- list(D1, D2i) 187 attr(result[[i]], "factors") <- 188 factors.aovlist(D2i, t.factor, strata = i, efactor = e.factor) 189 } 190 attr(result, "call") <- attr(object, "call") 191 attr(result, "e.factor") <- e.factor 192 attr(result, "t.factor") <- t.factor 193 class(result) <- c("aovprojlist", "listof") 194 result 195} 196 197terms.aovlist <- function(x, ...) 198{ 199 x <- attr(x, "terms") 200 terms(x, ...) 201} 202 203## wish of PR#13505 204as.data.frame.aovproj <- function(x, ...) as.data.frame(unclass(x), ...) 205