# Covariate Shift

A fundamental assumption in almost all [supervised learning](https://en.wikipedia.org/wiki/Supervised_learning) methods is that training and test samples are drawn from the same [probability distribution](https://en.wikipedia.org/wiki/Probability_distribution). However, in practice, this assumption is rarely satisfied and standard machine learning models may not work as well as anticipated. [Covariate shift](https://www.quora.com/What-is-Covariate-shift) refers to the situation where the probability distribution of covariates changes between training and test data. *Shikhar Gupta*, Master's student in Data Science, gives the following great visual representation of this in [How Dissimilar are my Train and Test Data](https://towardsdatascience.com/how-dis-similar-are-my-train-and-test-data-56af3923de9b) in the **Towards Data Science** blog.



It is clear that if covariate shift is not accounted for correctly, it can lead to poor generalization.

So why are we even talking about covariate shift? Well, one of the causes of covariate shift is [sample selection bias](https://en.wikipedia.org/wiki/Selection_bias). [Dataset Shift in Machine Learning](http://www.acad.bg/ebook/ml/The.MIT.Press.Dataset.Shift.in.Machine.Learning.Feb.2009.eBook-DDU.pdf) states that "Sample selection bias occurs when the training data points {$x_i$} (the sample) do not accurately represent the distribution of the test scenario (the population) due to a selection process for each item that is (usually implicitly) dependent on the target variable $y_i$." If you recall, during [feature construction](3.0-build-features.ipynb), we made a fleeting mention of the selection bias that was introduced by the way we sampled the training and test / validation sets. The training set was selected to consist entirely of deceased physicists and the test / validation set was selected to consist entirely of living physicists. The feature building process has already hinted that the physicists in these datasets have different characteristics. Remember the different feature columns created by one-hot encoding?

The data was purposely sampled in this way due to the realization that the Nobel Prize in Physics cannot be awarded posthumously. The selection bias is an inherent part of the problem, as one of the goals of this project is to try to predict the next Physics Nobel Laureates, who obviously must be alive. As a result, the selection bias is something that we have to live with. The aim here is to see if we can formally detect whether a covariate shift occurs between the training and test data. It is important to note that we will be using the validation set as a *proxy* for the test data, since the true performance of the model will be evaluated on the latter, which is meant to be unseen data. Naturally, this doesn't take into account that there may also be a covariate shift between the validation and the test set. However, hopefully the random sampling process that we employed to divide the living physicists into the validation and test sets will have mitigated this.


```python
import numpy as np
import pandas as pd
from sklearn.linear_model import LogisticRegression
from sklearn.metrics import make_scorer
from sklearn.metrics import matthews_corrcoef
from sklearn.model_selection import train_test_split
from sklearn.model_selection import GridSearchCV
from sklearn.utils import indices_to_mask

from src.features.features_utils import convert_categoricals_to_numerical
from src.models.metrics_utils import print_matthews_corrcoef
from src.stats.stats_utils import bootstrap_prediction
from src.stats.stats_utils import percentile_conf_int
from src.visualization.visualization_utils import plot_bootstrap_statistics
from src.visualization.visualization_utils import plot_logistic_regression_odds_ratio

%matplotlib inline
```

## Classifier Two-Sample Hypothesis Testing

A formal way of detecting whether there is a covariate shift between two samples (sets of identically and independently distributed examples) is to perform a [hypothesis test](https://en.wikipedia.org/wiki/Statistical_hypothesis_testing) known as a **two-sample test**. The goal of two-sample tests is to assess whether two samples, say $S_{train} \sim P$ and $S_{validation} \sim Q$, are drawn from the same probability distribution. Two-sample tests evaluate the difference between two distributions using the value of a [test statistic](https://en.wikipedia.org/wiki/Test_statistic) to either accept or reject the [null hypothesis](https://en.wikipedia.org/wiki/Null_hypothesis) $ H_0: P = Q$.

A rather elegant and informative way to perform a two-sample test is to train a binary classifier to distinguish between instances in $S_{train}$ and $S_{validation}$. The intuition is that if the null hypothesis is true, then the performance of such a binary classifier, as measured by an appropriate test statistic, will be approximately the same as random guessing. If the performance of such a classifier is better than chance-level, the null hypothesis is rejected in favor of the [alternative hypothesis](https://en.wikipedia.org/wiki/Alternative_hypothesis) $ H_1: P \neq Q$.

To test whether the null hypothesis $ H_0: P = Q$ is true, we will loosely follow the steps of *Lopez-Paz* and *Oquab* in [Revisiting Classifier Two-Sample Tests](https://arxiv.org/pdf/1610.06545). However, there are some notable differences to our approach here, which we will point out. Also, we will try not to over-burden the reader with heavy notation. We will be taking the following steps for both the original features and the features created from the topic modeling:

1. **Construct the dataset that is the union of the training and validation sets**, {$X = S_{train} \cup S_{validation}$, $y = y_{train} \cup y_{validation}$} where $X$ is the feature matrix and $y$ the target vector. Assign $y_{train} = 0$ and $y_{validation} = 1$ for all instances in the training and validation sets, respectively.
2. **Shuffle and split** $X$ **at random into disjoint training and test subsets** $X_{train}$ and $X_{test}$.
3. **Train a binary classifier** on $X_{train}$.
4. **Evaluate the performance of the binary classifier** by computing the Matthews Correlation Coefficient (MCC) as the test statistic on $X_{test}$: 
\begin{equation}
\hat{t} \equiv MCC = \frac{TP \times TN - FP \times FN}{{\sqrt{(TP + FP)(TP + FN)(TN + FP)(TN + FN)}}}
\end{equation}
where TP is the number of true positives, TN the number of true negatives, FP the number of false positives and FN the number of false negatives. This differs from the classification accuracy test statistic that is used in the paper. Accuracy is not an appropriate metric for the same reasons mentioned during our creation of the [baseline model](5.0-baseline-model.ipynb). There are many more instances in the training set than the validation set. MCC is a metric that will account for this imbalance of classes.
5. **Accept or reject the null hypothesis** by computing a 95% [bootstrap](https://en.wikipedia.org/wiki/Bootstrapping_(statistics)) [confidence interval](https://en.wikipedia.org/wiki/Confidence_interval) for $\hat{t} \equiv MCC$. If the confidence interval contains the value of $\hat{t} \equiv MCC$ corresponding to chance-level performance (i.e. $MCC = 0$) then accept the null hypothesis, otherwise reject it. Again this differs from the [p-value](https://en.wikipedia.org/wiki/P-value) approach taken in the paper; non-parametrics are needed here as we do not know the distribution of $\hat{t} \equiv MCC$. Anyway, confidence intervals are more informative than p-values as they elucidate the magnitude and precision of the estimated effect. Furthemore, [p-values and confidence intervals always agree about statistical significance](http://blog.minitab.com/blog/adventures-in-statistics-2/understanding-hypothesis-tests-confidence-intervals-and-confidence-levels), so the substitution of a confidence interval for a p-value is warranted.

First, let's construct the dataset that is the union of the training and validation data, making sure to convert the categorical fields to a numerical form that is suitable for building machine learning models. 


```python
train_features = pd.read_csv('../data/processed/train-features.csv')
train_features = convert_categoricals_to_numerical(train_features)
train_features.head()
```


```python
train_features_topics = pd.read_csv('../data/processed/train-features-topics.csv')
train_features_topics = convert_categoricals_to_numerical(train_features_topics)
train_features_topics.head()
```


```python
validation_features = pd.read_csv('../data/processed/validation-features.csv')
validation_features = convert_categoricals_to_numerical(validation_features)
validation_features.head()
```


```python
validation_features_topics = pd.read_csv('../data/processed/validation-features-topics.csv')
validation_features_topics = convert_categoricals_to_numerical(validation_features_topics)
validation_features_topics.head()
```


```python
X = train_features.append(validation_features)
assert(len(X) == len(train_features) + len(validation_features))
X.head()
```


```python
X_topics = train_features_topics.append(validation_features_topics)
assert(len(X_topics) == len(train_features_topics) + len(validation_features_topics))
X_topics.head()
```


```python
y = pd.Series(np.concatenate((np.zeros(len(train_features), dtype='int64'),
                              np.ones(len(validation_features), dtype='int64'))),
              index=train_features.index.append(validation_features.index))
assert(y.value_counts().equals(pd.Series([len(train_features), len(validation_features)], index=[0, 1])))
y.head()
```

Second, let's shuffle and split the dataset (in a stratified manner to maintain the class proportions) into disjoint training and test sets. We will use 80% of the data for training and 20% for testing.


```python
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.2, shuffle=True, stratify=y, random_state=0)
display(X_train.head())
display(X_test.head())
display(y_train.head())
y_test.head()
```


```python
X_train_topics = X_topics.loc[X_train.index, :]
X_test_topics = X_topics.loc[X_test.index, :]
display(X_train_topics.head())
X_test_topics.head()
```

Third, let's train a binary classifier on $X_{train}$. The choice of classifier is abitrary, so we will choose logistic regression, making sure to select the regularization parameter, $C$, and the regularization penalty, $L1$ or $L2$, via stratified 5-fold cross-validation. We will be using the MCC for scoring.


```python
def fit_logit_classifier(X, y):
    params = {'C': np.logspace(0, 1, 11), 'penalty': ['l1', 'l2']}
    mcc_scorer = make_scorer(matthews_corrcoef)
    clf = GridSearchCV(LogisticRegression(solver='liblinear'), param_grid=params, scoring=mcc_scorer,
                       cv=5, iid=False, return_train_score=True)
    clf.fit(X, y)
    return clf
```


```python
logit = fit_logit_classifier(X_train, y_train)
logit.best_params_, logit.best_estimator_
```


```python
logit_topics = fit_logit_classifier(X_train_topics, y_train)
logit_topics.best_params_, logit_topics.best_estimator_
```

Fourth, evaluate the performance of the classifier on $X_{test}$ using the MCC.


```python
print_matthews_corrcoef(
    matthews_corrcoef(y_train, logit.best_estimator_.predict(X_train)), 'Features', data_label='train')
test_mcc = matthews_corrcoef(y_test, logit.best_estimator_.predict(X_test))
print_matthews_corrcoef(test_mcc, 'Features', data_label='test')
print_matthews_corrcoef(
    matthews_corrcoef(y_train, logit_topics.best_estimator_.predict(X_train_topics)), 'Topics',
    data_label='train')
test_topics_mcc = matthews_corrcoef(y_test, logit_topics.best_estimator_.predict(X_test_topics))
print_matthews_corrcoef(test_topics_mcc, 'Topics', data_label='test')
```

Fifth, accept or reject the null hypothesis $ H_0: P = Q$ by computing a 95% bootstrap confidence interval.


```python
n_estimators = 1000
max_samples = 0.8
n_jobs = -1
alpha = 0.05
mccs = bootstrap_prediction(
    X, y, base_estimator=logit.best_estimator_, score_func=matthews_corrcoef, n_estimators=n_estimators,
    max_samples=max_samples, n_jobs=n_jobs, random_state=2)
conf_int = percentile_conf_int(mccs, alpha=alpha)
```


```python
stat_label = 'Matthews Correlation Coefficient (MCC)'
ax = plot_bootstrap_statistics(
    mccs, test_mcc, conf_int, alpha, 'test MCC:', stat_label,
    title='Features Bootstrap Matthews Correlation Coefficient (MCC) \nfor 1000 samples')
ax.set_xlim(0, 0.5)
ax.set_ylim(0, 300);
```


```python
mccs = bootstrap_prediction(
    X_topics, y, base_estimator=logit_topics.best_estimator_, score_func=matthews_corrcoef,
    n_estimators=n_estimators, max_samples=max_samples, n_jobs=n_jobs, random_state=3)
conf_int = percentile_conf_int(mccs, alpha=alpha)
```


```python
ax = plot_bootstrap_statistics(
    mccs, test_topics_mcc, conf_int, alpha, 'test MCC:', stat_label,
    title='Topics Bootstrap Matthews Correlation Coefficient (MCC) \nfor 1000 samples')
ax.set_xlim(0, 0.4)
ax.set_ylim(0, 300);
```

The figures illustrate that the distributions of the MCCs for the bootstrap samples are Gaussian-like for both the original features and the topics features. The MCC for the test set of the full dataset, along with the upper and lower values of the 95% confidence intervals, are shown. It is clear that for both sets of features, the confidence interval does not contain the value of chance-level performance ($MCC = 0$). Hence in both cases, there is sufficient evidence to reject the null hypothesis in favor of the alternate hypothesism, $ H_1: P \neq Q$. The conclusion is that the training and validation sets are drawn from different distributions. In other words, there is a covariate shift in both feature sets. We can see that the severity of the covariate shift has been reduced significantly by constructing features with topic modeling.

## Shifting Predictors

We can determine the predictors that exhibit a covariate shift by looking at the coefficients of the logistic regression models. Each coefficient represents the impact that the *presence* vs. *absence* of a predictor has on the [log odds ratio](https://en.wikipedia.org/wiki/Odds_ratio#Role_in_logistic_regression) of a physicist being from the validation set (as opposed to being from the training set). The change in [odds ratio](https://en.wikipedia.org/wiki/Odds_ratio) for each predictor can simply be computed by exponentiating its associated coefficient. 

Formally, a change in odds ratio of 1 for a particular predictor indicates that it is not shifting, whereas a value greater than 1 indicates a shift. As it is likely that a lot of predictors will have odds ratios of slightly over 1, we will loosely define a shifting predictor as one that has a change in odds ratio greater than 1.2. This will give us an idea of the predictors that contribute the most to the covariate shift in the data. These are plotted in the charts below for the two sets of predictors.


```python
ax = plot_logistic_regression_odds_ratio(
    logit.best_estimator_.coef_, top_n=1.2, columns=X.columns, title='Features covariate shifting',
    plotting_context='talk')
ax.figure.set_size_inches(20, 30)
```


```python
ax = plot_logistic_regression_odds_ratio(
    logit_topics.best_estimator_.coef_, top_n=1.2, columns=X_topics.columns,
    title='Topics covariate shifting')
```

The charts certainly make intuitive sense as they illustrate the following changes over time:

1. Increase in globalization of physics.
2. Concentration of physics in the hubs of North America and Europe.
3. Coming to prominence of major American institutions as places of study and work.
4. Broadening of the research fields of theoretical physics and astronomy.

Now that we know which features are shifting, what should be done about it? A simple solution could be to drop these features. However, this would result in some loss of information and we do not yet know if these features are important in predicting our target. Furthermore, this would raise the question as to what should be the minimum change in odds ratio for identifying drifting features. The value of 1.2 chosen above was rather *ad hoc* as it was intended only to illustrate the features that contribute the most to the covariate shift in the data. It is clear that a more principled approach is needed to deal with the covariate shift in the data.

5588c6ec970199986cf7d9fe5a16ba8640e5a613
ipynb
Jupyter Notebook
nobel_physics_prizes/notebooks/5.1-covariate-shift.ipynb
covuworie/nobel-physics-prizes
f89a32cd6eb9bbc9119a231bffee89b177ae847a
MIT
2019-08-21T05:35:42.000Z
2020-10-08T21:28:51.000Z
nobel_physics_prizes/notebooks/5.1-covariate-shift.ipynb
covuworie/nobel-physics-prizes
f89a32cd6eb9bbc9119a231bffee89b177ae847a
MIT
2018-09-01T23:15:59.000Z
2021-02-02T22:01:39.000Z
nobel_physics_prizes/notebooks/5.1-covariate-shift.ipynb
covuworie/nobel-physics-prizes
f89a32cd6eb9bbc9119a231bffee89b177ae847a
MIT
Qwen/Qwen-72B
1. YES
2. YES