parameter_restriction: Restrict the parameters space of OU and Brownian motion models.

Description

ou_diagH, ou_diagH_fixedtheta_diagSig, etc., restricts the OU model's parameters. For example, ou_diagH restricts the drift $H$ to diagonal matrix, and ou_diagH_fixedtheta_diagSig further restricts theta to be a constant and $\Sigma_x'$ to be diagonal. A Brownian motion model can be made by these restriction.

Usage

avail_restrictions
brn_diagSig(parfn)
ou_logdiagH(parfn)
dou_logdiagH(jacfn)
hou_logdiagH(hessfn)
ou_logdiagH_diagSig(parfn)
ou_logspdH_fixedtheta(parfn, theta)
ou_spdH_fixedSig(parfn, Sig)
ou_fixedH_diagSig(parfn, H)
dou_logdiagH_diagSig(jacfn)
dou_logspdH_fixedtheta(jacfn, theta)
dou_spdH_fixedSig(jacfn, Sig)
dou_fixedH_diagSig(jacfn, H)
hou_logdiagH_diagSig(hessfn)
hou_logspdH_fixedtheta(hessfn, jacfn, theta)
hou_spdH_fixedSig(hessfn, jacfn, Sig)
hou_spdH_fixedtheta_fixedSig(hessfn, jacfn, theta, Sig)
hou_fixedH_diagSig(hessfn, H)
nparams_ou_logdiagH(k)
nparams_brn(k)
nparams_ou_spdH_fixedSig(k)

Format

An object of class list of length 4.

Arguments

parfn: A function that maps from the user-parametrisation to the underlying Gaussian parameters. Each of them returns a vector of concatenated $(\Phi, w, V')$, where $V'$ is the lower triangular part of $V$, and accepts four arguments: a vector of parameters whose length is specified by the pardims argument to the glinv_gauss function, the branch length leading to the currently processing node, a vector of factors with three levels indicating which dimensions are missing or lost in the mother of the current node, and a vector of factors with the same three levels indicating missingness of the current node.
jacfn: A function that accepts the same arguments as parfn and returns the Jacobian of parfn.
hessfn: A function that accepts the same arguments as parfns and returns a list of three 3D arrays, named Phi, w, V respectively inside the list. ((hessfn)(...))$Phi[m,i,j] contains the cross second-order partial derivative of $\Phi_m$ (here we treat the matrix $\Phi$ as a column-major-flattened vector) with respect to the $i$-th and$j$-th parameters in the joint $(H,\theta,\Sigma_x')$ vector, and ((hessfn)(...))$w[m,i,j] and ((hessfn)(...))$V[m,i,j] analogously contains second-order derivative with respect to $w_m$ and $V'_m$.
H: A numerical vector containing the (flattened) fixed parameter $H$.
theta: A numerical vector containing the (flattened) fixed parameter $theta$.
Sig: A numerical vector containing the (flattened) fixed parameter $\Sigma_x'$.
k: An integer. The total number of dimensions of the multivariate traits.

Details

How reparametrisation and restriction works

In the simplest form, without any restriction or reparametrisation, the user typically needs to pass oupar, oujac, ouhess, all of which are simply functions which maps from the OU parameters $(H,\theta,\Sigma_x')$ to the Gaussian paramters $(\Phi_i,w_i,V'_i)$ for each node. For example:


        mod.full = glinv(tree, x0, my_data,
                         parfns  = oupar,
                         pardims = nparams_ou(k),
                         parjacs = oujac,
                         parhess = ouhess)

If one would like to restrict $H$ to only positively definite diagonal matrices, then the call should become


        mod.pddiag = glinv(tree, x0, my_data,
                           parfns  = ou_logdiagH(oupar),
                           pardims = nparams_ou_logdiagH(k),
                           parjacs = dou_logdiagH(oujac),
                           parhess = hou_logdiagH(ouhess))

Note that there is a naming convention that ou_* should be applied to `oupar`, dou_* to `oujac`, and hou_* to `ouhess`. d stands for `derivative' and h stands for `Hessian'.

In the above call, ou_logdiagH(oupar) accepts the oupar function as argument and returns a new function. This new function behaves the same way as oupar itself, except that it expects its first argument (which is the model parameters) to be of lower dimension, only consisting of $(h,\theta,\Sigma_x')$ where $h$ is the diagonal vector of $H$. The following example should be illustrative:


        f = ou_logdiagH(oupar)
        par.full = list(H     = matrix(c(3,0,0,2),2,2), # diagonal matrix
                        theta = c(4,5),
                        sig_x = c(1,0.1,1))
        par.restricted = list(H     = log(diag(par.full$H)),
                              theta = par.full$theta,
                              sig_x = par.full$sig_x)
        print(all.equal(f(unlist(par.restricted),1,NULL,NULL),
                        oupar(unlist(par.full),1,NULL,NULL)))
        # [1] TRUE

Pre-defined restrictions

The following table summarises all the pre-defined ou_* functions. See oupar for precise meaning of the $(H,\theta,\Sigma_x')$ mentioned below.

R function	Parameter Format after Restriction
`brn*`	$\Sigma_x'$. The Brownian motion. $H$ and $\theta$ are zero, thus missing.
`_diagH_`	$(h,\theta,\Sigma_x')$, with $h=diag(H)$, and H is a diagonal matrix
`_logdiagH_`	$(log(h),\theta,\Sigma_x')$, with $h=diag(H)$, and H is a diagonal matrix
`_symH_`	$(L,\theta,\Sigma_x')$, with $L$ being lower-triangular part of the symmetric matrix $H$
`_spdH_`	$(L,\theta,\Sigma_x')$, with $L$ being Cholesky factor of the S.P.D. matrix $H$
`_logspdH_`	$(L',\theta,\Sigma_x')$ where $L'$ equals $L$, except that on the diagonals $L'_i$ = $log L_i$
`_fixedH_`	$(\theta,\Sigma_x')$. $H$ is constant, hence missing
`_fixedtheta_`	$(H,\Sigma_x')$. $\theta$ is constant, hence missing
`_fixedSig_`	$(H,\theta)$. $\Sigma_x$ is constant, hence missing
`_diagSig_`	$(H,\theta,s)$ where $s=diag(\Sigma_x'$, with $\Sigma_x'$ being a diagonal matrix.

By Cholesky factor, we mean the only the non-zero part of the lower-triangular Cholesky factor. Restricting $\Sigma_x'$ to a diagonal matrix means that $\Sigma_x$ is also diagonal; and the variance of the Brownian motion is $log(diag(\Sigma_x'))$. In other words, the diagonal restriction is placed on $\Sigma_x'$, not $\Sigma_x$.

Finding a list of these restriction functions

One can use print(avail_restrictions) to see a list of all of these restriction function names.

Calling these restriction functions

All *ou_* or *brn* functions accepts the same arguemnts as ou_logdiagH, dou_logdiagH, hou_logdiagH, nparams_ou_logdiagH as shown in the Usage and Arguments section, except that:

If the reparametrisation contains any Cholesky decomposition (in other words, the function name contains spd or logspd) then in the Hessian-level reparameterisation function (named hou_*) an extra argument jacfn is required.
If the reparametrisation contains any fixed parameters, extra arguments H, theta, or Sig are required, depending what is fixed.

For example, in the Usage section, ou_logspdH_fixedtheta takes an extra argument theta because of (2), and hou_spdH_fixedSig takes extra argument two extra arguments because of both (1) and (2) are true.

R function	Parameter Format after Restriction
`brn*`	\(\Sigma_x'\). The Brownian motion. \(H\) and \(\theta\) are zero, thus missing.
`_diagH_`	\((h,\theta,\Sigma_x')\), with \(h=diag(H)\), and H is a diagonal matrix
`_logdiagH_`	\((log(h),\theta,\Sigma_x')\), with \(h=diag(H)\), and H is a diagonal matrix
`_symH_`	\((L,\theta,\Sigma_x')\), with \(L\) being lower-triangular part of the symmetric matrix \(H\)
`_spdH_`	\((L,\theta,\Sigma_x')\), with \(L\) being Cholesky factor of the S.P.D. matrix \(H\)
`_logspdH_`	\((L',\theta,\Sigma_x')\) where \(L'\) equals \(L\), except that on the diagonals \(L'_i\) = \(log L_i\)
`_fixedH_`	\((\theta,\Sigma_x')\). \(H\) is constant, hence missing
`_fixedtheta_`	\((H,\Sigma_x')\). \(\theta\) is constant, hence missing
`_fixedSig_`	\((H,\theta)\). \(\Sigma_x\) is constant, hence missing
`_diagSig_`	\((H,\theta,s)\) where \(s=diag(\Sigma_x'\), with \(\Sigma_x'\) being a diagonal matrix.