R/varcomp.R

In predictmeans: Predicted Means for Linear and Semiparametric Models

varcomp <- function(model, ci=TRUE, level=0.95) {
  if (inherits(model, "merMod")) {
    varcomp <- as.data.frame(VarCorr(model), order = "lower.tri")
    vcov_index <- -1:-length(fixef(model))
    random_vcov <- as.matrix(vcov_lmerMod(model, full = TRUE, ranpar = "var")[vcov_index, vcov_index])
    random_terms <- gsub("cov_", "", colnames(random_vcov))
    varcomp$SE <- sqrt(diag(random_vcov))
	if (ci){
    varcomp <- cbind(varcomp, confint(model, level=level)[1:length(random_terms), ])
    kk <- abs(varcomp$sdcor^2 - varcomp$vcov) > 0.0000000000001
    varcomp[, 7] <- ifelse(!kk, ifelse(varcomp[, 7] < 0, 0, varcomp[, 7]^2), varcomp[, 7]*varcomp$vcov/varcomp$sdcor)
    varcomp[, 8] <- ifelse(!kk, varcomp[, 8]^2, varcomp[, 8]*varcomp$vcov/varcomp$sdcor)	
    varcomp <- round(varcomp[,c(4,6:8)], 4)
	}else{
	  varcomp <- round(varcomp[,c(4,6)], 4)	
	}
	
    rownames(varcomp) <- random_terms
  }else if(inherits(model, "lme")){
    
    vcov=unlist(extract_varcomp(model))
    SE=sqrt(diag(varcomp_vcov(model)))
    
    # cbind(Std.error, vcov)
    confint_list <- try(intervals(model, which="var-cov", level=level), silent = TRUE)
    if (!inherits(confint_list, "try-error")){
      confint_matrix <- as.matrix(rbind(do.call("rbind", confint_list[[1]]), confint_list$sigma))
      confint_names <- c(as.vector(sapply(confint_list[[1]], rownames)), "residual")
      rownames(confint_matrix) <- confint_names
      confint_matrixN <- confint_matrix
      confint_matrixN["residual", ] <- confint_matrix["residual",]^2 
      sd_detect <- grepl("^sd", confint_names, fixed = FALSE)
      if (any(sd_detect)) {
        sd_names <- sub("\\)$", "", sub("^sd\\(", "", confint_names[sd_detect]))
        confint_matrixN[sd_detect, ] <- confint_matrix[sd_detect,]^2
      }
      
      cor_detect <- grepl("^cor", confint_names, fixed = FALSE) 
	  est. <- NULL
      if (any(cor_detect)) {
        cor_names <- sub("\\)$", "", sub("^cor\\(", "", confint_names[cor_detect]))
        sd_value <- confint_matrix[sd_detect, "est."]
        names(sd_value) <- sd_names
        sd_prod <- sapply(strsplit(cor_names, ","), function(x){sdv <- sd_value[x]; sdv[1]*sdv[2]})
        # print(sd_prod)
        # print(confint_matrixN[cor_detect, ])
        confint_matrixN[cor_detect, ] <- confint_matrix[cor_detect,]*sd_prod
      }
      rownames(confint_matrixN) <- as.character(round(confint_matrixN[,"est."], 4))
      confint_matrixN <- confint_matrixN[as.character(round(vcov, 4)),]
      varcomp <- cbind(vcov, SE, confint_matrixN)
	  varcomp <- subset(varcomp, select=-c(est.))
	  rownames(varcomp) <- sub("^Tau.", "", rownames(varcomp))
    }else{
      varcomp <- cbind(vcov, SE)      
    }
  }else{
    stop("The model must be a lme or lmer object!")
  }
  return(varcomp)  
}