Nothing
#' @title Get fit indices from objects
#' @description Get fit indices from objects for which a method exists.
#' @param x An object for which a method exists.
#' @param ... further arguments to be passed to or from other methods.
#' @return A data.frame.
#' @examples
#' \dontrun{
#' df <- iris[, 1, drop = FALSE]
#' names(df) <- "x"
#' res <- mx_mixture(model = "x ~ m{C}*1
#' x ~~ v{C}*x", classes = 1, data = df)
#' table_fit(res)
#' }
#' @export
get_fit <- function(x, ...) {
message("This function will be deprecated in the next version of tidySEM; use table_fit() instead.")
UseMethod("table_fit", x)
}
calc_fitindices <- function(model, ...){
UseMethod("calc_fitindices", model)
}
#' @export
calc_fitindices.MxModel <- function (model, type = NULL, ...){
sums <- .table_fit_mx(model, type = type)
if(isTRUE(type == "mixture")){
ll <- sums[["Minus2LogLikelihood"]]/-2
parameters <- sums[["estimatedParameters"]]
n <- sums[["numObs"]]
if(length(names(model@submodels)) < 2){
fits = c("Classes" = 1, "Entropy" = 1, "prob_min" = 1, "prob_max" = 1, "n_min" = 1, "n_max" = 1)
} else {
post_prob <- extract_postprob(model)
class <- apply(post_prob, 1, which.max)
class_tab <- table(factor(class, levels = 1:length(names(model@submodels))))
if (length(class_tab) == ncol(post_prob)) {
prop_n <- range(prop.table(class_tab))
}
else {
prop_n <- c(0, max(prop.table(class_tab)))
}
fits <- c("Classes" = length(names(model@submodels)),
ifelse(ncol(post_prob) == 1, 1, 1 + (1/(nrow(post_prob) *
log(ncol(post_prob)))) * (sum(rowSums(post_prob *
log(post_prob + 1e-12))))),
range(diag(classification_probs_mostlikely(post_prob, class))),
prop_n)
names(fits) <- c("Classes", "Entropy", "prob_min", "prob_max", "n_min", "n_max")
}
c(sums, fits)
} else {
sums
}
}
make_fitvector <- function(ll, parameters, n, postprob = NULL, fits = NULL){
out <- vector("numeric")
out["LogLik"] <- ll
out["parameters"] <- parameters
out["n"] <- n
out["AIC"] <- -2*ll + (2*parameters)
if(!is.null(fits)) out["AWE"] <- -2*(ll + unname(fits[1])) + 2*parameters*(3/2 + log(n))
out["BIC"] <- -2*ll + (parameters * log(n))
out["CAIC"] <- -2*ll + (parameters * (log(n)+1))
if(!is.null(fits)) out["CLC"] <- -2*ll + (2*unname(fits[1]))
out["KIC"] <- -2*ll + (3*(parameters + 1))
out["SABIC"] <- -2*ll + (parameters * log(((n+2)/24)))
if(!is.null(postprob)) out["ICL"] <- icl_default(postprob, out["BIC"])
c(out, fits)
}
#' @title Obtain latent class probabilities
#' @description Obtain latent class probabilities for an object for which a
#' method exists. See Details.
#' @details The following types are available:
#' \describe{
#' \item{"sum.posterior"}{A summary table of the posterior class
#' probabilities; this indicates what proportion of your data contributes to
#' each class.}
#' \item{"sum.mostlikely"}{A summary table of the most likely class
#' membership, based on the highest posterior class probability. Note that
#' this is subject to measurement error.}
#' \item{"mostlikely.class"}{If C is the true class of an observation, and N is
#' the most likely class based on the model, then this table shows the
#' probability P(N==i|C==j). The diagonal represents the probability that
#' observations in each class will be correctly classified.}
#' \item{"avg.mostlikely"}{Average posterior probabilities for each class, for
#' the subset of observations with most likely class of 1:k, where k is the
#' number of classes.}
#' \item{"individual"}{The posterior probability matrix, with dimensions n
#' (number of cases in the data) x k (number of classes).}
#' }
#' @param x An object for which a method exists.
#' @param type Character vector, indicating which types of probabilities to
#' extract. See Details.
#' @param ... Further arguments to be passed to or from other methods.
#' @return A data.frame.
#' @examples
#' \dontrun{
#' df <- iris[, 1, drop = FALSE]
#' names(df) <- "x"
#' res <- mx_mixture(model = "x ~ m{C}*1
#' x ~~ v{C}*x", classes = 1, data = df)
#' class_prob(res)
#' }
#' @export
class_prob <- function(x, type = c("sum.posterior", "sum.mostlikely", "mostlikely.class", "avg.mostlikely", "individual"), ...){
UseMethod("class_prob", x)
}
#' @method class_prob MxModel
#' @export
class_prob.MxModel <- function(x, type = c("sum.posterior", "sum.mostlikely", "mostlikely.class", "avg.mostlikely", "individual"), ...){
post_probs <- extract_postprob(x)
post_probs_pred <- cbind(post_probs, predicted = apply(post_probs, 1, which.max) )
out <- lapply(type, function(thetype){
switch(thetype,
"mostlikely.class" = classification_probs_mostlikely(post_probs),
"avg.mostlikely" = avgprobs_mostlikely(post_probs),
"sum.posterior" = sum_postprob(x),
"sum.mostlikely" = sum_mostlikely(x),
post_probs_pred)
})
names(out) <- type
out
}
sum_postprob <- function(model){
if(!inherits(model$expectation, "MxExpectationMixture")){
return(data.frame(class = "class1", count = model$data$numObs, proportion = 1))
}
weights_name <- model$expectation$weights
priors <- prop.table(model[[weights_name]]$values)
mix_names <- model$expectation$components
numObs <- model$data$numObs
data.frame(class = mix_names, count = as.vector(priors*numObs), proportion = as.vector(priors))
}
sum_mostlikely <- function(model){
if(!inherits(model$expectation, "MxExpectationMixture")){
return(return(data.frame(class = "class1", count = model$data$numObs, proportion = 1)))
}
post_prob <- extract_postprob(model)
classif <- table(apply(post_prob, 1, which.max))
weights_name <- model$expectation$weights
mix_names <- model$expectation$components
data.frame(class = mix_names, count = as.vector(classif), proportion = as.vector(prop.table(classif)))
}
extract_postprob <- function(model){
UseMethod("extract_postprob", model)
}
#' @export
extract_postprob.MxModel <- function(model){
if(!inherits(model$expectation, "MxExpectationMixture")){
return(matrix(1, nrow = model$data$numObs))
}
weights_name <- model$expectation$weights
priors <- model[[weights_name]]$values
mix_names <- model$expectation$components
liks <- sapply(mix_names, function(x){ attr(model@submodels[[x]], "fitfunction")$result })
#First calculate posteriors:
posteriors <- sweep(liks, 2, priors, "*")
#Then calculate participants' marginal likelihoods (probably not be necessary):
marglik <- rowSums(posteriors)
# Divide posteriors by marginal likelihood
posteriors <- posteriors/marglik
# then normalize them into probabilities:
posteriors/rowSums(posteriors)
}
#' @export
extract_postprob.mplus.model <- function(model){
model$savedata[grepl("^CPROB", names(model$savedata))]
}
classification_probs_mostlikely <- function (post_prob, class = NULL)
{
if (is.null(dim(post_prob)))
return(1)
if (is.null(class))
class <- apply(post_prob, 1, which.max)
avg_probs <- avgprobs_mostlikely(post_prob, class)
avg_probs[is.na(avg_probs)] <- 0
class_counts <- as.integer(table(ordered(class, levels = 1:ncol(post_prob)))) # Use ordered to ensure empty classes are included
tab <- avg_probs * class_counts
tab <- tab / matrix(colSums(avg_probs * class_counts), ncol = ncol(tab), nrow = nrow(tab), byrow = TRUE)
rownames(tab) <- paste0("assigned.", 1:nrow(tab))
colnames(tab) <- paste0("avgprob.", 1:nrow(tab))
return(t(tab))
}
avgprobs_mostlikely <- function (post_prob, class = NULL)
{
if (is.null(dim(post_prob)))
return(1)
if(is.null(class)) class <- apply(post_prob, 1, which.max)
tab <- t(sapply(1:ncol(post_prob), function(i) {
colMeans(post_prob[class == i, , drop = FALSE])
}))
rownames(tab) <- paste0("assigned.", 1:nrow(tab))
colnames(tab) <- paste0("meanprob.", colnames(tab))
return(tab)
}
# ICL method for mplus.model ----------------------------------------------
icl.mplus.model <- function(object, ...)
{
n <- object$summaries$Observations
z <- object$savedata[, grep("^CPROB", names(object$savedata))]
if(!is.null(dim(z))){
C <- z == apply(z, 1, max)
(-1*object$summaries$BIC) + 2*sum(C * log(apply(z, 1, function(x){x[which.max(x)]})+1e-12))
} else {
(-1*object$summaries$BIC) + 2*sum(z * log(z))
}
}
icl_default <- function(post_prob, BIC){
tryCatch({
if (!is.null(dim(post_prob))) {
C <- post_prob == apply(post_prob, 1, max)
(-1 * BIC) + 2 * sum(C * log(apply(post_prob,
1, function(x) {
x[which.max(x)]
}) + 1e-12))
}
else {
(-1 * BIC) + 2 * sum(post_prob * log(post_prob))
}
}, error = function(e){ NA })
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.