R/fixGMMAT.R

In devoges/fstat: fstat

#' fit_GMMAT
#'
#' @param Y numeric outcome vector
#' @param X numeric fixed effect design matrix -- do include intercept
#' @param Ss NAMED list of (n \times n) random effect covariance matrices
#'
#' @return a model fit with GMMAT::glmmkin
#' @export
fit_GMMAT <- function(Y,
                      X,
                      Ss) {
  Y <- c(Y)
  data <- model.matrix(Y ~ X) %>%
    as.data.frame() %>%
    remove_intercept() %>%
    as.data.frame() %>%
    add_column(Y = Y)
  
  n <- data %>% nrow
  p <- data %>% ncol
  
  m <- glmmkin(fixed = Y ~ .,
               data = data,
               family = gaussian(),
               kins = Ss,
               tol = .Machine$double.eps^0.5)
  m$theta %<>% `names<-`(c("Error", names(Ss)))
  
  Vs <- nlapply(names(Ss), function(n) Ss[[n]] * m$theta[n])
  m$V <- Reduce(`+`, Vs)
  diag(m$V) %<>% add(1)
  
  
  m$residuals %<>% c %>% unname %>% matrix(n, 1)
  
  m$RSS <- t(m$residuals) %mtimes% minv(m$V) %mtimes% m$residuals
  
  ##-- Calculate RML
  Vinv <- minv(m$V)
  P <- diag(n) - X %mtimes% minv(t(X) %mtimes% Vinv %mtimes% X) %mtimes% t(X) %mtimes% Vinv
  detV <- determinant(m$V, logarithm = TRUE)$modulus
  logYPVPY <- log(t(Y) %mtimes% t(P) %mtimes% Vinv %mtimes% P %mtimes% Y)
  rml <- -1 * detV +
    -1 * determinant(t(X) %mtimes% Vinv %mtimes% X, logarithm = TRUE)$modulus +
    -1 * (n - p) * logYPVPY
  m$rml <- rml[1]/2
  
  ml <- -1 * detV - n * logYPVPY
  m$ml <- ml[1]/2
  
  return(m)
}