# Load Turing.
using Turing
# Load RDatasets.
using RDatasets
# Load StatsPlots for visualizations and diagnostics.
using StatsPlots
# Functionality for splitting and normalizing the data.
using MLDataUtils: shuffleobs, splitobs, rescale!
# We need a softmax function which is provided by NNlib.
using NNlib: softmax
# Functionality for constructing arrays with identical elements efficiently.
using FillArrays
# Functionality for working with scaled identity matrices.
using LinearAlgebra
# Set a seed for reproducibility.
using Random
Random.seed!(0);Multinomial Logistic Regression
Multinomial logistic regression is an extension of logistic regression. Logistic regression is used to model problems in which there are exactly two possible discrete outcomes. Multinomial logistic regression is used to model problems in which there are two or more possible discrete outcomes.
In our example, we’ll be using the iris dataset. The iris multiclass problem aims to predict the species of a flower given measurements (in centimeters) of sepal length and width and petal length and width. There are three possible species: Iris setosa, Iris versicolor, and Iris virginica.
To start, let’s import all the libraries we’ll need.
Data Cleaning & Set Up
Now we’re going to import our dataset. Twenty rows of the dataset are shown below so you can get a good feel for what kind of data we have.
# Import the "iris" dataset.
data = RDatasets.dataset("datasets", "iris");
# Show twenty random rows.
data[rand(1:size(data, 1), 20), :]20×5 DataFrame
| Row | SepalLength | SepalWidth | PetalLength | PetalWidth | Species |
|---|---|---|---|---|---|
| Float64 | Float64 | Float64 | Float64 | Cat… | |
| 1 | 5.0 | 2.0 | 3.5 | 1.0 | versicolor |
| 2 | 5.4 | 3.7 | 1.5 | 0.2 | setosa |
| 3 | 7.2 | 3.0 | 5.8 | 1.6 | virginica |
| 4 | 4.8 | 3.0 | 1.4 | 0.1 | setosa |
| 5 | 5.7 | 2.8 | 4.1 | 1.3 | versicolor |
| 6 | 5.1 | 3.5 | 1.4 | 0.3 | setosa |
| 7 | 5.4 | 3.9 | 1.3 | 0.4 | setosa |
| 8 | 7.6 | 3.0 | 6.6 | 2.1 | virginica |
| 9 | 5.0 | 3.5 | 1.6 | 0.6 | setosa |
| 10 | 5.0 | 3.6 | 1.4 | 0.2 | setosa |
| 11 | 5.5 | 2.4 | 3.8 | 1.1 | versicolor |
| 12 | 6.1 | 2.6 | 5.6 | 1.4 | virginica |
| 13 | 4.4 | 3.0 | 1.3 | 0.2 | setosa |
| 14 | 7.0 | 3.2 | 4.7 | 1.4 | versicolor |
| 15 | 6.1 | 2.9 | 4.7 | 1.4 | versicolor |
| 16 | 7.7 | 3.0 | 6.1 | 2.3 | virginica |
| 17 | 6.4 | 2.7 | 5.3 | 1.9 | virginica |
| 18 | 5.1 | 3.3 | 1.7 | 0.5 | setosa |
| 19 | 6.7 | 3.1 | 4.7 | 1.5 | versicolor |
| 20 | 6.2 | 2.2 | 4.5 | 1.5 | versicolor |
In this data set, the outcome Species is currently coded as a string. We convert it to a numerical value by using indices 1, 2, and 3 to indicate species setosa, versicolor, and virginica, respectively.
# Recode the `Species` column.
species = ["setosa", "versicolor", "virginica"]
data[!, :Species_index] = indexin(data[!, :Species], species)
# Show twenty random rows of the new species columns
data[rand(1:size(data, 1), 20), [:Species, :Species_index]]20×2 DataFrame
| Row | Species | Species_index |
|---|---|---|
| Cat… | Union… | |
| 1 | setosa | 1 |
| 2 | versicolor | 2 |
| 3 | versicolor | 2 |
| 4 | setosa | 1 |
| 5 | versicolor | 2 |
| 6 | versicolor | 2 |
| 7 | versicolor | 2 |
| 8 | versicolor | 2 |
| 9 | virginica | 3 |
| 10 | versicolor | 2 |
| 11 | virginica | 3 |
| 12 | setosa | 1 |
| 13 | setosa | 1 |
| 14 | versicolor | 2 |
| 15 | setosa | 1 |
| 16 | setosa | 1 |
| 17 | virginica | 3 |
| 18 | versicolor | 2 |
| 19 | virginica | 3 |
| 20 | virginica | 3 |
After we’ve done that tidying, it’s time to split our dataset into training and testing sets, and separate the features and target from the data. Additionally, we must rescale our feature variables so that they are centered around zero by subtracting each column by the mean and dividing it by the standard deviation. Without this step, Turing’s sampler will have a hard time finding a place to start searching for parameter estimates.
# Split our dataset 50%/50% into training/test sets.
trainset, testset = splitobs(shuffleobs(data), 0.5)
# Define features and target.
features = [:SepalLength, :SepalWidth, :PetalLength, :PetalWidth]
target = :Species_index
# Turing requires data in matrix and vector form.
train_features = Matrix(trainset[!, features])
test_features = Matrix(testset[!, features])
train_target = trainset[!, target]
test_target = testset[!, target]
# Standardize the features.
μ, σ = rescale!(train_features; obsdim=1)
rescale!(test_features, μ, σ; obsdim=1);Model Declaration
Finally, we can define our model logistic_regression. It is a function that takes three arguments where
xis our set of independent variables;yis the element we want to predict;σis the standard deviation we want to assume for our priors.
We select the setosa species as the baseline class (the choice does not matter). Then we create the intercepts and vectors of coefficients for the other classes against that baseline. More concretely, we create scalar intercepts intercept_versicolor and intersept_virginica and coefficient vectors coefficients_versicolor and coefficients_virginica with four coefficients each for the features SepalLength, SepalWidth, PetalLength and PetalWidth. We assume a normal distribution with mean zero and standard deviation σ as prior for each scalar parameter. We want to find the posterior distribution of these, in total ten, parameters to be able to predict the species for any given set of features.
# Bayesian multinomial logistic regression
@model function logistic_regression(x, y, σ)
n = size(x, 1)
length(y) == n ||
throw(DimensionMismatch("number of observations in `x` and `y` is not equal"))
# Priors of intercepts and coefficients.
intercept_versicolor ~ Normal(0, σ)
intercept_virginica ~ Normal(0, σ)
coefficients_versicolor ~ MvNormal(Zeros(4), σ^2 * I)
coefficients_virginica ~ MvNormal(Zeros(4), σ^2 * I)
# Compute the likelihood of the observations.
values_versicolor = intercept_versicolor .+ x * coefficients_versicolor
values_virginica = intercept_virginica .+ x * coefficients_virginica
for i in 1:n
# the 0 corresponds to the base category `setosa`
v = softmax([0, values_versicolor[i], values_virginica[i]])
y[i] ~ Categorical(v)
end
end;Sampling
Now we can run our sampler. This time we’ll use NUTS to sample from our posterior.
setprogress!(false)m = logistic_regression(train_features, train_target, 1)
chain = sample(m, NUTS(), MCMCThreads(), 1_500, 3)Chains MCMC chain (1500×22×3 Array{Float64, 3}):
Iterations = 751:1:2250
Number of chains = 3
Samples per chain = 1500
Wall duration = 16.79 seconds
Compute duration = 14.34 seconds
parameters = intercept_versicolor, intercept_virginica, coefficients_versicolor[1], coefficients_versicolor[2], coefficients_versicolor[3], coefficients_versicolor[4], coefficients_virginica[1], coefficients_virginica[2], coefficients_virginica[3], coefficients_virginica[4]
internals = lp, n_steps, is_accept, acceptance_rate, log_density, hamiltonian_energy, hamiltonian_energy_error, max_hamiltonian_energy_error, tree_depth, numerical_error, step_size, nom_step_size
Summary Statistics
parameters mean std mcse ess_bulk ess_ ⋯
Symbol Float64 Float64 Float64 Float64 Flo ⋯
intercept_versicolor 1.0308 0.5157 0.0076 4680.1947 3003. ⋯
intercept_virginica -0.4564 0.6432 0.0100 4243.4040 3175. ⋯
coefficients_versicolor[1] 1.2846 0.6520 0.0097 4574.0259 3285. ⋯
coefficients_versicolor[2] -1.3241 0.5055 0.0075 4637.6511 3229. ⋯
coefficients_versicolor[3] 0.8528 0.7392 0.0106 4830.5701 3190. ⋯
coefficients_versicolor[4] 0.3110 0.7207 0.0101 5066.8048 3663. ⋯
coefficients_virginica[1] 0.8140 0.6780 0.0096 4955.9164 3222. ⋯
coefficients_virginica[2] -0.7314 0.6249 0.0091 4728.5727 3350. ⋯
coefficients_virginica[3] 2.4231 0.8547 0.0116 5456.6428 3286. ⋯
coefficients_virginica[4] 2.7416 0.7892 0.0107 5480.2993 3474. ⋯
3 columns omitted
Quantiles
parameters 2.5% 25.0% 50.0% 75.0% 97.5% ⋯
Symbol Float64 Float64 Float64 Float64 Float64 ⋯
intercept_versicolor 0.0589 0.6791 1.0240 1.3630 2.0560 ⋯
intercept_virginica -1.7448 -0.8843 -0.4433 -0.0215 0.7928 ⋯
coefficients_versicolor[1] 0.0313 0.8388 1.2861 1.7256 2.5933 ⋯
coefficients_versicolor[2] -2.3540 -1.6512 -1.3081 -0.9761 -0.3648 ⋯
coefficients_versicolor[3] -0.6153 0.3542 0.8571 1.3639 2.2709 ⋯
coefficients_versicolor[4] -1.1017 -0.1739 0.3128 0.7976 1.7228 ⋯
coefficients_virginica[1] -0.4806 0.3611 0.7983 1.2696 2.1476 ⋯
coefficients_virginica[2] -1.9594 -1.1606 -0.7299 -0.2977 0.4748 ⋯
coefficients_virginica[3] 0.8021 1.8299 2.4184 2.9937 4.1531 ⋯
coefficients_virginica[4] 1.1762 2.2105 2.7426 3.2652 4.2867 ⋯
Sampling With Multiple Threads
The sample() call above assumes that you have at least nchains threads available in your Julia instance. If you do not, the multiple chains will run sequentially, and you may notice a warning. For more information, see the Turing documentation on sampling multiple chains.
Since we ran multiple chains, we may as well do a spot check to make sure each chain converges around similar points.
plot(chain)Looks good!
We can also use the corner function from MCMCChains to show the distributions of the various parameters of our multinomial logistic regression. The corner function requires MCMCChains and StatsPlots.
# Only plotting the first 3 coefficients due to a bug in Plots.jl
corner(
chain,
MCMCChains.namesingroup(chain, :coefficients_versicolor)[1:3];
)# Only plotting the first 3 coefficients due to a bug in Plots.jl
corner(
chain,
MCMCChains.namesingroup(chain, :coefficients_virginica)[1:3];
)