diff --git a/.Rbuildignore b/.Rbuildignore index b5616ca..6ceb88c 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -3,3 +3,5 @@ ^\.travis\.yml$ .git cran-comments.md +^doc$ +^Meta$ diff --git a/.gitignore b/.gitignore index fc96c9a..06bacdc 100644 --- a/.gitignore +++ b/.gitignore @@ -41,4 +41,6 @@ vignettes/*.pdf *.o *.so *.cc -*.hpp \ No newline at end of file +*.hpp +doc +Meta diff --git a/DESCRIPTION b/DESCRIPTION index e0f772f..b72608b 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -37,7 +37,9 @@ LinkingTo: Rcpp (>= 0.12.16), RcppEigen (>= 0.3.3.4.0) SystemRequirements: GNU make +VignetteBuilder: knitr NeedsCompilation: yes RoxygenNote: 6.1.1 Suggests: - testthat + testthat, + knitr diff --git a/R/b_color.R b/R/b_color.R index c53daf4..bf566c9 100644 --- a/R/b_color.R +++ b/R/b_color.R @@ -92,6 +92,11 @@ b_color <- function(colors, control=NULL, suppress_warnings=TRUE) { + # multi core + if (chains > 1) { + options(mc.cores = parallel::detectCores()) + } + # prepare data n <- nrow(colors) diff --git a/R/b_linear.R b/R/b_linear.R index 39daeeb..b1496e8 100644 --- a/R/b_linear.R +++ b/R/b_linear.R @@ -55,6 +55,11 @@ b_linear <- function(x, control=NULL, suppress_warnings=TRUE) { + # multi core + if (chains > 1) { + options(mc.cores = parallel::detectCores()) + } + # prepare data n <- length(y) m <- length(unique(s)) diff --git a/R/b_reaction_time.R b/R/b_reaction_time.R index 5e9571c..33ad2f2 100644 --- a/R/b_reaction_time.R +++ b/R/b_reaction_time.R @@ -48,6 +48,11 @@ b_reaction_time <- function(t, control=NULL, suppress_warnings=TRUE) { + # multi core + if (chains > 1) { + options(mc.cores = parallel::detectCores()) + } + # prepare data n <- length(t) m <- length(unique(s)) diff --git a/R/b_success_rate.R b/R/b_success_rate.R index 780e8d4..764bd95 100644 --- a/R/b_success_rate.R +++ b/R/b_success_rate.R @@ -43,6 +43,11 @@ b_success_rate <- function(r, control=NULL, suppress_warnings=TRUE) { + # multi core + if (chains > 1) { + options(mc.cores = parallel::detectCores()) + } + # prepare data n <- length(r) m <- length(unique(s)) diff --git a/R/b_ttest.R b/R/b_ttest.R index 9de10e8..5870eda 100644 --- a/R/b_ttest.R +++ b/R/b_ttest.R @@ -42,6 +42,11 @@ b_ttest <- function(data, control=NULL, suppress_warnings=TRUE) { + # multi core + if (chains > 1) { + options(mc.cores = parallel::detectCores()) + } + # prepare data n <- length(data) diff --git a/R/color_class.R b/R/color_class.R index 7389180..93fffb0 100644 --- a/R/color_class.R +++ b/R/color_class.R @@ -58,6 +58,8 @@ #' #' plot_distributions_difference(`color_class`, hsv=`vector`): a visualization of the difference between the distribution of the first fit and a color defined with hsv components. You can also provide the rope and bins (number of bins in the histogram) parameters or visualize the comparison only through chosen color components (r, g, b, h, s, v) by using the pars parameter. #' +#' plot_hsv(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit thorough a circular visualization of hsv color components. +#' #' plot_fit_hsv(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit thorough a circular visualization of hsv color components. #' #' plot_means_hsv(`color_class`): a visualization of the difference between means of two fits through a circular visualization of hsv color components. You can also compare fit means with colors defined through rgb or hsv components (as points or as lines on the visualization). @@ -244,6 +246,9 @@ #' plot_distributions_difference(fit1, hsv=c(pi/2, 1, 1)) #' #' # plot the fitted distribution for hue against the hue data +#' plot_hsv(fit1) +#' +#' # plot the fitted distribution for hue against the hue data #' plot_fit_hsv(fit1) #' #' # visualize hue means of a single fit @@ -1989,6 +1994,28 @@ setMethod(f="plot_distributions_difference", signature(object="color_class"), de }) +#' @rdname color_class-plot_hsv +#' @exportMethod plot_hsv +setGeneric(name="plot_hsv", function(object) standardGeneric("plot_hsv")) + +#' @title plot_hsv +#' @description \code{plot_hsv} plots fitted model against the data. Use this function to explore the quality of your fit thorough a circular visualization of hsv color components. +#' @param object color_class object. +#' @rdname color_class-plot_hsv +#' @aliases plot_hsv_color +#' @return A ggplot visualization. +#' +#' @examples +#' # to use the function you first have to prepare the data and fit the model +#' # see class documentation for an example of the whole process +#' # along with an example of how to use this function +#' ?color_class +#' +setMethod(f="plot_hsv", signature(object="color_class"), definition=function(object) { + return(plot_fit_hsv(object)) +}) + + #' @rdname color_class-plot_fit_hsv #' @exportMethod plot_fit_hsv setGeneric(name="plot_fit_hsv", function(object) standardGeneric("plot_fit_hsv")) diff --git a/R/data.R b/R/data.R index 8589b91..01d9558 100644 --- a/R/data.R +++ b/R/data.R @@ -1,20 +1,147 @@ #' Datasets for bayes4psy examples -#' Small datasets for use in \pkg{rstanarm} examples and vignettes. +#' Example datasets for use in \pkg{rstanarm} examples and vignettes. #' -#' @name rstanarm-datasets +#' @name bayes4psy-datasets #' @aliases adaptation_level_small +#' #' @format #' \describe{ #' \item{\code{adaptation_level_small}}{ +#' Small dataset on subjects picking up weights and determining their weights from 1..10. +#' +#' 50 obs. of 3 variables +#' \itemize{ +#' \item \code{sequence} sequence index. +#' \item \code{weight} actual weight of the object. +#' \item \code{response} subject's estimation of weight. +#' } +#' } +#' \item{\code{adaptation_level}}{ #' Data on subjects picking up weights and determining their weights from 1..10. #' -#' 50 obs. of 5 variables +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 2900 obs. of 6 variables #' \itemize{ +#' \item \code{subject} subject index. +#' \item \code{group} group index. +#' \item \code{part} first or second part of the experiment. #' \item \code{sequence} sequence index. #' \item \code{weight} actual weight of the object. #' \item \code{response} subject's estimation of weight. #' } #' } +#' #' \item{\code{after_images_opponent_process}}{ +#' Colors predicted by the opponent process theory. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 6 obs. of 7 variables +#' \itemize{ +#' \item \code{stimuli} name of the color stimuli. +#' \item \code{r} value of the R component in the RGB model. +#' \item \code{g} value of the G component in the RGB model. +#' \item \code{b} value of the B component in the RGB model. +#' \item \code{h} value of the H component in the HSV model. +#' \item \code{s} value of the S component in the HSV model. +#' \item \code{v} value of the V component in the HSV model. +#' } +#' } +#' #' \item{\code{after_images_opponent_stimuli}}{ +#' Stimuli used in the after images experiment. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 6 obs. of 7 variables +#' \itemize{ +#' \item \code{r_s} value of the R component in the RGB model. +#' \item \code{g_s} value of the G component in the RGB model. +#' \item \code{b_s} value of the B component in the RGB model. +#' \item \code{stimuli} name of the color stimuli. +#' \item \code{h_s} value of the H component in the HSV model. +#' \item \code{s_s} value of the S component in the HSV model. +#' \item \code{v_s} value of the V component in the HSV model. +#' } +#' } +#' #' \item{\code{after_images_trichromatic}}{ +#' Colors predicted by the trichromatic theory. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 6 obs. of 7 variables +#' \itemize{ +#' \item \code{stimuli} name of the color stimuli. +#' \item \code{r} value of the R component in the RGB model. +#' \item \code{g} value of the G component in the RGB model. +#' \item \code{b} value of the B component in the RGB model. +#' \item \code{h} value of the H component in the HSV model. +#' \item \code{s} value of the S component in the HSV model. +#' \item \code{v} value of the V component in the HSV model. +#' } +#' } +#' #' \item{\code{after_images}}{ +#' Data gathered by the after images experiment. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 1311 obs. of 12 variables +#' \itemize{ +#' \item \code{subject} subject index. +#' \item \code{rt} reaction time. +#' \item \code{r} value of the R component in the RGB model of subject's response. +#' \item \code{g} value of the G component in the RGB model of subject's response. +#' \item \code{b} value of the B component in the RGB model of subject's response. +#' \item \code{stimuli} name of the color stimuli. +#' \item \code{r_s} value of the R component in the RGB model of the shown stimulus +#' \item \code{g_s} value of the G component in the RGB model of the shown stimulus +#' \item \code{b_s} value of the B component in the RGB model of the shown stimulus +#' \item \code{h_s} value of the H component in the HSV model of the shown stimulus +#' \item \code{s_s} value of the S component in the HSV model of the shown stimulus +#' \item \code{v_s} value of the V component in the HSV model of the shown stimulus +#' } +#' } +#' #' \item{\code{flanker}}{ +#' Data gathered by the flanker experiment. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 8256 obs. of 5 variables +#' \itemize{ +#' \item \code{subject} subject index. +#' \item \code{group} group index. +#' \item \code{congruencty} type of stimulus. +#' \item \code{result} was subject's reponse correct or wrong? +#' \item \code{rt} reaction time. +#' } +#' } +#' #' \item{\code{stroop_extended}}{ +#' All the data gathered by the Stroop experiment. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 41068 obs. of 5 variables +#' \itemize{ +#' \item \code{subject} subject ID. +#' \item \code{cond} type of condition. +#' \item \code{rt} reaction time. +#' \item \code{acc} was subject's reponse correct or wrong? +#' \item \code{age} age of subject. +#' } +#' } +#' #' \item{\code{stroop_simple}}{ +#' All the data gathered by the Stroop experiment. +#' +#' Source: Internal MBLab \url{www.mblab.si} repository. +#' +#' 61 obs. of 5 variables +#' \itemize{ +#' \item \code{subject} subject ID. +#' \item \code{reading_neutral} average response time for reading neutral stimuli. +#' \item \code{naming_neutral} average response time for naming neutral stimuli. +#' \item \code{reading_incongruent} average response time for reading incongruent stimuli. +#' \item \code{naming_incongruent} average response time for naming incongruent stimuli. +#' } +#' } #' } #' #' @examples diff --git a/bayes4psy.Rproj b/bayes4psy.Rproj index 744339d..2de124d 100644 --- a/bayes4psy.Rproj +++ b/bayes4psy.Rproj @@ -19,4 +19,4 @@ BuildType: Package PackageUseDevtools: Yes PackageInstallArgs: --no-multiarch --with-keep.source PackageCheckArgs: --no-multiarch -PackageRoxygenize: rd,collate,namespace +PackageRoxygenize: rd,collate,namespace,vignette diff --git a/data/adaptation_level.rda b/data/adaptation_level.rda new file mode 100644 index 0000000..72b63b0 Binary files /dev/null and b/data/adaptation_level.rda differ diff --git a/data/after_images.rda b/data/after_images.rda new file mode 100644 index 0000000..d9f344e Binary files /dev/null and b/data/after_images.rda differ diff --git a/data/after_images_opponent_process.rda b/data/after_images_opponent_process.rda new file mode 100644 index 0000000..f1aaf28 Binary files /dev/null and b/data/after_images_opponent_process.rda differ diff --git a/data/after_images_stimuli.rda b/data/after_images_stimuli.rda new file mode 100644 index 0000000..f7afe9f Binary files /dev/null and b/data/after_images_stimuli.rda differ diff --git a/data/after_images_trichromatic.rda b/data/after_images_trichromatic.rda new file mode 100644 index 0000000..4d7476b Binary files /dev/null and b/data/after_images_trichromatic.rda differ diff --git a/data/flanker.rda b/data/flanker.rda new file mode 100644 index 0000000..ccc2265 Binary files /dev/null and b/data/flanker.rda differ diff --git a/data/stroop_extended.rda b/data/stroop_extended.rda new file mode 100644 index 0000000..606c1b4 Binary files /dev/null and b/data/stroop_extended.rda differ diff --git a/data/stroop_simple.rda b/data/stroop_simple.rda new file mode 100644 index 0000000..75e188d Binary files /dev/null and b/data/stroop_simple.rda differ diff --git a/man/bayes4psy-datasets.Rd b/man/bayes4psy-datasets.Rd new file mode 100644 index 0000000..3194f8f --- /dev/null +++ b/man/bayes4psy-datasets.Rd @@ -0,0 +1,164 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\name{bayes4psy-datasets} +\alias{bayes4psy-datasets} +\alias{adaptation_level_small} +\title{Datasets for bayes4psy examples +Example datasets for use in \pkg{rstanarm} examples and vignettes.} +\format{\describe{ +\item{\code{adaptation_level_small}}{ +Small dataset on subjects picking up weights and determining their weights from 1..10. + +50 obs. of 3 variables +\itemize{ +\item \code{sequence} sequence index. +\item \code{weight} actual weight of the object. +\item \code{response} subject's estimation of weight. +} +} +\item{\code{adaptation_level}}{ +Data on subjects picking up weights and determining their weights from 1..10. + +Source: Internal MBLab \url{www.mblab.si} repository. + +2900 obs. of 6 variables +\itemize{ +\item \code{subject} subject index. +\item \code{group} group index. +\item \code{part} first or second part of the experiment. +\item \code{sequence} sequence index. +\item \code{weight} actual weight of the object. +\item \code{response} subject's estimation of weight. +} +} +#' \item{\code{after_images_opponent_process}}{ +Colors predicted by the opponent process theory. + +Source: Internal MBLab \url{www.mblab.si} repository. + +6 obs. of 7 variables +\itemize{ +\item \code{stimuli} name of the color stimuli. +\item \code{r} value of the R component in the RGB model. +\item \code{g} value of the G component in the RGB model. +\item \code{b} value of the B component in the RGB model. +\item \code{h} value of the H component in the HSV model. +\item \code{s} value of the S component in the HSV model. +\item \code{v} value of the V component in the HSV model. +} +} +#' \item{\code{after_images_opponent_stimuli}}{ +Stimuli used in the after images experiment. + +Source: Internal MBLab \url{www.mblab.si} repository. + +6 obs. of 7 variables +\itemize{ +\item \code{r_s} value of the R component in the RGB model. +\item \code{g_s} value of the G component in the RGB model. +\item \code{b_s} value of the B component in the RGB model. +\item \code{stimuli} name of the color stimuli. +\item \code{h_s} value of the H component in the HSV model. +\item \code{s_s} value of the S component in the HSV model. +\item \code{v_s} value of the V component in the HSV model. +} +} +#' \item{\code{after_images_trichromatic}}{ +Colors predicted by the trichromatic theory. + +Source: Internal MBLab \url{www.mblab.si} repository. + +6 obs. of 7 variables +\itemize{ +\item \code{stimuli} name of the color stimuli. +\item \code{r} value of the R component in the RGB model. +\item \code{g} value of the G component in the RGB model. +\item \code{b} value of the B component in the RGB model. +\item \code{h} value of the H component in the HSV model. +\item \code{s} value of the S component in the HSV model. +\item \code{v} value of the V component in the HSV model. +} +} +#' \item{\code{after_images}}{ +Data gathered by the after images experiment. + +Source: Internal MBLab \url{www.mblab.si} repository. + +1311 obs. of 12 variables +\itemize{ +\item \code{subject} subject index. +\item \code{rt} reaction time. +\item \code{r} value of the R component in the RGB model of subject's response. +\item \code{g} value of the G component in the RGB model of subject's response. +\item \code{b} value of the B component in the RGB model of subject's response. +\item \code{stimuli} name of the color stimuli. +\item \code{r_s} value of the R component in the RGB model of the shown stimulus +\item \code{g_s} value of the G component in the RGB model of the shown stimulus +\item \code{b_s} value of the B component in the RGB model of the shown stimulus +\item \code{h_s} value of the H component in the HSV model of the shown stimulus +\item \code{s_s} value of the S component in the HSV model of the shown stimulus +\item \code{v_s} value of the V component in the HSV model of the shown stimulus +} +} +#' \item{\code{flanker}}{ +Data gathered by the flanker experiment. + +Source: Internal MBLab \url{www.mblab.si} repository. + +8256 obs. of 5 variables +\itemize{ +\item \code{subject} subject index. +\item \code{group} group index. +\item \code{congruencty} type of stimulus. +\item \code{result} was subject's reponse correct or wrong? +\item \code{rt} reaction time. +} +} +#' \item{\code{stroop_extended}}{ +All the data gathered by the Stroop experiment. + +Source: Internal MBLab \url{www.mblab.si} repository. + +41068 obs. of 5 variables +\itemize{ +\item \code{subject} subject ID. +\item \code{cond} type of condition. +\item \code{rt} reaction time. +\item \code{acc} was subject's reponse correct or wrong? +\item \code{age} age of subject. +} +} +#' \item{\code{stroop_simple}}{ +All the data gathered by the Stroop experiment. + +Source: Internal MBLab \url{www.mblab.si} repository. + +61 obs. of 5 variables +\itemize{ +\item \code{subject} subject ID. +\item \code{reading_neutral} average response time for reading neutral stimuli. +\item \code{naming_neutral} average response time for naming neutral stimuli. +\item \code{reading_incongruent} average response time for reading incongruent stimuli. +\item \code{naming_incongruent} average response time for naming incongruent stimuli. +} +} +}} +\description{ +Datasets for bayes4psy examples +Example datasets for use in \pkg{rstanarm} examples and vignettes. +} +\examples{ + +# Example of Bayesian bootstraping on 'adaptation_level_small' dataset +# linear function of seqence vs. response +lm_statistic <- function(data) { + lm(sequence ~ response, data)$coef +} + +# load data +data <- adaptation_level_small + +# bootstrap +data_bootstrap <- b_bootstrap(data, lm_statistic, n1=1000, n2=1000) + +} diff --git a/man/color_class-class.Rd b/man/color_class-class.Rd index 823d3f7..e3ded3f 100644 --- a/man/color_class-class.Rd +++ b/man/color_class-class.Rd @@ -16,6 +16,12 @@ print(`color_class`): prints a more detailed summary of the fit show(`color_class`): prints a more detailed summary of the fit. +plot(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit. You can compare fit with underlying data only through chosen color components (r, g, b, h, s, v) by using the pars parameter. + +plot_fit(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit. You can compare fit with underlying data only through chosen color components (r, g, b, h, s, v) by using the pars parameter. + +plot_trace(`color_class`): traceplot for main fitted model parameters. + get_parameters(`color_class`): returns a dataframe with values of fitted parameters. compare_means(`color_class`, fit2=`color_class`): prints color difference between two fits. You can also provide the rope parameter or execute the comparison only through chosen color components (r, g, b, h, s, v) by using the pars parameter. @@ -58,15 +64,11 @@ plot_distributions_difference(`color_class`, rgb=`vector`): a visualization of t plot_distributions_difference(`color_class`, hsv=`vector`): a visualization of the difference between the distribution of the first fit and a color defined with hsv components. You can also provide the rope and bins (number of bins in the histogram) parameters or visualize the comparison only through chosen color components (r, g, b, h, s, v) by using the pars parameter. -plot_fit(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit. You can compare fit with underlying data only through chosen color components (r, g, b, h, s, v) by using the pars parameter. - plot_fit_hsv(`color_class`): plots fitted model against the data. Use this function to explore the quality of your fit thorough a circular visualization of hsv color components. plot_means_hsv(`color_class`): a visualization of the difference between means of two fits through a circular visualization of hsv color components. You can also compare fit means with colors defined through rgb or hsv components (as points or as lines on the visualization). plot_distributions_hsv(`color_class`): a visualization of distributions of one or two fits thorough a circular visualization of hsv color components. You can also compare fit means with colors defined through rgb or hsv components (as points or as lines on the visualization). - -plot_trace(`color_class`): traceplot for main fitted model parameters. } \section{Slots}{ @@ -150,6 +152,18 @@ summary(fit1) print(fit1) show(fit1) +# plot the fitted distribution against the data +plot(fit1) +plot_fit(fit1) + +# plot the fitted distribution against the data, +# specify only a subset of parameters for visualization +plot(fit1, pars=c("h", "s", "v")) +plot_fit(fit1, pars=c("h", "s", "v")) + +# traceplot of the fitted parameters +plot_trace(fit1) + # extract parameter values from the fit parameters <- get_parameters(fit1) @@ -241,13 +255,6 @@ plot_distributions_difference(fit1, rgb=c(255, 0, 0)) # and an hsv defined color plot_distributions_difference(fit1, hsv=c(pi/2, 1, 1)) -# plot the fitted distribution against the data -plot_fit(fit1) - -# plot the fitted distribution against the data, -# specify only a subset of parameters for visualization -plot_fit(fit1, pars=c("h", "s", "v")) - # plot the fitted distribution for hue against the hue data plot_fit_hsv(fit1) @@ -284,9 +291,6 @@ points <- list() points[[1]] <- c(pi, 1, 1) plot_distributions_hsv(fit1, fit2=fit2, points=points, lines=lines, hsv=TRUE) - -# traceplot of the fitted parameters -plot_trace(fit1) } } diff --git a/man/linear_class-class.Rd b/man/linear_class-class.Rd index c9a92ce..ee89d0b 100644 --- a/man/linear_class-class.Rd +++ b/man/linear_class-class.Rd @@ -16,6 +16,16 @@ print(`linear_class`): prints a more detailed summary of the fit show(`linear_class`): prints a more detailed summary of the fit. +plot(`linear_class`): plots fitted model against the data. Use this function to explore the quality of your fit. Fit will be plotted on the subject level. + +plot(`linear_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the subjects level (subjects=FALSE). + +plot_fit(`linear_class`): plots fitted model against the data. Use this function to explore the quality of your fit. Fit will be plotted on the subject level. + +plot_fit(`linear_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the subjects level (subjects=FALSE). + +plot_trace(`linear_class`): traceplot for main fitted model parameters. + get_parameters(`linear_class`): returns a dataframe with values of fitted parameters. get_subject_parameters(`linear_class`): returns a dataframe with values of fitted parameters for each subject in the hierarchical model. @@ -35,12 +45,6 @@ plot_distributions(`linear_class`): a visualization of the fitted distribution. plot_distributions(`linear_class`, fit2=`linear_class`): a visualization of two fitted distribution. plot_distributions_difference(`linear_class`, fit2=`linear_class`): a visualization of the difference between the distribution of the first group and the second group. You can plot only slope or intercept by using the par parameter. You can also provide the rope and bins (number of bins in the histogram) parameters. - -plot_fit(`linear_class`): plots fitted model against the data. Use this function to explore the quality of your fit. Fit will be plotted on the subject level. - -plot_fit(`linear_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the subjects level (subjects=FALSE). - -plot_trace(`linear_class`): traceplot for main fitted model parameters. } \section{Slots}{ @@ -88,6 +92,17 @@ summary(fit1) print(fit1) show(fit1) +# plot the fitted distribution against the data +plot(fit1) +plot_fit(fit1) + +# plot the fitted distribution against the data, +# plot on the top (group) level +plot(fit1, subjects=FALSE) +plot_fit(fit1, subjects=FALSE) + +# traceplot of the fitted parameters +plot_trace(fit1) # extract parameter values from the fit parameters <- get_parameters(fit1) @@ -143,16 +158,6 @@ plot_distributions_difference(fit1, fit2=fit2, rope_intercept=0.5, rope_slope=0. # visualize difference between distributions underlying two fits, plot slope only plot_distributions_difference(fit1, fit2=fit2, par="slope") - -# plot the fitted distribution against the data -plot_fit(fit1) - -# plot the fitted distribution against the data, -# plot on the top (group) level -plot_fit(fit1, subjects=FALSE) - -# traceplot of the fitted parameters -plot_trace(fit1) } } diff --git a/man/plot-color_class-missing-method.Rd b/man/plot-color_class-missing-method.Rd new file mode 100644 index 0000000..3d5c130 --- /dev/null +++ b/man/plot-color_class-missing-method.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/color_class.R +\docType{methods} +\name{plot,color_class,missing-method} +\alias{plot,color_class,missing-method} +\title{plot} +\usage{ +\S4method{plot}{color_class,missing}(x, y, ...) +} +\arguments{ +\item{...}{pars - components of comparison, a subset of (r, g, b, h, s, v).} + +\item{object}{color_class object.} +} +\description{ +\code{plot} plots fitted model against the data. Use this function to explore the quality of your fit. You can compare fit with underlying data only through chosen color components (r, g, b, h, s, v). +} +\examples{ +# to use the function you first have to prepare the data and fit the model +# see class documentation for an example of the whole process +# along with an example of how to use this function +?color_class + +} diff --git a/man/plot-linear_class-missing-method.Rd b/man/plot-linear_class-missing-method.Rd new file mode 100644 index 0000000..b5982fc --- /dev/null +++ b/man/plot-linear_class-missing-method.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/linear_class.R +\docType{methods} +\name{plot,linear_class,missing-method} +\alias{plot,linear_class,missing-method} +\title{plot} +\usage{ +\S4method{plot}{linear_class,missing}(x, y, ...) +} +\arguments{ +\item{...}{subjects - plot fits on a subject level (default = TRUE).} + +\item{object}{linear_class object.} +} +\description{ +\code{plot} plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the group level (subjects=FALSE). +} +\examples{ +# to use the function you first have to prepare the data and fit the model +# see class documentation for an example of the whole process +# along with an example of how to use this function +?linear_class + +} diff --git a/man/plot-reaction_time_class-missing-method.Rd b/man/plot-reaction_time_class-missing-method.Rd new file mode 100644 index 0000000..d9608aa --- /dev/null +++ b/man/plot-reaction_time_class-missing-method.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/reaction_time_class.R +\docType{methods} +\name{plot,reaction_time_class,missing-method} +\alias{plot,reaction_time_class,missing-method} +\title{plot} +\usage{ +\S4method{plot}{reaction_time_class,missing}(x, y, ...) +} +\arguments{ +\item{...}{subjects - plot fits on a subject level (default = TRUE).} + +\item{object}{reaction_time_class object.} +} +\description{ +\code{plot} plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subjects level (subjects=TRUE) or on the group level (subjects=FALSE). +} +\examples{ +# to use the function you first have to prepare the data and fit the model +# see class documentation for an example of the whole process +# along with an example of how to use this function +?reaction_time_class + +} diff --git a/man/plot-success_rate_class-missing-method.Rd b/man/plot-success_rate_class-missing-method.Rd new file mode 100644 index 0000000..5de911a --- /dev/null +++ b/man/plot-success_rate_class-missing-method.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/success_rate_class.R +\docType{methods} +\name{plot,success_rate_class,missing-method} +\alias{plot,success_rate_class,missing-method} +\title{plot} +\usage{ +\S4method{plot}{success_rate_class,missing}(x, y, ...) +} +\arguments{ +\item{...}{subjects - plot fits on a subject level (default = TRUE).} + +\item{object}{success_rate_class object.} +} +\description{ +\code{plot} plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subjects level (subjects=TRUE) or on the group level (subjects=FALSE). +} +\examples{ +# to use the function you first have to prepare the data and fit the model +# see class documentation for an example of the whole process +# along with an example of how to use this function +?success_rate_class + +} diff --git a/man/plot-ttest_class-missing-method.Rd b/man/plot-ttest_class-missing-method.Rd new file mode 100644 index 0000000..e2aa238 --- /dev/null +++ b/man/plot-ttest_class-missing-method.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ttest_class.R +\docType{methods} +\name{plot,ttest_class,missing-method} +\alias{plot,ttest_class,missing-method} +\title{plot} +\usage{ +\S4method{plot}{ttest_class,missing}(x) +} +\arguments{ +\item{object}{ttest_class object.} +} +\description{ +\code{plot} plots fitted model against the data. Use this function to explore the quality of your fit. +} +\examples{ +# to use the function you first have to prepare the data and fit the model +# see class documentation for an example of the whole process +# along with an example of how to use this function +?ttest_class + +} diff --git a/man/reaction_time_class-class.Rd b/man/reaction_time_class-class.Rd index 2eafd0e..afb5f69 100644 --- a/man/reaction_time_class-class.Rd +++ b/man/reaction_time_class-class.Rd @@ -16,6 +16,16 @@ print(`reaction_time_class`): prints a more detailed summary of the fit show(`reaction_time_class`): prints a more detailed summary of the fit. +plot(`reaction_time_class`): plots fitted model against the data. Use this function to explore the quality of your fit. + +plot(`reaction_time_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the group level (subjects=FALSE). + +plot_fit(`reaction_time_class`): plots fitted model against the data. Use this function to explore the quality of your fit. + +plot_fit(`reaction_time_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the group level (subjects=FALSE). + +plot_trace(`reaction_time_class`): traceplot for main fitted model parameters. + get_parameters(`reaction_time_class`): returns a dataframe with values of fitted parameters. get_subject_parameters(`reaction_time_class`): returns a dataframe with values of fitted parameters for each subject in the hierarchical model. @@ -47,12 +57,6 @@ plot_distributions(`reaction_time_class`, fits=`list`): a visualization of the d plot_distributions_difference(`reaction_time_class`, fit2=`reaction_time_class`): a visualization of the difference between the distribution of the first group and the second group. You can also provide the rope and bins (number of bins in the histogram) parameters. plot_distributions_difference(`reaction_time_class`, fits=`list`): a visualization of the difference between the distributions of multiple groups. You can also provide the rope and bins (number of bins in the histogram) parameters. - -plot_fit(`reaction_time_class`): plots fitted model against the data. Use this function to explore the quality of your fit. - -plot_fit(`reaction_time_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subject level (subjects=TRUE) or on the group level (subjects=FALSE). - -plot_trace(`reaction_time_class`): traceplot for main fitted model parameters. } \section{Slots}{ @@ -106,6 +110,18 @@ summary(fit1) print(fit1) show(fit1) +# plot the fitted distribution against the data +plot(fit1) +plot_fit(fit1) + +# plot the fitted distribution against the data, +# plot on the top (group) level +plot(fit1, subjects=FALSE) +plot_fit(fit1, subjects=FALSE) + +# traceplot of the fitted parameters +plot_trace(fit1) + # extract parameter values from the fit parameters <- get_parameters(fit1) @@ -168,16 +184,6 @@ plot_distributions_difference(fit1, fit2=fit2, rope=0.05) # visualize difference between distributions underlying multiple fits plot_distributions_difference(fit1, fits=fit_list) - -# plot the fitted distribution against the data -plot_fit(fit1) - -# plot the fitted distribution against the data, -# plot on the top (group) level -plot_fit(fit1, subjects=FALSE) - -# traceplot of the fitted parameters -plot_trace(fit1) } } diff --git a/man/rstanarm-datasets.Rd b/man/rstanarm-datasets.Rd deleted file mode 100644 index 3adf5b0..0000000 --- a/man/rstanarm-datasets.Rd +++ /dev/null @@ -1,38 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/data.R -\name{rstanarm-datasets} -\alias{rstanarm-datasets} -\alias{adaptation_level_small} -\title{Datasets for bayes4psy examples -Small datasets for use in \pkg{rstanarm} examples and vignettes.} -\format{\describe{ -\item{\code{adaptation_level_small}}{ -Data on subjects picking up weights and determining their weights from 1..10. - -50 obs. of 5 variables -\itemize{ -\item \code{sequence} sequence index. -\item \code{weight} actual weight of the object. -\item \code{response} subject's estimation of weight. -} -} -}} -\description{ -Datasets for bayes4psy examples -Small datasets for use in \pkg{rstanarm} examples and vignettes. -} -\examples{ - -# Example of Bayesian bootstraping on 'adaptation_level_small' dataset -# linear function of seqence vs. response -lm_statistic <- function(data) { - lm(sequence ~ response, data)$coef -} - -# load data -data <- adaptation_level_small - -# bootstrap -data_bootstrap <- b_bootstrap(data, lm_statistic, n1=1000, n2=1000) - -} diff --git a/man/success_rate_class-class.Rd b/man/success_rate_class-class.Rd index b0601ca..4f3e50b 100644 --- a/man/success_rate_class-class.Rd +++ b/man/success_rate_class-class.Rd @@ -16,8 +16,12 @@ print(`success_rate_class`): prints a more detailed summary of the fit show(`success_rate_class`): prints a more detailed summary of the fit. +plot(`success_rate_class`): plots fitted model against the data. Use this function to explore the quality of your fit. + plot(`success_rate_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subjects level (subjects=TRUE) or on the group level (subjects=FALSE). +plot_fit(`success_rate_class`): plots fitted model against the data. Use this function to explore the quality of your fit. + plot_fit(`success_rate_class`, subjects='boolean'): plots fitted model against the data. Use this function to explore the quality of your fit. You can plot on the subjects level (subjects=TRUE) or on the group level (subjects=FALSE). plot_trace(`success_rate_class`): traceplot for main fitted model parameters. @@ -102,6 +106,18 @@ summary(fit1) print(fit1) show(fit1) +# plot the fitted distribution against the data +plot(fit1) +plot_fit(fit1) + +# plot the fitted distribution against the data, +# plot on the top (group) level +plot(fit1, subjects=FALSE) +plot_fit(fit1, subjects=FALSE) + +# traceplot of the fitted parameters +plot_trace(fit1) + # extract parameter values from the fit parameters <- get_parameters(fit1) @@ -152,16 +168,6 @@ plot_distributions_difference(fit1, fit2=fit2, rope=0.05) # visualize difference between distributions underlying multiple fits plot_distributions_difference(fit1, fits=fit_list) - -# plot the fitted distribution against the data -plot_fit(fit1) - -# plot the fitted distribution against the data, -# plot on the top (group) level -plot_fit(fit1, subjects=FALSE) - -# traceplot of the fitted parameters -plot_trace(fit1) } } diff --git a/man/ttest_class-class.Rd b/man/ttest_class-class.Rd index 5a9b89d..080eb99 100644 --- a/man/ttest_class-class.Rd +++ b/man/ttest_class-class.Rd @@ -117,6 +117,13 @@ summary(fit1) print(fit1) show(fit1) +# plot the fitted distribution against the data +plot(fit1) +plot_fit(fit1) + +# traceplot of the fitted parameters +plot_trace(fit1) + # extract parameter values from the fit parameters <- get_parameters(fit1) @@ -201,12 +208,6 @@ plot_distributions_difference(fit1, mu=150, sigma=20) # visualize difference between distributions underlying multiple fits plot_distributions_difference(fit1, fits=fit_list) - -# plot the fitted distribution against the data -plot_fit(fit1) - -# traceplot of the fitted parameters -plot_trace(fit1) } } diff --git a/tests/testthat/test_bootstrap.R b/tests/testthat/test_bootstrap.R new file mode 100644 index 0000000..9304395 --- /dev/null +++ b/tests/testthat/test_bootstrap.R @@ -0,0 +1,25 @@ +library(bayes4psy) + +# set seed +seed <- 0 +set.seed(0) + +# set tolerance +tol <- 0.2 + +# linear function of seqence vs. response +lm_statistic <- function(data) { + lm(sequence ~ response, data)$coef +} + +# load data +data <- adaptation_level_small + +# bootstrap +data_bootstrap <- b_bootstrap(data, lm_statistic, n1=1000, n2=1000) + +# get_parameters +test_that("bootstrap", { + expect_equal(mean(data_bootstrap$response), -0.479, tolerance=tol) + expect_equal(mean(data_bootstrap$`(Intercept)`), 8.754, tolerance=tol) +}) diff --git a/vignettes/adaptation_level.Rmd b/vignettes/adaptation_level.Rmd new file mode 100644 index 0000000..ac7afdd --- /dev/null +++ b/vignettes/adaptation_level.Rmd @@ -0,0 +1,18 @@ +--- + title: "Analysis of the adaptation level experiment data using the Bayesian linear model" + author: "Jure Demšar, Erik Štrumbelj and Grega Repovš" + date: "`r Sys.Date()`" + output: + html_vignette: + toc: yes +--- + + + +# Introduction + +This vignette explains how to use the Bayesian linear model in the bayes4psy package. diff --git a/vignettes/flanker.Rmd b/vignettes/flanker.Rmd new file mode 100644 index 0000000..9ba58da --- /dev/null +++ b/vignettes/flanker.Rmd @@ -0,0 +1,208 @@ +--- + title: "Analysis of the flanker experiment data using the Bayesian reaction time model" + author: "Jure Demšar, Erik Štrumbelj and Grega Repovš" + date: "`r Sys.Date()`" + output: + html_vignette: + toc: yes +--- + + + +```{r, message=FALSE, warning=FALSE, echo=FALSE} +knitr::opts_chunk$set(fig.width=6, fig.height=4.5) +options(width=800) +``` + +# Introduction + +In the Eriksen flanker task participants are presented with an image of an odd number of arrows (usually five or seven). Their task is to indicate the orientation (left or right) of the middle arrow as quickly as possible whilst ignoring the flanking arrows on left and right. There are two types of stimuli in the task: in the **congruent** condition (e.g. ‘<<<<<<<‘) both the middle arrow and the flanking arrows point in the same direction, whereas in the **incongruent** condition (e.g. ‘<<<><<<‘) the middle arrow points to the opposite direction of the flanking arrows. + +As the participants have to consciously ignore and inhibit the misleading information provided by the flanking arrows in the incongruent condition, the performance in the incongruent condition is robustly worse than in the congruent condition, both in terms of longer reaction times as well as higher proportion of errors. The difference between reaction times and error rates in congruent and incongruent cases is a measure of subject's ability to focus his or her attention and inhibit distracting stimuli. + +In the illustration below we compare reaction times and error rates when solving the flanker task between the control group (healthy subjects) and the test group (subjects suffering from a certain medical condition). + +First, we load package **bayes4psy** and package **dplyr** for data wrangling. Second, we load the data and split them into control and test groups. For reaction time analysis we use only data where the response to the stimuli was correct: + +```{r, message=FALSE, warning=FALSE} +# libs +library(bayes4psy) +library(dplyr) + +# load data +data <- flanker + +# analyse only correct answers, load test and control data +control_rt <- data %>% filter(result == "correct" & + group == "control") + +test_rt <- data %>% filter(result == "correct" & + group == "test") +``` + +The model requires subjects to be indexed from $1$ to $n$. Control group subject indexes range from 22 to 45, so we cast it to 1 to 23. + +```{r, message=FALSE, warning=FALSE} +# control group subject indexes range is 22..45 cast it to 1..23 +# test group indexes are OK +control_rt$subject <- control_rt$subject - 21 +``` + +Now we are ready to fit the Bayesian reaction time model for both groups. The model function requires two parameters -- a vector of reaction times $t$ and the vector of subject indexes $s$. + +```{r, message=FALSE, warning=FALSE, results = 'hide'} +# fit +rt_control_fit <- b_reaction_time(t=control_rt$rt, + s=control_rt$subject) + +rt_test_fit <- b_reaction_time(t=test_rt$rt, + s=test_rt$subject) +``` + +Before we interpret the results, we check MCMC diagnostics and model fit. + +```{r, message=FALSE, warning=FALSE} +# plot trace +plot_trace(rt_control_fit) + +plot_trace(rt_test_fit) +``` + +The traceplots gives us no cause for concern regarding MCMC convergence and mixing. Note that the first 1000 iterations (shaded gray) are used for warmup (tuning of the MCMC algorithm) and are discarded. The next 1000 iterations are used for sampling. + + +```{r, message=FALSE, warning=FALSE} +# check fit (Rhat and n_eff) +print(rt_control_fit) + +#print(rt_test_fit) +``` + + +The output above provides further MCMC diagnostics, which again do not give us cause for concern (only output of one fit is provided for the sake of brevity). The convergence diagnostic **Rhat** is practically 1 for all parameters and there is little auto-correlation (possibly even some positive auto-correlation) -- effective sample sizes (**n\_eff**) are of the order of samples taken and Monte Carlo standard errors (**se\_mean**) are relatively small. + +What is a good-enough effective sample sizes depends on our goal. If we are interested in posterior quantities such as the more extreme percentiles, the effective sample sizes should be 10,000 or higher, if possible. If we are only interested in estimating the mean, 100 effective samples is in most cases enough for a practically negligible Monte Carlo error. + +We can increase the effective sample size by increasing the amount of MCMC iterations with the **iter** parameter. In our case we can achieve an effective sample size of 10,000 by setting **iter** to 4000. Because the MCMC diagnostics give us no reason for concern, we can leave the **warmup** parameter at its default value of 1000. + +```{r, eval=FALSE} +rt_control_fit <- b_reaction_time(t=control_rt$rt, + s=control_rt$subject, + iter=4000) + +rt_test_fit <- b_reaction_time(t=test_rt$rt, + s=test_rt$subject, + iter=4000) +``` + +Because we did not explicitly define any priors, flat (improper) priors were put on all of the model's parameters. In some cases, flat priors are a statement that we have no prior knowledge about the experiment results (in some sense). In general, even flat priors can express a preference for a certain region of parameter space. In practice, we will almost always have some prior information and we should incorporate it into the modelling process. + +Next, we check whether the model fits the data well by using the **plot** function. If we set the **subjects** parameter to **FALSE**, we will get a less detailed group level fit. The data are visualized as a blue region while the fit is visualized with a black line. In this case the model fits the underlying data well. + +```{r, message=FALSE, warning=FALSE} +# check fits +plot(rt_control_fit) + +plot(rt_test_fit) +``` + +We now use the **compare\_means** function to compare reaction times between healthy (control) and unhealthy (test) subjects. In the example below we use a rope (region of practical equivalence) interval of 0.01 s, meaning that differences smaller that 1/100 of a second are deemed as equal. The **compare\_means** function provides a user friendly output of the comparison and also returns the results in the form of a **data.frame**. + +```{r, message=FALSE, warning=FALSE} +# set rope (region of practical equivalence) interval to +/- 10ms +rope <- 0.01 + +# which mean is smaller/larger +rt_control_test <- compare_means(rt_control_fit, fit2=rt_test_fit, rope=rope) +``` + +The **compare\_means** function output contains probabilities that one group has shorter reaction times than the other, the probability that both groups are equal (if rope interval is provided) and the 95\% HDI (highest density interval) for the difference between groups. Based on the output we are quite certain (98\% +/- 0.5\%) that the healthy group's (**rt\_control\_fit**) expected reaction times are lower than the unhealthy group's (**rt\_test\_fit**). + +We can also visualize this difference by using the **plot\_means\_difference** function. The **plot\_means** function is an alternative for comparing **rt\_control\_fit** and **rt\_test\_fit** -- the function visualizes the parameters that determine the means of each model. + +```{r, message=FALSE, warning=FALSE} +# difference plot +plot_means_difference(rt_control_fit, fit2=rt_test_fit, rope=rope) +``` + +The graph above visualizes the difference between **rt\_control\_fit** and **rt\_test\_fit**. The histogram visualizes the distribution of the difference, vertical blue line denotes the mean, the black band at the bottom marks the 95\% HDI interval and the gray band marks the rope interval. Since the whole 95\% HDI of difference is negative and lies outside of the rope interval we can conclude that the statement that healthy subjects are faster than unhealthy ones is most likely correct. + +```{r, message=FALSE, warning=FALSE} +# visual comparsion of mean difference +plot_means(rt_control_fit, fit2=rt_test_fit) +``` + +Above you can see a visualization of means for **rt\_control\_fit** and **rt\_test\_fit**. Group 1 visualizes means for the healthy subjects and group 2 for the unhealthy subjects. + +Based on our analysis we can claim with high probability (98\% +/- 0.5\%) that healthy subjects have faster reaction times when solving the flanker task than unhealthy subjects. Next, we analyze if the same applies to success rates. + +The information about success of subject's is stored as correct/incorrect. However, the Bayesian success rate model requires binary inputs (0/1) so we have to transform the data. Also, like in the reaction time example, we have to correct the indexes of control group subjects. + +```{r, message=FALSE, warning=FALSE} +# map correct/incorrect/timeout to 1 (success) or 0 (fail) +data$result_numeric <- 0 +data[data$result == "correct", ]$result_numeric <- 1 + +# split to control and test groups +control_sr <- data %>% filter(group == "control") +test_sr <- data %>% filter(group == "test") + +# control group subject indexes range is 22..45 cast it to 1..23 +# test group indexes are OK +control_sr$subject <- control_sr$subject - 21 +``` + +Since the only prior information about the success rate of participants was the fact that success rate is located between 0 and 1, we used a beta distribution to put a uniform prior on the [0, 1] interval (we put a $\BE(1, 1)$ prior on the $p$ parameter). We execute the model fitting by running the **b\_success\_rate** function with appropriate input data. + +```{r, message=FALSE, warning=FALSE} +# priors +p_prior <- b_prior(family="beta", pars=c(1, 1)) + +# attach priors to relevant parameters +priors <- list(c("p", p_prior)) + +# fit +sr_control_fit <- b_success_rate(r=control_sr$result_numeric, + s=control_sr$subject, + priors=priors, + iter=4000) + +sr_test_fit <- b_success_rate(r=test_sr$result_numeric, + s=test_sr$subject, + priors=priors, + iter=4000) +``` + +The process for inspecting Bayesian fits is the same as above. When visually inspecting the quality of the fit (the **plot** function) we can set the **subjects** parameter to **FALSE**, which visualizes the fit on the group level. This offers a quicker, but less detailed method of inspection. Again one of the commands is commented out for the sake of brevity. + +```{r, message=FALSE, warning=FALSE} +# inspect control group fit +plot_trace(sr_control_fit) +plot(sr_control_fit, subjects=FALSE) +print(sr_control_fit) + +# inspect test group fit +plot_trace(sr_test_fit) +plot(sr_test_fit, subjects=FALSE) +#print(sr_test_fit) +``` + +Since diagnostic functions show no pressing issues and the fits look good we can proceed with the actual comparison between the two fitted models. We will again estimate the difference between two groups with **compare\_means**. + +```{r, message=FALSE, warning=FALSE} +# which one is higher +sr_control_test <- compare_means(sr_control_fit, fit2=sr_test_fit) +``` + +As we can see the success rate between the two groups is not that different. Since the probability that healthy group is more successful is only 53\% (+/- 1\%) and the 95\% HDI of the difference ([0.02, 0.02]) includes the 0 we cannot claim inequality. We can visualize this result by using the \code{plot\_means\_difference} function. + +```{r, message=FALSE, warning=FALSE} +# difference plot +plot_means_difference(sr_control_fit, fit2=sr_test_fit) +``` + +Above you can see the visualization of the difference between the **sr\_control\_fit** and the **sr\_test\_fit**. Histogram visualizes the distribution of the difference, vertical blue line denotes the mean difference and the black band at the bottom marks the 95\% HDI interval. Since the 95\% HDI of difference includes the value of 0 we cannot claim inequality. If we used a rope interval and the whole rope interval lied in the 95\% HDI interval we could claim equality.