library(magrittr)
[\begin{gathered} {\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{A_{32}}}&{{A_{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{A_{32}}}&{{A_{33}}} \end{array}} \right] = {\left[ {\begin{array}{{20}{c}} {{{A'}{00}}}&{{{A'}{10}}}&{{{A'}{20}}}&{{{A'}{30}}} \ 0&{{{A'}{11}}}&{{{A'}{21}}}&{{{A'}{31}}} \ 0&0&{{{A'}{22}}}&{{{A'}{32}}} \ 0&0&0&{{{A'}{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{A_{32}}}&{{A_{33}}} \end{array}} \right] \hfill \ = \left[ {\begin{array}{*{20}{c}} {{{A'}{00}}{A{00}} + {{A'}{10}}{A{10}} + {{A'}{20}}{A{20}} + {{A'}{30}}{A{30}}}&{}&{}&{} \ {{{A'}{11}}{A{10}} + {{A'}{21}}{A{20}} + {{A'}{31}}{A{30}}}&{{{A'}{11}}{A{11}} + {{A'}{21}}{A{21}} + {{A'}{31}}{A{31}}}&{}&{} \ {{{A'}{22}}{A{20}} + {{A'}{32}}{A{30}}}&{{{A'}{22}}{A{21}} + {{A'}{32}}{A{31}}}&{{{A'}{22}}{A{22}} + {{A'}{32}}{A{32}}}&{} \ {{{A'}{33}}{A{30}}}&{{{A'}{33}}{A{31}}}&{{{A'}{33}}{A{32}}}&{{{A'}{33}}{A{33}}} \end{array}} \right] \hfill \ \end{gathered} ]
[\begin{gathered} {\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right] = {\left[ {\begin{array}{{20}{c}} {{{A'}{00}}}&{{{A'}{10}}}&{{{A'}{20}}}&{{{A'}{30}}} \ 0&{{{A'}{11}}}&{{{A'}{21}}}&{{{A'}{31}}} \ 0&0&{{{A'}{22}}}&{{{0'}{32}}} \ 0&0&0&{{{A'}{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{A_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right] \hfill \ = \left[ {\begin{array}{*{20}{c}} {{{A'}{00}}{A{00}} + {{A'}{10}}{A{10}} + {{A'}{20}}{A{20}} + {{A'}{30}}{A{30}}}&{}&{}&{} \ {{{A'}{11}}{A{10}} + {{A'}{21}}{A{20}} + {{A'}{31}}{A{30}}}&{{{A'}{11}}{A{11}} + {{A'}{21}}{A{21}} + {{A'}{31}}{A{31}}}&{}&{} \ {{{A'}{22}}{A{20}} + {{0'}{32}}{A{30}}}&{{{A'}{22}}{A{21}} + {{0'}{32}}{A{31}}}&{{{A'}{22}}{A{22}} + {{0'}{32}}{0{32}}}&{} \ {{{A'}{33}}{A{30}}}&{{{A'}{33}}{A{31}}}&{{{A'}{33}}{0{32}}}&{{{A'}{33}}{A{33}}} \end{array}} \right] \hfill \ \end{gathered} ]
[\begin{gathered} {\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{0_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{0_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right] = {\left[ {\begin{array}{{20}{c}} {{{A'}{00}}}&{{{A'}{10}}}&{{{A'}{20}}}&{{{A'}{30}}} \ 0&{{{A'}{11}}}&{{{0'}{21}}}&{{{A'}{31}}} \ 0&0&{{{A'}{22}}}&{{{0'}{32}}} \ 0&0&0&{{{A'}{33}}} \end{array}} \right]^\prime }\left[ {\begin{array}{{20}{c}} {{A_{00}}}&0&0&0 \ {{A_{10}}}&{{A_{11}}}&0&0 \ {{A_{20}}}&{{0_{21}}}&{{A_{22}}}&0 \ {{A_{30}}}&{{A_{31}}}&{{0_{32}}}&{{A_{33}}} \end{array}} \right] \hfill \ = \left[ {\begin{array}{*{20}{c}} {{{A'}{00}}{A{00}} + {{A'}{10}}{A{10}} + {{A'}{20}}{A{20}} + {{A'}{30}}{A{30}}}&{}&{}&{} \ {{{A'}{11}}{A{10}} + {{0'}{21}}{A{20}} + {{A'}{31}}{A{30}}}&{{{A'}{11}}{A{11}} + {{0'}{21}}{0{21}} + {{A'}{31}}{A{31}}}&{}&{} \ {{{A'}{22}}{A{20}} + {{0'}{32}}{A{30}}}&{{{A'}{22}}{0{21}} + {{0'}{32}}{A{31}}}&{{{A'}{22}}{A{22}} + {{0'}{32}}{0{32}}}&{} \ {{{A'}{33}}{A{30}}}&{{{A'}{33}}{A{31}}}&{{{A'}{33}}{0{32}}}&{{{A'}{33}}{A{33}}} \end{array}} \right] \hfill \ \end{gathered} ]
[\begin{gathered} \left[ {\begin{array}{{20}{c}} {{a_{11}}}&{{a_{21}}}&{{a_{31}}}&{{a_{41}}}&{{a_{51}}} \ 0&{{a_{22}}}&{{a_{32}}}&{{a_{42}}}&{{a_{52}}} \ 0&0&{{a_{33}}}&{{a_{43}}}&{{a_{53}}} \ 0&0&0&{{a_{44}}}&{{a_{54}}} \ 0&0&0&0&{{a_{55}}} \end{array}} \right]\left[ {\begin{array}{{20}{c}} {{a_{11}}}&0&0&0&0 \ {{a_{21}}}&{{a_{22}}}&0&0&0 \ {{a_{31}}}&{{a_{32}}}&{{a_{33}}}&0&0 \ {{a_{41}}}&{{a_{42}}}&{{a_{43}}}&{{a_{44}}}&0 \ {{a_{51}}}&{{a_{52}}}&{{a_{53}}}&{{a_{54}}}&{{a_{55}}} \end{array}} \right] \hfill \ = \left[ {\begin{array}{*{20}{c}} {a_{11}^2 + a_{21}^2 + a_{31}^2 + a_{41}^2 + a_{51}^2}&{}&{}&{}&{} \ {{a_{22}}{a_{21}} + {a_{32}}{a_{31}} + {a_{42}}{a_{41}} + {a_{52}}{a_{51}}}&{a_{22}^2 + a_{32}^2 + a_{42}^2 + a_{52}^2}&{}&{}&{} \ {{a_{33}}{a_{31}} + {a_{43}}{a_{41}} + {a_{53}}{a_{51}}}&{{a_{33}}{a_{32}} + {a_{43}}{a_{42}} + {a_{53}}{a_{52}}}&{a_{33}^2 + a_{43}^2 + a_{53}^2}&{}&{} \ {{a_{44}}{a_{41}} + {a_{54}}{a_{51}}}&{{a_{44}}{a_{42}} + {a_{54}}{a_{52}}}&{{a_{44}}{a_{43}} + {a_{54}}{a_{53}}}&{a_{44}^2 + a_{54}^2}&{} \ {{a_{55}}{a_{51}}}&{{a_{55}}{a_{52}}}&{{a_{55}}{a_{53}}}&{{a_{55}}{a_{54}}}&{a_{55}^2} \end{array}} \right] \hfill \ \end{gathered} ]
Conjecture:
Consider a random effects model with a 3-level categorical factor, $F$, taking the values 0,1 and 2, and a continuous variable, $x$. The random effect assuming independence between the random effect associated with $F$ and that associated with $X$, can be written as:
$$ \delta = u_0 + f_1 u_1 + f_2 u_2 + x u_3 $$ where $f_1 = I_{F=1}$, $f_2 = I_{F=2}$, and
[\operatorname{Var} \left( {\begin{array}{{20}{c}} {{u_0}} \ {{u_1}} \ {{u_3}} \ {{u_4}} \end{array}} \right) = \left[ {\begin{array}{{20}{c}} {{g_{00}}}&{{g_{01}}}&{{g_{02}}}&{{g_{03}}} \ {{g_{10}}}&{{g_{11}}}&{{g_{12}}}&0 \ {{g_{20}}}&{{g_{21}}}&{{g_{22}}}&0 \ {{g_{30}}}&0&0&{{g_{33}}} \end{array}} \right]]
Questions to explore are verifying the coding and location invariance of the model.
[\begin{gathered} \operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{cov} (\delta |F = 2,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \end{gathered} ]
From this we see the following variance properties for $\delta$: [\begin{gathered} \operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \operatorname{cov} (\delta |F = 1,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \hfill \ \end{gathered} ]
The Log-Cholesky factorization used for the variance-covariance matrix of random effects in nlme
does not lend itself easily to selectively modeling some random effects as independent. This stems from the fact that the factorization is a 'right' factorization in which the variance matrix is expressed as:
$$
V = R'R
$$
where $R$ is an upper-triangular (or 'right-triangular') matrix.
Constraining off-diagonal upper-triangle elements of $R$ to be zero does not result in the corresponding element of $V$ being zero. That is, it does not result in a zero covariance. It produces, however, a zero conditional covariance.
TO DO: Demonstrate above.
On the other hand, a left-triangular factorization: $$ V = L'L $$ will have the property that 0 elements of $L$ in the lower triangle will result in 0 elements of $V$ in the corresponding position.
The chol
function in returns a right factorization. Note that the use of the pivot
argument to handle rank-deficient matrices produces a factor that is not triangular if a non-trivial pivot is applied. As a result, pivotting cannot be used to handle singular variances since the assumption of triangular parametrization would be violated.
(V <- diag(1:3) + 1) (f <- chol(V)) t(f) %*% f f[2,3] <- 0 f (v2 <- t(f) %*% f) # The conditional covariance is zero v2[2:3,2:3] - v2[2:3,1,drop=FALSE] %*% solve(v2[1,1],drop=FALSE) %*% v2[1,2:3,drop=FALSE] # How to do a left-factorization cholR <- function(x, ...) { # this can return a non-triangular factor, thus can't be used for our purposes ret <- chol(x, pivot = TRUE) pivot <- attr(ret, "pivot") ret[, order(pivot)] } cholL <- function(x, ...) { x[,] <- rev(x) ret <- chol(x) # not cholR (see above) ret[,] <- rev(ret) ret } # test rank deficient V V0 <- crossprod( cbind(1, 1:10, (1:10)^2, 11:20)) tryCatch(chol(V0), error = function(e) e) tryCatch(cholR(V0), error = function(e) e) lfac <- cholL(V) t(lfac) %*% lfac lfac[3,2] <- 0 t(lfac) %*% lfac
An inefficient, and probably numerically undesirable, way would construct the variance matrix and perform a left factorization: cholL( t(R) %*% R)
.
Instead, we use the qr
decomposition of R
, avoiding computing the variance matrix and the ensuing loss of precision.
(lfac <- cholL(V)) rfac <- qr.R(qr(lfac)) chol(V) t(rfac) %*% rfac - (t(lfac) %*% lfac) X <- cbind(1,1:10,(1:10)^2,11:20) qr.R(qr(X)) qr.R(qr(X, pivot = T)) L2R <- function(x, ...) { R <- qr.R(qr(x)) sign(diag(R)) * R } V cholL(V) %>% round(2) V - crossprod( cholL(V)) V - crossprod( chol(V)) V - crossprod( L2R(cholL(V))) R2L <- function(x, ...) { x[] <- rev(x) ret <- L2R(x, ...) ret[] <- rev(ret) ret } x <- matrix(rnorm(10000), 100) system.time( for(i in 1:1000) R2L(x) ) system.time( for(i in 1:1000) L2R(x) ) # Check that it works V Vr <- crossprod( matrix(rnorm(9),3)) svd(Vr)$d lc <- cholL(Vr) Vr - crossprod(lc) rc <- chol(Vr) Vr - crossprod(rc) lc2r <- L2R(lc) Vr - crossprod(lc2r) rc2l <- R2L(rc) Vr - crossprod(rc2l)
Seems that:
Note the following:
logChol_pd
in pdFactor.pdLogChol
returns a list of three components:crossprod(pdMatrix(object,factor=TRUE))
pdFactor
is the critical method and we need a pdFactor.pdInd
. The solve method transforms a pdInd object representing a variance matrix to a pdInd object representing the inverse of that matrix. If we let $L$ be a left-triangular factorization of $V = L'L$ then, if $F$ is the left-triangular factorization of $V^{-1}$, we have that $ F = R'^{-1}$ where $R$ is the right-triangular factor such that $V = R'R$.
Thus, solve.pdInd
would transform $L$ to $R'^{-1}$.
To take advantage of the structure of a pdInd
object, which is a factorization of a variance $V$ of the form:
[V = \left[ {\begin{array}{*{20}{c}} A&{B'} \ B&C \end{array}} \right]]
where $C$ is diagonal. In this case, $L$ has the form:
[L = \left[ {\begin{array}{*{20}{c}} D&0 \ E&F \end{array}} \right]]
conformably with the partitioning of $V$ with $F$ a diagonal matrix.
Now an important problem, in contrast with existing methods for pdMat objects, is that the left-triangular factor of the inverse cannot be expressed with the unconstrained parametrization for $V$.
If the result of solve
is only used through pdFactor
or pdMatrix
, then it is a simple thing to set a flag for the pdInd
object indicating that it represents the inverse. We can give pdFactor
and pdMatrix
the job of computing the correct right-triangular factor.
If the object created by solve.pdInd
is used in other ways then the problem might be more difficult.
One hopes that the approach of using solve
to create the representation of $V^{-1}$ arose simply from the fact that the pdMat
classes in the original version of nlme
all lent themselves to this approach.
Here are some relevant expressions:
[\begin{gathered} V = \left[ {\begin{array}{{20}{c}} A&{B'} \ B&C \end{array}} \right] \hfill \ L = \left[ {\begin{array}{{20}{c}} D&0 \ E&F \end{array}} \right] \hfill \ \end{gathered} ]
[\left[ {\begin{array}{{20}{c}} A&{B'} \ B&C \end{array}} \right] = \left[ {\begin{array}{{20}{c}} {D'}&{E'} \ 0&F \end{array}} \right]\left[ {\begin{array}{{20}{c}} D&0 \ E&F \end{array}} \right] = \left[ {\begin{array}{{20}{c}} {D'D + E'E}&{E'F} \ {FE}&{{F^2}} \end{array}} \right]]
[F = {C^{1/2}}]
[E = {F^{ - 1}}E]
[D'D = A - E'E]
[R = \left[ {\begin{array}{*{20}{c}} G&H \ 0&K \end{array}} \right]]
[V = R'R = \left[ {\begin{array}{{20}{c}} A&{B'} \ B&C \end{array}} \right] = \left[ {\begin{array}{{20}{c}} {G'}&0 \ {H'}&{K'} \end{array}} \right]\left[ {\begin{array}{{20}{c}} G&H \ 0&K \end{array}} \right] = \left[ {\begin{array}{{20}{c}} {G'G}&{G'H} \ {H'G}&{H'H + K'K} \end{array}} \right]]
The problem arises from the fact that, in the case of the right-triangular factorization, the diagonality of $C$ does not imply the diagonality of $K$, which is precisely why we need to use the left-triangular factorization in the first place.
Thus, when we invert and transpose $R$ to get the left-triangular factor of $V^{-1}$ we cannont use the parsimonious unrestricted parametrization of the pdInd
object.
To get around this, we will attempt to use the following strategy:
solve.pdInd
simply inverts the flag (from FALSE to TRUE or TRUE to FALSE).pdFactor.pdInd
returns a right-triangular factor for $V$ or for $V^{-1}$ depending on the value of the 'invert' attribute.Thus the burden of dealing with inversion is passed along to pdFactor.pdInd
.
library(nlme) methods(class="pdLogChol") methods(class="pdSymm") nlme:::coef.pdSymm methods(class="pdMat") nlme:::pdLogChol nlme:::pdSymm nlme:::pdMat pdFactorfull <- function (object) { Ncol <- round((-1 + sqrt(1 + 8 * length(object)))/2) .C(logChol_pd, Factor = double(Ncol * Ncol), as.integer(Ncol), as.double(object)) } environment(pdFactorfull) <- environment(pdFactor)
Vignettes are long form documentation commonly included in packages. Because they are part of the distribution of the package, they need to be as compact as possible. The html_vignette
output type provides a custom style sheet (and tweaks some options) to ensure that the resulting html is as small as possible. The html_vignette
format:
Note the various macros within the vignette
setion of the metadata block above. These are required in order to instruct R how to build the vignette. Note that you should change the title
field and the \VignetteIndexEntry
to match the title of your vignette.
The html_vignette
template includes a basic CSS theme. To override this theme you can specify your own CSS in the document metadata as follows:
output: rmarkdown::html_vignette: css: mystyles.css
The figure sizes have been customised so that you can easily put two images side-by-side.
plot(1:10) plot(10:1)
You can enable figure captions by fig_caption: yes
in YAML:
output: rmarkdown::html_vignette: fig_caption: yes
Then you can use the chunk option fig.cap = "Your figure caption."
in knitr.
You can write math expressions, e.g. $Y = X\beta + \epsilon$, footnotes^[A footnote here.], and tables, e.g. using knitr::kable()
.
knitr::kable(head(mtcars, 10))
Also a quote using >
:
"He who gives up [code] safety for [code] speed deserves neither." (via)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.