---
jupytext:
formats: ipynb,md:myst
text_representation:
extension: .md
format_name: myst
format_version: 0.13
jupytext_version: 1.11.5
kernelspec:
display_name: Python 3
language: python
name: python3
---
```{code-cell} ipython3
:tags: [remove-cell]
from ndslib.config import jupyter_startup
jupyter_startup()
```
(scikit-learn)=
# The scikit-learn package
Now that we have a grasp on some of the key concepts, we can start *doing* machine learning in Python. In this section, we'll introduce *scikit-learn* (often abbreviated `sklearn`), which is the primary package we'll be working with throughout this chapter. Scikit-learn is the most widely-used machine learning package in Python—and for that matter, probably in any programming language. Its popularity stems from its simple, elegant interface, [stellar documentation](https://scikit-learn.org/stable/documentation.html), and comprehensive support for many of the most widely used machine learning algorithms (the main domain scikit-learn doesn't cover is deep learning, which we'll discuss separately in {numref}`dl`). Scikit-learn provides well-organized, high-quality tools for virtually all aspects of the typical machine learning workflow, including data loading and preprocessing, feature extraction and feature selection, dimensionality reduction, model selection and evaluation, and so on. We'll touch on quite a few of these as we go along.
```{eval-rst}
.. index::
single: Scikit Learn
```
## The ABIDE-II dataset
To illustrate how scikit-learn works, we're going to need some data. Scikit-learn is built on top of numpy, so in theory, we could use numpy's random number generation routines to create suitable arrays for scikit-learn, just as we did earlier. But that would be kind of boring. We already understand the basics of numpy arrays at this point, so we can be a bit more ambitious here, and try to learn machine learning using a real neuroimaging dataset. For most of this chapter, we'll use a dataset drawn from the *Autism Brain Imaging Data Exchange II* (ABIDE II) project. ABIDE is an international consortium aimed at facilitating the study of autism spectrum disorder (ASD) by publicly releasing large-scale collections of structural neuroimaging data obtained from thousands of participants at dozens of research sites {cite}`di2017enhancing`. In this chapter, we'll use data from the second collection (hence the II in ABIDE II). To keep things simple, we're going to use a lightly preprocessed version of the ABIDE II dataset.
```{eval-rst}
.. index::
single: ABIDE-II
```
```{code-cell} ipython3
from ndslib.data import load_data
abide_data = load_data("abide2")
```
We can get a quick sense of the dataset's dimensions:
```{code-cell} ipython3
abide_data.shape
```
And here are the first five rows:
```{code-cell} ipython3
display(abide_data.head(5))
```
We can see that each row contains data for a single participant and each column represents a different variable. The first 5 columns contain key identifiers and phenotypic variables: respectively, these include the research `site` the subject comes from (our dataset combines data from 17 different sites); the `subject`'s unique ID, and their `age` (in years), `sex`, and diagnosis `group` (where 1 indicates an autism diagnosis and 2 indicates a control subject). The sixth column is not in the original dataset, and is a residualized version of the age variable that we'll talk more about later.s
The remaining 1440 columns contain 4 sets of structural brain features extracted using the widely-used FreeSurfer package {cite}`fischl2012freesurfer`. Respectively, these include measures of surface area (`fsArea`), volume (`fsVol`), cortical thickness (`fsCT`), and local gyrification (`fsLGI`). Each set of FreeSurfer features contains 360 variables, reflecting the 360 regions-of-interest in the Human Connectome Project's multi-modal parcellation {cite}`glasser2016mmp1`, which looks roughly like this (each color is a different ROI):
![](./figures/hcp_mmp1.png)
### Data in scikit-learn
Now that we have a suitable dataset, we need to get it to play well with scikit-learn. It turns out we don't have to do much! As we mentioned earlier, scikit-learn is built on top of numpy and expects numpy arrays as its inputs, but our `abide_data` variable is a pandas `DataFrame`. Fortunately, pandas `DataFrame` objects are built on top of the numpy array object, so it turn out that passing our data to scikit-learn is straightforward. The main thing we need to understand is how scikit-learn expects its inputs to be structured, and verify that our dataset respects those expectations.
#### Feature data: the `X` matrix
The most important thing to know is that all model-fitting routines in scikit-learn expect to receive a 2-dimensional numpy array—conventionally named `X`—as their mandatory first input. `X` is expected to contains samples (i.e., independent observations) in rows and features in columns (this should remind you of the "tidy data" mentioned in {numref}`pandas`). For example, suppose we want to build a machine learning model that predicts a person's chronological age from their structural brain data. Then our `X` data will be expected to contain each participant's data on a separate row, and the columns of that row will be the values extracted for the brain features.
```{eval-rst}
.. index::
single: Tidy Data
```
You'll probably notice that the expected structure of `X` is almost exactly the same as the structure of our ABIDE-II dataset. There too, we had participants in rows and features in columns. How convenient! The only thing that might concern us is that our ABIDE-II dataset currently contains not only the brain features extracted with Freesurfer, but also a few identifier/phenotype columns that we probably wouldn't want to pass to scikit-learn as part of the `X` matrix (e.g., because it wouldn't make sense to try to predict age from subject ID, and because `site` values are strings, whereas scikit-learn expects `X` to contain only numbers). For the sake of clarity, then, let's break up our `abide_data` dataset into two separate data frames: one containing brain variables (we'll call it `features`), and one containing phenotypic information (`phenotypes`). This isn't strictly necessary, as we can always select only a subset of columns from `abide_data` to use as `X`, but it will help keep our code clearer and tidier below. The `DataFrame.filter` method allows us to select variables by name. In this case, we put all features that match "fs"--which stands for freesurfer--in the `features` DataFrame:
```{code-cell} ipython3
features = abide_data.filter(like='fs')
```
And remember integer-based location indexing in pandas? Here we grab the first 6 columns:
```{code-cell} ipython3
phenotypes = abide_data.iloc[:, :6]
```
Let's verify that our `phenotypes` DataFrame now contains only the 5 phenotypic columns:
```{code-cell} ipython3
phenotypes.head(5)
```
#### Labels: the `y` vector
For unsupervised learning applications in scikit-learn, the feature data in `X` are all we need. But for supervised applications (i.e., classification and regression), where we're using `X` to try and recover some known ground truth, the features aren't enough; we also need labels—conventionally labeled `y`. Scikit-learn expects `y` to be a 1-dimensional array (or vector).
The variables we use as labels will vary across our examples, but let's start by assigning age to `y`. As mentioned above, Pandas objects ar built on top of a numpy array, and we can access this array by referring to the `values` attribute of the `Series` object:
```{code-cell} ipython3
y = phenotypes['age'].values
```
You can verify that looks right:
```{code-cell} ipython3
print(y)
```
## Regression example: brain-age prediction
Now we're ready to train some machine learning models! Let's start with a regression example. We'll use scikit-learn to try to predict the measured chronological age of subjects in our dataset from variation in their brain features. Predicting the biological age of a person from their neuroimaging data is actually a rather common use-case of machine learning. This is not just a neat party trick, though. There are reasons to think that a prediction of what is called "brain-age" could also be scientifically useful. That's because certain conditions are associated with specific kinds of errors in brain-age prediction. For further reading on this topic, refer to the papers mentioned below in {numref}`ml-sklearn-addtl-resources`.
### Estimators in scikit-learn: Basic usage
One of scikit-learn's most attractive features is its simple, unified interface for configuring estimators and fitting models. A good deal of the package consists of a very large set of `Estimator` classes you can use to do various forms of machine learning. In scikit-learn, an `Estimator` does exactly what the word *estimator* normally means in statistics: it implements a rule for estimating some parameter(s) of interest from data. In much of this chapter, we'll use the terms *model* and *estimator* interchangeably, though it would be more technically accurate to say that an estimator relies internally on a particular model to generate its estimates.
While scikit-learn includes hundreds of different estimators that generate predictions in very different ways, they all share a common *application programming interface* (API)—meaning, the user interacts with them in a very similar way. In particular, every estimator class implements a `fit()` method. When we call an estimator's `fit()` method, we're telling the estimator to take some training data and do something with it. For supervised learning estimators, which we'll focus on first, we have to pass two arguments to `fit()`: an `X` matrix containing the feature data we want to use to generate predictions, and a `y` vector containing the true labels or scores. The goal of the training process is to try to recover the labels in `y` from the feature data in `X`.
Once training is complete, we can call any supervised estimator's `predict()` method, which takes an `X` matrix as input and generates corresponding predictions for the `y` scores. The `X` matrix we provide to `predict()` can be exactly the same as the one we used to `fit()` the estimator, though for reasons that will become clear later, doing that is often a spectacularly bad idea.
Conceptually, we can think of `fit()` and `predict()` as mapping onto distinct training and application phases: in the training phase, our model learns how to make predictions, and in the application phase, we deterministically use the information our model has learned to actually make predictions.
```{eval-rst}
.. index::
single: Application Programming Interface
```
### Applying the `LinearRegression()` estimator
To see how this works in practice, let's try out one particular estimator in scikit-learn: "ordinary" least-squares (OLS) regression. We'll start out small, using just a handful of brain features (rather than all 1,440) to generate a prediction. Let's sample, say, 5 features at random from the full set we've stored in `features`.
```{eval-rst}
.. index::
single: Ordinary least-squares regression
```
We'll use the fact that pandas `DataFrame` objects have a helpful sample method
for randomly sampling rows/columns. Passing a `random_state` argument allows us
to pass a fixed random seed, giving us deterministic results. This is very is
helpful when you're writing a book and prefer the reader's results not to change
every time they re-run your code.
```{code-cell} ipython3
n_features = 5
X = features.sample(n_features, axis=1, random_state=100)
```
Next, we import and initialize our linear regression estimator. Scikit-learn organizes estimators and other tools into modules based on what they're used for and/or the kind of model they represent. Linear models, including `LinearRegression`, are generally found in the `sklearn.linear_model` module.
```{code-cell} ipython3
from sklearn.linear_model import LinearRegression
```
An important principle of Scikit Learn is that initialization of an estimator instance does not require any data to be passed. At initialization, only configuration parameters are passed. In particular, the `LinearRegression` estimator has relatively few configurable parameters, and for our purposes, the default settings are all fine (e.g., by default, `fit_intercept=True`, so scikit-learn will automatically add an intercept column to our predictor matrix). This means that we initialize this object without any inputs:
```{code-cell} ipython3
model = LinearRegression()
```
#### Fitting the model
Now we're ready to fit our data! As noted above, we do this by calling the `.fit()` method. This is true for every `Estimator` in scikit-learn.
```{code-cell} ipython3
model.fit(X, y)
```
Once we execute the above line, we have a fitted model. One thing we can do at this point is examine the estimated model parameters. The sklearn convention is to denote fitted parameters with a trailing underscore.
```{code-cell} ipython3
print("Estimated intercept:", model.intercept_)
print("Estimated coefficients:", model.coef_.round(2))
```
The coefficients vary dramatically in size. This isn't because 2 of them are much more important than the other 3; it's because the 4 sets of FreeSurfer features are on very different scales (e.g., the surface area features have much larger values than the cortical thickness features). Later in the tutorial, we'll explicitly standardize the features so they're all on the same scale. But for now it makes no difference, as the predictions we get out of a `LinearRegression` estimator are insensitive to scale.
#### Generating predictions
If we want to, we can use the parameter estimates extracted above to explicitly compute predictions for new observations. I.e., we would effectively be applying the following prediction equation (note that we're rounding to 2 digits for convenience; hence, the 4th variable drops out, because its coefficient is very close to 0):
$\hat{y} = 53.81 - 3.94x_1 + 0.02x_2 - 8.92x_3 - 0.01x_5$
But we don't need to do this manually; we can easily generate predicted scores for new `X` data using the `.predict()` method that all supervised scikit-learn estimators implement. For example, here are the predicted scores for our original `X` data:
```{code-cell} ipython3
y_pred = model.predict(X)
print(y_pred)
```
Let's visualize the predicted scores against the true scores. We'll make use of the `jointplot` plotting function in the Seaborn library we introduced in {numref}`viz`.
```{code-cell} ipython3
import seaborn as sns
g = sns.jointplot(x=y_pred, y=y, kind="reg").set_axis_labels("age (predicted)", "age (true)")
```
Two things jump out at us from this plot. First, the model appears to be capturing some signal in the data, in that the predicted age values clearly bear some non-trivial relationship to the true values. Second, the relationship appears to be non-linear, suggesting that our model might be misspecified. There are various ways we could potentially address this (as an exercise, try log-transforming age and repeating the estimation), but we won't worry about that here, as our goal is to learn machine learning in scikit-learn, not to produce publishable results.
The key takeaway is that, in just a few lines of code, we've initialized a linear regression model, fitted it to some data, and generated new predicted scores. This basic pattern is common to all supervised estimators in scikit-learn.
Just to underscore how little we had to do, here's the whole example again, in three lines, that:
1. Initialize the linear regression estimator.
2. Fit the model.
3. Generate predictions.
```{code-cell} ipython3
model = LinearRegression()
model.fit(X, y)
y_pred = model.predict(X)
```
Once we've fit our model, we're going to want to see some results. At this point, if you're used to a point-and-click package like SPSS, or a statistics-oriented language like R, you might expect to see a big dump of information that includes things like regression coefficients, standard errors, p-values, $R^2$ values, and so on. Well... we're not really going to talk about those things here. We *could* get that kind of information out of other Python packages, though (see {numref}`ml-sklearn-addtl-resources`). So, yes, we *could* do this kind of thing in Python. But it isn't really what machine learning (or scikit-learn) is about. Instead, the focus in machine learning is on prediction. Typically, we have some quantitative metric of predictive performance we care about, and to the degree that a model can produce better values on that metric, we incline to evaluate the model more favorably. This doesn't mean that we have to single-mindedly base our evaluation of model on *just* one quantity; in practice, many other considerations may come into play (e.g., computational efficiency, interpretability, etc.). The point is just that machine learning practitioners tend to care much more than traditional scientists do about what models can actually *do*, and are much less interested in drawing conclusions about the estimated values of the model's internal parameters. With that performance-oriented goal in mind, let's spend a bit of time thinking about how to measure the quality of a model's predictions.s
+++
### Quantifying performance in scikit-learn: the metrics module
There are many metrics we could use to quantify the quality of the predictions our model generates. Scikit-learn conveniently packages some of the most commonly used performance metrics in its `metrics` module. As is true of `Estimator` objects, all metric functions in scikit-learn have the same interface: we pass in the true scores and the model's predicted scores, respectively, and the function returns a number quantifying the relationship between the two.
#### The coefficient of determination ($R^2$)
```{eval-rst}
.. index::
single: Coefficient of determination
```
Let's see how this works using one of the most commonly used metrics in regression problems: the coefficient of determination, or $R^2$, which quantifies the proportion of variance in the outcome variable (e.g., age) explained by the fitted model. We start by importing the `r2_score` function from the `sklearn.metrics` module. Then, we generate predicted values so we can compare them with the ground truth values. Generally, scoring functions are called by passing an array of true scores and and an array of predicted scores as inputs, which is what we do here.
```{code-cell} ipython3
from sklearn.metrics import r2_score
y_pred = model.predict(X)
r2_score(y, y_pred)
```
Our fitted linear regression model explains around 20% of the variance in age. Whether this is good or bad is a matter of perspective, but at the very least it's clear that we can non-trivially predict people's ages from a handful of structural brain features.
```{admonition} Exercises
1. Replace `r2_score` in the above code with `mean_squared_error`, `mean_absolute_error`, or one of the other predefined metrics in the `metrics` module.
2. Write your own metric function. To be a valid metric function it should take the true scores and predicted scores (in that order) as the only arguments and compute a single number.
```
#### Built-in scoring
For convenience, supervised scikit-learn estimators have a `.score()` method we can use as an alternative to the above. Instead of generating predicted scores and then explicitly feeding them to a metric function like `r2_score`, we can call `.score()` directly on the estimator after the `fit()` step, and the prediction will be done implicitly. The only downside of this approach is that we lose the ability to specify which scoring metric to use; the choice is made for us. In the example below, we will initialize a linear regression, fit it with data and then, instead of generating predictions and comparing them to the data, we directly produce the score. The `LinearRegression` object uses $R^2$ for scoring.
```{code-cell} ipython3
est = LinearRegression()
est.fit(X, y)
print(est.score(X, y))
```
## Classification example: autism classification
```{eval-rst}
.. index::
single: Classification
```
Now let's look at classification. In this case, the target labels we're trying to predict are discrete. In general, one can always turn a regression problem into a classification problem by discretizing the data in some way. For example, we could binarize our continuous age variable around the median value, giving us two equal-sized YOUNG and OLD groups. You'll see this done a lot in the literature, but frankly, discretizing continuous data for classification purposes is almost always a really bad idea (if you're interested in learning more about that, refer to {cite}`MacCallum2002-du`), and we're mostly mentioning that it *can* be done in order to point out that it generally *shouldn't* be done.
Fortunately, no continuous variables were harmed in the making of this book. Instead of dichotomizing continuous labels, we'll use a different set of labels that are naturally discrete: diagnosis group. Recall that ABIDE II is a project primarily interested in understanding autism, and roughly half of the participants in the dataset are diagnosed with autism (one might reasonably argue that the pathologies *underlying* autism could be dimensional rather than discrete in nature, but the diagnoses themselves are definitely discrete). So let's see if we can predict autism diagnosis, rather than age, from structural brain features.
To get a very rough qualitative sense of how difficult a prediction problem this is likely to be, we can do a couple of things. First, we can look at the bivariate correlations between age and each feature in turn, sorting them by strength of association. The `.corrwith()` is a DataFrame method that allows us to correlate
each column with another column taken from a passed Series or DataFrame.
```{code-cell} ipython3
corrs = features.corrwith(phenotypes['group'])
print(corrs.sort_values().round(2))
```
We can immediately see that none of our 1,440 features are very strongly correlated with diagnosis group (the largest correlations are around r = 0.19). This doesn't mean that all is lost, however; even if each feature is individually only slightly predictive of diagnosis group individually, the full set of 1,440 features could still conceivably be very strongly predictive of diagnosis group in the aggregate.
One way to get a cursory sense of whether *that* might be true (i.e., whether combining features is likely to help separate autistic participants from controls) is to visualize diagnosis group as a function of a few brain features. We can't visualize very well in more than 3 dimensions, so let's pick the 3 most strongly correlated features and use those. To probe for separability of classes based on combinations of these variables, we plot three scatter plots of the relationships between them.
```{code-cell} ipython3
fs_vars = ['fsLGI_L_OFC_ROI', 'fsLGI_R_STSvp_ROI', 'fsCT_L_10d_ROI']
x, y, z = features[fs_vars].values.T
import matplotlib.pyplot as plt
fig, ax = plt.subplots(1, 3)
ax[0].scatter(x, y, c=phenotypes['group'])
ax[0].set(xlabel=fs_vars[0], ylabel=fs_vars[1])
ax[1].scatter(z, y, c=phenotypes['group'])
ax[1].set(xlabel=fs_vars[2], ylabel=fs_vars[1])
ax[2].scatter(x, z, c=phenotypes['group'])
ax[2].set(xlabel=fs_vars[0], ylabel=fs_vars[2])
fig.set_tight_layout("tight")
fig.set_size_inches([10,4])
```
Remember in {numref}`class`, when we mentioned that in real-world classification problems, classes are often very hard to separate cleanly, and distributions tend to overlap heavily in feature space? Well that's exactly what we're looking at here. The black and white dots are the two groups of participants, and it's pretty clear at a glance that there isn't going to be any hyperplane we could plunk down through that plot that would perfectly separate the two groups. Again though, we're still using only 3 of our features here, and there are 1,337 others. So let's see how we do when we scale things up to higher dimensions.
### Applying classifiers
Okay, onto the actual classification. How do we apply classification estimators in scikit-learn? It's a trick question! We've actually already covered it in the regression example. There's essentially no difference in the way we interact with regression and classification estimators; we just have to be careful not to pass discrete labels to regression estimators, or continuous labels to classification estimators. But mechanically, we construct and fit the models in the same way.
```{eval-rst}
.. index::
single: Gaussian Naive Bayes classifier
```
Let's try this out with a *gaussian Naive Bayes* (GNB) classifier. This is a
simple classification approach based on a "naive" application of Bayes' Theorem.
The naivete stems from the classifier's assumption that all features are
independent of one another once we condition on class assignment. This
assumption greatly simplifies analysis, though it's nearly always false. GNB is
a good classifier to use as a performance baseline, because it does surprisingly
well in many situations, and is extremely computationally efficient, so it
should be quick. Naive Bayes classifiers have no trouble handling large sets of
highly correlated features, and are also relatively resilient to overfitting
(we'll discuss overfitting in detail in {numref}`ml-overfitting`, so don't worry
if the term is new to you here). So we'll just throw everything we have at the
classifier, and try to predict diagnosis class from all 1,440 features.
```{code-cell} ipython3
from sklearn.naive_bayes import GaussianNB
gnb = GaussianNB()
y = phenotypes['group']
gnb = gnb.fit(features, y)
```
That's all there is to it! See how similar the code is to the regression example we did before?
### Evaluating classification performance
We already know how to obtain a performance estimate from the built-in scorer, so let's do that:
```{code-cell} ipython3
print(gnb.score(features, y))
```
In this case, the default scoring metric is just overall accuracy (i.e., the proportion of all samples that were assigned to the correct class). We see that the model classifies about 63% of samples correctly.
Is 63% good or bad? As always, it depends. It's obviously nowhere near good enough to be useful in practical applications; on the other hand, it's clearly better than randomly guessing. Or is it? That would only be the case if there is an equal number of participants that were diagnosed with autism as there are healthy controls. In some datasets, there is might be a *class imbalance*, in which case, randomly guessing might do better than 50% (for example, by guessing that everyone belongs to the more prevalent class). This is something that needs to be taken into account when evaluating classification results. Fortunately, we can compute some more nuanced scores.
#### Classification reports
The raw accuracy score is a good place to start when evaluating performance, but it often masks important nuances. We can get some more information using the `classification_report` utility, which breaks down classification performance into separate `precision`, `recall`, and `f1-score` metrics (we could also get each of these individually from the `sklearn.metrics` module). Precision is also knows as positive predictive value; it tells us the proportion of cases labeled as positive that truly are positive (e.g., the proportion of cases the classifier labels autistic that really are autistic). Recall (or sensitivity) tells us the proportion of true positive cases that were labeled as such by the classifier. The F1 score is the harmonic mean of recall and sensitivity, and thus tries to summarize both using a single quantity.
The `classification_report` function reports these quantities both for the dataset as a whole, and broken down separately for each class:
```{code-cell} ipython3
from sklearn.metrics import classification_report
y_pred = gnb.predict(features)
print(classification_report(y, y_pred))
```
Notice that while the precision in this case is almost identical for the two groups, the recall (or sensitivity) differs substantially: control subjects have a higher probability (78%) of being correctly labeled than do subjects with an autism diagnosis (44%). This suggests a bias in the model—namely, that the model is more likely to assign cases to the control group than the autism group. We can verify this by looking at the overall proportion of cases our trained model labels as controls:
```{code-cell} ipython3
(y_pred == 2).mean()
```
Sure enough, over two-thirds of cases are classified as controls, even though the ground truth is that only 54% of the sample is made up of controls.
Note that whether this bias in the trained model is a good or bad thing depends on one's context and goals. In the real world, there are often asymmetries in the costs associated with assigning different labels. For example, an HIV test should probably try very hard to minimize false negatives, because the cost of incorrectly telling someone they don't have HIV is likely to be higher than the cost of incorrectly telling them they *do* have HIV (seeing as follow-up tests will rapidly identify the latter error). On the other hand, if all we care about is maximizing overall classification accuracy for it's own sake (e.g., as a learning exercise), then we might be perfectly happy to accept a model with these kinds of class-wise biases, providing their introduction helps us improve our overall performance. We'll spend much more time talking about tradeoffs and context-dependence in {numref}`ml-selection`, when we talk about the bias-variance tradeoff.
## Clustering example: are there neural subtype of autism?
```{eval-rst}
.. index::
single: Clustering
```
Recall that clustering is a form of unsupervised learning where we seek to assign our observations to discrete clusters in the absence of any knowledge of the ground truth (if there is one). Clustering applications are common in psychiatry and psychiatric imaging, as researchers often want to determine whether, e.g., patients with a particular diagnosis appear to cluster into somewhat distinct subtypes—with potential implications for prognosis, treatment, etc.
We can ask this question with respect to the ABIDE data we have available. If we take only those subjects with an autism diagnosis, can we cluster subjects into discrete subtypes based on differences in the brain features we have available?
At this point it probably won't surprise you to hear that scikit-learn contains implementations of quite a few popular clustering algorithms. Clustering estimators in sklearn are located, as you might intuit, in the `sklearn.cluster` module. We'll focus our attention on what is arguably the most widely used clustering algorithm, namely, *k*-means. In *k*-means clustering, we assign samples to *k* discrete clusters in such a way as to minimize the distance from each observation to the centroid of the cluster it belongs to, and maximize the distance between the cluster centroids. We won't get into the details of the k-means algorithm here; instead, we'll simply demonstrate how we'd go about running a cluster analysis with scikit-learn. As you might expect by now, our code will look a lot like it did for regression and classification.
```{eval-rst}
.. index::
single: K-means clustering
```
In principle, we could use all 1,440 features if we wanted to (though computation would likely be slow). But the unsupervised nature of clustering means that evaluating clustering solutions in high dimensions can be quite difficult. So let's focus on clustering observations in just 2 dimensions, and if after that you want to crank up the number of features.
```{code-cell} ipython3
n_features = 2
```
Next, we will use a few lines of code that will get our `X` data. We start by selecting subjects. We select only the subjects with an autism diagnosis and then randomly sample brain features.
```{code-cell} ipython3
aut_grp = phenotypes['group'] == 1
dx_1 = features[aut_grp].sample(n_features, axis=1)
columns = dx_1.columns
```
The `sklearn.preprocessing` module contains a bunch of useful utilities. In this case we're going to standardize our columns (i.e., to transform them so they have a mean of 0 and variance of 1). Otherwise k-means will likely weight some features more than others.
```{code-cell} ipython3
from sklearn.preprocessing import scale
dx_1 = scale(dx_1)
```
The actual clustering is, as usual, just a couple of lines of code. We'll use `KMeans` here, and we would encourage you to experiment with others, but note that many of the implemented clustering algorithms will be very slow if you crank up the number of features. Here, we have to stipulate the number of clusters (k) in advance and this is set even before fitting to the data. Instead of fitting and then predicting, we do both of these operations in one shot, using the `fit_predict` method of the `KMeans` class.
```{code-cell} ipython3
from sklearn.cluster import KMeans
K = 4
km = KMeans(K)
clusters = km.fit_predict(dx_1)
```
Notice that the only real difference between the above code and our earlier regression/classification examples is that we're no longer passing in a `y` array of labels—because there isn't one! We're just trying to find meaningful structure in the data, with no access to any ground truth.
Also notice that we need to specify the number of clusters *k* ourselves. There are literally hundreds of methods people have developed to try to identify the "optimal" value of *k*, and some of them are included in the `sklearn.metrics` module. But for our purposes, let's just focus on the solution we get for $k = 4$.
Observe that the `clusters` array we get back by calling either `.predict()`, or (as in this case) `fit_predict` on a clustering estimator gives us the assigned cluster labels for each observation. Since we only used two features to do our clustering, we can easily visualize the results of the clustering in a two-dimensional scatter plot, where each point is colored based on the cluster that was assigned to it (this is what the `c` key-word argument to scatter means):
```{code-cell} ipython3
fig, ax = plt.subplots()
ax.scatter(dx_1[:, 0], dx_1[:, 1], c=clusters, edgecolor='k')
ax.scatter(*km.cluster_centers_.T, c=[0,1,2,3], edgecolor='k', linewidth=3, s=100)
ax.set_xlabel(columns[0])
ax.set_ylabel(columns[1])
g = ax.grid(None)
```
The clustering *looks* reasonable to the eye... but how far should we trust it? And how should we think about these cluster labels? Are they merely convenient descriptions of the data that facilitate exploration of differences between subgroups of participants, or do we think we're gaining deep insights into the "true" neural bases of different subtypes of autism? That's not a question we're going to try to answer here, but it's an important one to think about whenever one applies clustering methods to complex real-world data.
(ml-sklearn-addtl-resources)=
### Additional resources
As we mentioned in the text, in machine learning, we are usually more interesting in quantifying the predictive capacity of a particular model. However, we acknowledge that sometimes what you want to do is to fit and evaluate parametric models. This is possible to do in Python using a package called [StatsModels](https://www.statsmodels.org), which you can do . Another package that elegantly implements many of the statistical methods that behavioral scientists traditionally use (e.g., Analysis of Variance, or ANOVA) is [Pingouin](https://pingouin-stats.org/).
```{eval-rst}
.. index::
single: Statsmodels
```
```{eval-rst}
.. index::
single: Pingouin
```
To learn more about some of the concepts that we introduced here (e.g., Gaussian Naive Bayes and K-means clustering), we would particularly recommend the introductory statistical learning text [An Introduction to Statistical Learning](https://www.statlearning.com/). The book is an excellent introduction and is much more comprehensive than we could be in this chapter. Among it's many merits is also the fact that it is available to download for free from the book website.
There are many software packages that implement machine learning methods for analysis of neuroimaging data. One that stands out as particularly compatible with the Scikit Learn philosophy and API is the [nilearn](https://nilearn.github.io). One reason for this compatibility is that many of the contributors to the nilearn software are also developers of Scikit Learn.