In gmonette/Tcells: Fits Tcell dynamic models with mixed models
library(magrittr)
library(Tcells)
Triangular factorizations
In the 'lme' function in the 'nlme' package, the parametrization of the positive-definite variance matrix for random effects is achieved by modelling a factor of the inverse variance matrix: i.e. a matrix $A$ such that $\sigma^2 G^{-1} = A'A$.
The process is made easier if the model lends itself to an
uncontrained parametrization of $A$. For example, a model in which $G$ is any positive-definite diagonal matrix can be parametrized
by setting
[A = \left[ {\begin{array}{*{20}{c}}
{\sqrt {{g_{00}}} }&0&0&0 \
0&{\sqrt {{g_{11}}} }&0&0 \
0&0&{\sqrt {{g_{22}}} }&0 \
0&0&0&{\sqrt {{g_{33}}} }
\end{array}} \right]]
Using $phi_i = log( \sqrt(g_{ii}) ), i = 1,...,3$ as optimization parameters provides an unconstrained parametrization of the model.
Although constrained parametrizations can be used in 'lme', unconstrained parametrizations are, generally, more effective.
If $G$ can be any positive-definite matrix, one unconstrained parametrization is based on the Choleski decomposition: [{\sigma ^2}{G^{ - 1}} = R'R]
where $R$ is an upper-triangular matrix that can be chosen so that the diagonal elements are positive, since $G$ is positive definite.
The unconstrained parametrization uses the log of the diagonal elements of $R$ and the raw values of the upper-triangle of $R$.
Only a few models are provided in 'lme': a diagonal variance, a variance that is a multiple of the identity, full positive-definite matrices using a either a log Choleski parametrization or a 'log spectral' parametrization and, finally, block diagonal arrangements of the former.
Since most random effects models include a random intercept, the previous models do not include a sitution in which random effects are independent of each other but not of the random intercept. Assuming independence from the random intercept results in a model that lacks location invariance and that forces the location of minimal variance of the response surface to lie, arbitrarily, over the origin.
[\begin{aligned}
{\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] \
= \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] \
\end{aligned} ]
[\begin{aligned}
{\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] \
= \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] \
\end{aligned} ]
[\begin{aligned}
{\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] \
= \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] \
\end{aligned} ]
[\begin{aligned}
\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] \
= \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] \
\end{aligned} ]
What patterns of zeros in $L$ carry over to $V$?
Conjecture:
- Any pattern of zeros in which each column of zeros extends to the bottom row. Thus any subdiagonal rectangular block or a subdiagonal L-shape that extends to the bottom row will reproduce in the variance matrix.
- Two zeros (or blocks of zeros) in the first sub-diagonal (or first subdiagonal partition).
Models for categorical factors
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{aligned}
\operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{cov} (\delta |F = 2,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \
\end{aligned} ]
From this we see the following variance properties for $\delta$:
[\begin{aligned}
\operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \
\operatorname{cov} (\delta |F = 1,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \
\end{aligned} ]
Variance Factorizations
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
Is there an efficient way to turn a right factor into a left factor?
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)
Creating a pdClass based on pdLogChol
Seems that:
- we need pdFactor.pdInd which is called by pdMatrix etc.
- we need solve.pdInd which returns the 'pdInd' object of the inverse (we can take advantage of a few things here)
- we might not need coef.pdInd since coef.pdMat would work, but perhaps might as well
- we need, of course, pdInd and, I'm pretty sure, pdConstruct.pdInd
Note the following:
- solve.pdInd does not construct an inverse, just the pdInd representation of the inverse.
- So, where is the inverse right-factor constructed -- which I assume is needed for the ML? Perhaps pdFactor, which is exported by nlme.
- The C function
logChol_pd in pdFactor.pdLogChol returns a list of three components:
- The R factor as a vector, listed column by column
- The number of rows in R, as a single number
- The logChol coeficients: log of diagonal elements of R followed by off-diagonal elements
- A call to pdMatrix for a pdLogChol/pdSymm/pdMat object calls the 'pdMatrix.pdSymm' method that in turn calls 'pdFactor' if the argument 'factor' is TRUE. Otherwise, it calls pdMatrix.pdMat that calls
crossprod(pdMatrix(object,factor=TRUE))
- So, bottom line,
pdFactor is the critical method and we need a pdFactor.pdInd.
solve.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{aligned}
V = \left[ {\begin{array}{{20}{c}}
A&{B'} \
B&C
\end{array}} \right] \
L = \left[ {\begin{array}{{20}{c}}
D&0 \
E&F
\end{array}} \right] \
\end{aligned} ]
[\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:
- Add an 'invert' attribute to 'pdInd' object. Intially, it is set to FALSE.
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)
Examples
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:
- Never uses retina figures
- Has a smaller default figure size
- Uses a custom CSS stylesheet instead of the default Twitter Bootstrap style
Vignette Info
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.
Styles
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
Figures
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.
More Examples
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)
gmonette/Tcells documentation built on May 17, 2019, 7:25 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
library(magrittr) library(Tcells)
Triangular factorizations
In the 'lme' function in the 'nlme' package, the parametrization of the positive-definite variance matrix for random effects is achieved by modelling a factor of the inverse variance matrix: i.e. a matrix $A$ such that $\sigma^2 G^{-1} = A'A$.
The process is made easier if the model lends itself to an uncontrained parametrization of $A$. For example, a model in which $G$ is any positive-definite diagonal matrix can be parametrized by setting [A = \left[ {\begin{array}{*{20}{c}} {\sqrt {{g_{00}}} }&0&0&0 \ 0&{\sqrt {{g_{11}}} }&0&0 \ 0&0&{\sqrt {{g_{22}}} }&0 \ 0&0&0&{\sqrt {{g_{33}}} } \end{array}} \right]] Using $phi_i = log( \sqrt(g_{ii}) ), i = 1,...,3$ as optimization parameters provides an unconstrained parametrization of the model.
Although constrained parametrizations can be used in 'lme', unconstrained parametrizations are, generally, more effective.
If $G$ can be any positive-definite matrix, one unconstrained parametrization is based on the Choleski decomposition: [{\sigma ^2}{G^{ - 1}} = R'R] where $R$ is an upper-triangular matrix that can be chosen so that the diagonal elements are positive, since $G$ is positive definite.
The unconstrained parametrization uses the log of the diagonal elements of $R$ and the raw values of the upper-triangle of $R$.
Only a few models are provided in 'lme': a diagonal variance, a variance that is a multiple of the identity, full positive-definite matrices using a either a log Choleski parametrization or a 'log spectral' parametrization and, finally, block diagonal arrangements of the former.
Since most random effects models include a random intercept, the previous models do not include a sitution in which random effects are independent of each other but not of the random intercept. Assuming independence from the random intercept results in a model that lacks location invariance and that forces the location of minimal variance of the response surface to lie, arbitrarily, over the origin.
[\begin{aligned} {\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] \ = \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] \ \end{aligned} ]
[\begin{aligned} {\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] \ = \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] \ \end{aligned} ]
[\begin{aligned} {\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] \ = \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] \ \end{aligned} ]
[\begin{aligned} \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] \ = \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] \ \end{aligned} ]
What patterns of zeros in $L$ carry over to $V$?
Conjecture:
- Any pattern of zeros in which each column of zeros extends to the bottom row. Thus any subdiagonal rectangular block or a subdiagonal L-shape that extends to the bottom row will reproduce in the variance matrix.
- Two zeros (or blocks of zeros) in the first sub-diagonal (or first subdiagonal partition).
Models for categorical factors
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{aligned} \operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{cov} (\delta |F = 2,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \ \end{aligned} ]
From this we see the following variance properties for $\delta$: [\begin{aligned} \operatorname{var} (\delta |F = 0) = {g_{00}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{var} (\delta |F = 1) = {g_{00}} + 2{g_{01}} + {g_{11}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{var} (\delta |F = 2) = {g_{00}} + 2{g_{02}} + {g_{22}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{cov} (\delta |F = 0,\delta |F = 1) = {g_{00}} + {g_{01}} + 2{g_{03}}x + {g_{33}}{x^2} \ \operatorname{cov} (\delta |F = 1,\delta |F = 1) = {g_{00}} + {g_{01}} + {g_{20}} + {g_{21}} + 2{g_{03}}x + {g_{33}}{x^2} \ \end{aligned} ]
Variance Factorizations
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
Is there an efficient way to turn a right factor into a left factor?
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)
Creating a pdClass based on pdLogChol
Seems that:
- we need pdFactor.pdInd which is called by pdMatrix etc.
- we need solve.pdInd which returns the 'pdInd' object of the inverse (we can take advantage of a few things here)
- we might not need coef.pdInd since coef.pdMat would work, but perhaps might as well
- we need, of course, pdInd and, I'm pretty sure, pdConstruct.pdInd
Note the following:
- solve.pdInd does not construct an inverse, just the pdInd representation of the inverse.
- So, where is the inverse right-factor constructed -- which I assume is needed for the ML? Perhaps pdFactor, which is exported by nlme.
- The C function
logChol_pdinpdFactor.pdLogCholreturns a list of three components: - The R factor as a vector, listed column by column
- The number of rows in R, as a single number
- The logChol coeficients: log of diagonal elements of R followed by off-diagonal elements
- A call to pdMatrix for a pdLogChol/pdSymm/pdMat object calls the 'pdMatrix.pdSymm' method that in turn calls 'pdFactor' if the argument 'factor' is TRUE. Otherwise, it calls pdMatrix.pdMat that calls
crossprod(pdMatrix(object,factor=TRUE)) - So, bottom line,
pdFactoris the critical method and we need apdFactor.pdInd.
solve.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{aligned} V = \left[ {\begin{array}{{20}{c}} A&{B'} \ B&C \end{array}} \right] \ L = \left[ {\begin{array}{{20}{c}} D&0 \ E&F \end{array}} \right] \ \end{aligned} ]
[\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:
- Add an 'invert' attribute to 'pdInd' object. Intially, it is set to FALSE.
solve.pdIndsimply inverts the flag (from FALSE to TRUE or TRUE to FALSE).pdFactor.pdIndreturns 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)
Examples
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:
- Never uses retina figures
- Has a smaller default figure size
- Uses a custom CSS stylesheet instead of the default Twitter Bootstrap style
Vignette Info
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.
Styles
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
Figures
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.
More Examples
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)
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.