R/utilities.r

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

##
## Various support functions to assist the main functions
##

#' Create model matrices for \code{clme}
#'
#' @description
#' Parses formulas to creates model matrices for \code{clme}. 
#'
#' @param formula a formula defining a linear fixed or mixed effects model. The constrained effect(s) must come before any unconstrained covariates on the right-hand side of the expression. The first \code{ncon} terms will be assumed to be constrained. 
#' @param data data frame containing the variables in the model.
#' 
#' 
#' @note
#' The first term on the right-hand side of the formula should be the fixed effect 
#' with constrained coefficients. Random effects are represented with a vertical bar, 
#' so for example the random effect \code{U} would be included by 
#' \code{Y ~ X1 + (1|U)}.
#' 
#' The intercept is removed automatically. This is done to ensure that parameter 
#' estimates are of the means of interest, rather than being expressed as a mean 
#' with offsets.
#' 
#' @return
#' A list with the elements:
#' \tabular{rl}{
#'   Y       \tab response variable \cr
#'   X1      \tab design matrix for constrained effect \cr
#'   X2      \tab design matrix for covariates \cr
#'   P1      \tab number of constrained coefficients \cr
#'   U       \tab matrix of random effects \cr
#'   formula \tab the final formula call (automatically removes intercept) \cr
#'   dframe  \tab the dataframe containing the variables in the model \cr
#'   REidx   \tab an element to define random effect variance components \cr
#'   REnames \tab an element to define random effect variance components \cr
#' }
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' data( rat.blood )
#' model_terms_clme( mcv ~ time + temp + sex + (1|id) , data = rat.blood )
#' 
#' @importFrom lme4 lFormula
#' @export
#' 
model_terms_clme <- function( formula, data ){
  
  cc <- match.call()
  mf <- match.call( expand.dots = FALSE )
  cc$formula <- mf$formula <- update.formula( formula , . ~ . - 1 )
  
  m <- match( c("formula", "data"), names(mf), 0L )
  mf <- mf[ c(1L, m) ]
  # mf$drop.unused.levels <- TRUE
  
  ## Determine if random effects are present
  any_RE <- length( lme4::findbars(formula) )
  
  if( any_RE==0 ){
    ## ----- NO RANDOM EFFECTS -----
    mf[[1L]] <- quote(stats::model.frame)
    mf <- eval( mf, parent.frame() )
    mt <- attr(mf, "terms")
    
    if( !is.ordered( mf[,2] ) ){
      ## If this is stable, can remove the double xlev call and remove the else portion (same as below)
      warning( "Constrained effect is not an ordered factor, attempting to force ordering\n Assumed order: ", paste0( levels(mf[,2]), collapse=" < "), "\n See help(ordered) for more information")
      mf[,2] <- factor( mf[,2], ordered=TRUE )
      xlev <- levels( mf[,2] )
    } else{
      xlev <- levels( mf[,2] )
    }
    
    Y  <- model.response(mf, "numeric")
    X  <- model.matrix(mt, mf )
    P1 <- length( unique(mf[,2]) )
    X1 <- X[, 1:P1 , drop=FALSE]
    if( ncol(X) > P1 ){
      X2 <- X[,(P1+1):ncol(X), drop=FALSE]  
    } else{
      X2 <- NULL
    }
    U <- NULL
    
    dframe <- mf
    
    REnames <- NULL
    REidx   <- NULL
    
  } else{
    ## ----- YES RANDOM EFFECTS -----
    mf[[1L]] <- quote( lme4::lFormula )
    
    suppressMessages( 
      # clme_terms <- lme4::lFormula( formula, data=data ) 
      clme_terms <- eval( mf, parent.frame() )
    )
    
    if( !is.ordered( clme_terms$fr[,2] ) ){
      # stop( "Constrained effect is not an ordered factor")
      warning( "Constrained effect is not an ordered factor, attempting to force ordering\n Assumed order: ", paste0( levels(clme_terms$fr[,2]), collapse=" < "), "\n See help(ordered) for more information")
      clme_terms$fr[,2] <- factor( clme_terms$fr[,2], ordered=TRUE )
      xlev <- levels( clme_terms$fr[,2] )
    } else{
      xlev <- levels( clme_terms$fr[,2] )
    }
    
    Y  <- clme_terms$fr[,1]
    X  <- clme_terms$X
    P1 <- length( unique(clme_terms$fr[,2]) )
    X1 <- X[, 1:P1 , drop=FALSE]
    
    if( ncol(X) > P1 ){
      X2 <- X[,(P1+1):ncol(X), drop=FALSE]  
    } else{
      X2 <- NULL
    }
    
    U  <- t( as.matrix(clme_terms$reTrms$Zt) )
    
    dframe  <- clme_terms$fr
    REnames <- names(clme_terms$reTrms$flist)
    REidx   <- clme_terms$reTrms$Lind
    
  }
  
  ##############################
  return_obj <- list( Y=Y , X1=X1 , X2=X2 , P1=P1 , U=U , formula=cc$formula, 
                      dframe=dframe ,  REidx=REidx , REnames=REnames, xlev=xlev )
  
  return( return_obj )
  
  
}


##
## Some methods for class CLME
##

#' Constructor method for objects S3 class clme
#' 
#' @rdname as.clme
#' @export
#' 
is.clme <- function(x) inherits(x, "clme")


#' Constructor method for objects S3 class clme
#'
#' @description
#' Test if an object is of class \code{clme} or coerce an object to be such.
#' 
#' @rdname as.clme
#' 
#' @param x list with the elements corresponding to the output of \code{\link{clme}}.
#' @param ... space for additional arguments.
#' 
#' @return
#' Returns an object of the class \code{clme}.
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' data( rat.blood )
#' 
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' is.clme( clme.out )
#' as.clme( clme.out )
#' @export
#' 
as.clme <- function( x , ... ){
  
  if( is.clme(x) ){
    return(x)
  } else{
    
    err.flag  <- 0
    flagTheta <- flagSsq <- flagTsq <- flagCov <- flagW1 <- flagW2 <- flagP1 <- flagP2 <- flagConst <- ""
    
    if( !is.numeric(x$theta) ){
      err.flag  <- 1
      flagTheta <- " theta must be numeric \n"
      x$theta   <- numeric(0)
    }
    
    if( !is.numeric(x$ssq) ){
      err.flag <- 1
      flagSsq  <- " ssq must be numeric \n"
      x$ssq    <- numeric(0)
    } 
    
    if( !is.null(x$tsq) & !is.numeric(x$tsq) ){
      err.flag <- 1
      flagTsq  <- " if present, tau must be numeric \n"
      x$tsq    <- NULL
    }
    
    if( !is.matrix(x$cov.theta) || !is.numeric(x$cov.theta) ||
        nrow(x$cov.theta) != ncol(x$cov.theta) ||
        nrow(x$cov.theta) != length(x$theta)   ||
        sum(sum(abs(x$cov.theta - t(x$cov.theta)))) > sqrt(.Machine$double.eps) ){
      err.flag    <- 1
      flagCov     <- " cov.theta must be square, symmetric, numeric matrix with dimensions equal to length of theta\n"
      x$cov.theta <- matrix( numeric(0) , nrow=length(x$theta) , ncol=length(x$theta) )
    }
    
    if( !is.numeric(x$ts.glb) ){
      err.flag <- 1
      flagW1   <- " ts.glb must be numeric \n"
      x$ts.glb  <- numeric(0)
    } 
    
    if( !is.numeric(x$ts.ind) ){
      err.flag <- 1
      flagW2   <- " ts.ind must be numeric \n"
      x$ts.ind  <- numeric(0)
    }
    
    if( !is.numeric(x$p.value) || length(x$p.value) != length(x$ts.glb) ){
      err.flag   <- 1
      flagP1     <- " p.value must be numeric and of same length as ts.glb \n"
      x$p.value  <- numeric(0)
    } 
    
    if( !is.numeric(x$p.value.ind) || length(x$p.value.ind) != length(x$ts.ind) ){
      err.flag       <- 1
      flagP2         <- " p.value.ind must be numeric and of same length as ts.ind \n"
      x$p.value.ind  <- numeric(0)
    } 
    
    if( !is.list(x$constraints) ){
      err.flag        <- 1
      flagConst       <- " constraints must be list \n"
      x$constraints   <- list( A=matrix( numeric(0) ) )
    } else{
      cnames <- names(x$constraints)
      if( sum(cnames=="A") != 1 ){
        err.flag        <- 1
        flagConst       <- " constraints must contain element A\n"
        x$constraints$A <- matrix( numeric(0) , nrow=length(x$ts.ind ))
      }
    }
    
    
    if( err.flag==1 ){
      err.mssg <- paste( "coercing 'x' to class 'clme' produced errors: \n", 
                         flagTheta, flagSsq, flagTsq, flagCov, flagW1,
                         flagW2, flagP1, flagP2, flagConst, "output may not be valid." , sep = "")
      # warning(warn, sys.call(-1))
      warning( err.mssg )
    }
    
    class(x) <- "clme"
    return(x)
    
  }
}


################################################################################

#' Akaike information criterion
#'
#' @description
#' Calculates the Akaike and Bayesian information criterion for objects of class \code{clme}. 
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments.
#' @param k value multiplied by number of coefficients
#' 
#' @details
#' The log-likelihood is assumed to be the Normal distribution. The model uses residual bootstrap methodology, and Normality is neither required nor assumed. Therefore the log-likelihood and these information criterion may not be useful measures for comparing models.
#' For \code{k=2}, the function computes the AIC. To obtain BIC, set \eqn{k = log( n/(2*pi) )}; which the method \code{BIC.clme} does.
#' 
#' 
#' @return
#' Returns the information criterion (numeric).
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' 
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' AIC( clme.out )
#' AIC( clme.out, k=log( nobs(clme.out)/(2*pi) ) )
#' 
#' 
#' @method AIC clme
#' @export
#' 
AIC.clme <- function( object, ..., k=2 ){
  ## For BIC, set k = ln( n/(2*pi) )
  logl <- logLik.clme( object, ...)[1]
  kk   <- ncol(model.matrix.clme(object)) + length(object$tsq) + length(object$ssq) 
  # aic  <- 2*(kk - logl)
  aic  <- k*kk - 2*logl
  return(aic)
}

#' Akaike information criterion
#'
#' @description
#' Calculates the Akaike and Bayesian information criterion for objects of class \code{clme}. 
#' 
#' @rdname AIC.clme
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @method AIC summary.clme
#' @export
#' 
AIC.summary.clme <- function( object, ..., k=2 ){
  class(object) <- "clme"
  AIC( object, ..., k=k )
}



#' Bayesian information criterion
#'
#' @description
#' Calculates the Bayesian information criterion for objects of class \code{clme}. 
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments.
#' @param k value multiplied by number of coefficients
#' 
#' @details
#' The log-likelihood is assumed to be the Normal distribution. The model uses residual bootstrap methodology, and Normality is neither required nor assumed. Therefore the log-likelihood and these information criterion may not be useful measures for comparing models.
#' For \code{k=2}, the function computes the AIC. To obtain BIC, set \eqn{k = log( n/(2*pi) )}; which the method \code{BIC.clme} does.
#' 
#' 
#' @return
#' Returns the Bayesian information criterion (numeric).
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' 
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' BIC( clme.out )
#' BIC( clme.out, k=log( nobs(clme.out)/(2*pi) ) )
#' 
#' 
#' @method BIC clme
#' @export
#' 
BIC.clme <- function( object, ..., k=log(nobs(object)/(2*pi)) ){
  ## For BIC, set k = ln( n/(2*pi) )
  logl <- logLik( object, ...)[1]
  bic <- AIC( object, k=k )
  return(bic)
}

#' Bayesian information criterion
#'
#' @description
#' Calculates the Akaike and Bayesian information criterion for objects of class \code{clme}. 
#' 
#' @rdname BIC.clme
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @method BIC summary.clme
#' @export
#' 
BIC.summary.clme <- function( object, ..., k=log(nobs(object)/(2*pi)) ){
  class(object) <- "clme"
  BIC( object, ..., k=k )
}


#' Individual confidence intervals
#'
#' @description
#' Calculates confidence intervals for fixed effects parameter estimates in objects of class \code{clme}.
#' @rdname confint
#' 
#' @param object object of class \code{\link{clme}}.
#' @param parm parameter for which confidence intervals are computed (not used).
#' @param level nominal confidence level.
#' @param ... space for additional arguments.
#' 
#' @details
#' Confidence intervals are computed using Standard Normal critical values. 
#' Standard errors are taken from the covariance matrix of the unconstrained parameter estimates.
#' 
#' 
#' @return
#' Returns a matrix with two columns named lcl and ucl (lower and upper confidence limit).
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' confint( clme.out )
#' 
#' 
#' @method confint clme
#' @export
#' 
confint.clme <- function(object, parm, level=0.95, ...){
  ## More types of confidence intervals (e.g., bootstrap) may be added in the future.
  ## If so, default confidence interval will be the current methods
  
  ## Confidence intervals based on unconstrained variance-covariance matrix
  ## Actual covarage >= Nominal coverage
  
  cc     <- match.call()
  digits <- cc$digits
  if( is.null(digits) ){     digits <- 3 }
  if( !is.numeric(digits) ){ digits <- 3 }
  if( digits < 0  ){         digits <- 3 }
  
  
  alpha <- 1 - level
  theta <- fixef(object)
  cv    <- qnorm(1-alpha/2)
  varco <- vcov( object )
  lcl   <- as.numeric( format( round(theta - cv*sqrt(diag(varco)) , digits=digits)) )
  ucl   <- as.numeric( format( round(theta + cv*sqrt(diag(varco)) , digits=digits)) )
  
  ints <- cbind( lcl, ucl)
  
  ## Return intervals
  return( ints )
}

#' Individual confidence intervals
#'
#' @description
#' Calculates confidence intervals for fixed effects parameter estimates in objects of class \code{clme}.
#' @rdname confint
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @method confint summary.clme
#' @export
#' 
confint.summary.clme <- function(object, parm, level=0.95, ...){
  class(object) <- "clme"
  confint( object, parm, level, ... )
}


#' Extract fixed effects
#' 
#' @importFrom lme4 fixef
#' @method fixef clme
#' @export 
#' 
fixef.clme <- function( object, ...){ UseMethod("fixef") }


#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @importFrom lme4 fixef
#' @method fixef summary.clme
#' @export 
#' 
fixef.summary.clme <- function( object, ...){ 
  class(object) <- "clme"
  fixef(object, ...)
}


#' Extract fixed effects
#'
#' @description
#' Extracts the fixed effects estimates from objects of class \code{clme}. 
#' 
#' @rdname fixef.clme
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' @return
#' Returns a numeric vector.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' fixef( clme.out )
#' 
#' @importFrom lme4 fixef
#' @method fixef clme
#' @export
#' 
fixef.clme <- function( object, ... ){
  ## Print out the fixed effects
  if( is.clme(object) ){
    return( object$theta )
  } else{
    stop("'object' is not of class clme")
  }
}

#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @importFrom nlme fixed.effects
#' @export
#' 
fixed.effects <- function( object, ...){ UseMethod("fixed.effects") }

#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @importFrom nlme fixed.effects
#' @export
#' 
fixed.effects.summary.clme <- function( object, ...){
  class(object) <- "clme"
  fixef(object, ...)
}


#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' 
#' @importFrom nlme fixed.effects
#' @method fixed.effects clme
#' @export
#' 
fixed.effects.clme <- function( object , ... ){
  fixef.clme( object, ... )
}


#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @method coefficients clme
#' @export
#' 
coefficients.clme <- function( object, ... ){
  fixef.clme( object, ... )
}

#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @method coef clme
#' @export
#' 
coef.clme <- function( object, ... ){
  fixef.clme( object, ... )
}

#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @method coefficients summary.clme
#' @export
#' 
coefficients.summary.clme <- function( object, ... ){
  class(object) <- "clme"
  fixef.clme( object, ... )
}

#' Extract fixed effects
#' 
#' @rdname fixef.clme
#' @method coef summary.clme
#' @export
#' 
coef.summary.clme <- function( object, ... ){
  class(object) <- "clme"
  fixef.clme( object, ... )
}


#' Extract formula
#'
#' @description
#' Extracts the formula from objects of class \code{clme}. 
#' 
#' @param x object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' @details
#' The package \pkg{CLME} parametrizes the model with no intercept term. 
#' If an intercept was included, it will be removed automatically.
#' 
#' @return
#' Returns a formula object
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' formula( clme.out )
#' 
#' @method formula clme
#' @export
#' 
formula.clme <- function(x, ...){
  return( x$formula )
}


#' Log-likelihood
#'
#' @description
#' Computes the log-likelihood of the fitted model for objects of class \code{clme}. 
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' @details
#' The log-likelihood is computed using the Normal distribution. The model uses residual bootstrap 
#' methodology, and Normality is neither required nor assumed. Therefore the log-likelihood may 
#' not be a useful measure in the context of \pkg{CLME}.
#' 
#' @return
#' Numeric.
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' logLik( clme.out )
#' 
#' @method logLik clme
#' @export
#' 
logLik.clme <- function( object, ...){
  ## Residuals
  YY <- object$dframe[,1]
  XX <- model.matrix( object )
  TT <- fixef( object )
  RR <- YY - apply( XX , 1 , FUN=function(xx,tht){ sum(xx*tht) }, tht=TT )
  nn <- nobs(object)
  
  ## Covariance matrix (piecewise)
  ssq <- object$ssq
  Nks <- object$gfix
  ssqvec <- rep( ssq, Nks )
  RSiR <- c( t(RR) %*% (RR/ssqvec) )
  detS <- sum( log(ssqvec) )
  
  if( is.null(object$tsq) ){
    ## Fixed effects only
    RPiR   <- RSiR
    detPhi <- detS
  } else{
    ## Mixed Effects    
    UU     <- model.matrix( object , type="ranef" )
    tsq    <- object$tsq
    Qs     <- object$gran
    tsqvec <- rep( tsq, Qs )
    RSiU   <- matrix( apply( UU, 2, FUN=function(uu,rr){ sum(uu*rr) }, rr=(RR/ssqvec) ), nrow=1 )
    
    U1         <- apply( UU, 2, FUN=function(uu,sq){uu/sq}, sq=sqrt(ssqvec) )
    tusu       <- t(U1) %*% U1
    diag(tusu) <- diag(tusu) + 1/tsqvec
    tusui      <- solve( tusu )
    
    RPiR   <- RSiR - c( RSiU%*%(tusui%*%t(RSiU)) )
    detPhi <- log(det( tusu )) + sum( log(tsqvec) ) + detS
  }
  
  logL <- -0.5*( nn*log(2*pi) + detPhi + RPiR )
  
  return( logL )
  
}

#' Log-likelihood
#' 
#' @rdname logLik.clme
#' 
#' @seealso
#' \code{\link{logLik.clme}} 
#' 
#' @method logLik summary.clme
#' @export
#' 
logLik.summary.clme <- function( object, ...){
  class(object) <- "clme"
  logLik(object, ...)
}




#' Extract the model design matrix.
#'
#' @description
#' Extracts the fixed-effects design matrix from objects of class \code{clme}.
#' 
#' @param object an object of class \code{clme}.
#' @param type specify whether to return the fixed-effects or random-effects matrix.
#' @param ... space for additional arguments
#' 
#' 
#' @return
#' Returns a matrix.
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' \dontrun{
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' model.matrix( clme.out )
#' }
#' @method model.matrix clme
#' @export
#' 
model.matrix.clme <- function( object, type="fixef", ...){
  mmat <- model_terms_clme( object$formula, object$dframe )  
  if( type=="fixef" ){
    ## Return the fixed-effects matrix
    X1   <- mmat$X1
    X2   <- mmat$X2
    return( cbind(X1, X2) )  
  } else if( type=="ranef" ){
    ## Return the random-effects matrix
    return(mmat$U)
  }
}

#' Extract the model design matrix.
#' 
#' @rdname model.matrix.clme
#' 
#' @seealso
#' \code{\link{model.matrix.clme}} 
#' 
#' @method model.matrix summary.clme
#' @export
#' 
model.matrix.summary.clme <- function( object, ...){
  class(object) <- "clme"
  model.matrix(object, ...)
}

#' Number of observations
#'
#' @description
#' Obtains the number of observations used to fit an model for objects of class \code{clme}.
#' 
#' @param object an object of class \code{clme}.
#' @param ... space for additional arguments
#' 
#' 
#' @return
#' Numeric.
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' nobs( clme.out )
#' 
#' @method nobs clme
#' @export
#' 
nobs.clme <- function(object, ...){
  nrow( model.matrix.clme(object) )
}

#' Number of observations
#' 
#' @rdname nobs.clme
#' 
#' @seealso
#' \code{\link{nobs.clme}} 
#' 
#' @method nobs summary.clme
#' @export
#' 
nobs.summary.clme <- function( object, ...){
  class(object) <- "clme"
  nobs(object, ...)
}



#' Printout of fitted object.
#'
#' @description
#' Prints basic information on a fitted object of class \code{clme}.
#' 
#' @param x an object of class \code{clme}.
#' @param ... space for additional arguments
#' 
#' @return
#' Text printed to console.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' \dontrun{
#' data( rat.blood )
#' set.seed( 42 )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 10)
#' 
#' print( clme.out )
#' }
#' 
#' @method print clme
#' @export
#' 
print.clme <- function(x, ...){
  #cc     <- match.call()
  #digits <- cc$digits
  
  object <- x
  ## Print out residuals of specified type
  if( !is.clme(object) ){
    stop("'object' is not of class clme")
  }
  
  cat( "Linear mixed model subject to order restrictions\n")
  cat( "Formula: ")
  print( object$formula )  
  crit <- c(logLik.clme(object),
            AIC.clme(object),
            AIC( object, k=log(nobs.clme(object)/(2*pi)) ) )  
  
  critc <- format( crit , digits=5)
  
  cat( "\nlog-likelihood:", critc[1] )
  cat( "\nAIC:           ", critc[2] )
  cat( "\nBIC:           ", critc[3] )
  cat( "\n(log-likelihood, AIC, BIC computed under normality)")  
  
  cat( "\n\nFixed effect coefficients (theta): \n")
  print( fixef.clme(object)  )
  
  cat( "\nVariance components: \n")
  print( VarCorr.clme(object) )
  #cat( "\n\nModel based on", object$nsim, "bootstrap samples." )
}


#' Extract random effects
#' 
#' @param object object of class clme.
#' @param ... space for additional arguments
#' 
#' @rdname ranef.clme
#' @importFrom nlme ranef
#' @export
#' 
ranef.clme <- function( object, ...){ UseMethod("ranef") }

#' Extract random effects
#' 
#' @rdname ranef.clme
#' 
#' @importFrom nlme ranef
#' @export
#' 
ranef.summary.clme <- function( object, ...){
  class(object) <- "clme"
  ranef(object, ...)
}


#' Extract random effects
#'
#' @description
#' Extracts the random effects estimates from objects of class \code{clme}.
#' 
#' @rdname ranef.clme
#' 
#' @return
#' Returns a numeric vector.
#' 
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#' 
#' ranef( clme.out )
#' 
#' 
#' @importFrom nlme ranef
#' @method ranef clme
#' @export
#' 
ranef.clme <- function( object, ... ){
  ## Print out the random effects
  if( is.clme(object) ){
    return( object$random.effects )
  } else{
    stop("'object' is not of class clme")
  }
}

#' Extract random effects
#' 
#' @param object object of class clme.
#' @param ... space for additional arguments
#'  
#' @rdname ranef
#' @importFrom nlme random.effects
#' @export
#' 
random.effects <- function( object, ... ){ UseMethod("random.effects") }

#' Extract random effects
#' 
#' @rdname ranef
#' @importFrom nlme random.effects
#' @export
#' 
random.effects.summary.clme <- function( object, ...){
  class(object) <- "clme"
  ranef(object, ...)
}


#' Extract random effects
#' 
#' @rdname ranef.clme
#' @importFrom nlme random.effects
#' @method random.effects clme
#' @export
#' 
random.effects.clme <- function( object , ... ){
  ranef.clme( object, ... )
}



#' Various types of residuals 
#'
#' @description
#' Computes several types of residuals for objects of class \code{clme}. 
#' 
#' @param object object of class \code{\link{clme}}.
#' @param type type of residual (for mixed-effects models only).
#' @param ... space for additional arguments
#' 
#' @details
#' For fixed-effects models \eqn{Y = X\beta + \epsilon}{Y = X*b + e}, residuals are given as \deqn{\hat{e} = Y - X\hat{\beta}}{ ehat = Y - X*betahat}.
#' For mixed-effects models \eqn{Y = X\beta + + U\xi + \epsilon}{Y = X*b + U*xi + e}, three types of residuals are available.
#' \eqn{PA = Y - X\hat{\beta}}{ PA = Y - X*betahat}\\
#' \eqn{SS = U\hat{\xi}}{ SS = U*xihat}\\
#' \eqn{FM = Y - X\hat{\beta} - U\hat{\xi}}{ FM = Y - X*betahat - U*xihat}
#' 
#' @return
#' Returns a numeric matrix.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' \dontrun{
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'                  
#' residuals( clme.out, type='PA' )
#' }
#' @method residuals clme
#' @export
#' 
residuals.clme <- function( object, type="FM", ... ){
  ## Print out residuals of specified type
  if( is.clme(object) ){
    ridx <- which( c("PA", "SS", "FM")==type )
    if( ncol(object$residuals)<ridx ) ridx <- 1
    return( object$residuals[,ridx] )
  } else{
    stop("'object' is not of class clme")
  }
}

#' Various types of residuals 
#'
#' @rdname residuals.clme
#' 
#' @method residuals summary.clme
#' @export
#' 
residuals.summary.clme <- function( object, type="FM", ... ){
  class(object) <- "clme"
  residuals( object, type, ...)
}




#' Residual variance components
#'
#' @description
#' Extract residual variance components for objects of class \code{clme}.
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' 
#' @return
#' Numeric.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'          
#' sigma( clme.out )
#' 
#' 
#' @importFrom stats sigma
#' @method sigma clme
#' @export
#' 
# WAS IMPORTING sigma FROM lme4
sigma.clme <- function( object, ...){
  return( object$ssq )
}

#' Residual variance components
#'
#' @description
#' Extract residual variance components for objects of class \code{clme}.
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' 
#' @return
#' Numeric.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'                  
#' sigma( clme.out )
#' 
#' @method sigma summary.clme
#' @export
#' 
sigma.summary.clme <- function( object, ...){
  return( object$ssq )
}


#' Variance components
#' 
#' @param x object of class \code{\link{summary.clme}}.
#' @param sigma (unused at present).
#' @param rdig number of digits to round to (unused at present).
#'  
#' @rdname VarCorr
#' @export
#' 
VarCorr <- function( x, sigma, rdig ){ UseMethod("VarCorr") }


#' Variance components
#' 
#' @rdname VarCorr
#' @export
#' 
VarCorr.summary.clme <- function( x, sigma, rdig ){
  class(x) <- "clme"
  VarCorr(x, sigma=1, rdig=4)
}




#' Variance components.
#'
#' @description
#' Extracts variance components for objects of class \code{clme}. 
#' 
#' @rdname VarCorr
#' 
#' @return
#' Numeric.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'                  
#' VarCorr( clme.out )
#' 
#' 
#' 
#' @method VarCorr clme
#' @export
#' 
VarCorr.clme <- function(x, sigma, rdig ){
  # @importFrom lme4 VarCorr
  ## Print out variances or SDs
  ## Defines tiny class "varcorr_clme" to handle printing
  ## using the method: print.varcorr_clme
  if( !is.clme(x) ){
    stop("'x' is not of class clme")
  } else{
    varcomps <- matrix( sqrt(c(x$tsq, x$ssq )), ncol=1 )
    if( !is.null(x$tsq) & is.null(names(x$tsq)) ){
      names(x$tsq) <- paste0( "tau_", 1:length(x$tsq) )
    }
    rnames   <- c( "Source", names(x$tsq), names(x$ssq) )
    rownames(varcomps) <- rnames[-1]
    colnames(varcomps) <- "Std. Error"
    #class(varcomps)    <- "varcorr_clme"
    return( varcomps )
  }
}

## Leave this method out of alphabetical order so that
## is it right next to the VarCorr.clme method

#' Printout for variance components
#'
#' @description
#' Prints variance components of an objects of \code{clme}. 
#' 
#' @param object object of class \code{\link{clme}}.
#' @param rdig number of digits to round to.
#' @param ... space for additional arguments.
#' 
#' @return
#' Text printed to console.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' \dontrun{
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'                  
#' print.varcorr_clme( clme.out )
#' }
#' @importFrom stringr str_pad
#' 
#' 
print.varcorr_clme <- function(object, rdig=5, ...){
  
  rnames   <- c( "Source", rownames( object ) )
  rnames   <- str_pad(rnames, width=max(nchar(rnames)), side = "right", pad = " ")
  vars     <- format( object , digits=rdig )
  
  cat( rnames[1], "\t" , "Variance" )
  for( ii in 1:length(vars) ){
    cat( "\n", rnames[ii+1], "\t" , vars[ii] )
  }
  
  return(NULL)
  # @exportMethod print varcorr_clme
}



#' Variance-covariance matrix
#'
#' @description
#' Extracts variance-covariance matrix for objects of class \code{clme}.
#' 
#' @param object object of class \code{\link{clme}}.
#' @param ... space for additional arguments
#' 
#' 
#' @return
#' Numeric matrix.
#' 
#' @seealso
#' \code{\link{CLME-package}}
#' \code{\link{clme}}
#' 
#' @examples
#' 
#' data( rat.blood )
#' cons <- list(order = "simple", decreasing = FALSE, node = 1 )
#' clme.out <- clme(mcv ~ time + temp + sex + (1|id), data = rat.blood , 
#'                  constraints = cons, seed = 42, nsim = 0)
#'                  
#' vcov( clme.out )
#' 
#' @method vcov clme
#' @export
#' 
vcov.clme <- function(object, ...){
  ## Print out covariance matrix of theta
  if( is.clme(object) ){
    return( object$cov.theta )
  } else{
    stop("'object' is not of class clme")
  }  
}

#' Variance-covariance matrix
#'
#' @rdname vcov.clme
#' 
#' @method vcov summary.clme
#' @export
#' 
vcov.summary.clme <- function(object, ...){
  class(object) <- "clme"
  vcov(object, ...)
}


##
## Hidden functions to format characters / decimals
##



## Align the length of header with length of values
.align_table.clme <- function( tble, digits=4, ... ){
  cnames <- colnames( tble )
  for( ii in 1:length(cnames) ){
    ntitle <- nchar( cnames[ii] )
    maxc   <- max( nchar(tble[,ii]) )    
    if( ntitle > maxc ){
      tble[,ii] <- str_pad( tble[,ii], width=ntitle, side = "left", pad = " ")    
    }
    if( ntitle < maxc ){
      cnames[ii] <- str_pad( cnames[ii], width=(maxc+1), side = "right", pad = " ")    
    }
  }
  colnames(tble) <- cnames
  return( tble )
}