-
Notifications
You must be signed in to change notification settings - Fork 0
/
predictRTConf.R
167 lines (160 loc) · 8.65 KB
/
predictRTConf.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
#' Prediction of confidence rating and response time distribution for sequential sampling confidence models
#'
#' \code{predictConf} predicts the categorical response distribution of
#' decision and confidence ratings, \code{predictRT} computes the predicted
#' RT distribution (density) for the sequential sampling confidence model
#' specified by the argument \code{model}, given specific parameter constellations.
#' This function calls the respective functions for diffusion based
#' models (dynWEV and 2DSD: \code{\link{predictWEV}}) and race models (IRM, PCRM,
#' IRMt, and PCRMt: \code{\link{predictRM}}).
#'
#' @param paramDf a list or dataframe with one row. Column names should match the
#' names of the respective model parameters. For different stimulus
#' quality/mean drift rates, names should be `v1`, `v2`, `v3`,.... Different `s` parameters
#' are possible with `s1`, `s2`, `s3`... with equally many steps as for drift rates (same
#' for `sv` parameter in dynWEV and 2DSD).
#' Additionally, the confidence thresholds should be given by names with
#' `thetaUpper1`, `thetaUpper2`,..., `thetaLower1`,... or,
#' for symmetric thresholds only by `theta1`, `theta2`,....
#' @param model character scalar. One of "2DSD", "dynWEV", "IRM", "PCRM", "IRMt", or "PCRMt".
#' @param maxrt numeric. The maximum RT for the
#' integration/density computation. Default: 15 (for \code{predictConf} (integration)),
#' 9 (for \code{predictRT}).
#' @param subdivisions \code{integer} (default: 100).
#' For \code{predictConf} it is used as argument for the inner integral routine.
#' For \code{predictRT} it is the number of points for which the density is computed.
#' @param minrt numeric or NULL(default). The minimum rt for the density computation.
#' @param simult_conf logical, only relevant for dynWEV and 2DSD. Whether in the experiment
#' confidence was reported simultaneously with the decision, as then decision and confidence
#' judgment are assumed to have happened subsequent before response and computations are
#' different, when there is an observable interjudgment time (then `simult_conf` should be FALSE).
#' @param scaled logical. For \code{predictRT}. Whether the computed density
#' should be scaled to integrate to one (additional column `densscaled`). Otherwise the output
#' contains only the defective density (i.e. its integral is equal to the probability of a
#' response and not 1). If `TRUE`, the argument `DistConf` should be given, if available.
#' Default: `FALSE`.
#' @param DistConf `NULL` or `data.frame`. A `data.frame` or `matrix` with column
#' names, giving the distribution of response and rating choices for
#' different conditions and stimulus categories in the form of the output of
#' \code{predictConf}. It is only necessary, if `scaled=TRUE`, because these
#' probabilities are used for scaling. If `scaled=TRUE` and `DistConf=NULL`, it will be computed
#' with the function \code{predictRM_Conf}, which takes some time and the function will
#' throw a message. Default: `NULL`
#' @param stop.on.error logical. Argument directly passed on to integrate. Default is FALSE,
#' since the densities invoked may lead to slow convergence of the integrals (which are still
#' quite accurate) which causes R to throw an error.
#' @param .progress logical. If TRUE (default) a progress bar is drawn to the console.
#'
#' @return \code{predictConf} returns a `data.frame`/`tibble` with columns: `condition`, `stimulus`,
#' `response`, `rating`, `correct`, `p`, `info`, `err`. `p` is the predicted probability of a response
#' and `rating`, given the stimulus category and condition. `info` and `err` refer to the
#' respective outputs of the integration routine used for the computation.
#' \code{predictRT} returns a `data.frame`/`tibble` with columns: `condition`, `stimulus`,
#' `response`, `rating`, `correct`, `rt` and `dens` (and `densscaled`, if `scaled=TRUE`).
#'
#'
#' @details The function \code{predictConf} consists merely of an integration of
#' the reaction time density of the given model, \code{{d*model*}}, over the response
#' time in a reasonable interval (0 to `maxrt`). The function \code{predictRT} wraps
#' these density functions to a parameter set input and a data.frame output.
#' For the argument \code{paramDf}, the output of the fitting function \code{\link{fitRTConf}}
#' with the respective model may be used.
#'
#' @note Different parameters for different conditions are only allowed for drift rate,
#' \code{v}, drift rate variability, \code{sv} (in dynWEV and 2DSD), and process variability
#' `s`. All other parameters are used for all conditions.
#'
#' @author Sebastian Hellmann.
#'
#' @name predictRTConf
#' @importFrom stats integrate
#' @importFrom progress progress_bar
# @importFrom pracma integral
#' @aliases predictConf
#'
#' @examples
#' # Examples for "dynWEV" model (equivalent applicable for
#' # all other models (with different parameters!))
#'
#' # 1. Define some parameter set in a data.frame
#' paramDf <- data.frame(a=1.5,v1=0.2, v2=1, t0=0.1,z=0.52,
#' sz=0.3,sv=0.4, st0=0, tau=3, w=0.5,
#' theta1=1, svis=0.5, sigvis=0.8)
#'
#' # 2. Predict discrete Choice x Confidence distribution:
#' preds_Conf <- predictConf(paramDf, "dynWEV", maxrt = 25, simult_conf=TRUE)
#' head(preds_Conf)
#'
#'
#' # 3. Compute RT density
#' preds_RT <- predictRT(paramDf, "dynWEV") #(scaled=FALSE)
#' # same output with default rt-grid and without scaled density column:
#' preds_RT <- predictRT(paramDf, "dynWEV", maxrt=5, subdivisions=200,
#' minrt=paramDf$tau+paramDf$t0, simult_conf = TRUE,
#' scaled=TRUE, DistConf = preds_Conf)
#' head(preds_RT)
#' \donttest{
#' # produces a warning, if scaled=TRUE and DistConf missing
#' preds_RT <- predictRT(paramDf, "dynWEV",
#' scaled=TRUE)
#' }
#'
#' \donttest{
#' # Example of visualization
#' library(ggplot2)
#' preds_Conf$rating <- factor(preds_Conf$rating, labels=c("unsure", "sure"))
#' preds_RT$rating <- factor(preds_RT$rating, labels=c("unsure", "sure"))
#' ggplot(preds_Conf, aes(x=interaction(rating, response), y=p))+
#' geom_bar(stat="identity")+
#' facet_grid(cols=vars(stimulus), rows=vars(condition), labeller = "label_both")
#' ggplot(preds_RT, aes(x=rt, color=interaction(rating, response), y=densscaled))+
#' geom_line(stat="identity")+
#' facet_grid(cols=vars(stimulus), rows=vars(condition), labeller = "label_both")+
#' theme(legend.position = "bottom")+ ggtitle("Scaled Densities")
#' ggplot(aggregate(dens~rt+correct+rating+condition, preds_RT, mean),
#' aes(x=rt, color=rating, y=dens))+
#' geom_line(stat="identity")+
#' facet_grid(cols=vars(condition), rows=vars(correct), labeller = "label_both")+
#' theme(legend.position = "bottom")+ ggtitle("Non-Scaled Densities")
#' }
#' # Use PDFtoQuantiles to get predicted RT quantiles
#' head(PDFtoQuantiles(preds_RT, scaled = FALSE))
#'
#' @rdname predictRTConf
#' @export
predictConf <- function(paramDf, model=NULL,
maxrt=15, subdivisions = 100L, simult_conf = FALSE, stop.on.error=FALSE,
.progress=TRUE){
paramDf <- as.data.frame(paramDf)
paramDf <- paramDf[,!is.na(paramDf)]
if (is.null(model) && ("model" %in% names(paramDf))) model <- paramDf$model
#### Check model argument
if (grepl("RM", model)) {
res <- predictRM_Conf(paramDf, model, FALSE, maxrt, subdivisions, stop.on.error, .progress)
} else if (grepl("dynWEV|2DSD|dynaViTE", model)) {
res <- predictWEV_Conf(paramDf, model, maxrt, subdivisions,simult_conf, stop.on.error, .progress=.progress)
} else if (grepl("DDMConf", model)) {
res <- predictDDMConf_Conf(paramDf, maxrt, subdivisions, stop.on.error, .progress)
} else { stop("model not known.")}
return(res)
}
### Predict RT-distribution
#' @rdname predictRTConf
#' @export
predictRT <- function(paramDf, model=NULL,
maxrt=9, subdivisions = 100L, minrt=NULL, simult_conf=FALSE,
scaled = FALSE, DistConf=NULL,
.progress = TRUE) {
paramDf <- as.data.frame(paramDf)
paramDf <- paramDf[,!is.na(paramDf)]
if (is.null(model) && ("model" %in% names(paramDf))) model <- paramDf$model
#### Check model argument
if (grepl("RM", model)) {
res <- predictRM_RT(paramDf, model, FALSE, maxrt, subdivisions, minrt, scaled, DistConf, .progress)
} else if (grepl("dynWEV|2DSD|dynaViTE", model)) {
res <- predictWEV_RT(paramDf, model, maxrt, subdivisions, minrt, simult_conf, scaled, DistConf, .progress=.progress)
} else if (grepl("DDMConf", model)) {
res <- predictDDMConf_RT(paramDf, maxrt, subdivisions, minrt, scaled, DistConf, .progress)
} else { stop("model not known.")}
return(res)
}