R/fixlmer.R

In devoges/fstat: fstat

#' @useDynLib fstat, .registration = TRUE
#' @importFrom Rcpp sourceCpp
NULL

##-- found online at
##-- https://rstudio-pubs-static.s3.amazonaws.com/305732_7a0cfc5b535c41c8b56ffdc2e322a51d.html

myReTrms<-function(ListZ,ListVar=NULL){
  reTrms<-list()
  reTrms$Zt    <- Matrix::Matrix(t(Reduce('cbind',ListZ)),sparse=TRUE)
  reTrms$Zt <- as(reTrms$Zt, "dgCMatrix")
  reTrms$theta <- rep(1,length(ListZ))  # Initial Value of the covariance parameters
  reTrms$Lind  <- rep(1:length(ListZ),unlist(lapply(ListZ,ncol))) #an integer vector of indices determining the mapping of the elements of the theta vector to the "x" slot of Lambdat
  reTrms$Gp    <- as.integer(unname(cumsum(c(0,unlist(lapply(ListZ,ncol))))))
  reTrms$lower   <- rep(0,length(ListZ)) # lower bounds on the covariance parameters
  if (is.null(ListVar)) {
    reTrms$Lambdat <- Matrix(diag(rep(1,sum(unlist(lapply(ListZ,ncol))))),sparse=TRUE)
  }  else {
    reTrms$Lambdat <- bdiag(lapply(ListVar,chol))
  }
  reTrms$Ztlist <- lapply(ListZ,function(Z) as(Matrix(t(Z),sparse=T), "dgCMatrix"))
  reTrms$cnms   <-  as.list(names(ListZ)) ; names(reTrms$cnms)<- names(ListZ)
  # Flist is Not very clean (to say the least... )
  reTrms$flist <- lapply(ListZ,function(Z) {flist.factor<- as.factor(colnames(Z)[apply(Z,1,function(x) which(rmultinom(n=1,size=1,prob =abs(x)+0.1)==1) )]);
  levels(flist.factor)<-colnames(Z); return(flist.factor)}) #NULL # list of grouping factors used in the random-effects terms (used for computing initial variance ??)
  return(reTrms)
}

myranef<-function(model,condVar=TRUE){
  re.cond.mode<-tapply(model@u,mymod@pp$Lind,function(x) x)
  names(re.cond.mode)<- names(model@cnms)

  if (condVar) {
    Zt<-model@pp$Zt
    D  <- sigma(model)* t(model@pp$Lambdat) %*% model@pp$Lambdat
    Sigma<- t(Zt)%*%  D %*%Zt + sigma(model)*diag(rep(1,ncol(Zt)))
    var.cond <- D - Zt %*%solve(Sigma) %*% t(Zt)
    var.cond <- diag(var.cond)
    var.cond.mode <- tapply(var.cond,mymod@pp$Lind,function(x) x)
  }
  for (i in 1:length(re.cond.mode)) {
    re.cond.mode[[i]]<-data.frame(re.cond.mode[i])
    names(re.cond.mode[[i]])<-names(re.cond.mode[i])
    row.names(re.cond.mode[[i]]) <- levels(model@flist[[i]])

    if (condVar) attr(re.cond.mode[[i]],"postVar")=array(var.cond.mode[[i]],c(1,1,nrow(re.cond.mode[[i]])))
  }
  attr(re.cond.mode,"class") <- "ranef.mer"
  re.cond.mode
}

#' Fit lmer model with given fixed and random effect(s) design matrices
#'
#' @param Response n \times 1 vector of responses
#' @param X n \times p data.frame fixed effect design matrix
#' @param ListZ Named R list of size K, of n \times k_i data.frame random effect design matrices
#' @param REML TRUE/FALSE of whether to use REML (FALSE = ML)
#' @param eps control to pass to lmer
#'
#' @return a lmer model object (all functions available for lmer work on it)
#' @export
#'
#' @examples
#' X <- data.frame(matrix(1, 1000, 1))
#' Z1 <- data.frame(matrix(rnorm(2*1000), 1000, 2))
#' Z2 <- data.frame(matrix(rnorm(3*1000), 1000, 3))
#' colnames(X) <- c("(Intercept)")
#'
#' model <- new_lmer(Response = rnorm(1000),
#'                   X        = X,
#'                   ListZ    = list("Z1" = Z1,
#'                                   "Z2" = Z2))
#' summary(model)
fit_LMER <- function(Response,
                     X,
                     ListZ,
                     REML = TRUE){

  ListVar=NULL
  fr<-model.frame(Response~.,data.frame(Response=Response,X) )
  reTrms <- myReTrms(ListZ,ListVar)
  reTrms$Zt <- as(reTrms$Zt, "dgCMatrix")
  devfun <- lme4::mkLmerDevfun(fr, X, reTrms)
  opt <- lme4::optimizeLmer(devfun)
  return(lme4::mkMerMod(environment(devfun), opt, reTrms, fr = fr))
}