add exgaussian

README.md

@@ -9,7 +9,7 @@ This framework aims at moving away from a mere description of the data, to make
 ## Why Julia?
 
 [**Julia**](https://julialang.org/) - the new cool kid on the scientific block - is a modern programming language with many benefits when compared with R or Python.
-Importantly, it is currently the only language in which we can fit all the cognitive models under a Bayesian framework using a unified interface like [**Turing**](https://turing.ml/) and [**SSM**](https://github.com/itsdfish/SequentialSamplingModels.jl).
+Importantly, it is currently the only language in which we can fit all the cognitive models under a Bayesian framework using a unified interface like [**Turing**](https://turing.ml/) and [**SequentialSamplingModels**](https://github.com/itsdfish/SequentialSamplingModels.jl).
 
 ## Why Bayesian?
 

content/.quarto/cites/index.json

@@ -1 +1 @@
-{"index.qmd":[],"1_introduction.qmd":[],"2_predictors.qmd":[],"5_individual.qmd":[],"references.qmd":[],"4_1_Normal.qmd":["wagenmakers2008diffusion","theriault2024check","lo2015transform","schramm2019reaction"],"3_scales.qmd":[],"4_rt.qmd":["wagenmakers2008diffusion","heathcote2012linear","theriault2024check","lo2015transform","schramm2019reaction"]}
+{"references.qmd":[],"4_rt.qmd":["wagenmakers2008diffusion","heathcote2012linear","theriault2024check","lo2015transform","schramm2019reaction","balota2011moving","matzke2009psychological","hohle1965inferred","kieffaber2006switch","matzke2009psychological"],"2_predictors.qmd":[],"5_individual.qmd":[],"3_scales.qmd":[],"4_1_Normal.qmd":["wagenmakers2008diffusion","theriault2024check","lo2015transform","schramm2019reaction"],"1_introduction.qmd":[],"index.qmd":[]}

content/.quarto/xref/1a47137c

@@ -1 +1 @@
-{"options":{"chapters":true},"entries":[],"headings":[]}
+{"headings":[],"entries":[],"options":{"chapters":true}}

content/.quarto/xref/26afb962

@@ -1 +1 @@
-{"entries":[],"options":{"chapters":true},"headings":["categorical-predictors-condition-group","ordered-predictors-likert-scales","interactions","non-linear-relationships-polynomial-gams"]}
+{"options":{"chapters":true},"headings":["categorical-predictors-condition-group","ordered-predictors-likert-scales","interactions","non-linear-relationships-polynomial-gams"],"entries":[]}

content/.quarto/xref/6cbe9151

@@ -1 +1 @@
-{"entries":[],"options":{"chapters":true},"headings":["preface","why-julia","why-bayesian","the-plan","looking-for-coauthors"]}
+{"options":{"chapters":true},"entries":[],"headings":["preface","why-julia","why-bayesian","the-plan","looking-for-coauthors"]}

content/.quarto/xref/efe17597

@@ -1 +1 @@
-{"headings":["the-data","descriptive-models","gaussian-aka-linear-model","model-specification","posterior-predictive-check","scaled-gaussian-model","solution-1-directional-effect-of-condition","solution-2-avoid-exploring-negative-variance-values","the-problem-with-linear-models","shifted-lognormal-model","model","more-conditional-parameters","exgaussian-model","wald-model","generative-models-ddm","other-models-lba-lnr","including-random-effects","additional-resources"],"options":{"chapters":true},"entries":[]}
+{"entries":[],"options":{"chapters":true},"headings":["the-data","descriptive-models","gaussian-aka-linear-model","model-specification","posterior-predictive-check","scaled-gaussian-model","solution-1-directional-effect-of-condition","solution-2-avoid-exploring-negative-variance-values","the-problem-with-linear-models","shifted-lognormal-model","prior-on-minimum-rt","model-specification-1","interpretation","exgaussian-model","model-specification-2","wald-model","generative-models-ddm","other-models-lba-lnr","including-random-effects","additional-resources"]}

content/4_rt.qmd

@@ -85,10 +85,10 @@ In order to fit a Linear Model for RTs, we need to set a prior on all these para
 @model function model_Gaussian(rt; condition=nothing)
 
     # Set priors on variance, intercept and effect of condition
-    σ ~ truncated(Normal(0, 1); lower=0)
+    σ ~ truncated(Normal(0, 0.5); lower=0)
 
     μ_intercept ~ truncated(Normal(0, 1); lower=0)
-    μ_condition ~ Normal(0, 0.5)
+    μ_condition ~ Normal(0, 0.3)
 
     for i in 1:length(rt)
         μ = μ_intercept + μ_condition * condition[i]
@@ -161,9 +161,9 @@ And such effect would lead to a sigma $\sigma$ of **0.14 - 0.15 = -0.01**, which
 
 #### Solution 1: Directional Effect of Condition
 
-One possible (but not recommended) solution is to simply make it impossible for the effect of condition to be negative. 
+One possible (but not recommended) solution is to simply make it impossible for the effect of condition to be negative by *Truncating* the prior to a lower bound of 0. 
 This can work in our case, because we know that the comparison condition is likely to have a higher variance than the reference condition (the intercept) - and if it wasn't the case, we could have changed the reference factor.
-
+However, this is not a good practice as we are enforcing a very strong a priori specific direction of the effect, which is not always justified.
 
 ```{julia}
 #| code-fold: false
@@ -173,9 +173,9 @@ This can work in our case, because we know that the comparison condition is like
 
     # Priors
     μ_intercept ~ truncated(Normal(0, 1); lower=0)
-    μ_condition ~ Normal(0, 0.5)
+    μ_condition ~ Normal(0, 0.3)
 
-    σ_intercept ~ truncated(Normal(0, 1); lower=0)  # Same prior as previously
+    σ_intercept ~ truncated(Normal(0, 0.5); lower=0)  # Same prior as previously
     σ_condition ~ truncated(Normal(0, 0.1); lower=0)  # Enforce positivity
 
     for i in 1:length(rt)
@@ -211,9 +211,9 @@ This can be done by adding a conditional statement when sigma $\sigma$ is negati
 
     # Priors
     μ_intercept ~ truncated(Normal(0, 1); lower=0)
-    μ_condition ~ Normal(0, 0.5)
+    μ_condition ~ Normal(0, 0.3)
 
-    σ_intercept ~ truncated(Normal(0, 1); lower=0)
+    σ_intercept ~ truncated(Normal(0, 0.5); lower=0)
     σ_condition ~ Normal(0, 0.1)
 
     for i in 1:length(rt)
@@ -273,128 +273,113 @@ Instead, rather than applying arbitrary data transformation, it would be better
 ### Shifted LogNormal Model
 
 One of the obvious candidate alternative to the log-transformation would be to use a model with a Log-transformed Normal distribution.
+A LogNormal distribution is a distribution of a random variable whose logarithm is normally distributed.
+A **Shifted** LogNormal model introduces a shift (a delay) parameter *tau* $\tau$ that corresponds to the minimum "starting time" of the response process.
+
+We need to set a prior for this parameter, which is usually truncated between 0 (to exclude negative minimum times) and the minimum RT of the data (the logic being that the minimum delay for response must be lower than the faster response actually observed).
+
+While $Uniform(0, min(RT))$ is a common choice of prior, it is not ideal as it implies that all values between 0 and the minimum RT are equally likely, which is not the case.
+Indeed, psychology research has shown that such minimum response time for Humans is often betwen 100 and 250 ms. 
+Moreover, in our case, we explicitly removed all RTs below 180 ms, suggesting that the minimum response time is more likely to approach 180 ms than 0 ms.
 
-New parameter, $\tau$ (Tau for delay), which corresponds to the "starting time". We need to set a prior for this parameter, which is usually truncated between 0 and the minimum RT of the data (the logic being that the minimum delay for response must be lower than the faster response actually observed).
+#### Prior on Minimum RT
 
+Instead of a $Uniform$ prior, we will use a $Gamma(1.1, 11)$ distribution (truncated at min. RT), as this particular parameterization reflects the low probability of very low minimum RTs (near 0) and a steadily increasing probability for increasing times.  
 ```{julia}
-xaxis = range(0, 1, 1000)
-lines(xaxis, pdf.(Gamma(1.1, 11), xaxis))
+xaxis = range(0, 0.3, 1000)
+fig = lines(xaxis, pdf.(Gamma(1.1, 11), xaxis); color=:blue, label="Gamma(1.1, 11)")
+vlines!([minimum(df.RT)]; color="red", linestyle=:dash, label="Min. RT = 0.18 s")
+axislegend()
+fig
 ```
 
-#### Model
+
+#### Model Specification
 
 ```{julia}
 #| code-fold: false
+#| output: false
 
-@model function model_lognormal(rt; min_rt=minimum(df.RT), condition=nothing)
+@model function model_LogNormal(rt; min_rt=minimum(df.RT), condition=nothing)
 
     # Priors 
-    σ ~ truncated(Normal(0, 0.5); lower=0)
     τ ~ truncated(Gamma(1.1, 11); upper=min_rt)
 
-    μ_intercept ~ Normal(0, 2)
-    μ_condition ~ Normal(0, 0.5)
+    μ_intercept ~ Normal(0, exp(1))  # On the log-scale: exp(μ) to get value in seconds
+    μ_condition ~ Normal(0, exp(0.3))
+
+    σ_intercept ~ truncated(Normal(0, 0.5); lower=0)
+    σ_condition ~ Normal(0, 0.1)
 
     for i in 1:length(rt)
         μ = μ_intercept + μ_condition * condition[i]
+        σ = σ_intercept + σ_condition * condition[i]
+        if σ < 0  # Avoid negative variance values
+            Turing.@addlogprob! -Inf
+            return nothing
+        end
         rt[i] ~ ShiftedLogNormal(μ, σ, τ)
     end
 end
 
-model = model_lognormal(df.RT; condition=df.Accuracy)
-chain_lognormal = sample(model, NUTS(), 400)
-
-# Summary (95% CI)
-quantile(chain_lognormal; q=[0.025, 0.975])
+model = model_LogNormal(df.RT; condition=df.Accuracy)
+chain_LogNormal = sample(model, NUTS(), 400)
 ```
 
+#### Interpretation
+
 ```{julia}
-#| output: false
+#| code-fold: false
 
-pred = predict(model_lognormal([(missing) for i in 1:length(df.RT)]; condition=df.Accuracy), chain_lognormal)
-pred = Array(pred)
+hpd(chain_LogNormal; alpha=0.05)
 ```
 
+
 ```{julia}
-#| fig-width: 10
-#| fig-height: 7
+pred = predict(model_LogNormal([(missing) for i in 1:length(df.RT)]; condition=df.Accuracy), chain_LogNormal)
+pred = Array(pred)
 
 fig = plot_distribution(df, "Predictions made by Shifted LogNormal Model")
-for i in 1:length(chain_lognormal)
+for i in 1:length(chain_LogNormal)
     lines!(Makie.KernelDensity.kde(pred[:, i]), color=ifelse(df.Accuracy[i] == 1, "#388E3C", "#D32F2F"), alpha=0.1)
 end
 fig
 ```
 
-#### More Conditional Parameters
-
-```{julia}
-xaxis = range(-6, 2, 1000)
-fig = Figure()
-ax = Axis(fig[1, 1])
-lines!(xaxis, pdf.(-Weibull(2, 2.5)+1, xaxis); color="red")
-lines!(xaxis, pdf.(-Weibull(2.5, 2)+1, xaxis); color="orange")
-lines!(xaxis, pdf.(-Weibull(2.5, 2.5)+1, xaxis); color="blue")
-lines!(xaxis, pdf.(-Weibull(2.5, 3)+1, xaxis); color="green")
-ax2 = Axis(fig[1, 2])
-lines!(exp.(xaxis), pdf.(-Weibull(2, 2.5)+1, xaxis); color="red")
-lines!(exp.(xaxis), pdf.(-Weibull(2.5, 2)+1, xaxis); color="orange")
-lines!(exp.(xaxis), pdf.(-Weibull(2.5, 2.5)+1, xaxis); color="blue")
-lines!(exp.(xaxis), pdf.(-Weibull(2.5, 3)+1, xaxis); color="green")
-fig
-```
-
-```{julia}
-#| code-fold: false
+This model provides a much better fit to the data, and confirms that the `Accuracy` condition is associated with higher RTs and higher variability (i.e., a larger distribution width).
 
-@model function model_lognormal2(rt; min_rt=minimum(df.RT), condition=nothing)
+### ExGaussian Model
 
-    # Priors 
-    τ ~ truncated(Gamma(1.1, 11); upper=min_rt)
+Another popular model to describe RTs uses the **ExGaussian**  distribution, i.e., the *Exponentially-modified Gaussian* distribution [@balota2011moving; @matzke2009psychological].
 
-    μ_intercept ~ Normal(0, 2)
-    μ_condition ~ Normal(0, 0.5)
+This distribution is a convolution of normal and exponential distributions and has three parameters, namely *mu* $\mu$ and *sigma* $\sigma$ - the mean and standard deviation of the Gaussian distribution - and *tau* $\tau$ - the exponential component of the distribution (note that although denoted by the same letter, it does not correspond directly to a shift of the distribution). 
+Intuitively, these parameters reflect the centrality, the width and the tail dominance, respectively.
 
-    σ_intercept ~ -Weibull(2.5, 3) + 1
-    σ_condition ~ Normal(0, 0.01)
+![](media/rt_exgaussian.gif)
 
-    for i in 1:length(rt)
-        μ = μ_intercept + μ_condition * condition[i]
-        σ = σ_intercept + σ_condition * condition[i]
-        rt[i] ~ ShiftedLogNormal(μ, exp(σ), τ)
-    end
-end
 
-model = model_lognormal2(df.RT; condition=df.Accuracy)
-chain_lognormal2 = sample(model, NUTS(), 400)
+Beyond the descriptive value of these types of models, some have tried to interpret their parameters in terms of **cognitive mechanisms**, arguing for instance that changes in the Gaussian components ($\mu$ and $\sigma$) reflect changes in attentional processes [e.g., "the time required for organization and execution of the motor response"; @hohle1965inferred], whereas changes in the exponential component ($\tau$) reflect changes in intentional (i.e., decision-related) processes [@kieffaber2006switch]. 
+However, @matzke2009psychological demonstrate that there is likely no direct correspondence between ex-Gaussian parameters and cognitive mechanisms, and underline their value primarily as **descriptive tools**, rather than models of cognition *per se*.
 
-# Summary (95% CI)
-quantile(chain_lognormal2; q=[0.025, 0.975])
-```
+Descriptively, the three parameters can be interpreted as:
 
-```{julia}
-#| output: false
+- **Mu** $\mu$: The location / centrality of the RTs. Would correspond to the mean in a symmetrical distribution.
+- **Sigma** $\sigma$: The variability and dispersion of the RTs. Akin to the standard deviation in normal distributions.
+- **Tau** $\tau$: Tail weight / skewness of the distribution.
 
-pred = predict(model_lognormal2([(missing) for i in 1:length(df.RT)]; condition=df.Accuracy), chain_lognormal2)
-pred = Array(pred)
-```
-
-```{julia}
-#| fig-width: 10
-#| fig-height: 7
-
-fig = plot_distribution(df, "Predictions made by Shifted LogNormal Model")
-for i in 1:length(chain_lognormal2)
-    lines!(Makie.KernelDensity.kde(pred[:, i]), color=ifelse(df.Accuracy[i] == 1, "#388E3C", "#D32F2F"), alpha=0.1)
-end
-fig
-```
+::: {.callout-important}
+Note that these parameters are not independent with respect to distribution characteristics, such as the empirical mean and SD. 
+Below is an example of different distributions with the same location (*mu* $\mu$) and dispersion (*sigma* $\sigma$) parameters.
+Although only the tail weight parameter (*tau* $\tau$) is changed, the whole distribution appears to shift is centre of mass. 
+Hence, one should be careful note to interpret the values of *mu* $\mu$ directly as the "mean" or the distribution peak and *sigma* $\sigma$ as the SD or the "width".
+:::
 
+![](media/rt_exgaussian2.gif)
 
-### ExGaussian Model
+#### Model Specification
 
+TODO.
 
-- [**Ex-Gaussian models in R: A Tutorial**](https://dominiquemakowski.github.io/easyRT/articles/exgaussian.html)
 
 ### Wald Model
 

content/index.qmd

@@ -3,7 +3,7 @@
 ## Why Julia?
 
 [**Julia**](https://julialang.org/) - the new cool kid on the scientific block - is a modern programming language with many benefits when compared with R or Python.
-Importantly, it is currently the only language in which we can fit all the cognitive models under a Bayesian framework using a unified interface like [**Turing**](https://turing.ml/) and [**SSM**](https://github.com/itsdfish/SequentialSamplingModels.jl).
+Importantly, it is currently the only language in which we can fit all the cognitive models under a Bayesian framework using a unified interface like [**Turing**](https://turing.ml/) and [**SequentialSamplingModels**](https://github.com/itsdfish/SequentialSamplingModels.jl).
 
 ## Why Bayesian?
 

content/media/animations_rt.jl

@@ -71,28 +71,93 @@ end
 frames = range(0, 1, length=60)
 record(make_animation, fig, "rt_normal.gif", frames; framerate=30)
 
-# LogNormal =====================================================================================
+# ExGaussian =====================================================================================
 
 # Parameters
-# μ = Observable(0.0)
-# σ = Observable(0.2)
-
-# x = range(-0.1, 0.8, length=100)
-# # y = @lift(pdf.(Normal($μ, $σ), x))
-
-# # Initialize the figure
-# fig = Figure()
-# ax = Axis(
-#     fig[1, 1],
-#     # title=@lift("LogNormal(μ = $(round($μ, digits = 1)), σ =  $(round($σ, digits = 2)))"),
-#     xlabel="RT (s)",
-#     ylabel="Distribution",
-#     yticksvisible=false,
-#     xticksvisible=false,
-#     yticklabelsvisible=false,
-# )
-
-# density!(ax, df.RT, bandwidth=0.01, npoints=1000, color=:grey)
-# lines!(x, pdf.(LogNormal(), x), linestyle=:dot, color=:red)  # Best fitting line
-# lines!(x, y, linewidth=4, color=:orange)
-# fig
+μ = Observable(0.0)
+σ = Observable(0.4)
+τ = Observable(0.1)
+
+x = range(-0.1, 2, length=1000)
+y = @lift(pdf.(ExGaussian($μ, $σ, $τ), x))
+
+# Initialize the figure
+fig = Figure()
+ax = Axis(
+    fig[1, 1],
+    title=@lift("ExGaussian(μ = $(round($μ, digits = 1)), σ =  $(round($σ, digits = 2)), τ = $(round($τ, digits = 2)))"),
+    xlabel="RT (s)",
+    ylabel="Distribution",
+    yticksvisible=false,
+    xticksvisible=false,
+    yticklabelsvisible=false,
+)
+density!(ax, df.RT, npoints=1000, color=:grey)
+lines!(x, pdf.(ExGaussian(0.4, 0.06, 0.2), x), linestyle=:dot, color=:red)  # Best fitting line
+lines!(x, y, linewidth=4, color=:orange)
+fig
+
+function make_animation(frame)
+    if frame < 0.2
+        μ[] = change_param(frame; frame_range=(0, 0.2), param_range=(0, 0.4))
+    end
+    if frame >= 0.2 && frame < 0.4
+        σ[] = change_param(frame; frame_range=(0.2, 0.4), param_range=(0.4, 0.1))
+    end
+    if frame >= 0.4 && frame < 0.6
+        τ[] = change_param(frame; frame_range=(0.4, 0.6), param_range=(0.1, 0.4))
+    end
+    # Return to normal
+    if frame >= 0.7
+        μ[] = change_param(frame; frame_range=(0.7, 1), param_range=(0.4, 0))
+        σ[] = change_param(frame; frame_range=(0.7, 1), param_range=(0.1, 0.4))
+        τ[] = change_param(frame; frame_range=(0.7, 1), param_range=(0.4, 0.1))
+    end
+end
+
+# animation settings
+frames = range(0, 1, length=60)
+record(make_animation, fig, "rt_exgaussian.gif", frames; framerate=30)
+
+
+# ExGaussian 2 =====================================================================================
+# Parameters
+μ = Observable(0.3)
+σ = Observable(0.2)
+τ = Observable(0.006)
+
+x = range(-0.4, 2, length=1000)
+y = @lift(pdf.(ExGaussian($μ, $σ, $τ), x))
+
+m = Observable(mean(rand(ExGaussian(0.3, 0.2, 0.006), 1000)))
+
+# Initialize the figure
+fig = Figure()
+ax = Axis(
+    fig[1, 1],
+    title=@lift("ExGaussian(μ = $(round($μ, digits = 1)), σ =  $(round($σ, digits = 2)), τ = $(round($τ, digits = 3)))"),
+    xlabel="RT (s)",
+    ylabel="Distribution",
+    yticksvisible=false,
+    xticksvisible=false,
+    yticklabelsvisible=false,
+)
+lines!(x, y, linewidth=4, color=:orange)
+vlines!(m, color=:red, linestyle=:dot, label="Average RT")
+leg = axislegend(position=:rt)
+fig
+
+function make_animation(frame)
+    if frame < 0.5
+        τ[] = change_param(frame; frame_range=(0, 0.5), param_range=(0.006, 0.4))
+    end
+    # Return to normal
+    if frame >= 0.5
+        τ[] = change_param(frame; frame_range=(0.5, 1), param_range=(0.4, 0.006))
+    end
+    m[] = mean(rand(ExGaussian(0.3, 0.2, τ[]), 100_000))
+end
+
+# animation settings
+frames = range(0, 1, length=60)
+record(make_animation, fig, "rt_exgaussian2.gif", frames; framerate=30)