R/clme.em.r

In jelsema/CLME: Constrained Inference for Linear Mixed Effects Models

#' @title Constrained EM algorithm for linear fixed or mixed effects models.
#'
#' @description \code{clme_em} is the general function, it will call the others.
#' These Expectation-maximization (EM) algorithms estimate model parameters and 
#' compute a test statistic.
#'
#' @rdname clme_em
#'
#' @param Y \eqn{N \times 1}{Nx1} vector of response data.
#' @param X1 \eqn{N \times p_1}{Nxp1} design matrix.
#' @param X2 optional \eqn{N \times p_2}{Nxp2} matrix of covariates.
#' @param U optional \eqn{N \times c}{Nxc} matrix of random effects.
#' @param Nks optional \eqn{K \times 1}{Kx1} vector of group sizes. 
#' @param Qs optional \eqn{Q \times 1}{Qx1} vector of group sizes for random effects.
#' @param constraints list containing the constraints. See Details. 
#' @param mq.phi optional MINQUE estimates of variance parameters.
#' @param tsf function to calculate the test statistic. 
#' @param tsf.ind function to calculate the test statistic for individual constrats. See Details for further information.
#' @param mySolver solver to use in isotonization (passed to \code{activeSet}). 
#' @param em.eps criterion for convergence for the EM algorithm. 
#' @param em.iter maximum number of iterations permitted for the EM algorithm. 
#' @param all_pair logical, whether all pairwise comparisons should be considered (constraints will be ignored).
#' @param dvar fixed values to replace bootstrap variance of 0.
#' @param verbose if \code{TRUE}, function prints messages on progress of the EM algorithm.
#' @param ... space for additional arguments.
#'
#' @details 
#' Argument \code{constraints} is a list including at least the elements \code{A}, \code{B}, and \code{Anull}. This argument can be generated by function \code{\link{create.constraints}}.
#'
#' @return
#' The function returns a list with the elements:
#' \itemize{
#' \item{\code{theta}}{ coefficient estimates.}
#' \item{\code{theta.null}}{ vector of coefficient estimates under the null hypothesis.}
#' \item{\code{ssq}}{ estimate of residual variance term(s).}
#' \item{\code{tsq}}{ estimate of variance components for any random effects.}
#' \item{\code{cov.theta}}{ covariance matrix of the unconstrained coefficients. }
#' \item{\code{ts.glb}}{ test statistic for the global hypothesis.}
#' \item{\code{ts.ind}}{ test statistics for each of the constraints.}
#' \item{\code{mySolver}}{ the solver used for isotonization.}
#' }
#' 
#' @note
#' There are few error catches in these functions. If only the EM estimates are desired, 
#' users are recommended to run \code{\link{clme}} setting \code{nsim=0}.
#' 
#' By default, homogeneous variances are assumed for the residuals and (if included) 
#' random effects. Heterogeneity can be induced using the arguments \code{Nks} and \code{Qs},
#'  which refer to the vectors \eqn{ (n_{1}, n_{2}, \ldots, n_{k}) }{(n1, n2 ,... , nk)} and
#'   \eqn{ (c_{1}, c_{2}, \ldots, c_{q}) }{(c1, c2 ,... , cq)}, respectively. See 
#'   \code{\link{CLME-package}} for further explanation the model and these values.
#' 
#' See \code{\link{w.stat}} and \code{\link{lrt.stat}} for more details on using custom 
#' test statistics.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' \code{\link{create.constraints}}
#' \code{\link{lrt.stat}}
#' \code{\link{w.stat}}
#' 
#' @examples
#' data( rat.blood )
#' 
#' model_mats <- model_terms_clme( mcv ~ time + temp + sex + (1|id), data = rat.blood )
#' 
#' Y  <- model_mats$Y
#' X1 <- model_mats$X1
#' X2 <- model_mats$X2
#' U  <- model_mats$U
#' 
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' 
#' clme.out <- clme_em(Y = Y, X1 = X1, X2 = X2, U = U, constraints = cons)
#' 
#' @importFrom isotone activeSet
#' @importFrom isotone gpava
#' @export
#' 
clme_em   <- function( Y, X1, X2 = NULL, U = NULL, Nks = nrow(X1),
                     Qs = ncol(U), constraints, mq.phi = NULL, tsf = lrt.stat, 
                     tsf.ind = w.stat.ind, mySolver="LS", em.iter = 500, 
                     em.eps = 0.0001, all_pair = FALSE, dvar = NULL, verbose = FALSE, ... ){
  
  
  ##
  ## Development plans: 
  ## - make function inputs more flexible like clme()
  ## - this may make passing the arguments need to be more detailed, but
  ##   passing arguments TO clme.em (e.g. from clme() ) would be less detailed.
  ##
  em_call  <- as.list( environment() )  
  dots     <- as.list(substitute(list(...)))[-1L]
  new_call <- append( em_call, dots )
  
  
  if( is.null(U) ){
    em_results <- do.call( "clme_em_fixed" , new_call )
  } else{
    em_results <- do.call( "clme_em_mixed" , new_call )
  }
  
  
  
  
  
  return( em_results )

}