R/helpers.R

In ClaireKelley/WeMix: Weighted Mixed-Effects Models using Multilevel Pseudo Maximum Likelihood Estimation

# This helper funciton returns just the coefficient from the WeMix results object.
#' @export
coef.WeMixResults <- function(object, ...) {
  object$coef
}

# This helper function prints the coefficient from the WeMix results object.
#' @export
print.WeMixResults <- function(x, ...) {
  print(x$coef)
}

# This helper funciton creates a coefficient matrix from WeMix results.
#' @importFrom stats cov2cor
#' @export
summary.WeMixResults <- function(object, ...) {
  x0 <- object
  
  #for adaptive quad models 
  if(object$is_adaptive){
    VC <- makeSandwich(object) # calculates the sandwhich estimator of variance 
    object$vc <- VC$VC #sets variances returned by model to variances calculated by sandwich estimator
    se <- VC$se[1:length(x0$coef)]
    re_se <- VC$se[-(1:length(x0$coef))] #not used currently 
    x0$varVC <- list(NULL,VC$VC[names(x0$vars),names(x0$vars)]) #This is hard coded for 2 levels, with a NULL because the are no refs at level 1
    # build the var mat output, first adding SEvcov
    varDF <- x0$varDF
    varDF$SEvcov <- NA
    names(re_se) <- gsub(":",".",names(re_se),fixed=TRUE) # change : to .
    for(i in 1:nrow(varDF)) {
      if(varDF$fullGroup[i] %in% names(re_se)) {
        varDF$SEvcov[i] <- re_se[varDF$fullGroup[i]] 
      }
    }
  } else{
    se <- object$SE
    re_se <- rep(0,length(x0$vars))
    # build the var mat output from this
    varDF <- x0$varDF
  }
  object$coef <- as.matrix(data.frame(Estimate=x0$coef, "Std. Error"=se, "t value"=x0$coef/se))
  # make.names breaks/fixes colnames with spaces, un-break/-fix them.
  colnames(object$coef)[2:3] <- c("Std. Error", "t value")
  rownames(object$coef) <- names(x0$coef)
  # build the var mat output
  varsmat <- varDF[is.na(x0$varDF$var2),c("level","grp", "var1","vcov","SEvcov")]
  varsmat$st <- sqrt(varsmat$vcov)
  # add variance estimates
  colnames(varsmat) <- c("Level", "Group", "Name", "Variance", "Std. Error", "Std.Dev.")
  for(li in 2:x0$levels) {
    vc <- as.matrix(x0$varVC[[li]])
    cr <- cov2cor(vc)
    if(ncol(vc)>1) {
      for(i in 2:ncol(cr)) {
        for(j in 1:(i-1)){
          varsmat[varsmat$Level==li & varsmat$Name==rownames(cr)[i],paste0("Corr",j)] <- cr[i,j]
        }
      }
    }
  }
  object$varsmat <- varsmat
  object$vars <- cbind(x0$vars, re_se, sqrt(x0$vars))
  class(object) <- "summaryWeMixResults"
  return(object)
}

# This helper funciton creates a coefficient matrix from WeMix results.
#' @export
#' @importFrom stats printCoefmat
print.summaryWeMixResults <- function(x, ...) {
  cat("Call:\n")
  print(x$call)
  cat("\nVariance terms:\n")
  varsmat <- x$varsmat
  i <- 1
  while(paste0("Corr",i) %in% colnames(varsmat)) {
    # round correlations to two digits
    varsmat[[paste0("Corr",i)]] <- as.character(round(varsmat[[paste0("Corr",i)]],2))
    i <- i + 1
  }
  print(varsmat, na.print="", row.names=FALSE, digits=3)
  cat("Groups:\n")
  print(x$ngroups)
  cat("\nFixed Effects:\n")
  printCoefmat(x$coef, ...)
  cat("\nlnl=",format(x$lnl,nsmall=2),"\n")
  if(length(x$ICC) > 0) {
    cat("Intraclass Correlation=",format(x$ICC, nsmall=3, digits=3),"\n")
  }
}


#' This function calculates robust standard errors using the sandwich estimator of the standard errors  
#' following Rabe-Hesketh & Skrondal 2006. For linear models, robust standard errors are returned from the main estimation function. 
#' This function is only used for post-hoc estimation for non-linear models. 
#' @importFrom numDeriv jacobian
#' @keywords internal
makeSandwich <- function(fittedModel) {
  #make same estimator based on 

  par <- c(fittedModel$coef, fittedModel$vars)
  FisherInfInverse <- solve(fittedModel$invHessian)
  
  if (fittedModel$robustSE ==  TRUE){
    L <- jacobian(fittedModel$lnlf, par, top=FALSE)
    nGroups <- nrow(L)
    nParams <- length(par)
    J <- matrix(0,ncol=nParams, nrow=nParams)
    for (i in 1:nGroups){
      J <- J + L[i,] %*% t(L[i,])
    }
    J <- (nGroups/(nGroups-1))*J
    SE <- FisherInfInverse%*%J%*%FisherInfInverse
    colnames(SE) <- rownames(SE)<-names(par)
    se <- sqrt(diag(SE)) 
  
    #return NA for obs with variance below threshold 
    se[which(log(fittedModel$vars) < -3.6) + length(fittedModel$coef) ] <- NA
    diag(SE) <- se^2 
    
    names(se) <- colnames(SE)
  } else {
    
    SE <- -1* FisherInfInverse
    colnames(SE) <- rownames(SE)<-names(par)
    se <- sqrt(diag(SE))
    names(se) <- names(par)
  }
  return(list(VC=SE,se=se)) 
 

}