From 6d657a45f1b42c71cc6b5d43322d218e540ebc39 Mon Sep 17 00:00:00 2001 From: Jose Storopoli Date: Wed, 10 Aug 2022 18:31:59 -0300 Subject: [PATCH] Ordinal regression (#55) * added esoph dataset * renamed everything to 01-10 and added template of ordinal * added CairoMakie and AoG * changed ordering of the other content * draft ordinal-regression * added CategoricalArrays and Bijectors * final draft ordinal-regression * fix date in ordinal_reg --- Project.toml | 4 + README.md | 27 +- _layout/pgwrap.html | 25 +- _literate/{1_why_Julia.jl => 01_why_Julia.jl} | 0 .../{2_bayes_stats.jl => 02_bayes_stats.jl} | 0 _literate/{3_prob_dist.jl => 03_prob_dist.jl} | 0 _literate/{4_Turing.jl => 04_Turing.jl} | 0 _literate/{5_MCMC.jl => 05_MCMC.jl} | 0 .../{6_linear_reg.jl => 06_linear_reg.jl} | 0 .../{7_logistic_reg.jl => 07_logistic_reg.jl} | 0 _literate/08_ordinal_reg.jl | 450 ++++++++++++++++++ _literate/{8_count_reg.jl => 09_count_reg.jl} | 2 +- .../{9_robust_reg.jl => 10_robust_reg.jl} | 2 +- ...evel_models.jl => 11_multilevel_models.jl} | 0 ...1_Turing_tricks.jl => 12_Turing_tricks.jl} | 0 .../{12_epi_models.jl => 13_epi_models.jl} | 0 datasets/esoph.csv | 89 ++++ index.md | 27 +- pages/{1_why_Julia.md => 01_why_Julia.md} | 2 +- pages/{2_bayes_stats.md => 02_bayes_stats.md} | 2 +- pages/{3_prob_dist.md => 03_prob_dist.md} | 2 +- pages/{4_Turing.md => 04_Turing.md} | 2 +- pages/{5_MCMC.md => 05_MCMC.md} | 2 +- pages/{6_linear_reg.md => 06_linear_reg.md} | 2 +- .../{7_logistic_reg.md => 07_logistic_reg.md} | 2 +- pages/08_ordinal_reg.md | 9 + pages/{8_count_reg.md => 09_count_reg.md} | 2 +- pages/{9_robust_reg.md => 10_robust_reg.md} | 2 +- ...evel_models.md => 11_multilevel_models.md} | 2 +- ...1_Turing_tricks.md => 12_Turing_tricks.md} | 2 +- pages/{12_epi_models.md => 13_epi_models.md} | 2 +- 31 files changed, 607 insertions(+), 52 deletions(-) rename _literate/{1_why_Julia.jl => 01_why_Julia.jl} (100%) rename _literate/{2_bayes_stats.jl => 02_bayes_stats.jl} (100%) rename _literate/{3_prob_dist.jl => 03_prob_dist.jl} (100%) rename _literate/{4_Turing.jl => 04_Turing.jl} (100%) rename _literate/{5_MCMC.jl => 05_MCMC.jl} (100%) rename _literate/{6_linear_reg.jl => 06_linear_reg.jl} (100%) rename _literate/{7_logistic_reg.jl => 07_logistic_reg.jl} (100%) create mode 100644 _literate/08_ordinal_reg.jl rename _literate/{8_count_reg.jl => 09_count_reg.jl} (99%) rename _literate/{9_robust_reg.jl => 10_robust_reg.jl} (99%) rename _literate/{10_multilevel_models.jl => 11_multilevel_models.jl} (100%) rename _literate/{11_Turing_tricks.jl => 12_Turing_tricks.jl} (100%) rename _literate/{12_epi_models.jl => 13_epi_models.jl} (100%) create mode 100644 datasets/esoph.csv rename pages/{1_why_Julia.md => 01_why_Julia.md} (77%) rename pages/{2_bayes_stats.md => 02_bayes_stats.md} (76%) rename pages/{3_prob_dist.md => 03_prob_dist.md} (77%) rename pages/{4_Turing.md => 04_Turing.md} (78%) rename pages/{5_MCMC.md => 05_MCMC.md} (79%) rename pages/{6_linear_reg.md => 06_linear_reg.md} (76%) rename pages/{7_logistic_reg.md => 07_logistic_reg.md} (76%) create mode 100644 pages/08_ordinal_reg.md rename pages/{8_count_reg.md => 09_count_reg.md} (77%) rename pages/{9_robust_reg.md => 10_robust_reg.md} (76%) rename pages/{10_multilevel_models.md => 11_multilevel_models.md} (73%) rename pages/{11_Turing_tricks.md => 12_Turing_tricks.md} (75%) rename pages/{12_epi_models.md => 13_epi_models.md} (76%) diff --git a/Project.toml b/Project.toml index 918e0b2d..3e37b559 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,10 @@ [deps] +AlgebraOfGraphics = "cbdf2221-f076-402e-a563-3d30da359d67" BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf" +Bijectors = "76274a88-744f-5084-9051-94815aaf08c4" CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b" +CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0" +CategoricalArrays = "324d7699-5711-5eae-9e2f-1d82baa6b597" Chain = "8be319e6-bccf-4806-a6f7-6fae938471bc" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" DifferentialEquations = "0c46a032-eb83-5123-abaf-570d42b7fbaa" diff --git a/README.md b/README.md index 6495f35e..0ce5d480 100644 --- a/README.md +++ b/README.md @@ -89,18 +89,19 @@ To configure a local environment: ## Tutorials -1. [**Why Julia?**](https://storopoli.github.io/Bayesian-Julia/pages/1_why_Julia/) -2. [**What is Bayesian Statistics?**](https://storopoli.github.io/Bayesian-Julia/pages/2_bayes_stats/) -3. [**Common Probability Distributions**](https://storopoli.github.io/Bayesian-Julia/pages/3_prob_dist/) -4. [**How to use Turing**](https://storopoli.github.io/Bayesian-Julia/pages/4_Turing/) -5. [**Markov Chain Monte Carlo (MCMC)**](https://storopoli.github.io/Bayesian-Julia/pages/5_MCMC/) -6. [**Bayesian Linear Regression**](https://storopoli.github.io/Bayesian-Julia/pages/6_linear_reg/) -7. [**Bayesian Logistic Regression**](https://storopoli.github.io/Bayesian-Julia/pages/7_logistic_reg/) -8. [**Bayesian Regression with Count Data**](https://storopoli.github.io/Bayesian-Julia/pages/8_count_reg/) -9. [**Robust Bayesian Regression**](https://storopoli.github.io/Bayesian-Julia/pages/9_robust_reg/) -10. [**Multilevel Models (a.k.a. Hierarchical Models)**](https://storopoli.github.io/Bayesian-Julia/pages/10_multilevel_models/) -11. [**Computational Tricks with Turing (Non-Centered Parametrization and QR Decomposition)**](https://storopoli.github.io/Bayesian-Julia/pages/11_Turing_tricks/) -12. [**Epidemiological Models using ODE Solvers in Turing**](https://storopoli.github.io/Bayesian-Julia/pages/12_epi_models/) +1. [**Why Julia?**](https://storopoli.github.io/Bayesian-Julia/pages/01_why_Julia/) +2. [**What is Bayesian Statistics?**](https://storopoli.github.io/Bayesian-Julia/pages/02_bayes_stats/) +3. [**Common Probability Distributions**](https://storopoli.github.io/Bayesian-Julia/pages/03_prob_dist/) +4. [**How to use Turing**](https://storopoli.github.io/Bayesian-Julia/pages/04_Turing/) +5. [**Markov Chain Monte Carlo (MCMC)**](https://storopoli.github.io/Bayesian-Julia/pages/05_MCMC/) +6. [**Bayesian Linear Regression**](https://storopoli.github.io/Bayesian-Julia/pages/06_linear_reg/) +7. [**Bayesian Logistic Regression**](https://storopoli.github.io/Bayesian-Julia/pages/07_logistic_reg/) +8. [**Bayesian Ordinal Regression**](https://storopoli.github.io/Bayesian-Julia/pages/08_ordinal_reg/) +9. [**Bayesian Regression with Count Data**](https://storopoli.github.io/Bayesian-Julia/pages/09_count_reg/) +10. [**Robust Bayesian Regression**](https://storopoli.github.io/Bayesian-Julia/pages/10_robust_reg/) +11. [**Multilevel Models (a.k.a. Hierarchical Models)**](https://storopoli.github.io/Bayesian-Julia/pages/11_multilevel_models/) +12. [**Computational Tricks with Turing (Non-Centered Parametrization and QR Decomposition)**](https://storopoli.github.io/Bayesian-Julia/pages/12_Turing_tricks/) +13. [**Epidemiological Models using ODE Solvers in Turing**](https://storopoli.github.io/Bayesian-Julia/pages/13_epi_models/) ## Datasets @@ -193,4 +194,4 @@ The references are divided in **books**, **papers**, **software**, and **dataset This content is licensed under [Creative Commons Attribution-ShareAlike 4.0 Internacional](http://creativecommons.org/licenses/by-sa/4.0/). -[![CC BY-SA 4.0](https://licensebuttons.net/l/by-sa/4.0/88x31.png)](http://creativecommons.org/licenses/by-sa/4.0/) \ No newline at end of file +[![CC BY-SA 4.0](https://licensebuttons.net/l/by-sa/4.0/88x31.png)](http://creativecommons.org/licenses/by-sa/4.0/) diff --git a/_layout/pgwrap.html b/_layout/pgwrap.html index 5bcbc2f8..73ef7995 100644 --- a/_layout/pgwrap.html +++ b/_layout/pgwrap.html @@ -11,18 +11,19 @@
diff --git a/_literate/1_why_Julia.jl b/_literate/01_why_Julia.jl similarity index 100% rename from _literate/1_why_Julia.jl rename to _literate/01_why_Julia.jl diff --git a/_literate/2_bayes_stats.jl b/_literate/02_bayes_stats.jl similarity index 100% rename from _literate/2_bayes_stats.jl rename to _literate/02_bayes_stats.jl diff --git a/_literate/3_prob_dist.jl b/_literate/03_prob_dist.jl similarity index 100% rename from _literate/3_prob_dist.jl rename to _literate/03_prob_dist.jl diff --git a/_literate/4_Turing.jl b/_literate/04_Turing.jl similarity index 100% rename from _literate/4_Turing.jl rename to _literate/04_Turing.jl diff --git a/_literate/5_MCMC.jl b/_literate/05_MCMC.jl similarity index 100% rename from _literate/5_MCMC.jl rename to _literate/05_MCMC.jl diff --git a/_literate/6_linear_reg.jl b/_literate/06_linear_reg.jl similarity index 100% rename from _literate/6_linear_reg.jl rename to _literate/06_linear_reg.jl diff --git a/_literate/7_logistic_reg.jl b/_literate/07_logistic_reg.jl similarity index 100% rename from _literate/7_logistic_reg.jl rename to _literate/07_logistic_reg.jl diff --git a/_literate/08_ordinal_reg.jl b/_literate/08_ordinal_reg.jl new file mode 100644 index 00000000..d8d9fc05 --- /dev/null +++ b/_literate/08_ordinal_reg.jl @@ -0,0 +1,450 @@ +# # Bayesian Ordinal Regression + +# Leaving the universe of linear models, we start to venture into generalized linear models (GLM). The second is +# **ordinal regression**. + +# A ordinal regression behaves exactly like a linear model: it makes a prediction simply by computing a weighted +# sum of the independent variables $\mathbf{X}$ by the estimated coefficients $\boldsymbol{\beta}$, +# but now we have not only one intercept but several intercepts $\alpha_k$ for $k \in K$. + +# We use ordinal regression when our dependent variable is ordinal. +# That means it has different that have a "natural order"**. +# Most important, the distance between values is not the same. +# For example, imagine a pain score scale that goes from 1 to 10. +# The distance between 1 and 2 is different from the distance 9 to 10. +# Another example is opinion pools with their ubiquously disagree-agree range +# of plausible values. +# These are also known as Likert scale variables. +# The distance between "disagree" to "not agree or disagree" is different +# than the distance between "agree" and "strongly agree". + +# This assumption is what we call the "metric" assumption, +# also called as the equidistant assumption. +# Almost always when we model an ordinal dependent variable this assumption is violated. +# Thus, we cannot blindly employ linear regression here. + +# ## How to deal with Ordered Discrete Dependent Variable? + +# So, how we deal with ordered discrete responses in our dependent variable? +# This is similar with the previous logistic regression approach. +# We have to do a **non-linear transformation of the dependent variable**. + +# ### Cumulative Distribution Function (CDF) + +# In the case of **ordinal regression**, we need to first transform the dependent variable into a **cumulative scale**. +# We need to calculate the **cumulative distribution function** (CDF) of our dependent variable: + +# $$P(Y \leq y) = \sum^y_{i=y_{\text{min}}} P(Y = i)$$ + +# The **CDF is a monotonic increasing function** that depicts the **probability of our random variable $Y$ taking values less than a certain value**. +# In our case, the discrete ordinal case, these values can be represented as positive integers ranging from 1 to the length of possible values. +# For instance, a 6-categorical ordered response variable will have 6 possible values, and all their CDFs will lie between 0 and 1. +# Furthermore, their sum will be 1; since it represents the total probability of the variable taking any of the possible values, i.e. 100%. + +# ### Log-cumulative-odds + +# That is still not enough, we need to apply the **logit function to the CDF**: + +# $$\mathrm{logit}(x) = \mathrm{logistic}^{-1}(x) = \ln\left(\frac{x}{1 -x}\right)$$ + +# where $\ln$ is the natural logarithm function. + +# The logit is the **inverse of the logistic transformation**, +# it takes as a input any number between 0 and 1 +# (where a probability is the perfect candidate) and outputs a real number, +# which we call the **log-odds**. + +# Since we are taking the log-odds of the CDF, we can call this complete transformation as +# **log-odds of the CDF**, or **log-cumulative-odds**. + +# ### $K-1$ Intercepts + +# Now, the next question is: what do we do with the log-cumulative-odds? + +# **We need the log-cumulative-odds because it allows us to construct different intercepts for the possible values our ordinal dependent variable**. +# We create an unique intercept for each possible outcome $k \in K$. + +# Notice that the highest probable value of $X$ will always have a log-cumulative-odds of $\infty$, since for $p=1$: + +# $$\ln \frac{p}{1-p} = \ln \frac{1}{1-1} = \ln 0 = \infty$$ + +# Thus, we only need $K-1$ intercepts for a $K$ possible depedent variables' response values. +# These are known as **cut points**. + +# Each intercept implies a CDF for each value $K$. +# This allows us to **violate the equidistant assumption** absent in most ordinal variables. + +# Each intercept implies a log-cumulative-odds for each $k \in K$. +# We also need to **undo the cumulative nature of the $K-1$ intercepts**. +# We can accomplish this by first converting a **log-cumulative-odds back to a cumulative probability**. +# This is done by undoing the logit transformation and applying the logistic function: + +# $$\mathrm{logit}^{-1}(x) = \mathrm{logistic}(x) = \frac{1}{1 + e^{-x}}$$ + +# Then, finally, we can remove the **cumulative from "cumulative probability"** by +# **subtraction of each of the $k$ cut points by their previous $k-1$ cut point**: + +# $$P(Y=k) = P(Y \leq k) - P(Y \leq k-1)$$ + +# where $Y$ is the depedent variable and $k \in K$ are the cut points for each intercept. + +# Let me show you an example with some syntethic data. + +using DataFrames +using CairoMakie +using AlgebraOfGraphics +using Distributions +using StatsFuns: logit + +# Here we have a discrete variable `x` with 6 possible ordered values as response. +# The values range from 1 to 6 having probability, respectively: +# 10%, 15%, 33%, 25%, 10%, and 7%; +# represented with the `probs` vector. +# The underlying distribution is represented by a +# `Categorical` distribution from `Distributions.jl`, +# which takes a vector of probabilities as parameters (`probs`). + +# For each value we are calculating: + +# 1. **P**robability **M**ass **F**unction with the `pdf` function +# 2. **C**umulative **D**istribution **F**unction with the `cdf` function +# 3. **Log-cumulative-odds** with the `logit` transformation of the CDF + +# In the plot below there are 3 subplots: +# - Upper corner: histogram of `x` +# - Left lower corner: CDF of `x` +# - Right lower corner: log-cumulative-odds of `x` + +let + probs = [0.10, 0.15, 0.33, 0.25, 0.10, 0.07] + dist = Categorical(probs) + x = 1:length(probs) + x_pmf=pdf.(dist, x) + x_cdf=cdf.(dist, x) + x_logodds_cdf=logit.(x_cdf) + df = DataFrame(; + x, + x_pmf, + x_cdf, + x_logodds_cdf) + labels = ["CDF", "Log-cumulative-odds"] + fig = Figure() + plt1 = data(df) * + mapping(:x, :x_pmf) * + visual(BarPlot) + plt2 = data(df) * + mapping(:x, + [:x_cdf, :x_logodds_cdf]; + col=dims(1) => renamer(labels)) * + visual(ScatterLines) + axis=(; xticks=1:6) + draw!(fig[1, 2:3], plt1; axis) + draw!(fig[2, 1:4], plt2; + axis, + facet=(; linkyaxes=:none)) + fig + save(joinpath(@OUTPUT, "logodds.svg"), fig); # hide +end + +# \fig{logodds} +# \center{*Ordinal Dependent Variable*} \\ + +# As we can see, we have $K-1$ (in our case $6-1=5$) intercept values in log-cumulative-odds. +# You can carly see that these values they violate the **equidistant assumption** +# for metric response values. +# The spacing between the cut points are not the same, they vary. + +# ## Adding Coefficients $\boldsymbol{\beta}$ + +# Ok, the $K-1$ intercepts $\boldsymbol{\alpha}$ are done. +# Now let's add coefficients to act as covariate effects in our ordinal regression model. + +# We transformed everything into log-odds scale so that we can add effects +# (coefficients multiplying a covariate) or basal rates (intercepts) together. +# For each $k \in K-1$, we calculate: + +# $$\phi_k = \alpha_k + \beta_i x_i$$ + +# where $\alpha_k$ is the log-cumulative-odds for the $k \in K-1$ intercepts, +# $\beta_i$ is the coefficient for the $i$th covariate $x$. +# Finally, $\phi_k$ represents the linear predictor for the $k$th intercept. + +# Observe that the coefficient $\beta$ is being added to a log-cumulative-odds, +# such that, it will be expressed in a log-cumulative-odds also. + +# We can express it in matrix form: + +# $$\boldsymbol{\phi} = \boldsymbol{\alpha} + \mathbf{X} \cdot \boldsymbol{\beta}$$ + +# where $\boldsymbol{\phi}$, $\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$ are vectors +# and $\mathbf{X}$ is the data matrix where each row is an observation and each column a covariate. + +# This still obeys the ordered constraint on the dependent variable possible values. + +# #### How to Interpret Coefficient $\boldsymbol{\beta}$? + +# Now, suppose we have found our ideal value for our $\boldsymbol{\beta}$. +# **How we would interpret our $\boldsymbol{\beta}$ estimated values?** + +# First, to recap, $\boldsymbol{\beta}$ measures the strength of association between +# the covariate $\mathbf{x}$ and depedent variable $\mathbf{y}$. +# But, $\boldsymbol{\beta}$ is nested inside a non-linear transformation called +# logistic function: + +# $$\mathrm{logistic}(\boldsymbol{\beta}) = \frac{1}{1 + e^{-\boldsymbol{\beta}}}$$ + +# So, our first step is to **undo the logistic function**. +# There is a function that is called **logit function** that is the **inverse** +# of the logistic function: + +# $$\mathrm{logistic}^{-1}(x) = \mathrm{logit}(x) = \ln\left(\frac{x}{1 -x}\right)$$ + +# where $\ln$ is the natural logarithm function. + +# If we analyze closely the logit function we can find that +# inside the $\ln$ there is a disguised odds in the $\frac{x}{1 -x}$. +# Hence, our $\boldsymbol{\beta}$ can be interpreted as the +# **logarithm of the odds**, or in short form: the **log-odds**. + +# We already saw how odds, log-odds, and probability are related +# in the previous logistic regression tutorial. +# So you might want to go back there to get the full explanation. + +# The log-odds are the key to interpret coefficient $\boldsymbol{\beta}$**. +# Any positive value of $\beta$ means that there exists a positive association between $x$ and $y$, while any negative values indicate a negative association between $x$ and $y$. +# Values close to 0 demonstrates the lack of association between $x$ and $y$. + +# ## Likelihood + +# We have almost everything we need for our ordinal regression. +# We are only missing a final touch. +# Currently our **logistic function outputs a vector of probabilities** that sums to 1. + +# All of the intercepts $\alpha_k$ along with the coefficients $\beta_i$ are in +# log-cumulative-odds scale. +# If we apply the logistic function to the linear predictors $\phi_k$ we get $K-1$ +# probabilities: one for each $\phi_k$ + +# We need a **likelihood that can handle a vector of probabilities and outputs a single +# discrete value**. +# The **categorical distribution is the perfect candidate**. + +# ## Bayesian Ordinal Regression + +# Now we have all the components for our Bayesian ordinal regression specification: + +# $$ +# \begin{aligned} +# \mathbf{y} &\sim \text{Categorical}(\mathbf{p}) \\ +# \mathbf{p} &= \text{Logistic}(\boldsymbol{\phi}) \\ +# \boldsymbol{\phi} &= \boldsymbol{\alpha} + \mathbf{X} \cdot \boldsymbol{\beta} \\ +# \alpha_1 &= \text{CDF}(y_1) \\ +# \alpha_k &= \text{CDF}(y_k) - \text{CDF}(y_{k-1}) \quad \text{for} \quad 1 < k < K-1 \\ +# \alpha_{K-1} &= 1 - \text{CDF}(y_{K-1}) +# \end{aligned} +# $$ + +# where: + +# * $\mathbf{y}$ -- ordered discrete dependent variable. +# * $\mathbf{p}$ -- probability vector of length $K$. +# * $K$ -- number of possible values $\mathbf{y}$ can take, i.e. number of ordered discrete values. +# * $\boldsymbol{\phi}$ -- log-cumulative-odds, i.e. cut points considering the intercepts and covariates effect. +# * $\alpha_k$ -- intercept in log-cumulative-odds for each $k \in K-1$. +# * $\mathbf{X}$ -- covariate data matrix. +# * $\boldsymbol{\beta}$ -- coefficient vector of the same length as the number of columns in $\mathbf{X}$. +# * $\mathrm{logistic}$ -- logistic function. +# * $\mathrm{CDF}$ -- **c**umulative **d**istribution **f**unction. + +# What remains is to specify the model parameters' prior distributions: + +# * Prior Distribution of $\boldsymbol{\alpha}$ -- Knowledge we possess regarding the model's intercepts. +# * Prior Distribution of $\boldsymbol{\beta}$ -- Knowledge we possess regarding the model's independent variables' coefficients. + +# Our goal is to instantiate an ordinal regression with the observed data ($\mathbf{y}$ and $\mathbf{X}$) and find the posterior +# distribution of our model's parameters of interest ($\boldsymbol{\alpha}$ and $\boldsymbol{\beta}$). This means to find the full posterior +# distribution of: + +# $$ P(\boldsymbol{\theta} \mid \mathbf{y}) = P(\boldsymbol{\alpha}, \boldsymbol{\beta} \mid \mathbf{y}) $$ + +# Note that contrary to the linear regression, which used a Gaussian/normal likelihood function, we don't have an error +# parameter $\sigma$ in our ordinal regression. This is due to the Categorical distribution not having +# a "scale" parameter such as the $\sigma$ parameter in the Gaussian/normal distribution. + +# This is easily accomplished with Turing: + +using Turing +using Bijectors +using LazyArrays +using LinearAlgebra +using Random:seed! +seed!(123) +setprogress!(false) # hide + +@model function ordreg(X, y; predictors=size(X, 2), ncateg=maximum(y)) + #priors + cutpoints ~ Bijectors.ordered(filldist(TDist(3) * 5, ncateg - 1)) + β ~ filldist(TDist(3) * 2.5, predictors) + + #likelihood + y ~ arraydist([OrderedLogistic(X[i, :] ⋅ β, cutpoints) for i in 1:length(y)]) +end; + +# First, let's deal with the new stuff in our model: the **`Bijectors.ordered`**. +# As I've said in the [4. **How to use Turing**](/pages/04_Turing/), +# Turing has a rich ecossystem of packages. +# Bijectors implements a set of functions for transforming constrained random variables +# (e.g. simplexes, intervals) to Euclidean space. +# Here we are defining `cutpoints` as a `ncateg - 1` vector of Student-$t$ distributions +# with mean 0, standard deviation 5 and degrees of freedom $\nu = 3$. +# Remember that we only need $K-1$ cutpoints for all of our $K$ intercepts. +# And we are also contraining it to be an ordered vector with `Bijectors.ordered`, +# such that for all cutpoints $c_i$ we have $c_1 < c_2 < ... c_{k-1}$. + +# As before, we are giving $\boldsymbol{\beta}$ a very weakly informative priors of a +# Student-$t$ distribution centered on 0 with variance 1 and degrees of freedom $\nu = 3$. +# That wide-tailed $t$ distribution will cover all possible values for our coefficients. +# Remember the Student-$t$ also has support over all the real number line $\in (-\infty, +\infty)$. Also the `filldist()` is a nice Turing's function which takes any univariate or multivariate distribution and returns another distribution that repeats the input distribution. + +# Finally, in the likelihood, +# Turing's `arraydist()` function wraps an array of distributions returning a new distribution sampling from the individual +# distributions. +# And we use some indexing inside an array literal. + +# ## Example - Esoph + +# For our example, I will use a famous dataset called `esoph` (Breslow & Day, 1980), +# which is data from a case-control study of (o)esophageal cancer in Ille-et-Vilaine, France. +# It has records for 88 age/alcohol/tobacco combinations: + +# * `agegp`: Age group +# * `1`: 25-34 years +# * `2`: 35-44 years +# * `3`: 45-54 years +# * `4`: 55-64 years +# * `5`: 65-74 years +# * `6`: 75+ years +# * `alcgp`: Alcohol consumption +# * `1`: 0-39 g/day +# * `2`: 40-79 g/day +# * `3`: 80-119 g/day +# * `4`: 120+ g/day +# * `tobgp`: Tobacco consumption +# * `1`: 0-9 g/day +# * `2`: 10-19 g/day +# * `3`: 20-29 g/day +# * `4`: 30+ g/day +# * `ncases`: Number of cases +# * `ncontrols`: Number of controls + +# Ok let's read our data with `CSV.jl` and output into a `DataFrame` from `DataFrames.jl`: +using DataFrames, CSV, HTTP + +url = "https://raw.githubusercontent.com/storopoli/Bayesian-Julia/master/datasets/esoph.csv" +esoph = CSV.read(HTTP.get(url).body, DataFrame) +describe(esoph) + +# As you can see from the `describe()` output 58% of the respondents switched wells and 42% percent of respondents +# somehow are engaged in community organizations. The average years of education of the household's head is approximate +# 5 years and ranges from 0 (no education at all) to 17 years. The distance to safe drinking water is measured in meters +# and averages 48m ranging from less than 1m to 339m. Regarding arsenic levels I cannot comment because the only thing I +# know that it is toxic and you probably would never want to have your well contaminated with it. Here, we believe that all +# of those variables somehow influence the probability of a respondent switch to a safe well. + +# Now let's us instantiate our model with the data. +# But here I need to do some data wrangling to create the data matrix `X`. +# Specifically, I need to convert all of the categorical variables to integer values, +# representing the ordinal values of both our independent and also dependent variables: + +using CategoricalArrays + +transform!( + esoph, + :agegp => + x -> categorical( + x; levels=["25-34", "35-44", "45-54", "55-64", "65-74", "75+"], ordered=true + ), + :alcgp => + x -> categorical(x; levels=["0-39g/day", "40-79", "80-119", "120+"], ordered=true), + :tobgp => + x -> categorical(x; levels=["0-9g/day", "10-19", "20-29", "30+"], ordered=true); + renamecols=false, +) +transform!(esoph, [:agegp, :alcgp, :tobgp] .=> ByRow(levelcode); renamecols=false) + +X = Matrix(select(esoph, [:agegp, :alcgp])) +y = esoph[:, :tobgp] +model = ordreg(X, y); + +# And, finally, we will sample from the Turing model. We will be using the default `NUTS()` sampler with `2_000` samples, with +# 4 Markov chains using multiple threads `MCMCThreads()`: + +chain = sample(model, NUTS(), MCMCThreads(), 2_000, 4) +summarystats(chain) + +# We had no problem with the Markov chains as all the `rhat` are well below `1.01` (or above `0.99`). +# Note that the coefficients are in log-odds scale. They are the natural log of the odds[^logit], and odds is +# defined as: + +# $$ \text{odds} = \frac{p}{1-p} $$ + +# where $p$ is a probability. So log-odds is defined as: + +# $$ \log(\text{odds}) = \log \left( \frac{p}{1-x} \right) $$ + +# So in order to get odds from a log-odds we must undo the log operation with a exponentiation. +# This translates to: + +# $$ \text{odds} = \exp ( \log ( \text{odds} )) $$ + +# We can do this with a transformation in a `DataFrame` constructed from a `Chains` object: + +using Chain + +@chain quantile(chain) begin + DataFrame + select(_, + :parameters, + names(_, r"%") .=> ByRow(exp), + renamecols=false) +end + +# Our interpretation of odds is the same as in betting games. Anything below 1 signals a unlikely probability that $y$ will be $1$. +# And anything above 1 increases the probability of $y$ being $1$, while 1 itself is a neutral odds for $y$ being either $1$ or $0$. +# Since I am not a gambling man, let's talk about probabilities. So I will create a function called `logodds2prob()` that converts +# log-odds to probabilities: + +function logodds2prob(logodds::Float64) + return exp(logodds) / (1 + exp(logodds)) +end + +@chain quantile(chain) begin + DataFrame + select(_, + :parameters, + names(_, r"%") .=> ByRow(logodds2prob), + renamecols=false) +end + +# There you go, much better now. Let's analyze our results. +# The `cutpoints` is the basal rate of the probability of our dependent variable +# having values less than a certain value. +# For example the cutpoint for having values less than `2` which its code represents +# the tobacco comsumption of 10-19 g/day has a median value of 20%. + +# Now let's take a look at our coefficients +# All coefficients whose 95% credible intervals captures the value $\frac{1}{2} = 0.5$ tells +# that the effect on the propensity of tobacco comsumption is inconclusive. +# It is pretty much similar to a 95% credible interval that captures the 0 in +# the linear regression coefficients. + +# That's how you interpret 95% credible intervals from a `quantile()` output of a ordinal regression `Chains` +# object converted from log-odds to probability. + +# ## Footnotes + +# [^logit]: actually the [logit](https://en.wikipedia.org/wiki/Logit) function or the log-odds is the logarithm of the odds $\frac{p}{1-p}$ where $p$ is a probability. + +# ## References + +# Breslow, N. E. & Day, N. E. (1980). **Statistical Methods in Cancer Research. Volume 1: The Analysis of Case-Control Studies**. IARC Lyon / Oxford University Press. diff --git a/_literate/8_count_reg.jl b/_literate/09_count_reg.jl similarity index 99% rename from _literate/8_count_reg.jl rename to _literate/09_count_reg.jl index 19fa9323..77f3321f 100644 --- a/_literate/8_count_reg.jl +++ b/_literate/09_count_reg.jl @@ -1,6 +1,6 @@ # # Bayesian Regression with Count Data -# Leaving the universe of linear models, we start to venture into generalized linear models (GLM). The second of these is +# Leaving the universe of linear models, we start to venture into generalized linear models (GLM). The third of these is # **regression with count data** (also called Poisson regression). # A regression with count data behaves exactly like a linear model: it makes a prediction simply by computing a weighted diff --git a/_literate/9_robust_reg.jl b/_literate/10_robust_reg.jl similarity index 99% rename from _literate/9_robust_reg.jl rename to _literate/10_robust_reg.jl index f98be50e..068a8c2d 100644 --- a/_literate/9_robust_reg.jl +++ b/_literate/10_robust_reg.jl @@ -1,6 +1,6 @@ # # Robust Bayesian Regression -# Leaving the universe of linear models, we start to venture into generalized linear models (GLM). The third of these is +# Leaving the universe of linear models, we start to venture into generalized linear models (GLM). The fourth of these is # **robust regression**. # A regression with count data behaves exactly like a linear model: it makes a prediction simply by computing a weighted diff --git a/_literate/10_multilevel_models.jl b/_literate/11_multilevel_models.jl similarity index 100% rename from _literate/10_multilevel_models.jl rename to _literate/11_multilevel_models.jl diff --git a/_literate/11_Turing_tricks.jl b/_literate/12_Turing_tricks.jl similarity index 100% rename from _literate/11_Turing_tricks.jl rename to _literate/12_Turing_tricks.jl diff --git a/_literate/12_epi_models.jl b/_literate/13_epi_models.jl similarity index 100% rename from _literate/12_epi_models.jl rename to _literate/13_epi_models.jl diff --git a/datasets/esoph.csv b/datasets/esoph.csv new file mode 100644 index 00000000..7e43ff57 --- /dev/null +++ b/datasets/esoph.csv @@ -0,0 +1,89 @@ +"agegp","alcgp","tobgp","ncases","ncontrols" +"25-34","0-39g/day","0-9g/day",0,40 +"25-34","0-39g/day","10-19",0,10 +"25-34","0-39g/day","20-29",0,6 +"25-34","0-39g/day","30+",0,5 +"25-34","40-79","0-9g/day",0,27 +"25-34","40-79","10-19",0,7 +"25-34","40-79","20-29",0,4 +"25-34","40-79","30+",0,7 +"25-34","80-119","0-9g/day",0,2 +"25-34","80-119","10-19",0,1 +"25-34","80-119","30+",0,2 +"25-34","120+","0-9g/day",0,1 +"25-34","120+","10-19",1,0 +"25-34","120+","20-29",0,1 +"25-34","120+","30+",0,2 +"35-44","0-39g/day","0-9g/day",0,60 +"35-44","0-39g/day","10-19",1,13 +"35-44","0-39g/day","20-29",0,7 +"35-44","0-39g/day","30+",0,8 +"35-44","40-79","0-9g/day",0,35 +"35-44","40-79","10-19",3,20 +"35-44","40-79","20-29",1,13 +"35-44","40-79","30+",0,8 +"35-44","80-119","0-9g/day",0,11 +"35-44","80-119","10-19",0,6 +"35-44","80-119","20-29",0,2 +"35-44","80-119","30+",0,1 +"35-44","120+","0-9g/day",2,1 +"35-44","120+","10-19",0,3 +"35-44","120+","20-29",2,2 +"45-54","0-39g/day","0-9g/day",1,45 +"45-54","0-39g/day","10-19",0,18 +"45-54","0-39g/day","20-29",0,10 +"45-54","0-39g/day","30+",0,4 +"45-54","40-79","0-9g/day",6,32 +"45-54","40-79","10-19",4,17 +"45-54","40-79","20-29",5,10 +"45-54","40-79","30+",5,2 +"45-54","80-119","0-9g/day",3,13 +"45-54","80-119","10-19",6,8 +"45-54","80-119","20-29",1,4 +"45-54","80-119","30+",2,2 +"45-54","120+","0-9g/day",4,0 +"45-54","120+","10-19",3,1 +"45-54","120+","20-29",2,1 +"45-54","120+","30+",4,0 +"55-64","0-39g/day","0-9g/day",2,47 +"55-64","0-39g/day","10-19",3,19 +"55-64","0-39g/day","20-29",3,9 +"55-64","0-39g/day","30+",4,2 +"55-64","40-79","0-9g/day",9,31 +"55-64","40-79","10-19",6,15 +"55-64","40-79","20-29",4,13 +"55-64","40-79","30+",3,3 +"55-64","80-119","0-9g/day",9,9 +"55-64","80-119","10-19",8,7 +"55-64","80-119","20-29",3,3 +"55-64","80-119","30+",4,0 +"55-64","120+","0-9g/day",5,5 +"55-64","120+","10-19",6,1 +"55-64","120+","20-29",2,1 +"55-64","120+","30+",5,1 +"65-74","0-39g/day","0-9g/day",5,43 +"65-74","0-39g/day","10-19",4,10 +"65-74","0-39g/day","20-29",2,5 +"65-74","0-39g/day","30+",0,2 +"65-74","40-79","0-9g/day",17,17 +"65-74","40-79","10-19",3,7 +"65-74","40-79","20-29",5,4 +"65-74","80-119","0-9g/day",6,7 +"65-74","80-119","10-19",4,8 +"65-74","80-119","20-29",2,1 +"65-74","80-119","30+",1,0 +"65-74","120+","0-9g/day",3,1 +"65-74","120+","10-19",1,1 +"65-74","120+","20-29",1,0 +"65-74","120+","30+",1,0 +"75+","0-39g/day","0-9g/day",1,17 +"75+","0-39g/day","10-19",2,4 +"75+","0-39g/day","30+",1,2 +"75+","40-79","0-9g/day",2,3 +"75+","40-79","10-19",1,2 +"75+","40-79","20-29",0,3 +"75+","40-79","30+",1,0 +"75+","80-119","0-9g/day",1,0 +"75+","80-119","10-19",1,0 +"75+","120+","0-9g/day",2,0 +"75+","120+","10-19",1,0 diff --git a/index.md b/index.md index 055ff346..5cdbd040 100644 --- a/index.md +++ b/index.md @@ -70,18 +70,19 @@ Pkg.instantiate() ## Tutorials -1. [**Why Julia?**](/pages/1_why_Julia/) -2. [**What is Bayesian Statistics?**](/pages/2_bayes_stats/) -3. [**Common Probability Distributions**](/pages/3_prob_dist/) -4. [**How to use Turing**](/pages/4_Turing/) -5. [**Markov Chain Monte Carlo (MCMC)**](/pages/5_MCMC/) -6. [**Bayesian Linear Regression**](/pages/6_linear_reg/) -7. [**Bayesian Logistic Regression**](/pages/7_logistic_reg/) -8. [**Bayesian Regression with Count Data**](/pages/8_count_reg/) -9. [**Robust Bayesian Regression**](/pages/9_robust_reg/) -10. [**Multilevel Models (a.k.a. Hierarchical Models)**](/pages/10_multilevel_models/) -11. [**Computational Tricks with Turing (Non-Centered Parametrization and QR Decomposition)**](/pages/11_Turing_tricks/) -12. [**Epidemiological Models using ODE Solvers in Turing**](/pages/12_epi_models/) +1. [**Why Julia?**](/pages/01_why_Julia/) +2. [**What is Bayesian Statistics?**](/pages/02_bayes_stats/) +3. [**Common Probability Distributions**](/pages/03_prob_dist/) +4. [**How to use Turing**](/pages/04_Turing/) +5. [**Markov Chain Monte Carlo (MCMC)**](/pages/05_MCMC/) +6. [**Bayesian Linear Regression**](/pages/06_linear_reg/) +7. [**Bayesian Logistic Regression**](/pages/07_logistic_reg/) +8. [**Bayesian Ordinal Regression**](/pages/08_ordinal_reg/) +9. [**Bayesian Regression with Count Data**](/pages/09_count_reg/) +10. [**Robust Bayesian Regression**](/pages/10_robust_reg/) +11. [**Multilevel Models (a.k.a. Hierarchical Models)**](/pages/11_multilevel_models/) +12. [**Computational Tricks with Turing (Non-Centered Parametrization and QR Decomposition)**](/pages/12_Turing_tricks/) +13. [**Epidemiological Models using ODE Solvers in Turing**](/pages/13_epi_models/) ## Datasets @@ -199,4 +200,4 @@ list = ["$(p.name) $(p.version)" for p in deps] sort!(list) println(join(list, '\n')) ``` -\output{packages} \ No newline at end of file +\output{packages} diff --git a/pages/1_why_Julia.md b/pages/01_why_Julia.md similarity index 77% rename from pages/1_why_Julia.md rename to pages/01_why_Julia.md index 54680ea0..60c98848 100644 --- a/pages/1_why_Julia.md +++ b/pages/01_why_Julia.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/1_why_Julia.jl} +\literate{/_literate/01_why_Julia.jl} diff --git a/pages/2_bayes_stats.md b/pages/02_bayes_stats.md similarity index 76% rename from pages/2_bayes_stats.md rename to pages/02_bayes_stats.md index 95306693..df073c9d 100644 --- a/pages/2_bayes_stats.md +++ b/pages/02_bayes_stats.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/2_bayes_stats.jl} +\literate{/_literate/02_bayes_stats.jl} diff --git a/pages/3_prob_dist.md b/pages/03_prob_dist.md similarity index 77% rename from pages/3_prob_dist.md rename to pages/03_prob_dist.md index 688ad5e0..1b2de11d 100644 --- a/pages/3_prob_dist.md +++ b/pages/03_prob_dist.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/3_prob_dist.jl} +\literate{/_literate/03_prob_dist.jl} diff --git a/pages/4_Turing.md b/pages/04_Turing.md similarity index 78% rename from pages/4_Turing.md rename to pages/04_Turing.md index bf19202d..cc9d8e95 100644 --- a/pages/4_Turing.md +++ b/pages/04_Turing.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/4_Turing.jl} +\literate{/_literate/04_Turing.jl} diff --git a/pages/5_MCMC.md b/pages/05_MCMC.md similarity index 79% rename from pages/5_MCMC.md rename to pages/05_MCMC.md index 6f0d0d60..1ef807b4 100644 --- a/pages/5_MCMC.md +++ b/pages/05_MCMC.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/5_MCMC.jl} +\literate{/_literate/05_MCMC.jl} diff --git a/pages/6_linear_reg.md b/pages/06_linear_reg.md similarity index 76% rename from pages/6_linear_reg.md rename to pages/06_linear_reg.md index 970dec17..97ad8293 100644 --- a/pages/6_linear_reg.md +++ b/pages/06_linear_reg.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/6_linear_reg.jl} +\literate{/_literate/06_linear_reg.jl} diff --git a/pages/7_logistic_reg.md b/pages/07_logistic_reg.md similarity index 76% rename from pages/7_logistic_reg.md rename to pages/07_logistic_reg.md index a365b89a..ba63b681 100644 --- a/pages/7_logistic_reg.md +++ b/pages/07_logistic_reg.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/7_logistic_reg.jl} +\literate{/_literate/07_logistic_reg.jl} diff --git a/pages/08_ordinal_reg.md b/pages/08_ordinal_reg.md new file mode 100644 index 00000000..c611c27f --- /dev/null +++ b/pages/08_ordinal_reg.md @@ -0,0 +1,9 @@ +@def tags = ["syntax", "code"] +@def date = Date("2022-08-10") +@def showall = true +@def hascode = true +@def mintoclevel = 2 + +\toc + +\literate{/_literate/08_ordinal_reg.jl} diff --git a/pages/8_count_reg.md b/pages/09_count_reg.md similarity index 77% rename from pages/8_count_reg.md rename to pages/09_count_reg.md index e3fac76b..bc200273 100644 --- a/pages/8_count_reg.md +++ b/pages/09_count_reg.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/8_count_reg.jl} +\literate{/_literate/09_count_reg.jl} diff --git a/pages/9_robust_reg.md b/pages/10_robust_reg.md similarity index 76% rename from pages/9_robust_reg.md rename to pages/10_robust_reg.md index 200ae586..5cbeedce 100644 --- a/pages/9_robust_reg.md +++ b/pages/10_robust_reg.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/9_robust_reg.jl} +\literate{/_literate/10_robust_reg.jl} diff --git a/pages/10_multilevel_models.md b/pages/11_multilevel_models.md similarity index 73% rename from pages/10_multilevel_models.md rename to pages/11_multilevel_models.md index 77b9d403..94654b4f 100644 --- a/pages/10_multilevel_models.md +++ b/pages/11_multilevel_models.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/10_multilevel_models.jl} +\literate{/_literate/11_multilevel_models.jl} diff --git a/pages/11_Turing_tricks.md b/pages/12_Turing_tricks.md similarity index 75% rename from pages/11_Turing_tricks.md rename to pages/12_Turing_tricks.md index b7fd6742..738836a5 100644 --- a/pages/11_Turing_tricks.md +++ b/pages/12_Turing_tricks.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/11_Turing_tricks.jl} +\literate{/_literate/12_Turing_tricks.jl} diff --git a/pages/12_epi_models.md b/pages/13_epi_models.md similarity index 76% rename from pages/12_epi_models.md rename to pages/13_epi_models.md index 045066eb..4d89c7e1 100644 --- a/pages/12_epi_models.md +++ b/pages/13_epi_models.md @@ -6,4 +6,4 @@ \toc -\literate{/_literate/12_epi_models.jl} +\literate{/_literate/13_epi_models.jl}