Form component of GAMM covariance matrix

This is a service routine for gamm. Given, \(V\), an estimated covariance matrix obtained using extract.lme.cov2 this routine forms a matrix square root of \( X^TV^{-1}X\) as efficiently as possible, given the structure of \(V\) (usually sparse).

formXtViX(V,X)

Arguments

V

A data covariance matrix list returned from extract.lme.cov2

X

A model matrix.

Details

The covariance matrix returned by extract.lme.cov2 may be in a packed and re-ordered format, since it is usually sparse. Hence a special service routine is required to form the required products involving this matrix.

Value

A matrix, R such that crossprod(R) gives \( X^TV^{-1}X\).

References

For lme see:

Pinheiro J.C. and Bates, D.M. (2000) Mixed effects Models in S and S-PLUS. Springer

For details of how GAMMs are set up for estimation using lme see:

Wood, S.N. (2006) Low rank scale invariant tensor product smooths for Generalized Additive Mixed Models. Biometrics 62(4):1025-1036

Author

Simon N. Wood simon.wood@r-project.org

See also

gamm, extract.lme.cov2

Examples

require(mgcv)
library(nlme)
data(ergoStool)
b <- lme(effort ~ Type, data=ergoStool, random=~1|Subject)
V1 <- extract.lme.cov(b, ergoStool)
V2 <- extract.lme.cov2(b, ergoStool)
X <- model.matrix(b, data=ergoStool)
crossprod(formXtViX(V2, X))
#>             (Intercept)    TypeT2    TypeT3    TypeT4
#> (Intercept)    4.330827  1.082707  1.082707  1.082707
#> TypeT2         1.082707  5.846203 -1.587832 -1.587832
#> TypeT3         1.082707 -1.587832  5.846203 -1.587832
#> TypeT4         1.082707 -1.587832 -1.587832  5.846203
t(X)
#>             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
#> (Intercept) 1 1 1 1 1 1 1 1 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1
#> TypeT2      0 1 0 0 0 1 0 0 0  1  0  0  0  1  0  0  0  1  0  0  0  1  0  0  0
#> TypeT3      0 0 1 0 0 0 1 0 0  0  1  0  0  0  1  0  0  0  1  0  0  0  1  0  0
#> TypeT4      0 0 0 1 0 0 0 1 0  0  0  1  0  0  0  1  0  0  0  1  0  0  0  1  0
#>             26 27 28 29 30 31 32 33 34 35 36
#> (Intercept)  1  1  1  1  1  1  1  1  1  1  1
#> TypeT2       1  0  0  0  1  0  0  0  1  0  0
#> TypeT3       0  1  0  0  0  1  0  0  0  1  0
#> TypeT4       0  0  1  0  0  0  1  0  0  0  1
#> attr(,"assign")
#> [1] 0 1 1 1
#> attr(,"contrasts")
#> attr(,"contrasts")$Type
#> [1] "contr.treatment"
#>