if(getOption("test_mplus")) { oldwd <- getwd() testdir <- file.path(tempdir(), "comparefitstats") dir.create(testdir) setwd(testdir) on.exit({ setwd(oldwd) unlink(testdir, recursive = TRUE) }) library(tibble) iris_df <- iris names(iris_df) <- gsub("\\.", "_", names(iris_df)) test_that("fit stats are the same for both mplus and mclust", { tab_mclust <- tibble::tribble( ~ stat, ~ val, "Model", 1, "Classes", 2, "LogLik", -271.2943387, "parameters", 7, "n", 32, "AIC", 556.5886774, "AWE", 610.4725311, "BIC", 566.8488288, "CAIC", 573.8488288, "CLC", 544.2251264, "KIC", 566.5886774, "SABIC", 545.0268243, "ICL", -569.7289846, "Entropy", 0.818224466, "prob_min", 0.95731999, "prob_max", 0.960443609, "n_min", 0.46875, "n_max", 0.53125, "BLRT_val", 22.53901766, "BLRT_p", 0.00990099 ) # for benchmarking # createMixtures(classes = 2, filename_stem = "mtcars-mpg-hp", rdata = mtcars[, c("mpg", "hp")]) # runModels(target = "mtcars-mpg-hp_2_class.inp") # o <- readModels("mtcars-mpg-hp_2_class.out") tab_benchmark <- tibble::tribble( # these are from benchmarking code above ~ stat, ~ val, "Model", 1, "Classes", 2, "LogLik", -271.292, "parameters", 7, "n", 32, "AIC", 556.584, "BIC", 566.844, "Entropy", 0.819, "BLRT_val", 22.544 ) tmp <- capture_output({ m_cars_mplus <- estimate_profiles( mtcars[, c("mpg", "hp")], n_profiles = 2, models = 1, package = "MplusAutomation" ) }) m_cars_mclust <- estimate_profiles(mtcars[, c("mpg", "hp")], n_profiles = 2, models = 1) m_cars_mx <- estimate_profiles(mtcars[, c("mpg", "hp")], n_profiles = 2, models = 1, package = "openmx") expect_equal(as.vector(m_cars_mplus$model_1_class_2$fit), tab_mclust$val, tolerance = .0001) expect_equal(as.vector(m_cars_mclust$model_1_class_2$fit), tab_mclust$val, tolerance = .0001) expect_equal( as.vector(m_cars_mplus$model_1_class_2$fit), as.vector(m_cars_mclust$model_1_class_2$fit), tolerance = .0001 ) expect_equal( as.vector(m_cars_mplus$model_1_class_2$fit[1:18]), as.vector(m_cars_mx$model_1_class_2$fit[1:18]), tolerance = .001 ) expect_equal( as.vector(m_cars_mclust$model_1_class_2$fit[1:18]), as.vector(m_cars_mx$model_1_class_2$fit[1:18]), tolerance = .001 ) kept_cases <- complete.cases(pisaUSA15[, c("broad_interest", "enjoyment", "self_efficacy")]) tmp <- capture_output({ m_pisa_mplus <- estimate_profiles( pisaUSA15[kept_cases, c("broad_interest", "enjoyment", "self_efficacy")], n_profiles = 3, models = 1, package = "MplusAutomation", keepfiles = TRUE ) }) m_pisa_mclust <- estimate_profiles(pisaUSA15[kept_cases, c("broad_interest", "enjoyment", "self_efficacy")], n_profiles = 3, models = 1) # THis test is NOT equal, probably due to local maximum? expect_equal( as.vector(m_pisa_mplus$model_1_class_3$fit), as.vector(m_pisa_mclust$model_1_class_3$fit), tolerance = 3 ) m_cars_mplus_benchmark_stats <- tibble( stat = names(m_cars_mplus$model_1_class_2$fit), val = m_cars_mplus$model_1_class_2$fit ) m_cars_mclust_benchmark_stats <- tibble( stat = names(m_cars_mclust$model_1_class_2$fit), val = m_cars_mclust$model_1_class_2$fit ) m_cars_mplus_benchmark_stats <- m_cars_mplus_benchmark_stats[m_cars_mplus_benchmark_stats$stat %in% tab_benchmark$stat,] m_cars_mclust_benchmark_stats <- m_cars_mclust_benchmark_stats[m_cars_mclust_benchmark_stats$stat %in% tab_benchmark$stat,] expect_equivalent(m_cars_mplus_benchmark_stats$val, tab_benchmark$val, tolerance = .0001) expect_equivalent(m_cars_mclust_benchmark_stats$val, tab_benchmark$val, tolerance = .0001) }) test_that("fit stats are equal to those from models estimated externally in Mplus", { tmp <- capture_output({ # Model 1 suppressWarnings({ m_cars_mplus_model1 <- estimate_profiles( iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 1, package = "MplusAutomation" ) m_cars_mclust_model1 <- estimate_profiles(iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 1) # Model 2 m_cars_mplus_model2 <- estimate_profiles( iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 2, package = "MplusAutomation" ) m_cars_mclust_model2 <- estimate_profiles(iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 2) # Model 3 m_cars_mplus_model3 <- estimate_profiles( iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 3, package = "MplusAutomation" ) m_cars_mclust_model3 <- estimate_profiles(iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 3) # Model 6 m_cars_mplus_model6 <- estimate_profiles( iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 6, package = "MplusAutomation" ) m_cars_mclust_model6 <- estimate_profiles(iris_df[, c("Sepal_Length", "Sepal_Width", "Petal_Length", "Petal_Width")], n_profiles = 2, models = 6) }) }) iris_log_lik <- tribble(~ model, ~ logLik, 1, -488.915, 2, -386.185, 3, -296.448, 6, -214.355) iris_log_lik$mplus <- c( m_cars_mplus_model1$model_1_class_2$fit[3], m_cars_mplus_model2$model_2_class_2$fit[3], m_cars_mplus_model3$model_3_class_2$fit[3], m_cars_mplus_model6$model_6_class_2$fit[3] ) iris_log_lik$mclust <- c( m_cars_mclust_model1$model_1_class_2$model$loglik, m_cars_mclust_model2$model_2_class_2$model$loglik, m_cars_mclust_model3$model_3_class_2$model$loglik, m_cars_mclust_model6$model_6_class_2$model$loglik ) expect_equivalent(iris_log_lik$logLik, iris_log_lik$mplus, tolerance = .0001) expect_equivalent(iris_log_lik$logLik, iris_log_lik$mclust, tolerance = .0001) expect_equivalent(iris_log_lik$mplus, iris_log_lik$mclust, tolerance = .0001) }) }