fit.re {RestoreNet}

Fit the random-effects model

Description

This function builds the design matrix of the random-effects model and returns the fitted values and the corresponding statistics.

Usage

fit.re(
  theta_0,
  Y,
  rct.lst,
  maxit = 10000,
  factr = 1e+07,
  pgtol = 1e-08,
  lmm = 100,
  maxemit = 100,
  eps = 1e-05,
  trace = TRUE,
  verbose = TRUE
)

Arguments

theta_0

A p-dimensional vector parameter as the initial guess for the inference.

Y

A 3-dimensional array whose dimensions are the time, the cell type and the clone respectively.

rct.lst

list of biochemical reactions. A differentiation move from cell type "A" to cell type "B" must be coded as "A->B" Duplication of cell "A" must be coded as "A->1" Death of cell "A" must be coded as "A->0"

maxit

maximum number of iterations for the optimization step. This argument is passed to optim() function. Details on "maxit" can be found in "optim()" documentation page.

factr

controls the convergence of the "L-BFGS-B" method. Convergence occurs when the reduction in the objective is within this factor of the machine tolerance. Default is 1e7, that is a tolerance of about 1e-8. This argument is passed to optim() function.

pgtol

helps control the convergence of the "L-BFGS-B" method. It is a tolerance on the projected gradient in the current search direction. This defaults to zero, when the check is suppressed. This argument is passed to optim() function.

lmm

is an integer giving the number of BFGS updates retained in the "L-BFGS-B" method, It defaults to 5. This argument is passed to optim() function.

maxemit

maximum number of iterations for the expectation-maximization algorithm.

eps

relative error for the value x and the objective function f(x) that has to be optimized in the expectation-maximization algorithm.

trace

Non-negative integer. If positive, tracing information on the progress of the optimization is produced. This parameter is also passed to the optim() function. Higher values may produce more tracing information: for method "L-BFGS-B" there are six levels of tracing. (To understand exactly what these do see the source code: higher levels give more detail.)

verbose

(defaults to TRUE) Logical value. If TRUE, then information messages on the progress of the algorithm are printed to the console.

Value

A 3-length list. First element is the output returned by "optim()" function (see "optim()" documentation for details) along with the conditional expectation E[u \vert y] and variance V[u \vert y] of the latent states u given the observed states y from the last step of the expectation-maximization algorithm. Second element is a vector of statistics associated to the fitted random-effects model:

nPar	number of parameters of the base(null) model
cll	value of the conditional log-likelihood, in this case just the log-likelihood
mll	value of the marginal log-likelihood, in this case just the log-likelihood
cAIC	conditional Akaike Information Criterion (cAIC), in this case simply the AIC.
mAIC	marginal Akaike Information Criterion (mAIC), in this case simply the AIC.
Chi2	value of the \chi^2 statistic (y - M\theta)'S^{-1}(y - M\theta).
p-value	p-value of the \chi^2 test for the null hypothesis that Chi2 follows a \chi^2 distribution with n - nPar degrees of freedom.
KLdiv	Kullback-Leibler divergence of the random-effects model from the null model.
KLdiv/N	Rescaled Kullback-Leibler divergence of the random-effects model from the null model.
BhattDist_nullCond	Bhattacharyya distance between the random-effects model and the null model.
BhattDist_nullCond/N	Rescaled Bhattacharyya distance between the random-effects model and the null model.

The third element, called "design", is a list including:

M	A n \times K dimensional (design) matrix.
M_bdiag	An \times Jp dimensional block-diagonal design matrix.
V	A p \times K dimensional net-effect matrix.

Examples

rcts <- c("A->1", "B->1", "C->1", "D->1",
          "A->0", "B->0", "C->0", "D->0",
          "A->B", "A->C", "C->D") ## set of reactions
ctps <- head(LETTERS,4)
nC <- 3 ## number of clones
S <- 10 ## trajectory length
tau <- 1 ## for tau-leaping algorithm
u_1 <- c(.2, .15, .17, .09*5,
         .001, .007, .004, .002,
         .13, .15, .08)
u_2 <- c(.2, .15, .17, .09,
         .001, .007, .004, .002,
         .13, .15, .08)
u_3 <- c(.2, .15, .17*3, .09,
         .001, .007, .004, .002,
         .13, .15, .08)
theta_allcls <- cbind(u_1, u_2, u_3) ## clone-specific parameters
rownames(theta_allcls) <- rcts
s20 <- 1 ## additional noise
Y <- array(data = NA,
           dim = c(S + 1, length(ctps), nC),
           dimnames = list(seq(from = 0, to = S*tau, by = tau),
                           ctps,
                           1:nC)) ## empty array to store simulations
Y0 <- c(100,0,0,0) ## initial state
names(Y0) <- ctps
for (cl in 1:nC) { ## loop over clones
  Y[,,cl] <- get.sim.tl(Yt = Y0,
                        theta = theta_allcls[,cl],
                        S = S,
                        s2 = s20,
                        tau = tau,
                        rct.lst = rcts,
                        verbose = TRUE)
}
null.res <- fit.null(Y = Y,
                     rct.lst = rcts,
                     maxit = 0, ## needs to be increased (>=100) for real applications
                     lmm = 0, ## needs to be increased (>=5) for real applications
) ## null model fitting

re.res <- fit.re(theta_0 = null.res$fit$par,
                 Y = Y,
                 rct.lst = rcts,
                 maxit = 0, ## needs to be increased (>=100) for real applications
                 lmm = 0, ## needs to be increased (>=5) for real applications
                 maxemit = 1 ## needs to be increased (>= 100) for real applications
) ## random-effects model fitting