# Get the first two-order smoothing moments from a fitted Kalman Reaction Network. This function returns the first two-order smoothing moments from a Kalman Reaction Network previously fitted on a clonal tracking dataset. ```r get.sMoments(res.fit, X = NULL, cell.cols = NULL) ``` ## Arguments - `res.fit`: A list returned by get.fit() containing the information of a fitted Kalman Reaction Network. - `X`: Stochastic process. A 3-dimensional array whose dimensions are the time, the cell type and the clone respectively. - `cell.cols`: Color legend for the cell types. Defaults to NULL, in which case no color legend for the cell types is provided. ## Returns No return value. ## Examples ```r rcts <- c("HSC->T", ## reactions "HSC->M", "T->0", "M->0") cnstr <- c("theta\\[\\'HSC->T\\'\\]=(theta\\[\\'T->0\\'\\])", "theta\\[\\'HSC->M\\'\\]=(theta\\[\\'M->0\\'\\])") latsts <- "HSC" ## latent cell types ctps <- unique(setdiff(c(sapply(rcts, function(r){ ## all cell types as.vector(unlist(strsplit(r, split = "->", fixed = TRUE))) }, simplify = "array")), c("0", "1"))) Y0 <- Y_CT$WAS[,setdiff(ctps,"HSC"),] topClones <- 2 Y0 <- Y0[,,names(head(sort(apply(Y0!=0, 3, sum), decreasing = TRUE), topClones)),drop=FALSE] ## cluster parameters: cl <- parallel::makeCluster(2, type = "PSOCK") ## initial condition: X0 <- rep(0, length(ctps)) names(X0) <- ctps X0["HSC"] <- 1 ## mean vector and covariance matrix of X0: m_0 <- replicate(dim(Y0)[3], X0, simplify = "array") colnames(m_0) <- dimnames(Y0)[[3]] P_0 <- Matrix::Diagonal(length(ctps) * dim(Y0)[3], 10) rownames(P_0) <- colnames(P_0) <- rep(dimnames(Y0)[[3]], each = length(ctps)) ## fit Karen on data: res.fit <- get.fit(rct.lst = rcts, constr.lst = cnstr, latSts.lst = latsts, ct.lst = ctps, Y = Y0, m0 = m_0, P0 = P_0, cl = cl, list(nLQR = 1, lmm = 0, ## needs to be >=5 for real applications pgtol = 0, relErrfct = 1e-5, tol = 1e-3, maxit = 0, ## needs to be increased for real applications maxitEM = 1, ## needs to be increased for real applications trace = 1, verbose = TRUE, FORCEP = FALSE)) parallel::stopCluster(cl) oldpar <- par(no.readonly = TRUE) par(mar = c(5,5,2,2), mfrow = c(1,3)) get.sMoments(res.fit) par(oldpar) ```

get.sMoments function