1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
|
#' The ResidualMatrix class
#'
#' The ResidualMatrix class supports delayed calculation of the residuals from a linear model fit.
#' This serves as a light-weight representation of what would otherwise be a large dense matrix in memory.
#' It also enables efficient matrix multiplication based on features of the the original matrix (e.g., sparsity).
#'
#' @section Construction:
#' \code{ResidualMatrix(x, design=NULL, keep=NULL)} returns a ResidualMatrix object, given:
#' \itemize{
#' \item \code{x}, a matrix-like object.
#' This can alternatively be a ResidualMatrixSeed, in which case \code{design} and \code{keep} are ignored.
#' \item \code{design}, a numeric matrix containing the experimental design,
#' to be used for linear model fitting on each \emph{column} of \code{x}.
#' This defaults to an intercept-only matrix.
#' \item \code{keep}, an integer vector specifying the columns of \code{design} to \emph{not} regress out.
#' By default, all columns of \code{design} are regressed out.
#' \item \code{restrict}, an integer or logical vector specifying the rows of \code{x} to use for model fitting.
#' If \code{NULL}, all rows of \code{x} are used.
#' }
#'
#' When \code{keep=NULL}, the ResidualMatrix contains values equivalent to \code{lm.fit(x=design, y=x)$residuals}.
#'
#' @section Methods:
#' In the following code chunks, \code{x} is a ResidualMatrix object:
#' \itemize{
#' \item \code{x[i, j, .., drop=FALSE]} will return a ResidualMatrix object for the specified row and column subsets,
#' or a numeric vector if either \code{i} or \code{j} are of length 1.
#' \item \code{t(x)} will return a ResidualMatrix object with transposed contents.
#' \item \code{dimnames(x) <- value} will return a ResidualMatrix object where the rows and columns are renamed by \code{value},
#' a list of two character vectors (or \code{NULL}).
#' }
#'
#' \code{\link{colSums}(x)}, \code{\link{colMeans}(x)}, \code{\link{rowSums}(x)} and \code{\link{rowMeans}(x)}
#' will return the relevant statistics for a ResidualMatrix \code{x}.
#'
#' \code{\%*\%}, \code{\link{crossprod}} and \code{\link{tcrossprod}} can also be applied
#' where one or both of the arguments are ResidualMatrix objects.
#'
#' ResidualMatrix objects are derived from \linkS4class{DelayedMatrix} objects and support all of valid operations on the latter.
#' All operations not listed here will use the underlying \pkg{DelayedArray} machinery.
#' Unary or binary operations will generally create a new DelayedMatrix instance containing a \linkS4class{ResidualMatrixSeed}.
#'
#' @author
#' Aaron Lun
#'
#' @examples
#' design <- model.matrix(~gl(5, 50))
#'
#' library(Matrix)
#' y0 <- rsparsematrix(nrow(design), 200, 0.1)
#' y <- ResidualMatrix(y0, design)
#' y
#'
#' # For comparison:
#' fit <- lm.fit(x=design, y=as.matrix(y0))
#' DelayedArray(fit$residuals)
#'
#' # Keeping some of the factors:
#' y2 <- ResidualMatrix(y0, design, keep=1:2)
#' y2
#' DelayedArray(fit$residuals + design[,1:2] %*% fit$coefficients[1:2,])
#'
#' # Matrix multiplication:
#' crossprod(y)
#' tcrossprod(y)
#' y %*% rnorm(200)
#'
#' @aliases
#' ResidualMatrix
#' ResidualMatrix-class
#' dimnames<-,ResidualMatrix,ANY-method
#' t,ResidualMatrix-method
#' [,ResidualMatrix,ANY,ANY,ANY-method
#' colSums,ResidualMatrix-method
#' rowSums,ResidualMatrix-method
#' colMeans,ResidualMatrix-method
#' rowMeans,ResidualMatrix-method
#' %*%,ANY,ResidualMatrix-method
#' %*%,ResidualMatrix,ANY-method
#' %*%,ResidualMatrix,ResidualMatrix-method
#' crossprod,ResidualMatrix,missing-method
#' crossprod,ResidualMatrix,ANY-method
#' crossprod,ANY,ResidualMatrix-method
#' crossprod,ResidualMatrix,ResidualMatrix-method
#' tcrossprod,ResidualMatrix,missing-method
#' tcrossprod,ResidualMatrix,ANY-method
#' tcrossprod,ANY,ResidualMatrix-method
#' tcrossprod,ResidualMatrix,ResidualMatrix-method
#'
#' @name ResidualMatrix-class
#' @docType class
NULL
#' @export
#' @importFrom DelayedArray DelayedArray
ResidualMatrix <- function(x, design=NULL, keep=NULL, restrict=NULL) {
DelayedArray(ResidualMatrixSeed(x, design, keep, restrict))
}
#' @export
#' @importFrom DelayedArray DelayedArray new_DelayedArray
setMethod("DelayedArray", "ResidualMatrixSeed",
function(seed) new_DelayedArray(seed, Class="ResidualMatrix")
)
###################################
# Overridden utilities from DelayedArray, for efficiency.
#' @export
#' @importFrom DelayedArray DelayedArray seed
setReplaceMethod("dimnames", "ResidualMatrix", function(x, value) {
DelayedArray(rename_ResidualMatrixSeed(seed(x), value))
})
#' @export
#' @importFrom Matrix t
#' @importFrom DelayedArray DelayedArray seed
setMethod("t", "ResidualMatrix", function(x) {
DelayedArray(transpose_ResidualMatrixSeed(seed(x)))
})
#' @export
#' @importFrom DelayedArray DelayedArray seed
setMethod("[", "ResidualMatrix", function(x, i, j, ..., drop=TRUE) {
if (missing(i)) i <- NULL
if (missing(j)) j <- NULL
rseed <- seed(x)
if (was_transposed <- is_transposed(rseed)) {
rseed <- transpose_ResidualMatrixSeed(rseed)
tmp <- i
i <- j
j <- tmp
}
rseed <- subset_ResidualMatrixSeed(rseed, i, j)
if (was_transposed) {
rseed <- transpose_ResidualMatrixSeed(rseed)
}
out <- DelayedArray(rseed)
if (drop && any(dim(out)==1L)) {
return(drop(out))
}
out
})
|
