Generative model to classify parkinson subject

Generative model is part of machine learning, and is usually used when you have easy access to your data and you want to infer the distribution on it. We will take a look at one generative model widely known : Naive Bayes.
Thomas Bayes is a well known mathematician, with work focused in probability. Bayes theory [1] is central in machine learning and has numerous applications in automatic classification.

tl;dr¶

1. Generated data from PPMI cohort initiative
2. Choice of model and parameter estimation on the training set (testing set to check accuracy).
3. Evaluation of the likelihood for all the classes on the validation set.
4. Posterior estimation and MAP for class prediction.
5. Classsifier accuracy study

1. Theory¶

Consider two event $a$ and $b$, the Bayesian rule define the probability $p$ that $a$ occurs, knowing the result of the event $b$ :

It can be interpreted as is :
if we know a priori the probabilities of the events_ $a$, $b$, and the result of the event $b$ (knowing result of $a$), then we can deduct the result $a$ (knowing result of $b$).

This rule will allow us to infer things about the income of events, with some aprioris. The idea behing generative models is to learn the probability $p(b,a) = p(a) . p(b|a)$, which will later be used to classify the data $a$ with $p(a|b)$ through bayesian rule.

2. Classifying parkinson patient with naive bayes classifier¶

2.1 Parkinson disease biomarkers¶

Parkinson disorder is a neurodegenerative disease, characterized by the loss of cells creating the dopamine . Dopamine is an organic chemical ensuring the communication between the neurons in the brain. Early diagnosis of parkinson is crucial because to initiate neuroprotective therapies and enhancing their success rate. Here we will use biomarkers using single-photon emission computed tomography (SPECT) from one region of the brain (striatal binding ratio from right putamen). Then, we will try to classify healthy subject and parkinson, following the paper from prashanth et al.

Given the SBR of the righ putamen $\mathbf{x}$, and the outcome $\Omega$ for a patient to have parkinson $\omega_p$ or being healthy $\omega_h$, $\omega_i = \{\omega_p, \omega_h\}$, the bayesian rule can be rewritten as : \begin{equation} P(\Omega = \omega_i\mid \mathbf{x}) = \frac{p(\mathbf{x}\mid\Omega = \omega_i) \times p(\Omega = \omega_i)}{p(\mathbf{x})} \end{equation}

$p(\omega_i)$ and $p(X)$ are the marginal probabilities, and $p(\mathbf{x}|\omega_i).p(\omega_i)$ is known as the likelihood (refer to post \url{}) which is used to create the generative model.

So let's load and visualize the data, for the rest of the post we suppose that the sample $\mathbf{x}$ are independent and identically distributed. To be in agreement with the data use agreement, "fake" data will be sampled from a model I generated using the trial see appendix A for more details.

2.2 Data modelisation¶

Given the data, you want to model the distribution of the two conditions. In machine learning, we usually refer the probability as a function called pdf (probability density function). Because we hypothetized that the data are coming from a gaussian, we can calculate its model parameters ($\mu$, $\sigma$) (usually you want to calculate the likelihood of your data to verify if it's truly gaussian, or use a test like Student's t-test). To estimate the parameters, we will divide the dataset into ($1/3$) training set (that we will use to tune our model), and ($2/3$) validation set (used later to evaluate accuracy).

2.3 Prior estimations (with learning set)¶

Now that we have our models $p(x_i | \omega_i)$, we can evaluate the likelihood of the validation data for our classes and make the predictions $p(x_i | \Omega)$. So we just fix the data and evaluate all the possible probabilities for the two classes. Usually it is not recommended to use the learning set to check the prediction accuracy but instead a testing set, this to avoid overfitting your model. Here for simplicity of this post, we will just use the training set.

The apriori probabilities will just be the number of parkinson subjects $n_{\omega_0}=3000$ divided by the total number of samples $n$. The marginal probability $p(X)$ is just the probability to have a data (for example the probability of a successfull scanning), we will assume that all subjects were successfully scanned.

Let's try with a SBR right putamen of 1.3.

A common error is to assume that the likelihood is a pdf but it is not the case, the total area is not equal to 1. You can see that for $\mathbf{x}=1.3$, the likelihoods are really close so you cannot really make any prediction there.

2.4 Posterior evaluation (with validation set)¶

Now that we calculated the likelihood evaluated for the data, we can evaluate the posterior probability (predict what condition the subject has), given the marginals and the likelihood (sometimes, it can be simpler and faster to evaluate just the likelihoods to do predictions).

After we calculated the posterior for every data, the MAP (maximum aposteriori prediction) can be used to find to which class our measure comes from. If $p(\omega_p|\mathbf{x}) > p(\omega_h|\mathbf{x})$ that means that the subject is more likely to have early parkinson.

2.5 Precision of the classifier¶

There are many tools to compare the precision of a classifier with another. I will focus on two widely used : the confusion matrix and the ROC curve.

The confusion matrix is usefull to compare how our classifier react to one data or another. It can be used for many classes but in our case we have just two conditions, either the subject has early parkinson $\omega_i = 1$, or is healthy $\omega_i = 0$.

Another usefull tool is called receiver operating characteristic (ROC) curve, it was used in during World War II to analyse the ability of radars to detect ennemies.

The idea is to visualize the true positive rate (TPR), and false positive rate (FNR) for some specific data points. A perfect is an horizontal line passing through (0,1) to (1,1). For this example, we will take uniformly cut points inside all the validation data. Then we can calculate step by step the TPR and FPR, by taking all the data inferior to the specific cut-point. In our case, a true positive is an early parkinson subject predicted as parkinson, a false positive is an healthy patient predicted as parkinson. If you want more details on ROC curve, check this : https://twitter.com/cecilejanssens/status/1104134423673479169/photo/1

I advice you to use ROC curve implementation from scikit learn, which is more precise, and less prone to floating point arithmetic error.

Conclusion¶

We can see that the bayesian classifier can perform quite well to predict healthy patient (precision healthy > 95%) but there are still improvements to be made to detect early parkinson condition. Maybe the feature is not enough separable ? Of course the biggest challenge for this type of classifiers is the need to know the distribution of the data. That's why many works simply bypass the data feature extraction and modelling (because the data is too complicated or highly-dimensionnal) by using a complex classifier and tuning its parameters (like deep NN, random forest...) or non-parametric generative model (gaussian processes, kernel density estimation..). This with the hope that the high dimensionnality of the model can help to learn the latent distribution of the data itself.

Appendix A¶

Because I am not allowed to share the data from the PPMI cohort, we will just estimate the probability density function representing the two conditions. We could fit a simple non-parametric generative model with gaussian kernel estimation on the SBR on right putamen, using scikit-learn (https://scikit-learn.org/stable/modules/density.html). After estimating the distribution of the data, we can sample the data from these two distributions as above.

import pickle
import numpy as np
import matplotlib.pyplot as plt
from sklearn.neighbors import KernelDensity

# fix random state
seed = 0
np.random.seed(seed)

# Loading of data from PPMI database
with open('YOUR_DIR/parkinson_data.pickle', 'rb') as handle:
    data = pickle.load(handle)

nb_classes = len(np.unique(data["label"]))
n_samples = len(data["label"])
idx = np.random.permutation(n_samples)
right_putamen = data["SBR"][idx, 2].reshape(-1, 1)
labels = data["label"][idx] == 1

model_parkinson = KernelDensity(kernel='gaussian', bandwidth=0.1).fit(right_putamen[labels])
model_healthy = KernelDensity(kernel='gaussian', bandwidth=0.1).fit(right_putamen[~labels])
data = {"parkinson":model_parkinson, "healthy":model_healthy}
with open('./data/parkinson/parkinson_data.pickle', 'wb') as handle:
    pickle.dump(data, handle, protocol=pickle.HIGHEST_PROTOCOL)

inputs_parkinson = np.linspace(np.min(right_putamen[labels]), np.max(right_putamen[labels]), 1000).reshape(-1, 1)
inputs_healthy = np.linspace(np.min(right_putamen[~labels]), np.max(right_putamen[~labels]), 1000).reshape(-1, 1)
pdf_parkinson = model_parkinson.score_samples(inputs_parkinson)
pdf_healthy = model_healthy.score_samples(inputs_healthy)

plt.figure()
plt.scatter(inputs_parkinson, pdf_parkinson, marker='+' ,color='darkred', linewidth=2, label='Parkinson')
plt.scatter(inputs_healthy, pdf_healthy, marker='+' ,color='darkgreen', linewidth=2, label='Healthy')
plt.legend()
plt.title('pdf SBR')
plt.xlabel('SBR right putamen')
plt.ylabel('pdf')
plt.show()

References ¶