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