version 0.2-0

cran · Mar 11, 2019 · e8fe110 · e8fe110
1 parent 5361ceb
commit e8fe110
Show file tree

Hide file tree

Showing 16 changed files with 222 additions and 100 deletions.
diff --git a/DESCRIPTION b/DESCRIPTION
@@ -1,6 +1,6 @@
 Package: MPTmultiverse
 Title: Multiverse Analysis of Multinomial Processing Tree Models
-Version: 0.1
+Version: 0.2-0
 Description: 
     Statistical or cognitive modeling usually requires a number of more or less 
     arbitrary choices creating one specific path through a 'garden of forking paths'. 
@@ -30,18 +30,19 @@ URL: https://github.com/mpt-network/MPTmultiverse
 BugReports: https://github.com/mpt-network/MPTmultiverse/issues
 Depends: R (>= 2.11.1),
 Imports: parallel, magrittr, tidyr, dplyr, tibble, rlang, reshape2,
-        ggplot2, MPTinR, TreeBUGS, runjags, coda, purrr, broom, readr
+        ggplot2, MPTinR, TreeBUGS, runjags, coda, purrr, readr,
+        limSolve
 Suggests: knitr, rmarkdown, testthat
 LazyData: yes
 VignetteBuilder: knitr
-RoxygenNote: 6.1.0
+RoxygenNote: 6.1.1
 License: GPL-2
 NeedsCompilation: no
-Packaged: 2018-11-11 22:31:44 UTC; henrik
+Packaged: 2019-03-11 21:24:50 UTC; henrik
 Author: Henrik Singmann [aut, cre] (<https://orcid.org/0000-0002-4842-3657>),
   Daniel W. Heck [aut],
   Marius Barth [aut],
   Frederik Aust [ctb] (<https://orcid.org/0000-0003-4900-788X>)
 Maintainer: Henrik Singmann <singmann@gmail.com>
 Repository: CRAN
-Date/Publication: 2018-11-19 18:20:11 UTC
+Date/Publication: 2019-03-11 23:52:42 UTC
diff --git a/MD5 b/MD5
@@ -1,19 +1,20 @@
-7635fa884787afbb4dc3192fd4521976 *DESCRIPTION
-2149fa0da9ce2823a90580d233e4c21e *NAMESPACE
-da8c321c996700f376f6a3d69f7a0dbb *R/check_results.R
+0aef9496de27a3a45758b4d3da803151 *DESCRIPTION
+358c7c50c5ff4ec4854c1f7003fd0d78 *NAMESPACE
+c58b5b86bc3755c9cdaa14f15e36b8ec *R/check_results.R
 762033a6f12bdd6e72ce504f4a153dd3 *R/fit_mpt.R
-268d2b32c4dda924ddaa7e0cb570479b *R/get_eqn_categories.R
-af7946c208cd196e88af53373a4f33c4 *R/make_results_row.R
-0df8672d644e4d8c6453ae5d4b8eb88e *R/mpt_options.R
-8a1667a7cf0820e062aba0a874e52272 *R/mptinr.R
+a2de0d9664ddbf8e4fa1e7af30b05f3b *R/get_eqn_categories.R
+d8604ad9aa9ca4ecfbb5af1f6747c42a *R/make_results_row.R
+ffc2311b2c3894ea2fedd163f2bb996b *R/mpt_options.R
+8e8fac68ec6a93596de3cfee0a4d0125 *R/mptinr.R
 c0860dc3bb02743f1ff358a5c9afed03 *R/onload.R
-49dc6af309539032e9cb10f5da3faf9e *R/plot_multiverseMPT.R
+5bdbe313d4f00a95ac9e6359030c6836 *R/plot_multiverseMPT.R
 faf7d231083e8038c413991ca81d3b2c *R/prep_data_fitting.R
-7c851b5ee634302bb8742988a71b65b8 *R/treebugs.R
+9dcca834565312ea9744c6dc927241ad *R/treebugs.R
+7ba90e1c6fef8537e2ed9d137a801f96 *R/write_results.R
 cf1d73fdd11829c6652ca98c8def0c50 *README.md
-4ec29cb38de0b055cf1edf21e1e10fd9 *build/vignette.rds
+6dd1b953436a21b77ba666b4c357e285 *build/vignette.rds
 5ddc745fc0e241ffe1a2801a5fdf854a *inst/doc/introduction-bayen_kuhlmann_2011.R
-117c3dcc8768cfdc1c51eb766478067e *inst/doc/introduction-bayen_kuhlmann_2011.html
+8aa59c2a33a86224ab4fe840fedd51bb *inst/doc/introduction-bayen_kuhlmann_2011.html
 8fb0a3da73fb660d23e6e7fb0b8f1dd8 *inst/doc/introduction-bayen_kuhlmann_2011.rmd
 378c931503b45c257da3b87e4c5e7490 *inst/extdata/2HTSM_Submodel4.eqn
 31de4fabf5cb90e987b2d418a49fed02 *inst/extdata/Kuhlmann_dl7.csv
@@ -25,13 +26,14 @@ d9ebea39a1adb10102dbe3cc044395dc *inst/extdata/prospective_memory.eqn
 a7e9b12163013673deb6c74abeb4cd55 *man/fit_mpt.Rd
 9a3023bfc66e5b420f7df718499fb2fa *man/get_pb_output.Rd
 8f3811a251eb5a19d2bbfeea686f1d85 *man/make_results_row.Rd
-1f46692603cd31496c0614cf11b4fdef *man/mpt_options.Rd
+91e20cb185f30a2236eeb41c8108fb7d *man/mpt_options.Rd
 70d444459d24aac210c2b85acb17d520 *man/plot.multiverseMPT.Rd
 a859ce4d222c2a86d9a3cc7a04364d4a *man/write_check_results.Rd
+65b8904cc88c4eea226b24f46b9be466 *man/write_results.Rd
 dd4bfb25654cdb367a4c4485a712e332 *tests/testthat.R
 9c4f3d2af60ba67b298f959e6e07735f *tests/testthat/test-data-structures.R
 16a5cbce2a4d7eceb5a9ba2d0cdac561 *tests/testthat/test-identifiability_checks.R
-a5854ab6a762c544f0f64035ddbf0c08 *tests/testthat/test-mptinr.R
+bd1e0f4d37b5f4d79272abc5d149d85d *tests/testthat/test-mptinr.R
 15b9a91654c2a73e11ff9b74ca2b83e8 *tests/testthat/test_treebugs.R
 7e3b16937fce40bab6396676351497c7 *vignettes/NOT_USED/1_bayen_kuhlmann_2011.html
 378c931503b45c257da3b87e4c5e7490 *vignettes/NOT_USED/2HTSM_Submodel4.eqn

diff --git a/NAMESPACE b/NAMESPACE
@@ -5,6 +5,7 @@ export(check_results)
 export(fit_mpt)
 export(mpt_options)
 export(write_check_results)
+export(write_results)
 importFrom(graphics,plot)
 importFrom(magrittr,"%>%")
 importFrom(rlang,.data)

diff --git a/R/check_results.R b/R/check_results.R
@@ -31,37 +31,38 @@ check_results <- function(results) {
 
   cat("## MPTinR: no pooling\n")
 
+  mpt_no_pool <- c("asymptotic", "PB/MLE", "NPB/MLE")
+  mpt_no_pool <- mpt_no_pool[mpt_no_pool %in% results$method]
   tryCatch({
-    for(meth in c("asymptotic", "PB/MLE", "NPB/MLE")){
+    for(meth in mpt_no_pool){
 
-      conv_mptinr_no <- results %>% 
-        dplyr::filter(.data$package == "MPTinR" & .data$pooling == "no" & .data$method == meth) %>% 
-        dplyr::select("convergence") %>% 
-        tidyr::unnest()
+      # conv_mptinr_no <- results %>% 
+      #   dplyr::filter(.data$package == "MPTinR" & .data$pooling == "no" & .data$method == meth) %>% 
+      #   dplyr::select("convergence") %>% 
+      #   tidyr::unnest()
 
       not_id <- results %>% 
         dplyr::filter(.data$package == "MPTinR" & .data$pooling == "no" & .data$method == meth) %>% 
         dplyr::select("est_indiv") %>% 
         tidyr::unnest() %>% 
-        dplyr::group_by(.data$condition) %>% 
+        dplyr::group_by(.data$condition, .data$core) %>% 
         dplyr::summarise(proportion = mean(!.data$identifiable))
 
       not_id2 <- results %>% 
         dplyr::filter(.data$package == "MPTinR" & .data$pooling == "no" & .data$method == meth) %>% 
         dplyr::select("est_indiv") %>% 
         tidyr::unnest() %>%
         dplyr::filter(!.data$identifiable) %>% 
-        dplyr::group_by(.data$condition) %>% 
-        dplyr::summarise(not_identified = list(broom::tidy(table(.data$parameter)))) %>% 
-        tidyr::unnest(.data$not_identified) %>% 
-        suppressWarnings()
-
+        dplyr::group_by(.data$condition, .data$core, .data$parameter) %>% 
+        dplyr::count() %>% 
+        dplyr::ungroup()
+
       if (any(not_id$proportion > 0)) {
         cat("Based on", meth, "method, proportion of participants with non-identified parameters:\n")
-        cat(format(not_id)[-c(1,3)], "", sep = "\n")
+        cat(format(not_id, n = Inf)[-c(1,3)], "", sep = "\n")
 
         cat("Based on", meth, "CIs, table of non-identified parameters:\n")
-        cat(format(not_id2)[-c(1,3)], sep = "\n")
+        cat(format(not_id2, n = Inf)[-c(1,3)], sep = "\n")
 
       } else {
         cat("Based on", meth, "CIs, all parameters of all participants seem to be identifiable.\n")
@@ -87,9 +88,12 @@ check_results <- function(results) {
     comp_prob <- (conv_mptinr_comp$convergence != 0) | 
       (conv_mptinr_comp$rank.fisher != conv_mptinr_comp$n.parameters)
 
-    if (any(comp_prob)) {
+    if (any(comp_prob, na.rm = TRUE)) {
       cat("Convergence problems:\n")
       cat(format(conv_mptinr_comp[comp_prob,])[-c(1,3)], "", sep = "\n")
+    } else if (any(is.na(comp_prob))) {
+      cat("Convergence problems:\n")
+      cat(format(conv_mptinr_comp[is.na(comp_prob),])[-c(1,3)], "", sep = "\n")
     } else {
       cat("No convergence problems.\n")
     }

diff --git a/R/get_eqn_categories.R b/R/get_eqn_categories.R
@@ -17,7 +17,7 @@ get_eqn_categories <- function (model.filename)
   tmp.ordered <- tmp.in[order(tmp.in$V1), , drop = FALSE]
   tmp.spl <- split(tmp.ordered, factor(tmp.ordered$V1))
   tmp.spl <- lapply(tmp.spl, function(d.f) d.f[order(d.f[, 2]), ])
-  unlist(lapply(tmp.spl, function(x) unique(x$V2)))
+  as.character(unlist(lapply(tmp.spl, function(x) unique(x$V2))))
   # model <- lapply(tmp.spl, parse.eqn)
   # names(model) <- NULL
   # model
@@ -38,5 +38,5 @@ get_eqn_trees <- function(model_file) {
   n_cat <- lapply(X = tmp_split, FUN = function(x) {
     length(unique(x$V2))
   })
-  as.character(mapply(FUN = rep, each = n_cat, x = names(tmp_split)))
+  as.character(unlist(mapply(FUN = rep, each = n_cat, x = names(tmp_split))))
 }
diff --git a/R/make_results_row.R b/R/make_results_row.R
@@ -38,7 +38,6 @@ make_results_row <- function(
   package,
   method,
   data,
-  # parameters,
   id,
   condition,
   core = NULL  # character vector specifying which are core parameters
@@ -94,8 +93,8 @@ make_results_row <- function(
     colnames(est_group)[ncol(est_group)] <- paste0("ci_", getOption("MPTmultiverse")$ci_size[i])
   }
 
-
-  # group comparisons
+  # ----------------------------------------------------------------------------
+  # test_between: group comparisons
   if (length(conditions) > 1) {
 
     pairs <- utils::combn(
@@ -128,8 +127,22 @@ make_results_row <- function(
     }
     test_between <- dplyr::bind_rows(tmp_test_between) 
   } else {
-    test_between <- tibble::tibble()
+    # Return a zero-row tibble if no between-Ss condition is analyzed ----
+    test_between <- tibble::tibble(
+      parameter = character(0)
+      , core = logical(0)
+      , condition1 = character(0)
+      , condition2 = character(0)
+      , est_diff = numeric(0)
+      , se = numeric(0)
+      , p = numeric(0)
+    )
+    CI <- getOption("MPTmultiverse")$ci_size
+    for (i in seq_along(CI)) {
+      test_between[[paste0("ci_", CI[i])]] <- numeric(0)
+    }
   }
+
   ## est_covariate <- ##MISSING
 
   if (method == "trait"){

diff --git a/R/mpt_options.R b/R/mpt_options.R
@@ -21,7 +21,7 @@
 #'   \item{\code{max_ci_indiv}: }{Numeric. Used for excluding individual parameter estimates in the bootstrap approaches. If the range of the CI (i.e., distance between minimum and maximum) is larger than this value, the estimate is excluded from the group-level estimates.}
 #'   \item{\code{silent_jags}: }{Logical. Whether to suppress JAGS output.}
 # ' TODO  \item{\code{catch_warnings}: }{Logical. Whether to store warnings and errors as additional columns in the output.}
-#'   \item{\code{save_models}: }{Logical.}
+#'   \item{\code{save_models}: }{Logical. Default is \code{FALSE} which does not save the individual MCMC samples in \code{.RData} files. Instead only summairzes are retained in \code{results} object.}
 #' }
 #'   
 #' 
@@ -127,7 +127,7 @@ set_default_options <- function() {
     , ci_size = c(.025, .1, .9, .975)
     , max_ci_indiv = .99
     , n.CPU = parallel::detectCores()
-    , save_models = TRUE
+    , save_models = FALSE
   )
 }
 

diff --git a/R/mptinr.R b/R/mptinr.R
@@ -57,8 +57,9 @@ mpt_mptinr <- function(
     )
   }
 
-  for(i in 1:length(res)) {
-    res[[i]]$test_homogeneity[[1]] <- homogeneity_tests
+  for(i in seq_len(length(res))) {
+    if (nrow(res[[i]]) > 0) 
+      res[[i]]$test_homogeneity[[1]] <- homogeneity_tests
   }
 
   # return
@@ -333,7 +334,7 @@ mpt_mptinr_no <- function(
 
 
   if ("pb" %in% bootstrap) {
-    res[["pb_no"]] <- get_pb_results(dataset = dataset
+    try(res[["pb_no"]] <- get_pb_results(dataset = dataset
                                      , prepared = prepared
                                      , model = model
                                      , id = id
@@ -342,10 +343,10 @@ mpt_mptinr_no <- function(
                                      , fit_mptinr = fit_mptinr
                                      , additional_time = additional_time
                                      , convergence = convergence
-                                     , core = core)
+                                     , core = core))
   }
   if ("npb" %in% bootstrap) {
-    res[["npb_no"]] <- get_pb_results(dataset = dataset
+    try(res[["npb_no"]] <- get_pb_results(dataset = dataset
                                      , prepared = prepared
                                      , model = model
                                      , id = id
@@ -354,7 +355,7 @@ mpt_mptinr_no <- function(
                                      , fit_mptinr = fit_mptinr
                                      , additional_time = additional_time
                                      , convergence = convergence
-                                     , core = core)
+                                     , core = core))
   }
   # return
   dplyr::bind_rows(res)
@@ -740,7 +741,10 @@ mpt_mptinr_complete <- function(dataset,
       fit_mptinr_tmp$parameters[
         res$est_group[[1]][res$est_group[[1]]$condition == prepared$conditions[i], ]$parameter, "estimates"]
 
-    par_se <- sqrt(diag(solve(fit_mptinr_tmp$hessian[[1]])))
+    par_se <- rep(NA_real_, length(rownames(fit_mptinr_tmp$parameters)))
+    par_se <- tryCatch(sqrt(diag(solve(fit_mptinr_tmp$hessian[[1]]))),
+                       error = function(x) 
+                         sqrt(diag(limSolve::Solve(fit_mptinr_tmp$hessian[[1]]))))
     names(par_se) <- rownames(fit_mptinr_tmp$parameters)
 
     res$est_group[[1]][

diff --git a/R/plot_multiverseMPT.R b/R/plot_multiverseMPT.R
@@ -105,7 +105,14 @@
 
 plot.multiverseMPT <- function(x, which = "est", save = FALSE, ...){
 
-  shapes <- c(16, 18, 15, 1, 0, 8, 11, 12, 4, 6)
+  args <- list(...)
+
+  if(is.null(args$shapes)) {
+    shapes <- seq_len(nrow(x))
+  } else {
+    shapes <- args$shapes
+  }
+
 
   results <- x
   prefix <- paste0(gsub("\\.eqn", "", results$model[1]), "_", 

diff --git a/R/treebugs.R b/R/treebugs.R
@@ -222,7 +222,9 @@ mpt_treebugs <- function (
         for(p in parameters){
           test_between <- TreeBUGS::betweenSubjectMPT(treebugs_fit[[i]], treebugs_fit[[j]], 
                                             par1 = p, stat = "x-y")
-          test_summ <- TreeBUGS::summarizeMCMC(test_between$mcmc, probs = CI_SIZE)
+          test_summ <- TreeBUGS::summarizeMCMC(test_between$mcmc, 
+                                               probs = CI_SIZE, 
+                                               batchSize = 2)
           bayesp <- mean(do.call("rbind", test_between$mcmc) <= 0)
 
           sel_row <- 
@@ -269,7 +271,8 @@ mpt_treebugs <- function (
 
     for (i in seq_along(conditions)){
       mcmc <- treebugs_fit[[i]]$runjags$mcmc[,sel_rho]
-      rho_summ <- TreeBUGS::summarizeMCMC(mcmc, probs = CI_SIZE)
+      rho_summ <- TreeBUGS::summarizeMCMC(mcmc, probs = CI_SIZE,
+                                          batchSize = 2)
       bayesp <- colMeans(do.call("rbind", mcmc) <= 0)
       res <- data.frame(par_mat,
                         rho_summ[,c("Mean", "SD")], 

diff --git a/R/write_results.R b/R/write_results.R
@@ -0,0 +1,37 @@
+#' Write Results of Multiverse Analysis to csv-Files
+#' 
+#' Exports the results to csv format. 
+#' 
+#' @param results An object of class multiverseMPT.
+#' @param path a path where to save the files (e.g., \code{"C:/results/modelX_dataY_"})
+#' @importFrom rlang .data
+#' @export
+write_results <- function(results, path = "MPTmultiverse_"
+                          # what = c("est_group", "est_indiv", "est_rho", 
+                          #          "test_between", "gof", "gof_group", "gof_indiv", 
+                          #          "fungibility", "test_homogeneity")
+                          ){
+
+  # TODO: allow to specify which columns should be unnested and exported
+  #       => requires to work with expressions/tidyverse-issues
+
+  readr::write_csv(tidyr::unnest(results, .data$est_group), 
+                   path = paste0(path,"est_group.csv"))
+  readr::write_csv(tidyr::unnest(results, .data$est_indiv), 
+                   path = paste0(path,"est_indiv.csv"))
+  readr::write_csv(tidyr::unnest(results, .data$est_rho), 
+                   paste0(path,"est_rho.csv"))
+
+  readr::write_csv(tidyr::unnest(results, .data$test_between), 
+                   paste0(path,"test_between.csv"))
+
+  readr::write_csv(tidyr::unnest(results, .data$gof), 
+                   paste0(path,"gof.csv"))
+  readr::write_csv(tidyr::unnest(results, .data$gof_indiv), 
+                   paste0(path,"gof_indiv.csv"))
+  readr::write_csv(tidyr::unnest(results, .data$gof_group), 
+                   paste0(path,"gof_group.csv"))
+
+  readr::write_csv(tidyr::unnest(results, .data$fungibility), 
+                   paste0(path,"fungibility.csv"))
+}
diff --git a/build/vignette.rds b/build/vignette.rds