The Matrix Sweep Operator

The swp function “sweeps” a matrix on the rows and columns given in index to produce a new matrix with those rows and columns “partialled out” by orthogonalization. This was defined as a fundamental statistical operation in multivariate methods by Beaton (1964) and expanded by Dempster (1969). It is closely related to orthogonal projection, but applied to a cross-products or covariance matrix, rather than to data.

Usage

swp(M, index)

Arguments

M

a numeric matrix

index

a numeric vector indicating the rows/columns to be swept. The entries must be less than or equal to the number or rows or columns in M. If missing, the function sweeps on all rows/columns 1:min(dim(M)).

Value

the matrix M with rows and columns in indices swept.

Details

If M is the partitioned matrix $$\left[ \begin{array}{cc} \mathbf {R} & \mathbf {S} \\ \mathbf {T} & \mathbf {U} \end{array} \right]$$ where \(R\) is \(q \times q\) then swp(M, 1:q) gives $$\left[ \begin{array}{cc} \mathbf {R}^{-1} & \mathbf {R}^{-1}\mathbf {S} \\ -\mathbf {TR}^{-1} & \mathbf {U}-\mathbf {TR}^{-1}\mathbf {S} \\ \end{array} \right]$$

References

Beaton, A. E. (1964), The Use of Special Matrix Operations in Statistical Calculus, Princeton, NJ: Educational Testing Service.

Dempster, A. P. (1969) Elements of Continuous Multivariate Analysis. Addison-Wesley, Reading, Mass.

See also

Proj, QR

Examples

data(therapy)
mod3 <- lm(therapy ~ perstest + IE + sex, data=therapy)
X <- model.matrix(mod3)
XY <- cbind(X, therapy=therapy$therapy)
XY
#>        (Intercept) perstest IE sexM therapy
#> John             1       26  3    1      32
#> Susan            1       24  4    0      40
#> Mary             1       22  8    0      44
#> Paul             1       33  4    1      44
#> Jenny            1       27  6    0      48
#> Rick             1       36  4    1      52
#> Cathy            1       30 10    0      56
#> Robert           1       38  4    1      56
#> Lisa             1       30 12    0      60
#> Tina             1       34 15    0      68
M <- crossprod(XY)
swp(M, 1)
#>             (Intercept) perstest  IE  sexM therapy
#> (Intercept)        -0.1       30   7   0.4      50
#> perstest           30.0      250  10  13.0     300
#> IE                  7.0       10 152 -13.0     304
#> sexM                0.4       13 -13   2.4     -16
#> therapy            50.0      300 304 -16.0    1000
swp(M, 1:2)
#>             (Intercept) perstest     IE    sexM therapy
#> (Intercept)       -3.70    0.120   5.80  -1.160    14.0
#> perstest           0.12   -0.004   0.04   0.052     1.2
#> IE                 5.80    0.040 151.60 -13.520   292.0
#> sexM              -1.16    0.052 -13.52   1.724   -31.6
#> therapy           14.00    1.200 292.00 -31.600   640.0