R/summarySimResult.R

In simsem: SIMulated Structural Equation Modeling

# A collection of summary*, get*, and set* functions for the result object


# summaryParam: This function will summarize the obtained parameter estimates
# and standard error.

summaryParam <- function(object, alpha = 0.05, std = FALSE, detail = FALSE,
                         improper = TRUE, digits = NULL, matchParam = FALSE) {
  object <- clean(object, improper = improper)
  usedCoef <- object@coef
  usedSe <- object@se
  if (std) {
    if (length(object@stdCoef) == 0L)
      stop("The standardized coefficients cannot be summarized because there",
           " are no standardized coefficients in the object")
    usedCoef <- object@stdCoef
    usedSe <- object@stdSe
  }
  coef <- colMeans(usedCoef, na.rm = TRUE)
  real.se <- sapply(usedCoef, sd, na.rm = TRUE)
  result <- cbind(coef, real.se)
  colnames(result) <- c("Estimate Average", "Estimate SD")
  ## 26 Feb 2017, Terry moved above if-block, as suggested by Minkko
    estimated.se <- colMeans(usedSe, na.rm = TRUE)
    estimated.se[estimated.se == 0] <- NA
    crit <- qnorm(1 - alpha/2)
  if (!all(is.na(usedSe))) {
    z <- usedCoef / usedSe
    sig <- abs(z) > crit
    pow <- apply(sig, 2, mean, na.rm = TRUE)
    result <- cbind(result, "Average SE" = estimated.se, "Power (Not equal 0)" = pow)
  }
  if (!std && length(object@stdCoef) != 0) {
    stdCoef <- colMeans(object@stdCoef, na.rm = TRUE)
    stdRealSE <- sapply(object@stdCoef, sd, na.rm = TRUE)
    stdEstSE <- colMeans(object@stdSe, na.rm = TRUE)
    leftover <- setdiff(rownames(result), names(stdCoef))
    if (length(leftover) > 0) {
      temp <- rep(NA, length(leftover))
      names(temp) <- leftover
      stdCoef <- c(stdCoef, temp)
      stdRealSE <- c(stdRealSE, temp)
      stdEstSE <- c(stdEstSE, temp)
    }
    resultStd <- cbind("Std Est" = stdCoef, "Std Est SD" = stdRealSE,
                       "Std Ave SE" = stdEstSE)
    result <- cbind(result, resultStd[rownames(result),])
  }

  paramExist <- !(all(dim(object@paramValue) == 1) && is.na(object@paramValue))
  stdParamExist <- !(all(dim(object@stdParamValue) == 1) && is.na(object@stdParamValue))

  if ((!std & paramExist) | (std & stdParamExist)) {
    paramValue <- object@paramValue
    if (std) paramValue <- object@stdParamValue
    targetVar <- match(colnames(usedCoef), colnames(paramValue))
    targetVar <- targetVar[!is.na(targetVar)]
    paramValue <- paramValue[, targetVar]
    if (matchParam) result <- result[colnames(paramValue),]
    if ((nrow(result) == ncol(paramValue)) && all(rownames(result) ==
                                                  colnames(paramValue))) {
      nRep <- object@nRep
      nParam <- nrow(result)
      if (nrow(paramValue) == 1)
        paramValue <- matrix(unlist(rep(paramValue, nRep)), nRep, nParam,
                             byrow = T)
      biasParam <- usedCoef[,rownames(result)] - paramValue
      lowerBound <- object@cilower
      upperBound <- object@ciupper
      if (std) {
        lowerBound <- usedCoef - crit * usedSe
        upperBound <- usedCoef + crit * usedSe
      }
      selectci <- colnames(lowerBound) %in% rownames(result)
      lowerBound <- lowerBound[,selectci]
      upperBound <- upperBound[,selectci]

      noci <- setdiff(rownames(result), colnames(lowerBound))
      if (length(noci) > 0) {
        if (length(selectci) > 0) warning("Some CIs are Wald CI and others",
                                          " are calculated inside the simulation.")
        selectCoef <- usedCoef[,noci]
        selectSE <- usedSe[,noci]

        lowerBound <- cbind(lowerBound, selectCoef - crit * selectSE)
        upperBound <- cbind(upperBound, selectCoef + crit * selectSE)
      }
      cover <- (paramValue > lowerBound) & (paramValue < upperBound)
      average.param <- apply(paramValue, 2, mean, na.rm = TRUE)
      sd.param <- apply(paramValue, 2, sd, na.rm = TRUE)
      average.bias <- apply(biasParam, 2, mean, na.rm = TRUE)
      perc.cover <- apply(cover, 2, mean, na.rm = TRUE)
      sd.bias <- apply(biasParam, 2, sd, na.rm = TRUE)
      perc.cover[estimated.se[rownames(result)] == 0] <- NA
      result2 <- cbind(average.param, sd.param, average.bias, sd.bias, perc.cover)
      colnames(result2) <- c("Average Param", "SD Param", "Average Bias", "SD Bias",
                             "Coverage")
      if (nrow(object@paramValue) == 1)
        result2 <- result2[, c(1, 3, 5)]
      result <- cbind(result, result2)
      if (detail) {
        relative.bias <- biasParam/paramValue
        relBias <- apply(relative.bias, 2, mean, na.rm = TRUE)
        relBias[is.nan(relBias)] <- NA
        std.bias <- NULL
        relative.bias.se <- NULL
        if (nrow(object@paramValue) == 1) {
          std.bias <- average.bias/real.se
          relative.bias.se <- (estimated.se - real.se) / real.se
        } else {
          std.bias <- average.bias/sd.bias
          relative.bias.se <- (estimated.se - sd.bias) / sd.bias
        }
        width <- upperBound - lowerBound
        average.width <- apply(width, 2, mean, na.rm = TRUE)
        sd.width <- apply(width, 2, sd, na.rm = TRUE)
        belowLowerBound <- paramValue < lowerBound
        aboveUpperBound <- paramValue > upperBound
        perc.lower <- apply(belowLowerBound, 2, mean, na.rm = TRUE)
        perc.upper <- apply(aboveUpperBound, 2, mean, na.rm = TRUE)
        result3 <- cbind(relBias, std.bias, relative.bias.se, perc.lower,
                         perc.upper, average.width, sd.width)
        colnames(result3) <- c("Rel Bias", "Std Bias", "Rel SE Bias",
                               "Not Cover Below", "Not Cover Above",
                               "Average CI Width", "SD CI Width")
        result <- cbind(result, result3)
      }
    }
  }

  if (nrow(object@FMI1) > 1 & ncol(object@FMI1) >= 1) {
    FMI1 <- object@FMI1
    average.FMI1 <- apply(FMI1, 2, mean, na.rm = TRUE)
    sd.FMI1 <- apply(FMI1, 2, sd, na.rm = TRUE)
    resultFMI <- cbind(average.FMI1, sd.FMI1)
    colnames(resultFMI) <- c("Average FMI1", "SD FMI1")
    targetParam <- matchLavaanName(rownames(result), rownames(resultFMI))
    targetParam <- targetParam[!is.na(targetParam)]
    result <- cbind(result, resultFMI[targetParam,])
  }

  if (nrow(object@FMI2) > 1 & ncol(object@FMI2) >= 1) {
    FMI2 <- object@FMI2
    average.FMI2 <- apply(FMI2, 2, mean, na.rm = TRUE)
    sd.FMI2 <- apply(FMI2, 2, sd, na.rm = TRUE)

    resultFMI <- cbind(average.FMI2, sd.FMI2)
    colnames(resultFMI) <- c("Average FMI2", "SD FMI2")
    targetParam <- matchLavaanName(rownames(result), rownames(resultFMI))
    targetParam <- targetParam[!is.na(targetParam)]
    result <- cbind(result, resultFMI[targetParam,])
  }

  if (length(unique(object@n)) > 1) {
    corCoefN <- cor(cbind(usedCoef, object@n),
                    use = "pairwise.complete.obs")[colnames(usedCoef), "object@n"]
    corSeN <- cor(cbind(usedSe, object@n),
                  use = "pairwise.complete.obs")[colnames(usedSe), "object@n"]
    result <- cbind(result, r_coef.n = corCoefN, r_se.n = corSeN)
  }
  if (length(unique(object@pmMCAR)) > 1) {
    corCoefMCAR <- cor(cbind(usedCoef, object@pmMCAR),
                       use = "pairwise.complete.obs")[colnames(usedCoef), "object@pmMCAR"]
    corSeMCAR <- cor(cbind(usedSe, object@pmMCAR),
                     use = "pairwise.complete.obs")[colnames(usedSe), "object@pmMCAR"]
    result <- cbind(result, r_coef.pmMCAR = corCoefMCAR, r_se.pmMCAR = corSeMCAR)
  }
  if (length(unique(object@pmMAR)) > 1) {
    corCoefMAR <- cor(cbind(usedCoef, object@pmMAR),
                      use = "pairwise.complete.obs")[colnames(usedCoef), "object@pmMAR"]
    corSeMAR <- cor(cbind(usedSe, object@pmMAR),
                    use = "pairwise.complete.obs")[colnames(usedSe), "object@pmMAR"]
    result <- cbind(result, r_coef.pmMAR = corCoefMAR, r_se.pmMAR = corSeMAR)
  }
  if (!is.null(digits)) {
    result <- round(result, digits)
  }
  return(as.data.frame(result))
}

matchLavaanName <- function(name1, name2) {
  result1 <- match(name1, name2)
  result2 <- match(name1, switchLavaanName(name2))
  result1[is.na(result1)] <- result2[is.na(result1)]
  result1
}

switchLavaanName <- function(name) {
  nameX <- strsplit(name, "\\.")
  nameGroup <- sapply(nameX, "[", 1)
  nameParam <- sapply(nameX, "[", 2)
  target <- grep("~~", nameParam)
  varTarget <- strsplit(nameParam[target], "~~")
  nameParam[target] <- paste0(sapply(varTarget, "[", 2), "~~", sapply(varTarget, "[", 1))
  nameX <- paste0(nameGroup, ".", nameParam)
  nameX
}

# summaryMisspec: This function will summarize the obtained fit indices and
# generate a data frame.

summaryMisspec <- function(object, improper = TRUE) {
  object <- clean(object, improper = improper)
  if (all(dim(object@misspecValue) == 0)) {
    stop("This object does not have any model misspecification.")
  } else {
    misspecAverage <- colMeans(object@misspecValue, na.rm = TRUE)
    misspecSE <- sapply(object@misspecValue, sd, na.rm = TRUE)
    popAverage <- colMeans(object@popFit, na.rm = TRUE)
    popSE <- sapply(object@popFit, sd, na.rm = TRUE)
    mis <- data.frame(mean = misspecAverage, sd = misspecSE)
    pop <- data.frame(mean = popAverage, sd = popSE)
    return(rbind(pop, mis))
  }
}

# summaryPopulation: Summarize population values behind data generation model

summaryPopulation <- function(object, std = FALSE, improper = TRUE) {
  object <- clean(object, improper = improper)
  paramValue <- object@paramValue
  if (std) {
    if (all(dim(object@stdParamValue) == 1) && is.na(object@stdParamValue))
      stop("The standardized parameters cannot be summarized because there are",
           " no standardized parameters in the object")
    paramValue <- object@stdParamValue
  }
  nRep <- nrow(paramValue)
  nParam <- ncol(paramValue)
  result <- NULL
  if (nrow(paramValue) == 1) {
    result <- matrix(paramValue, nrow = 1)
    rownames(result) <- "Population Value"
    colnames(result) <- colnames(paramValue)
  } else {
    average.param <- apply(paramValue, 2, mean, na.rm = TRUE)
    sd.param <- apply(paramValue, 2, sd, na.rm = TRUE)
    result <- rbind(average.param, sd.param)
    rownames(result) <- c("Average", "SD")
  }
  return(result)
}

# summaryFit: This function will summarize the obtained fit indices and
# generate a data frame.

summaryFit <- function(object, alpha = NULL, improper = TRUE, usedFit = NULL) {
  cleanObj <- clean(object, improper = improper)
  usedFit <- cleanUsedFit(usedFit, colnames(object@fit))
  condition <- c(length(unique(object@pmMCAR)) > 1, length(unique(object@pmMAR)) >
                   1, length(unique(object@n)) > 1)
  if (any(condition)) {
    if (is.null(alpha))
      alpha <- 0.05
    values <- list()
    ifelse(condition[3], values[[3]] <- round(seq(min(object@n), max(object@n),
                                                  length.out = 5)), values[[3]] <- NA)
    ifelse(condition[1], values[[1]] <- seq(min(object@pmMCAR), max(object@pmMCAR),
                                            length.out = 5), values[[1]] <- NA)
    ifelse(condition[2], values[[2]] <- seq(min(object@pmMAR), max(object@pmMAR),
                                            length.out = 5), values[[2]] <- NA)
    m <- do.call(expand.grid, values)
    FUN <- function(vec, obj, alpha, usedFit) {
      as.numeric(getCutoff(obj, alpha, revDirec = FALSE, usedFit = usedFit,
                           nVal = vec[3], pmMCARval = vec[1], pmMARval = vec[2]))
    }
    cutoffs <- sapply(as.data.frame(t(m)), FUN, obj = object, alpha = alpha,
                      usedFit = usedFit)
    mSelect <- as.matrix(m[, condition])
    colnames(mSelect) <- c("%MCAR", "%MAR", "N")[condition]
    rownames(cutoffs) <- usedFit
    colnames(cutoffs) <- NULL
    result <- data.frame(mSelect, t(cutoffs))
    rownames(result) <- NULL
  } else {
    if (is.null(alpha))
      alpha <- c(0.1, 0.05, 0.01, 0.001)
    cutoffs <- sapply(alpha, getCutoff, object = cleanObj, usedFit = usedFit)
    cutoffs <- matrix(as.numeric(cutoffs), length(usedFit), length(alpha))
    if (ncol(as.matrix(cutoffs)) == 1) {
      cutoffs <- t(cutoffs)
      # rownames(cutoffs) <- usedFit
    }

    fit <- as.data.frame(cleanObj@fit[, usedFit])
    meanfit <- apply(fit, 2, mean, na.rm = TRUE)
    sdfit <- apply(fit, 2, sd, na.rm = TRUE)
    if (length(alpha) == 1) cutoffs <- t(cutoffs)
    result <- cbind(cutoffs, meanfit, sdfit)
    colnames(result) <- c(alpha, "Mean", "SD")
    rownames(result) <- usedFit
    names(dimnames(result)) <- c("Fit Indices", "Alpha")
    # print(as.data.frame(cutoffs))
  }

  return(result)
}

# summaryMisspec: This function will summarize the obtained fit indices and
# generate a data frame.

summaryConverge <- function(object, std = FALSE, improper = TRUE) {
  result <- list()
  converged <- object@converged == 0
  numnonconverged <- sum(!converged)
  if (numnonconverged == 0)
    stop("You are good! All replications were converged!")
  result <- c(result, list(Converged = c(num.converged = sum(converged),
                                         num.nonconverged = numnonconverged)))
  reasons <- c("Nonconvergent" = sum(object@converged %in% 1:2),
               "Improper SE" = sum(object@converged == 3),
               "Nonpositive definite matrix" = sum(object@converged == 4),
               # "Improper Variance" = sum(object@converged == 4),
               # "Improper Correlation" = sum(object@converged == 5),
               # "Not-positive-definite model-implied covariance matrix of latent variables" = sum(object@converged == 6),
               "Optimal estimates were not guaranteed" = sum(object@converged == 7))
  reasons <- as.matrix(reasons)
  colnames(reasons) <- "count"
  result <- c(result, list("Nonconvergent Reasons" = reasons))
  n <- object@n
  pmMCAR <- object@pmMCAR
  pmMAR <- object@pmMAR
  paramValue <- object@paramValue
  if (std) {
    if (all(dim(object@stdParamValue) == 1) && is.na(object@stdParamValue))
      stop("The standardized parameters cannot be summarized because there are",
           " no standardized parameters in the object")
    paramValue <- object@stdParamValue
  }
  misspecValue <- object@misspecValue
  popFit <- object@popFit
  nonconverged <- !converged
  improprep <- rep(FALSE, length(converged))
  if (improper) {
    nonconverged <- object@converged %in% 1:2
    improprep <- object@converged %in% 3:7
  }
  if (length(unique(n)) > 1) {
    temp1 <- n[converged]
    temp2 <- n[nonconverged]
    temp3 <- n[improprep]
    resultTemp <- c(mean.convergence = mean(temp1), sd.convergence = sd(temp1))
    resultDiff <- NULL
    if (length(temp2) > 0) {
      resultTemp <- c(resultTemp, mean.nonconvergence = mean(temp2), sd.nonconvergence = sd(temp2))
      resultDiff <- c(resultDiff, diff.non.mean = mean(temp2) - mean(temp1), diff.non.sd = sd(temp2) - sd(temp1))
    }
    if (length(temp3) > 0) {
      resultTemp <- c(resultTemp, mean.improper = mean(temp3),
                      sd.improper = sd(temp3))
      resultDiff <- c(resultDiff, diff.improper.mean = mean(temp3) - mean(temp1),
                      diff.improper.sd = sd(temp3) - sd(temp1))
    }
    result <- c(result, list(n = c(resultTemp, resultDiff)))
  }
  if (length(unique(pmMCAR)) > 1) {
    temp1 <- pmMCAR[converged]
    temp2 <- pmMCAR[nonconverged]
    temp3 <- pmMCAR[improprep]
    resultTemp <- c(mean.convergence = mean(temp1), sd.convergence = sd(temp1))
    resultDiff <- NULL
    if (length(temp2) > 0) {
      resultTemp <- c(resultTemp, mean.nonconvergence = mean(temp2),
                      sd.nonconvergence = sd(temp2))
      resultDiff <- c(resultDiff, diff.non.mean = mean(temp2) - mean(temp1),
                      diff.non.sd = sd(temp2) - sd(temp1))
    }
    if (length(temp3) > 0) {
      resultTemp <- c(resultTemp, mean.improper = mean(temp3),
                      sd.improper = sd(temp3))
      resultDiff <- c(resultDiff, diff.improper.mean = mean(temp3) - mean(temp1),
                      diff.improper.sd = sd(temp3) - sd(temp1))
    }
    result <- c(result, list(pmMCAR = c(resultTemp, resultDiff)))
  }
  if (length(unique(pmMAR)) > 1) {
    temp1 <- pmMAR[converged]
    temp2 <- pmMAR[nonconverged]
    temp3 <- pmMAR[improprep]
    resultTemp <- c(mean.convergence = mean(temp1), sd.convergence = sd(temp1))
    resultDiff <- NULL
    if (length(temp2) > 0) {
      resultTemp <- c(resultTemp, mean.nonconvergence = mean(temp2),
                      sd.nonconvergence = sd(temp2))
      resultDiff <- c(resultDiff, diff.non.mean = mean(temp2) - mean(temp1),
                      diff.non.sd = sd(temp2) - sd(temp1))
    }
    if (length(temp3) > 0) {
      resultTemp <- c(resultTemp, mean.improper = mean(temp3),
                      sd.improper = sd(temp3))
      resultDiff <- c(resultDiff, diff.improper.mean = mean(temp3) - mean(temp1),
                      diff.improper.sd = sd(temp3) - sd(temp1))
    }
    result <- c(result, list(pmMAR = c(resultTemp, resultDiff)))
  }
  if (nrow(paramValue) > 1) {
    temp1 <- paramValue[converged, , drop = FALSE]
    temp2 <- paramValue[nonconverged, , drop = FALSE]
    temp3 <- paramValue[improprep, , drop = FALSE]
    resultTemp <- cbind(mean.convergence = apply(temp1, 2, mean),
                        sd.convergence = apply(temp1, 2, sd))
    resultDiff <- NULL
    if (nrow(temp2) > 0) {
      resultTemp <- cbind(resultTemp, mean.nonconvergence = apply(temp2, 2, mean),
                          sd.nonconvergence = apply(temp2, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.non.mean = apply(temp2, 2, mean) - apply(temp1, 2, mean),
                          diff.non.sd = apply(temp2, 2, sd) - apply(temp1, 2, sd))
    }
    if (nrow(temp3) > 0) {
      resultTemp <- cbind(resultTemp, mean.improper = apply(temp3, 2, mean),
                          sd.improper = apply(temp3, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.improper.mean = apply(temp3, 2, mean) - apply(temp1, 2, mean),
                          diff.improper.sd = apply(temp3, 2, sd) - apply(temp1, 2, sd))
    }
    result <- c(result, list(paramValue = cbind(resultTemp, resultDiff)))
  }
  if (!all(dim(misspecValue) == 0) && nrow(misspecValue) > 1) {
    temp1 <- misspecValue[converged, ]
    temp2 <- misspecValue[nonconverged, ]
    temp3 <- misspecValue[improprep, ]
    resultTemp <- cbind(mean.convergence = apply(temp1, 2, mean),
                        sd.convergence = apply(temp1, 2, sd))
    resultDiff <- NULL
    if (nrow(temp2) > 0) {
      resultTemp <- cbind(resultTemp, mean.nonconvergence = apply(temp2, 2, mean),
                          sd.nonconvergence = apply(temp2, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.non.mean = apply(temp2, 2, mean) - apply(temp1, 2, mean),
                          diff.non.sd = apply(temp2, 2, sd) - apply(temp1, 2, sd))
    }
    if (nrow(temp3) > 0) {
      resultTemp <- cbind(resultTemp, mean.improper = apply(temp3, 2, mean),
                          sd.improper = apply(temp3, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.improper.mean = apply(temp3, 2, mean) - apply(temp1, 2, mean),
                          diff.improper.sd = apply(temp3, 2, sd) - apply(temp1, 2, sd))
    }
    result <- c(result, list(misspecValue = cbind(resultTemp, resultDiff)))
  }
  if (!all(dim(popFit) == 0) && nrow(popFit) > 1) {
    temp1 <- popFit[converged, ]
    temp2 <- popFit[nonconverged, ]
    temp3 <- popFit[improprep, ]
    resultTemp <- cbind(mean.convergence = apply(temp1, 2, mean),
                        sd.convergence = apply(temp1, 2, sd))
    resultDiff <- NULL
    if (nrow(temp2) > 0) {
      resultTemp <- cbind(resultTemp, mean.nonconvergence = apply(temp2, 2, mean),
                          sd.nonconvergence = apply(temp2, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.non.mean = apply(temp2, 2, mean) - apply(temp1, 2, mean),
                          diff.non.sd = apply(temp2, 2, sd) - apply(temp1, 2, sd))
    }
    if (nrow(temp3) > 0) {
      resultTemp <- cbind(resultTemp, mean.improper = apply(temp3, 2, mean),
                          sd.improper = apply(temp3, 2, sd))
      resultDiff <- cbind(resultDiff,
                          diff.improper.mean = apply(temp3, 2, mean) - apply(temp1, 2, mean),
                          diff.improper.sd = apply(temp3, 2, sd) - apply(temp1, 2, sd))
    }
    result <- c(result, list(popFit = cbind(resultTemp, resultDiff)))
  }
  return(result)
}

# setPopulation: Set population parameter values

setPopulation <- function(target, population) {
  popParam <- NULL
  if (is(population, "SimSem")) {
    psl <- generate(population, n = 20, params = TRUE)$psl
    paramSet <- lapply(psl, "[[", 1)
    generatedgen <- population@dgen
    if (!is.list(generatedgen[[1]])) {
      generatedgen <- list(generatedgen)
    }
    for (i in seq_along(paramSet)) {
      popParam <- c(popParam, parsePopulation(generatedgen[[i]],
                                              paramSet[[i]], group = i))
    }
    extraParamIndex <- population@pt$op %in% c(">", "<", "==", ":=")
    extraParamName <- NULL
    if (any(extraParamIndex)) {
      extraparam <- collapseExtraParam(paramSet, population@dgen, fill = TRUE,
                                       con = population@con)
      extraParamName <- renameExtraParam(population@con$lhs, population@con$op,
                                         population@con$rhs)
      popParam[extraParamIndex] <- extraparam
    }
    index <- ((population@pt$free != 0) & !(duplicated(population@pt$free))) | extraParamIndex
    popParam <- popParam[index]
    tempFit4Names <- lavaan::lavaan(population@pt,
                                    sample.nobs = rep(200, max(population@pt$group)))
    names(popParam) <- c(names(coef(tempFit4Names)), extraParamName)
  } else if (is(population, "lavaan")) {
    pt <- parTable(population)
    popParam <- pt$est
    names(popParam) <- lavaan::lav_partable_labels(pt)
  } else if (is(population, "character")) {
    pt <- parTable(lavaan::lavaan(population))
    popParam <- pt$est
    names(popParam) <- lavaan::lav_partable_labels(pt)
  } else if (is(population, "list")) {
    pt <- population
    popParam <- pt$est
    names(popParam) <- lavaan::lav_partable_labels(pt)
  } else {
    stop("population must be a simsem object, lavaan object, lavaan script",
         " for data generation, or parameter table")
  }
  if (is.vector(popParam)) {
    popParam <- popParam[!is.na(popParam)]
    tmp <- names(popParam)
    popParam <- data.frame(t(matrix(popParam)))
    colnames(popParam) <- tmp
  } else if (is.matrix(popParam)) {
    popParam <- data.frame(popParam)
  }
  target@paramValue <- popParam
  return(target)
}

# getPopulation: Description: Extract the population value from an object

getPopulation <- function(object, std = FALSE, improper = TRUE, nonconverged = FALSE) {
  toextract <- "param"
  if (std) toextract <- "stdparam"
  inspect(object, toextract, improper = improper, nonconverged = nonconverged)
}

# getExtraOutput: Extract the extra output that users set in the 'outfun' argument

getExtraOutput <- function(object, improper = TRUE, nonconverged = FALSE) {
  targetRep <- 0
  if (improper) targetRep <- c(targetRep, 3:7)
  if (nonconverged) targetRep <- c(targetRep, 1:2)
  if (length(object@extraOut) == 0) {
    stop("This simulation result does not contain any extra results")
  } else {
    return(object@extraOut[object@converged %in% targetRep])
  }
}

setMethod("coef", "SimResult",
          function(object, improper = TRUE, nonconverged = FALSE) {
            inspect(object, "coef", improper = improper, nonconverged = nonconverged)
          })

setGeneric("inspect",
           function(object, what="coef", improper = TRUE, nonconverged = FALSE)
             standardGeneric("inspect"))
setMethod("inspect", "SimResult", function(object, what = "coef",
                                           improper = TRUE, nonconverged = FALSE) {
  targetRep <- 0
  if (improper) targetRep <- c(targetRep, 3:7)
  if (nonconverged) targetRep <- c(targetRep, 1:2)
  targetRep <- object@converged %in% targetRep

  if (length(what) > 1) {
    stop("`what' arguments contains multiple arguments; only one is allowed")
  }

  # be case insensitive
  what <- tolower(what)

  if (what == "type" ||
      what == "model.type" ||
      what == "modeltype") {
    return(object@modelType)
  } else if (what == "nrep" ||
             what == "numrep") {
    return(object@nRep)
  } else if (what == "param" ||
             what == "paramvalue" ||
             what == "parameter" ||
             what == "parameter.value" ||
             what == "paramvalues" ||
             what == "parameter.values" ||
             what == "parameters") {
    target <- object@paramValue
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "stdparam" ||
             what == "stdparamvalue" ||
             what == "stdparameter" ||
             what == "stdparameter.value" ||
             what == "std.param" ||
             what == "std.param.value" ||
             what == "std.param.values" ||
             what == "std.parameters" ||
             what == "stdparamvalues" ||
             what == "std.parameter.values" ||
             what == "standardized.parameters") {
    target <- object@stdParamValue
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "se" ||
             what == "std.err" ||
             what == "standard.errors") {
    return(object@se[targetRep, , drop=FALSE])
  } else if (what == "misspec" ||
             what == "misspecification" ||
             what == "misspecified" ||
             what == "misspecified.param") {
    target <- object@misspecValue
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "coef" ||
             what == "coefficients" ||
             what == "parameter.estimates" ||
             what == "estimates" ||
             what == "x" ||
             what == "est") {
    return(object@coef[targetRep, , drop = FALSE])
  } else if (what == "popfit" ||
             what == "popmisfit" ||
             what == "popmis" ||
             what == "misfit") {
    target <- object@popFit
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "fmi" ||
             what == "fmi1") {
    target <- object@FMI1
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "fmi2") {
    target <- object@FMI2
    if (nrow(target) > 1) target <- target[targetRep, , drop = FALSE]
    return(target)
  } else if (what == "fit" ||
             what == "fitmeasures" ||
             what == "fit.measures" ||
             what == "fit.indices") {
    return(object@fit[targetRep, , drop=FALSE])
  } else if (what == "std" ||
             what == "std.coef" ||
             what == "standardized" ||
             what == "standardizedsolution" ||
             what == "standardized.solution") {
    return(object@stdCoef[targetRep, , drop=FALSE])
  } else if (what == "stdse" ||
             what == "std.se" ||
             what == "standardizedse" ||
             what == "standardized.se") {
    return(object@stdSe[targetRep, , drop=FALSE])
  } else if (what == "cilower" ||
             what == "ci.lower" ||
             what == "lowerci" ||
             what == "lower.ci") {
    return(object@cilower[targetRep, , drop=FALSE])
  } else if (what == "ciupper" ||
             what == "ci.upper" ||
             what == "upperci" ||
             what == "upper.ci") {
    return(object@ciupper[targetRep, , drop=FALSE])
  } else if (what == "ciwidth" ||
             what == "ci.width" ||
             what == "widthci" ||
             what == "width.ci" ||
             what == "width") {
    lower <- as.matrix(object@cilower)
    upper <- as.matrix(object@ciupper)
    if (nrow(lower) <= 1) stop("There are no lower and upper bounds of",
                               " confidence intervals saved in this simulation.")
    return(upper[targetRep, , drop = FALSE] - lower[targetRep, , drop = FALSE])
  } else if (what == "seed") {
    return(summarySeed(object))
  } else if (what == "n" ||
             what == "samplesize" ||
             what == "groupsize" ||
             what == "groupn" ||
             what == "ngroup" ||
             what == "sample.size") {
    return(object@nobs[targetRep, , drop=FALSE])
  } else if (what == "ntotal" ||
             what == "totaln" ||
             what == "totalsamplesize" ||
             what == "total.sample.size") {
    return(object@n[targetRep])
  } else if (what == "mcar" ||
             what == "pmmcar" ||
             what == "pmcar") {
    target <- object@pmMCAR
    if (length(target) > 1) target <- target[targetRep]
    return(target)
  } else if (what == "mar" ||
             what == "pmmar" ||
             what == "pmar") {
    target <- object@pmMAR
    if (length(target) > 1) target <- target[targetRep]
    return(target)
  } else if (what == "extra" ||
             what == "extraout" ||
             what == "misc") {
    return(getExtraOutput(object, improper = improper, nonconverged = nonconverged))
  } else if (what == "time"  ||
             what == "timing") {
    return(summaryTime(object))
  } else if (what == "converged") {
    lab <- c("converged", "nonconverged", "nonconvergedMI", "improperSE", "improperVariance", "improperCorrelation", "improperCovLv", "nonOptimal")
    lab <- lab[sort(unique(object@converged)) + 1]
    out <- factor(object@converged, labels = lab)
    return(out)
  } else {
    stop("unknown `what' argument in inspect function: `", what, "'")
  }

})

# summaryPopulation: Summarize population values behind data generation model

summaryTime <- function(object, units = "seconds") {
  timing <- object@timing
  timing1 <- timing[-which(names(timing) %in% c("StartTime", "EndTime"))]
  units <- tolower(units)
  if (units %in% c("secs", "seconds", "second", "sec", "s")) {
    # Do nothing
  } else if (units %in% c("min", "mins", "minute", "minutes", "m")) {
    timing1 <- lapply(timing1, "/", 60)
  } else if (units %in% c("hours", "hour", "h")) {
    timing1 <- lapply(timing1, "/", 3600)
  } else if (units %in% c("days", "day", "d")) {
    timing1 <- lapply(timing1, "/", 86400)
  } else {
    stop("unknown `units' argument in summaryTime function: `", units, "'")
  }

  #format(.POSIXct(dt,tz="GMT"), "%H:%M:%S")
  cat("============ Wall Time ============\n")
  t0.txt <- sprintf("  %-72s", "1. Error Checking and setting up data-generation and analysis template:")
  t1.txt <- sprintf("  %10.3f", timing1$SimulationParams)
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-72s", "2. Set combinations of n, pmMCAR, and pmMAR:")
  t1.txt <- sprintf("  %10.3f", timing1$RandomSimParams)
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-72s", "3. Setting up simulation conditions for each replication:")
  t1.txt <- sprintf("  %10.3f", timing1$SimConditions)
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-72s", "4. Total time elapsed running all replications:")
  t1.txt <- sprintf("  %10.3f", timing1$RunReplications)
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-72s", "5. Combining outputs from different replications:")
  t1.txt <- sprintf("  %10.3f", timing1$CombineResults)
  cat(t0.txt, t1.txt, "\n", "\n", sep="")

  cat("============ Average Time in Each Replication ============\n")
  t0.txt <- sprintf("  %-46s", "1. Data Generation:")
  t1.txt <- sprintf("  %10.3f", timing1$InReps[1])
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-46s", "2. Impose Missing Values:")
  t1.txt <- sprintf("  %10.3f", timing1$InReps[2])
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-46s", "3. User-defined Data-Transformation Function:")
  t1.txt <- sprintf("  %10.3f", timing1$InReps[3])
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-46s", "4. Main Data Analysis:")
  t1.txt <- sprintf("  %10.3f", timing1$InReps[4])
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-46s", "5. Extracting Outputs:")
  t1.txt <- sprintf("  %10.3f", timing1$InReps[5])
  cat(t0.txt, t1.txt, "\n", "\n", sep="")

  cat("============ Summary ============\n")
  t0.txt <- sprintf("  %-28s", "Start Time:")
  t1.txt <- format(timing$StartTime, "  %Y-%m-%d %H:%M:%S")
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-28s", "End Time:")
  t1.txt <- format(timing$EndTime, "  %Y-%m-%d %H:%M:%S")
  cat(t0.txt, t1.txt, "\n", sep="")
  totaltime <- Reduce("+", timing1[-which("InReps" == names(timing1))])
  t0.txt <- sprintf("  %-28s", "Wall (Actual) Time:")
  t1.txt <- sprintf("  %10.3f", totaltime)
  cat(t0.txt, t1.txt, "\n", sep="")
  inreptime <- max(sum(timing1$InReps * object@nRep), timing1$RunReplications)
  systemtime <- totaltime - timing1$RunReplications + inreptime
  t0.txt <- sprintf("  %-28s", "System (Processors) Time:")
  t1.txt <- sprintf("  %10.3f", systemtime)
  cat(t0.txt, t1.txt, "\n", sep="")
  t0.txt <- sprintf("  %-28s", "Units:")
  cat(t0.txt, "  ", units, "\n", sep="")
}

summarySeed <- function(object) {
  seed <- object@seed
  list("Seed number" = object@seed[1], "L'Ecuyer seed of the last replication" = object@seed[-1])
}