miloR: Differential neighbourhood abundance testing on a graph

Documented in fitGeneticPLGlmm fitPLGlmm

# Generated by using Rcpp::compileAttributes() -> do not edit by hand
# Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393

#' GLMM parameter estimation using pseudo-likelihood with a custom covariance matrix
#'
#' Iteratively estimate GLMM fixed and random effect parameters, and variance
#' component parameters using Fisher scoring based on the Pseudo-likelihood
#' approximation to a Normal loglihood. This function incorporates a user-defined
#' covariance matrix, e.g. a kinship matrix for genetic analyses.
#'
#' @param Z mat - sparse matrix that maps random effect variable levels to
#' observations
#' @param X mat - sparse matrix that maps fixed effect variables to
#' observations
#' @param K mat - sparse matrix that defines the known covariance patterns between
#' individual observations. For example, a kinship matrix will then adjust for the
#' known/estimated genetic relationships between observations.
#' @param muvec vec vector of estimated phenotype means
#' @param offsets vec vector of model offsets
#' @param curr_theta vec vector of initial parameter estimates
#' @param curr_beta vec vector of initial beta estimates
#' @param curr_u vec of initial u estimates
#' @param curr_sigma vec of initial sigma estimates
#' @param curr_G mat c X c matrix of variance components
#' @param y vec of observed counts
#' @param u_indices List a List, each element contains the indices of Z relevant
#' to each RE and all its levels
#' @param theta_conv double Convergence tolerance for paramter estimates
#' @param rlevels List containing mapping of RE variables to individual
#' levels
#' @param curr_disp double Dispersion parameter estimate
#' @param REML bool - use REML for variance component estimation
#' @param maxit int maximum number of iterations if theta_conv is FALSE
#' @param solver string which solver to use - either HE (Haseman-Elston regression) or Fisher scoring
#' @param vardist string which variance form to use NB = negative binomial, P=Poisson [not yet implemented]/
#'
#' @details Fit a NB-GLMM to the counts provided in \emph{y}. The model uses an iterative approach that
#' switches between the joint fixed and random effect parameter inference, and the variance component
#' estimation. A pseudo-likelihood approach is adopted to minimise the log-likelihood of the model
#' given the parameter estimates. The fixed and random effect parameters are estimated using
#' Hendersons mixed model equations, and the variance component parameters are then estimated with
#' the specified solver, i.e. Fisher scoring, Haseman-Elston or constrained Haseman-Elston regression. As
#' the domain of the variance components is [0, +\code{Inf}], any negative variance component estimates will
#' trigger the switch to the HE-NNLS solver until the model converges.
#'
#' @return A \code{list} containing the following elements (note: return types are dictated by Rcpp, so the R
#' types are described here):
#' \describe{
#' \item{\code{FE}:}{\code{numeric} vector of fixed effect parameter estimates.}
#' \item{\code{RE}:}{\code{list} of the same length as the number of random effect variables. Each slot contains the best
#' linear unbiased predictors (BLUPs) for the levels of the corresponding RE variable.}
#' \item{\code{Sigma:}}{\code{numeric} vector of variance component estimates, 1 per random effect variable. For this model the
#' last variance component corresponds to the input \emph{K} matrix.}
#' \item{\code{converged:}}{\code{logical} scalar of whether the model has reached the convergence tolerance or not.}
#' \item{\code{Iters:}}{\code{numeric} scalar with the number of iterations that the model ran for. Is strictly <= \code{max.iter}.}
#' \item{\code{Dispersion:}}{\code{numeric} scalar of the dispersion estimate computed off-line}
#' \item{\code{Hessian:}}{\code{matrix} of 2nd derivative elements from the fixed and random effect parameter inference.}
#' \item{\code{SE:}}{\code{matrix} of standard error estimates, derived from the hessian, i.e. the square roots of the diagonal elements.}
#' \item{\code{t:}}{\code{numeric} vector containing the compute t-score for each fixed effect variable.}
#' \item{\code{COEFF:}}{\code{matrix} containing the coefficient matrix from the mixed model equations.}
#' \item{\code{P:}}{\code{matrix} containing the elements of the REML projection matrix.}
#' \item{\code{Vpartial:}}{\code{list} containing the partial derivatives of the (pseudo)variance matrix with respect to each variance
#' component.}
#' \item{\code{Ginv:}}{\code{matrix} of the inverse variance components broadcast to the full Z matrix.}
#' \item{\code{Vsinv:}}{\code{matrix} of the inverse pseudovariance.}
#' \item{\code{Winv:}}{\code{matrix} of the inverse elements of W = D^-1 V D^-1}
#' \item{\code{VCOV:}}{\code{matrix} of the variance-covariance for all model fixed and random effect variable parameter estimates.
#' This is required to compute the degrees of freedom for the fixed effect parameter inference.}
#' \item{\code{CONVLIST:}}{\code{list} of \code{list} containing the parameter estimates and differences between current and previous
#' iteration estimates at each model iteration. These are included for each fixed effect, random effect and variance component parameter.
#' The list elements for each iteration are: \emph{ThetaDiff}, \emph{SigmaDiff}, \emph{beta}, \emph{u}, \emph{sigma}.}
#' }
#'
#' @author Mike Morgan
#'
#' @examples
#' NULL
#'
#' @name fitGeneticPLGlmm
#'
fitGeneticPLGlmm <- function(Z, X, K, muvec, offsets, curr_beta, curr_theta, curr_u, curr_sigma, curr_G, y, u_indices, theta_conv, rlevels, curr_disp, REML, maxit, solver, vardist) {
    .Call('_miloR_fitGeneticPLGlmm', PACKAGE = 'miloR', Z, X, K, muvec, offsets, curr_beta, curr_theta, curr_u, curr_sigma, curr_G, y, u_indices, theta_conv, rlevels, curr_disp, REML, maxit, solver, vardist)
}

#' GLMM parameter estimation using pseudo-likelihood
#'
#' Iteratively estimate GLMM fixed and random effect parameters, and variance
#' component parameters using Fisher scoring based on the Pseudo-likelihood
#' approximation to a Normal loglihood.
#'
#' @param Z mat - sparse matrix that maps random effect variable levels to
#' observations
#' @param X mat - sparse matrix that maps fixed effect variables to
#' observations
#' @param muvec vec vector of estimated phenotype means
#' @param offsets vec vector of model offsets
#' @param curr_theta vec vector of initial parameter estimates
#' @param curr_beta vec vector of initial beta estimates
#' @param curr_u vec of initial u estimates
#' @param curr_sigma vec of initial sigma estimates
#' @param curr_G mat c X c matrix of variance components
#' @param y vec of observed counts
#' @param u_indices List a List, each element contains the indices of Z relevant
#' to each RE and all its levels
#' @param theta_conv double Convergence tolerance for paramter estimates
#' @param rlevels List containing mapping of RE variables to individual
#' levels
#' @param curr_disp double Dispersion parameter estimate
#' @param REML bool - use REML for variance component estimation
#' @param maxit int maximum number of iterations if theta_conv is FALSE
#' @param solver string which solver to use - either HE (Haseman-Elston regression) or Fisher scoring
#' @param vardist string which variance form to use NB = negative binomial, P=Poisson [not yet implemented.]
#'
#' @details Fit a NB-GLMM to the counts provided in \emph{y}. The model uses an iterative approach that
#' switches between the joint fixed and random effect parameter inference, and the variance component
#' estimation. A pseudo-likelihood approach is adopted to minimise the log-likelihood of the model
#' given the parameter estimates. The fixed and random effect parameters are estimated using
#' Hendersons mixed model equations, and the variance component parameters are then estimated with
#' the specified solver, i.e. Fisher scoring, Haseman-Elston or constrained Haseman-Elston regression. As
#' the domain of the variance components is [0, +\code{Inf}], any negative variance component estimates will
#' trigger the switch to the HE-NNLS solver until the model converges.
#'
#' @return A \code{list} containing the following elements (note: return types are dictated by Rcpp, so the R
#' types are described here):
#' \describe{
#' \item{\code{FE}:}{\code{numeric} vector of fixed effect parameter estimates.}
#' \item{\code{RE}:}{\code{list} of the same length as the number of random effect variables. Each slot contains the best
#' linear unbiased predictors (BLUPs) for the levels of the corresponding RE variable.}
#' \item{\code{Sigma:}}{\code{numeric} vector of variance component estimates, 1 per random effect variable.}
#' \item{\code{converged:}}{\code{logical} scalar of whether the model has reached the convergence tolerance or not.}
#' \item{\code{Iters:}}{\code{numeric} scalar with the number of iterations that the model ran for. Is strictly <= \code{max.iter}.}
#' \item{\code{Dispersion:}}{\code{numeric} scalar of the dispersion estimate computed off-line}
#' \item{\code{Hessian:}}{\code{matrix} of 2nd derivative elements from the fixed and random effect parameter inference.}
#' \item{\code{SE:}}{\code{matrix} of standard error estimates, derived from the hessian, i.e. the square roots of the diagonal elements.}
#' \item{\code{t:}}{\code{numeric} vector containing the compute t-score for each fixed effect variable.}
#' \item{\code{COEFF:}}{\code{matrix} containing the coefficient matrix from the mixed model equations.}
#' \item{\code{P:}}{\code{matrix} containing the elements of the REML projection matrix.}
#' \item{\code{Vpartial:}}{\code{list} containing the partial derivatives of the (pseudo)variance matrix with respect to each variance
#' component.}
#' \item{\code{Ginv:}}{\code{matrix} of the inverse variance components broadcast to the full Z matrix.}
#' \item{\code{Vsinv:}}{\code{matrix} of the inverse pseudovariance.}
#' \item{\code{Winv:}}{\code{matrix} of the inverse elements of W = D^-1 V D^-1}
#' \item{\code{VCOV:}}{\code{matrix} of the variance-covariance for all model fixed and random effect variable parameter estimates.
#' This is required to compute the degrees of freedom for the fixed effect parameter inference.}
#' \item{\code{CONVLIST:}}{\code{list} of \code{list} containing the parameter estimates and differences between current and previous
#' iteration estimates at each model iteration. These are included for each fixed effect, random effect and variance component parameter.
#' The list elements for each iteration are: \emph{ThetaDiff}, \emph{SigmaDiff}, \emph{beta}, \emph{u}, \emph{sigma}.}
#' }
#'
#' @author Mike Morgan
#'
#' @examples
#' NULL
#'
#' @name fitPLGlmm
fitPLGlmm <- function(Z, X, muvec, offsets, curr_beta, curr_theta, curr_u, curr_sigma, curr_G, y, u_indices, theta_conv, rlevels, curr_disp, REML, maxit, solver, vardist) {
    .Call('_miloR_fitPLGlmm', PACKAGE = 'miloR', Z, X, muvec, offsets, curr_beta, curr_theta, curr_u, curr_sigma, curr_G, y, u_indices, theta_conv, rlevels, curr_disp, REML, maxit, solver, vardist)
}