diff --git a/.nojekyll b/.nojekyll
index b6f4951..4dd38eb 100644
--- a/.nojekyll
+++ b/.nojekyll
@@ -1 +1 @@
-2e7347ed
\ No newline at end of file
+7bab3c81
\ No newline at end of file
diff --git a/1_introduction.html b/1_introduction.html
index 5871718..e8c972e 100644
--- a/1_introduction.html
+++ b/1_introduction.html
@@ -2,7 +2,7 @@
Family: gaussian
+ Links: mu = identity; sigma = identity
+Formula: mpg ~ mo(gear)
+ Data: mtcars (Number of observations: 32)
+ Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
+ total post-warmup draws = 4000
+
+Regression Coefficients:
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+Intercept 16.51 1.31 14.01 19.19 1.00 2106 2213
+mogear 3.88 1.02 1.81 5.86 1.00 2160 2315
+
+Monotonic Simplex Parameters:
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+mogear1[1] 0.84 0.13 0.52 0.99 1.00 2698 1862
+mogear1[2] 0.16 0.13 0.01 0.48 1.00 2698 1862
+
+Further Distributional Parameters:
+ Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
+sigma 5.02 0.68 3.90 6.58 1.00 2850 2485
+
+Draws were sampled using sampling(NUTS). For each parameter, Bulk_ESS
+and Tail_ESS are effective sample size measures, and Rhat is the potential
+scale reduction factor on split chains (at convergence, Rhat = 1).
+
#| code-fold: false
+
+# using Downloads, CSV, DataFrames, Random
+# using Turing, Distributions
+# using GLMakie
+# using LinearAlgebra
+
+# # Define the monotonic function
+# function monotonic(scale::Vector{Float64}, i::Int)
+# if i == 0
+# return 0.0
+# else
+# return length(scale) * sum(scale[1:i])
+# end
+# end
+
+# # Test
+# using RDatasets
+# data = dataset("datasets", "mtcars")
+
+# # Response variable
+# y = data[!, :MPG]
+# x = data[!, :Gear]
+
+# # Number of observations
+# N = length(y)
+
+# # Length of simplex
+# Jmo = maximum(x)
+
+# # Prior concentration for the simplex (using a uniform prior for simplicity)
+# con_simo_1 = ones(Jmo)
+
+# # Turing Model Specification
+# @model function monotonic_model(y, x, N, Jmo, con_simo_1)
+# # Parameters
+# Intercept ~ Normal(19.2, 5.4)
+# simo_1 ~ Dirichlet(con_simo_1)
+# sigma ~ truncated(Normal(0, 5.4), 0, Inf)
+# bsp ~ Normal(0, 1)
+
+# # Linear predictor
+# mu = Vector{Float64}(undef, N)
+# for n in 1:N
+# mu[n] = Intercept + bsp * monotonic(simo_1, x[n])
+# end
+
+# # Likelihood
+# y ~ MvNormal(mu, sigma)
+# end
+
+# # Run the model
+# model = monotonic_model(y, x, N, Jmo, con_simo_1)
+# chain = sample(model, NUTS(), 1000)
+
+# # Display the results
+# display(chain)
+
2.4 Non-linear Relationships
-
-
usingDownloads, CSV, DataFrames, Random
-usingTuring, Distributions
-usingGLMakie
-
-Random.seed!(123) # For reproducibility
-
-df = CSV.read(Downloads.download("https://raw.githubusercontent.com/DominiqueMakowski/CognitiveModels/main/data/nonlinear.csv"), DataFrame)
-
-# Show 10 first rows
-scatter(df.Age, df.SexualDrive, color=:black)
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
-
-
-
-
+
It is a common misconception that “linear models” can only model straight relationships between variables. In fact, the “linear” part refers to the relationship between the parameters (between the outcome varialbe and the predictors).
+
#| code-fold: false
+
+using Downloads, CSV, DataFrames, Random
+using Turing, Distributions
+using GLMakie
+
+Random.seed!(123) # For reproducibility
+
+df = CSV.read(Downloads.download("https://raw.githubusercontent.com/DominiqueMakowski/CognitiveModels/main/data/nonlinear.csv"), DataFrame)
+
+# Show 10 first rows
+scatter(df.Age, df.SexualDrive, color=:black)
Let’s take the data from Makowski et al. (2023) that contains data from participants that underwent the Mini-IPIP6 personality test and the PID-5 BF questionnaire for “maladaptive” personality. We will focus on the “Disinhibition” trait from the PID-5 BF questionnaire. Note that altough it is usually computed as the average of items a 4-point Likert scales [0-3], this study used analog slides to obtain more finer-grained scores.
-
+
usingDownloads, CSV, DataFrames, RandomusingTuring, Distributions, StatsFunsusingGLMakie
-
-Random.seed!(123) # For reproducibility
-
-df = CSV.read(Downloads.download("https://raw.githubusercontent.com/DominiqueMakowski/CognitiveModels/main/data/makowski2023.csv"), DataFrame)
-
-# Show 10 first rows
-first(df, 10)
-
-# Plot the distribution of "Disinhibition"
-hist(df.Disinhibition, normalization =:pdf, color=:darkred, bins=40)
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
+usingSubjectiveScalesModels
+
+Random.seed!(123) # For reproducibility
+
+df = CSV.read(Downloads.download("https://raw.githubusercontent.com/DominiqueMakowski/CognitiveModels/main/data/makowski2023.csv"), DataFrame)
+
+# Show 10 first rows
+first(df, 10)
+
+# Plot the distribution of "Disinhibition"
+hist(df.Disinhibition, normalization =:pdf, color=:darkred, bins=40)
-
+
We will then fit a simple Gaussian model (an “intercept-only” linear model) that estimates the mean and the standard deviation of our variable of interest.
-
-
@modelfunctionmodel_Gaussian(x)
-
-# Priors
- σ ~truncated(Normal(0, 1); lower=0) # Strictly positive half normal distribution
- μ ~Normal(0, 3)
-
-# Iterate through every observation
-for i in1:length(x)
-# Likelihood family
- x[i] ~Normal(μ, σ)
-end
-end
-
-# Fit the model with the data
-fit_Gaussian =model_Gaussian(df.Disinhibition)
-# Sample results using MCMC
-chain_Gaussian =sample(fit_Gaussian, NUTS(), 400)
+
+
@modelfunctionmodel_Gaussian(x)
+
+# Priors
+ σ ~truncated(Normal(0, 1); lower=0) # Strictly positive half normal distribution
+ μ ~Normal(0, 3)
+
+# Iterate through every observation
+for i in1:length(x)
+# Likelihood family
+ x[i] ~Normal(μ, σ)
+end
+end
+
+# Fit the model with the data
+fit_Gaussian =model_Gaussian(df.Disinhibition)
+# Sample results using MCMC
+chain_Gaussian =sample(fit_Gaussian, NUTS(), 400)
Let see if the model managed to recover the mean and standard deviation of the data:
-
+
Code
-
println("Mean of the data: $(round(mean(df.Disinhibition); digits=3)) vs. mean from the model: $(round(mean(chain_Gaussian[:μ]); digits=3))")
-println("SD of the data: $(round(std(df.Disinhibition); digits=3)) vs. SD from the model: $(round(mean(chain_Gaussian[:σ]); digits=3))")
+
println("Mean of the data: $(round(mean(df.Disinhibition); digits=3)) vs. mean from the model: $(round(mean(chain_Gaussian[:μ]); digits=3))")
+println("SD of the data: $(round(std(df.Disinhibition); digits=3)) vs. SD from the model: $(round(mean(chain_Gaussian[:σ]); digits=3))")
-
Mean of the data: 1.064 vs. mean from the model: 1.066
-SD of the data: 0.616 vs. SD from the model: 0.618
+
Mean of the data: 1.064 vs. mean from the model: 1.063
+SD of the data: 0.616 vs. SD from the model: 0.619
Impressive! The model managed to almost perfectly recover the mean and standard deviation of the data. That means we must have a good model, right? Not so fast!
Linear models are by definition designed at recovering the mean of the outcome variables (and its SD, assuming it is inavriant across groups). That does not mean that they can capture the full complexity of the data.
Let us then jump straight into generating predictions from the model and plotting the results against the actual data to see how well the model fits the data (a procedure called the posterior predictive check).
-
+
Code
-
pred =predict(model_Gaussian([(missing) for i in1:length(df.Disinhibition)]), chain_Gaussian)
-pred =Array(pred)
+
pred =predict(model_Gaussian([(missing) for i in1:length(df.Disinhibition)]), chain_Gaussian)
+pred =Array(pred)
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
-
+
+
+
+
+
+
+
+
+Code Tip
+
+
+
+
The plotting functions from the Makie package can take a tuple of a color and an alpha value (e.g., color=(:red, 0.3)) to set the transparency of the plotted elements.
-
As we can see, the model assumes that the data is normally distributed, in a way that allows for negative values and values above 3, which are not possible because our variable is, by design, bounded between 0 and 3 (at it is the result of the mean of variables between 0 and 3). The linear model might thus not be the best choice for this type of data.
+
+
As we can see, the model assumes that the data is normally distributed, in a way that allows for negative values and values above 3, which are not possible as our variable is - by design - bounded between 0 and 3 (at it is the result of the mean of variables between 0 and 3). The linear model might thus not be the best choice for this type of data.
3.2 Rescaling
-
Continuous variables can be trivially rescaled, which is often done to improve the interpretability of the results. For instance, a z-score is a rescaled variable with a mean of 0 and a standard deviation of 1. Importantly, rescaling variables does not change the variable’s distribution or the absolute relationship between variables. It does not alter the fundamental conclusions from an analysis, but can help with the interpretation of the results (or computational performance).
-
Two common rescalings are standardization (mean=0, SD=1) and normalization (range 0-1). The benefits of standardization are the interpretation in terms of deviation (which can be compared across variables). The benefits of normalization are the interpretation in terms of “proportion” (percentage): a value of \(0.5\) (i.e., \(50\%\)) means that the value is in the middle of the range. The latter is particularly useful for bounded variables, as it redefines their bounds to 0 and 1 and allows for a more intuitive interpretation (e.g., “Condition B led to an increase of 20% of the rating on that scale”).
+
Continuous variables can be trivially rescaled, which is often done to improve the interpretability of the results. For instance, a z-score is a rescaled variable with a mean of 0 and a standard deviation of 1 that allows for an interpretation in terms of deviation from the mean. Importantly, rescaling variables does not change the variable’s distribution or the absolute relationship between variables. In other words, it does not alter the fundamental conclusions from an analysis, but can help with the interpretation of the results (or computational performance).
+
Two common rescalings are standardization (mean=0, SD=1) and normalization (to a 0-1 range). The benefits of standardization are the interpretation in terms of the magnitude of deviation relative to its own sample (which can be compared across variables). The benefits of normalization are the interpretation in terms of “proportion” (percentage): a value of \(0.5\) (i.e., \(50\%\)) means that the value is in the middle of the range. The latter is particularly useful for bounded variables, as it redefines their bounds to 0 and 1 and allows for a more intuitive interpretation (e.g., “Condition B led to an increase of 20% of the rating on that scale”).
Let’s rescale the “Disinhibition” variable to a 0-1 range:
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
-
+
3.3 Modified Beta Distribution
-
One good potential alternative to linear models for bounded variables is to use a Beta distribution instead of a Gaussian distribution, as the Beta distribution is bounded between 0 and 1 (not including them). Moreover, Beta distributions are powerful and can model a wide range of shapes, including normal-like distributions, but also uniformly spread data and data clustered at one or both ends.
+
One good potential alternative to linear models for bounded variables is to use a Beta distribution instead of a Gaussian distribution, as the Beta distribution is naturally bounded between 0 and 1 (not including them). Moreover, Beta distributions are powerful and can model a wide variety of shapes, including normal-like distributions, but also uniformly spread data and data clustered at one or both ends.
The Beta distribution is typically defined by two parameters, alpha\(\alpha\) and beta\(\beta\), which are the shape parameters. Unfortunately, these parameters are not very intuitive, and so we often use a “reparametrization” of the Beta distribution to define it by its mean mu\(\mu\) and “precision” phi\(\phi\) (referring to the narrowness of the distribution). This is particularly convenient in the context of regressions, as these parameters are more interpretable and can be directly linked to the predictors.
-
Here is the code to redefine a Beta distribution based on the mean mu\(\mu\) and precision phi\(\phi\), that converts them to the shape parameters alpha\(\alpha\) and beta\(\beta\).
TODO: Replace if it is supported somewhere (e.g., Distributions.jl)
-
-
-
Note that for this modified Beta distribution, the mean mu\(\mu\) and precision phi\(\phi\) can impact the distribution in suprising ways.
-
+
We will use a reparametrization of the Beta distribution based on the mean mu\(\mu\) and precision phi\(\phi\), that converts them to the shape parameters alpha\(\alpha\) and beta\(\beta\). You can find details about this reparametrization in the documentation of the BetaPhi2() function of the SubjectiveScalesModels.jl package.
+
3.4 Beta Models
Note that we have a suited distribution for our bounded variable, we can now fit a Beta model to the rescaled variable. However, there is one important issue to address: the Beta distribution is not defined at exactly 0 and 1, and we currently rescaled our variable to be between 0 and 1, possibly including them.
One common trick is to actually rescale our variable to be within \([0, 1]\) by nudging the zeros and ones to be just above and below, respectively. For this, one can use the function eps(), which returns the smallest possible number. For instance, one can rescale the variable to be in the range [eps(), 1 - eps()], equivalent to \([0.000...1, 0.999...]\).
For the priors, we will use a Beta distribution \(Beta(1.25, 1.25)\) for mu\(\mu\) that is naturally bounded at \(]0, 1[\), peaks at 0.5, and assign less plausibility to extreme values. A Gamma distribution \(Gamma(1.5, 15)\) for phi\(\phi\) is a good choice for the precision, as it is naturally bounded at \(]0, +\infty[\). That said, as we demonstrate below, it is often more convenient to express phi\(\phi\) on the log scale (making it more convenient to express priors).
-
+
For the priors, we will use a Beta distribution \(Beta(1.25, 1.25)\) for mu\(\mu\) that is naturally bounded at \(]0, 1[\), peaks at 0.5, and assign less plausibility to extreme values. A Gamma distribution \(Gamma(1.5, 15)\) for phi\(\phi\) is a good choice for the precision, as it is naturally bounded at \(]0, +\infty[\).
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
-
-
+
+
We can now fit a Beta model to the rescaled variable.
Note that it is also common to express mu\(\mu\) on the logit scale. In other words, the prior on mu\(\mu\) can be specified using any unbounded distributions (e.g., \(Normal(0, 1)\)), which is convenient to set effect coefficients. The resulting value is passed through a logistic function that transforms any values to the [0, 1] range (suited for mu\(\mu\)). We will demonstrate this below.
+
Note that it is actually convenient to express the parameters with constraints such as a specific range or sign (e.g., positive) on transformed scales (e.g., the log scale) rather than the original scale. We will demonstrate this better way of parametrizing a model below.
Let us see make a posterior predictive check to see how well the model fits the data.
-
+
Code
-
pred =predict(model_Beta([(missing) for i in1:nrow(df)]), chain_Beta)
-pred =Array(pred)
-
-fig =hist(df.Disinhibition3, normalization =:pdf, color=:darkred, bins=40)
-for i in1:length(chain_Beta)
-density!(pred[i, :], color=(:black, 0), strokecolor = (:black, 0.1), strokewidth =3, boundary=(0, 1))
-end
-fig
+
pred =predict(model_Beta([(missing) for i in1:nrow(df)]), posteriors_Beta)
+pred =Array(pred)
+
+bins =range(0, 1, length=40)
+
+fig =hist(df.Disinhibition3, normalization =:pdf, color=:darkred, bins=bins)
+for i in1:length(posteriors_Beta)
+hist!(pred[i, :], normalization =:pdf, bins=bins, color=(:black, 0.01))
+end
+fig
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
+
+
+
+
+
+
+
+
+
+
+Code Tip
-
-
+
+
The distribution of bounded data is typically difficult to visualize. Density plots will tend to be very distorted at the edges (due to the Gaussian kernel used), and histograms will be dependent on the binning. One option is to specify the bin edges in a consistent way.
-
It is… quite terrible. Why?
+
+
As we can see, the model produced a distribution that appears to be collapsed to the left, not reflecting the actual data. The model fit appears quite terrible. Why?
3.5 Excluding Extreme Observations
-
One of the main issues is that, as you can se from the histogram, there is a high number of observations clumped at zero. This creates a bimodal distribution which makes standard unimodal distributions fail to capture the data (note this issues is not addressed by linear models, which estimates will get biased by this configuration away from the actual “mean” of the variable).
-
One simple, although not ideal, solution is to exclude extreme values (zeros or ones). Beyond the statistical sanitization benefits, one could argue that these “floor” and “ceiling” effects might correspond to a different cognitive process (this will be important in the later part of this chapter).
+
One of the main issues is that, as you can se from the histogram, there is a high number of observations clumped at zero. This is a very common and often overlooked issue in psychological science, where participants tend to use the extreme values of the scale (0 and 1) more often than the intermediate values. This creates a bimodal distribution which makes standard unimodal distributions fail to capture the data (note this issue is also present with linear models, which estimates will get biased by this configuration away from the actual “mean” of the variable).
+
One simple, although not ideal, solution could be to exclude extreme values (zeros or ones). Beyond the statistical sanitization benefits, one could argue that these “floor” and “ceiling” effects might correspond to a different cognitive process (this will be important in the latter part of this chapter).
This time, we will add another trick to make the model more robust (note that this is a general improvement that we are introducing here but that is not related to the current issue at hand of the extreme values). The current parameter mu\(\mu\) is defined on the \(]0, 1[\) range. Although this is not an issue in our model where we don’t have any predictors, these types of bounded parameters can be a bit problematic in the context of regressions, where the effect of predictors can push the parameter outside of its bounds. For example, imagine that the algorithm pics a value of \(0.45\) for mu\(\mu\) from the prior, and then picks a value of \(+0.30\) for the effect of a potential predictor (e.g., an experimental condition). This would result in a value of \(0.75\), which is outside the range of possible, and would make the model fail to converge.
-
One common solution (at the heard of the so-called logistic models) is to express mu\(\mu\) on the logit scale using the logistic function (available in the StatsFuns package). The logistic function is a simple transformation that maps any value (from the unbounded range \(]-\infty, +\infty[\))) to the 0, 1 range. We can now specify the prior for mu\(\mu\) on the logit scale as a normal distribution (e.g., \(Normal(0, 3)\)) without worrying about the estimates going out of bounds.
-
+
One common solution (at the heard of the so-called logistic models) is to express mu\(\mu\) on the logit scale, which transforms from a bounded [0, 1] range it to an unbounded \(]-\infty, +\infty[\) range. The priors on mu\(\mu\) can then be specified using any unbounded distributions (e.g., \(Normal(0, 3)\)), without worrying about the estimates going out of bounds. The parameter is then transformed back to the \(]0, 1[\) range using the logistic function (available in the StatsFuns package) before being evaluated.
+
The logistic function is a simple transformation that maps any value from the unbounded range \(]-\infty, +\infty[\)) back to the 0, 1 range.
+
Code
-
xaxis =range(-10, 10, length=1000)
-
-fig =Figure()
-ax1 =Axis(fig[1:2, 1], xlabel="Value of μ on the logit scale", ylabel="Actual value of μ", title="Logistic function")
-lines!(ax1, xaxis, logistic.(xaxis), color=:red, linewidth=2)
-ax2 =Axis(fig[1, 2], xlabel="Value of μ on the logit scale", ylabel="Plausibility", title="Prior for μ ~ Normal(0, 3)", yticksvisible=false, yticklabelsvisible=false,)
-lines!(ax2, xaxis, pdf.(Normal(0, 3), xaxis), color=:blue, linewidth=2)
-ax3 =Axis(fig[2, 2], xlabel="Value of μ after logistic transformation", ylabel="Plausibility", yticksvisible=false, yticklabelsvisible=false,)
-lines!(ax3, logistic.(xaxis), pdf.(Normal(0, 3), xaxis), color=:green, linewidth=2)
-fig
+
xaxis =range(-10, 10, length=1000)
+
+fig =Figure()
+ax1 =Axis(fig[1:2, 1], xlabel="Value of μ on the logit scale", ylabel="Actual value of μ", title="Logistic function")
+lines!(ax1, xaxis, logistic.(xaxis), color=:red, linewidth=2)
+ax2 =Axis(fig[1, 2], xlabel="Value of μ on the logit scale", ylabel="Plausibility", title="Prior for μ ~ Normal(0, 3)", yticksvisible=false, yticklabelsvisible=false,)
+lines!(ax2, xaxis, pdf.(Normal(0, 3), xaxis), color=:blue, linewidth=2)
+ax3 =Axis(fig[2, 2], xlabel="Value of μ after logistic transformation", ylabel="Plausibility", yticksvisible=false, yticklabelsvisible=false,)
+lines!(ax3, logistic.(xaxis), pdf.(Normal(0, 3), xaxis), color=:green, linewidth=2)
+fig
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
Because precision phi\(\phi\) is a positive parameter that can go from 0 to big values, with a particular threshold value at 1 (where the distribution becomes flat), it is often more convenient to express it on the log scale. This means that the value will be exponentiated before being used in the Beta distribution. Similarly to above, changing the link function does not change the model per se, but it makes it more convenient to set priors and often makes the model more efficient.
-
Expressing phi\(\phi\) on the log scale is convenient because we can set a \(Normal\) prior centered around zero, which will result (after the exponential transformation) in a distribution peaking at 1.
-
+
Similarly, because precision phi\(\phi\) is a positive parameter that can go from 0 to big values, with a particular threshold value at 1 (where the distribution becomes flat), it is often more convenient to express it on the log scale to be able to set priors on effects without worrying about potential combinations of the parameters that would result in negative values of phi\(\phi\), which would make the model fail to converge.
+
Having the parameter on the log scale simply means that its value will be exponentiated before being used in the Beta distribution (the log transformation being the inverse transformation of the exponential function). Expressing phi\(\phi\) on the log scale is convenient because we can set a \(Normal\) prior centered around zero, which will result (after the exponential transformation) in a distribution peaking at 1 (which is a good prior for the precision parameter).
+
Importantly, changing the link function of parameters does not change the model per se, but it makes it more convenient to set priors and often makes the model more efficient.
+
Code
-
xaxis =range(-10, 10, length=1000)
-
-fig =Figure()
-ax1 =Axis(fig[1:2, 1], xlabel="Value of ϕ on the log scale", ylabel="Actual value of ϕ", title="Exponential function")
-lines!(ax1, xaxis, exp.(xaxis), color=:red, linewidth=2)
-xlims!(-10, 10)
-ax2 =Axis(fig[1, 2], xlabel="Value of ϕ on the log scale", ylabel="Plausibility", title="Prior for ϕ ~ Normal(0, 2)", yticksvisible=false, yticklabelsvisible=false,)
-lines!(ax2, xaxis, pdf.(Normal(0, 2), xaxis), color=:blue, linewidth=2)
-ax3 =Axis(fig[2, 2], xlabel="Value of ϕ after exponential transformation", ylabel="Plausibility", yticksvisible=false, yticklabelsvisible=false,)
-lines!(ax3, exp.(xaxis), pdf.(Normal(0, 2), xaxis), color=:green, linewidth=2)
-xlims!(-1, 20)
-fig
+
xaxis =range(-10, 10, length=1000)
+
+fig =Figure()
+ax1 =Axis(fig[1:2, 1], xlabel="Value of ϕ on the log scale", ylabel="Actual value of ϕ", title="Exponential function")
+lines!(ax1, xaxis, exp.(xaxis), color=:red, linewidth=2)
+xlims!(-10, 10)
+ax2 =Axis(fig[1, 2], xlabel="Value of ϕ on the log scale", ylabel="Plausibility", title="Prior for ϕ ~ Normal(0, 2)", yticksvisible=false, yticklabelsvisible=false,)
+lines!(ax2, xaxis, pdf.(Normal(0, 2), xaxis), color=:blue, linewidth=2)
+ax3 =Axis(fig[2, 2], xlabel="Value of ϕ after exponential transformation", ylabel="Plausibility", yticksvisible=false, yticklabelsvisible=false,)
+lines!(ax3, exp.(xaxis), pdf.(Normal(0, 2), xaxis), color=:green, linewidth=2)
+xlims!(-1, 20)
+fig
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
+
+
+
+
+
Let us now refit the Beta model using these link functions.
+
+
@modelfunctionmodel_Beta(x)
+ μ ~Normal(0, 3) # On the logit scale
+ ϕ ~Normal(0, 2) # On the log scale
+
+for i in1:length(x)
+ x[i] ~BetaPhi2(logistic(μ), exp(ϕ))
+end
+end
+
+# Refit
+fit_Beta =model_Beta(var_noextreme)
+posteriors_Beta =sample(fit_Beta, NUTS(), 500)
+
The only caveat is that we need to transform the parameters back to the original scale when interpreting the results.
@modelfunctionmodel_Beta(x)
- μ ~Normal(0, 3) # On the logit scale
- ϕ ~Normal(0, 2) # On the log scale
-
-for i in1:length(x)
- x[i] ~BetaMuPhi(logistic(μ), exp(ϕ))
-end
-end
-
-# Refit
-fit_Beta =model_Beta(var_noextreme)
-chain_Beta =sample(fit_Beta, NUTS(), 500)
3.5.3 Conclusion
Let us now make a posterior predictive check to see how well the model fits the data.
-
+
Code
-
pred =predict(model_Beta([(missing) for i in1:length(var_noextreme)]), chain_Beta)
-pred =Array(pred)
-
-fig =hist(var_noextreme, normalization =:pdf, color=:darkred, bins=40)
-for i in1:length(chain_Beta)
-density!(pred[i, :], color=(:black, 0), strokecolor = (:black, 0.1), strokewidth =3, boundary=(0, 1))
-end
-fig
+
pred =predict(model_Beta([(missing) for i in1:length(var_noextreme)]), posteriors_Beta)
+pred =Array(pred)
+
+bins =range(0, 1, length=40)
+
+fig =hist(var_noextreme, normalization =:pdf, color=:darkred, bins=bins)
+for i in1:length(posteriors_Beta)
+hist!(pred[i, :], normalization =:pdf, bins=bins, color=(:black, 0.01))
+end
+fig
-
-
┌ Warning: Found `resolution` in the theme when creating a `Scene`. The `resolution` keyword for `Scene`s and `Figure`s has been deprecated. Use `Figure(; size = ...` or `Scene(; size = ...)` instead, which better reflects that this is a unitless size and not a pixel resolution. The key could also come from `set_theme!` calls or related theming functions.
-└ @ Makie C:\Users\domma\.julia\packages\Makie\VRavR\src\scenes.jl:220
-
-
-
+
+
-
Removing the extreme values improved the fit. However, it is not a perfect solution.
+
Although removing the extreme values did improve the fit, it is not a perfect solution, as it neglects the importance of these values and the potential cognitive processes behind them.