create plotting function

R/fitConfModels.R

@@ -42,11 +42,10 @@
 #' @param n.cores `integer`. Number of cores used for parallelization. If NULL (default), the available
 #' number of cores -1 will be used.
 
-#' @return Gives data frame with one row for each combination of model and
-#' participant. Columns include a model and participant column,
-#' one column for each estimated parameter for the different models (parameters
-#' that are not present in a specific model (row) but in other models are
-#' filled with NAs.
+#' @return Gives `data.frame` with one row for each combination of model and
+#' participant. There are different columns for the model, the participant ID, and one
+#' one column for each estimated model parameter (parameters
+#' not present in a specific model are filled with NAs).
 #' Additional information  about the fit is provided in additional columns:
 #' - `negLogLik` (negative log-likelihood of the best-fitting set of parameters),
 #' - `k` (number of parameters),

R/plotConfModelFit.R

@@ -0,0 +1,143 @@
+#' @title Plot the prediction of fitted parameters of one model of confidence over the corresponding real data
+#'
+#' @description The `plotConfModelFit` function plots the predicted distribution of discrimination responses
+#' and confidence ratings created from a `data.frame` of parameters obtaind from \code{\link{fitConfModels}}
+#' and overlays the predicted distribution over the data to which the model parameters were fitted.
+#'
+#' @param data  a `data.frame` where each row is one trial, containing following
+#' variables:
+#' * \code{diffCond} (optional; different levels of discriminability,
+#'    should be a factor with levels ordered from hardest to easiest),
+#' * \code{rating} (discrete confidence judgments, should be a factor with levels
+#'    ordered from lowest confidence to highest confidence;
+#'    otherwise will be transformed to factor with a warning),
+#' * \code{stimulus} (stimulus category in a binary choice task,
+#'    should be a factor with two levels, otherwise it will be transformed to
+#'    a factor with a warning),
+#' * \code{correct} (encoding whether the response was correct; should  be 0 for
+#'    incorrect responses and 1 for correct responses)
+#' * \code{participant} (some group ID, most often a participant identifier;
+#'    the models given in the second argument are fitted to each subset of `data`
+#'    determined by the different values of this column)
+#'
+#' @param fitted_pars a `data.frame` with one row for each participant and model parameters in different columns.
+#' fitted_pars also may contain a column called  `model`specifying the model to be visualized.
+#' If there is no model column in data or if there are multiple models in fitted_pars,
+#' it is necessary to specify the model argument.
+#'
+#' @param model `character`. See \code{\link{fitConf}} for all available models
+#'
+#' @return a `data.frame` with one row for each combination of model and
+#' participant. There are different columns for the model, the participant ID, and one
+#' one column for each estimated model parameter (parameters
+#' not present in a specific model are filled with NAs)
+#'
+#' #' @examples
+#' # 1. Select two subjects from the masked orientation discrimination experiment
+#' data <- subset(MaskOri, participant %in% c(1:2))
+#' head(data)
+#'
+#' # 2. Fit some models to each subject of the masked orientation discrimination experiment
+#' \donttest{
+#'   # Fitting several models to several subjects takes quite some time
+#'   # (about 10 minutes per model fit per participant on a 2.8GHz processor
+#'   # with the default values of nInits and nRestart).
+#'   # If you want to fit more than just two subjects,
+#'   # we strongly recommend setting .parallel=TRUE
+#'   Fits <- fitConfModels(data, models = "ITGc", .parallel = FALSE)
+#' }
+#' # 3. Plot the predicted probabilies based on model and fitted parameter over the observed relative frequencies.
+#'
+#' \donttest{
+#'   # Fitting several models to several subjects takes quite some time
+#'   # (about 10 minutes per model fit per participant on a 2.8GHz processor
+#'   # with the default values of nInits and nRestart).
+#'   # If you want to fit more than just two subjects,
+#'   # we strongly recommend setting .parallel=TRUE
+#'   myPlottedFit <- plotConfModelFit(data, Fits)
+#'   myPlottedFit
+#' }
+#' @import ggplot2
+#' @importFrom plyr ddply transform summarise
+#' @importFrom Rmisc summarySEwithin
+#'
+#' @export
+plotConfModelFit <- function(data, fitted_pars, model = NULL){
+  if(is.null(model)){
+    if("model" %in% colnames(fitted_pars)){
+      if(length(unique(fitted_pars$model))==1){
+        model <- unique(fitted_pars$model)
+      } else {
+        stop("Please use the model argument to specify which model should be used")
+      }
+    }else {
+      stop("Please specify which model should be used")
+    }
+  }
+  if (is.null(data$diffCond)) data$diffCond <- factor(1)
+  if (!is.factor(data$diffCond)) {
+    data$diffCond <- factor(data$diffCond)
+  }
+  if(length(unique(data$stimulus)) != 2) {
+    stop("There must be exacltly two different possible values of stimulus")
+  }
+
+  if (!is.factor(data$stimulus)) {
+    data$stimulus <- factor(data$stimulus)
+  }
+  if (!is.factor(data$rating)) {
+    data$rating <- factor(data$rating)
+  }
+  if(!all(data$correct %in% c(0,1))) stop("correct should be 1 or 0")
+
+  myColor <- switch(model, 'GN' = 1, 'IG' = 2, 'ITGc'  = 3, 'ITGcm' = 4, 'logN' = 5,
+                    'logWEV' = 6,'PDA' = 7,  'WEV' = 8, 'SDT' = 9) # models are color coded
+
+  # 1. First aggregate on the level of subjects
+
+  AggDist <-
+    plyr::ddply(data,
+          ~  diffCond * rating * stimulus * correct * participant, #,
+          plyr::summarise, p = length(rating),  .drop=FALSE)
+
+  AggDist <- plyr::ddply(AggDist, ~ diffCond * stimulus,
+                         transform, N = sum(p))
+  AggDist$p <- AggDist$p / AggDist$N
+
+
+  AggDist$rating <- as.numeric(AggDist$rating)
+  AggDist$correct <-
+    factor(AggDist$correct)
+  levels(AggDist$correct)  <-
+    c("incorrect", "correct")
+
+   # 2. aggregate across subjects
+  AggDist <-
+    Rmisc::summarySEwithin(AggDist , measurevar = "p",
+                           withinvars = c("diffCond", "correct", "rating", "stimulus"),
+                           idvar = "participant",
+                           na.rm = TRUE, .drop = TRUE)
+  AggDist$rating <- as.numeric(AggDist$rating)
+  levels(AggDist$stimulus) <- c("S = -1", "S = 1")
+  levels(AggDist$diffCond) <- paste("K =", as.numeric(levels(AggDist$diffCond)))
+
+  # 3) create a plot with the observed data
+  PlotObsVsPred  <-
+    ggplot(AggDist,
+           aes(x=rating, y = p)) +
+    facet_grid(diffCond ~ stimulus+correct) +
+    geom_bar(stat="identity",
+             fill ="white", color = "black") +
+    geom_errorbar(aes(ymin=p-se, ymax=p+se), width=0) +
+    xlab("Confidence rating") +
+    ylab("probability") +
+    theme(strip.text.y = element_text(angle=0)) +
+    theme_minimal()
+
+
+  # create a plot from the observed data
+
+
+  PlotObsVsPred
+}
+

README.rmd

@@ -1,6 +1,6 @@
 ---
 title: "README"
-author: "Manuel Rausch, >Sascha Meyen, and Sebastian Hellmann"
+author: "Manuel Rausch, Sascha Meyen, and Sebastian Hellmann"
 date: "2024-10-13"
 output:
   md_document:
@@ -190,10 +190,10 @@ $$RMI = \frac{meta-I}{\max_{\text{accuracy}}\{meta-I\}}$$
 
 As these measures are prone to estimation bias, the package offers a simple
 bias reduction mechanism in which the observed frequencies of
-stimulus-response combinations are taken as the underyling probability
+stimulus-response combinations are taken as the underlying probability
 distribution. From this, Monte-Carlo simulations are conducted to estimate
-and subtract the bias in these measures. Note that there provably is no way
-to completely remove this bias.
+and subtract the bias from these measures. Note that there is probably no way
+to remove this bias completely.
 
 ## Installation
 
@@ -324,6 +324,7 @@ the data (notably, this is not the case though in the demo data set), the functi
 
 ```{r}        
 MetaDs <- fitMetaDprime(data = MaskOri, model="ML", .parallel = TRUE)
+head(MetaDs)
 ```
 
 Information theoretic measures of metacognition can be obtained by the function `estimateMetaI`. It expects the same kind of data.frame as `fitMetaDprime` 
@@ -332,7 +333,7 @@ The preferred way to estimate these measures is with bias reduction, but this ma
 
 ```{r}        
 metaIMeasures <- estimateMetaI(data = MaskOri, bias_reduction = TRUE)
-metaIMeasures 
+head(metaIMeasures)
 ```
 
 For the impatient and for testing purposes, bias_reduction can be turned off to increase computation speed: 

TestScript.R

@@ -356,18 +356,28 @@ Plot_recov_metaDprime_F
 
 # 4. meta-I and co
 
-MetaInfoMeasures  <-
-  MaskOri %>%
-  mutate(y = fct_recode(stimulus, "-1" = "0", "1" = "90")) %>%
-  mutate(r = factor(ifelse(response=="0", -1, 1) * as.numeric(rating))) %>%
-  group_by(participant) %>% #group_modify(~ estimate_meta_I(.x))
-  dplyr::summarise(meta_I = estimate_meta_I(pick(y,r)),
-                   meta_Ir1 = estimate_meta_Ir1(pick(y,r)),
-                   meta_Ir1_acc = estimate_meta_Ir1_acc(pick(y,r)),
-                   meta_Ir2  = estimate_meta_Ir2(pick(y,r)),
-                   RMI = estimate_RMI(pick(y,r))
-                   )  # need dplyr::summarise because plyr::summarise is used somewhere else
+MetaDs <- fitMetaDprime(data = MaskOri, model="ML", .parallel = TRUE)
+MetaInfoMeasures  <- estimateMetaI(data = MaskOri, bias_reduction = F)
+
+
+merge(MetaDs %>% select(participant, Ratio),
+      MetaInfoMeasures %>%  select(participant, meta_Ir1)) %>%
+  ggplot(aes(x=Ratio, y=meta_Ir1)) +
+  geom_point() + geom_smooth(method="lm", se=F)+
+  theme_minimal()
+
+
+# 5) Plotting fits
 
+PlotFitSDT <- plotConfModelFit(MaskOri, fitted_pars, model="SDT")
+PlotFitGN <- plotConfModelFit(MaskOri, fitted_pars, model="GN")
+PlotFitLogN <- plotConfModelFit(MaskOri, fitted_pars, model="logN")
+PlotFitWEV <- plotConfModelFit(MaskOri, fitted_pars, model="WEV")
+PlotFitLogWEV <- plotConfModelFit(MaskOri, fitted_pars, model="logWEV")
+PlotFitITGcm <- plotConfModelFit(MaskOri, fitted_pars, model="ITGcm")
+PlotFitITGc <- plotConfModelFit(MaskOri, fitted_pars, model="ITGc")
+PlotFitIG <- plotConfModelFit(MaskOri, fitted_pars, model="IG")
+PlotFitPDA <- plotConfModelFit(MaskOri, fitted_pars, model="PDA")
 
 
 save(fitted_pars, PlotFitsBICWeights,
@@ -385,3 +395,4 @@ save(fitted_pars, PlotFitsBICWeights,
      recov_metaDprime_F, Plot_recov_metaDprime_F,
      MetaInfoMeasures,
      file = "TestResults.RData")
+PlotFitSDT <- plotConfModelFit(MaskOri, fitatted_pars, model="SDT")