Repository: scikit-learn/scikit-learn
Pull request: DOC Improve documentation regarding some pitfalls in interpretation
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at
the contribution guidelines: https://github.com/scikit-learn/scikit-learn/blob/main/CONTRIBUTING.md
-->

#### Reference Issues/PRs
Fixes #19413
<!--
Example: Fixes #1234. See also #3456.
Please use keywords (e.g., Fixes) to create link to the issues or pull requests
you resolved, so that they will automatically be closed when your pull request
is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->


#### What does this implement/fix? Explain your changes.
In the [Common pitfalls in interpretation of coefficients of linear models](https://scikit-learn.org/stable/auto_examples/inspection/plot_linear_model_coefficient_interpretation.html):
* add warning on wrongly giving a causal interpretation to  coefficients.
* add warning on how interpretation from a model may not apply to the Data Generating Process 
   when the model is poor or sample used is not representative of the population.

Add a tutorial to show--via simulation--that coefficients can be biased in the presence of unobserved
confounders. 

#### Any other comments?


<!--
Please be aware that we are a loose team of volunteers so patience is
necessary; assistance handling other issues is very welcome. We value
all user contributions, no matter how minor they are. If we are slow to
review, either the pull request needs some benchmarking, tinkering,
convincing, etc. or more likely the reviewers are simply busy. In either
case, we ask for your understanding during the review process.
For more information, see our FAQ on this topic:
http://scikit-learn.org/dev/faq.html#why-is-my-pull-request-not-getting-any-attention.

Thanks for contributing!
-->


Problem statement:
Improve documentation regarding some pitfalls in interpretation
We have now a couple of tutorials in the `inspection` section of our documentation. We try to document the pitfalls regarding interpreting the output of the inspection methods. However, it seems that we could improve our discussions several aspects:

* we should emphasize that the inspection tools will always give us some insights regarding our model. However, if our model is not a good one (e.g. bad statistical performance) an interpretation trying to put into perspective our model with the generative model that created the data is potential wrong.
* we do not warn about the non-causal interpretation. Indeed, most of the tools used do not allow to make causal interpretation and we should be assertive regarding this thing.

One possibility would be to recall these facts in the different tutorial by adding some warnings. It could also be interesting to add new tutorials:

* for the first point, we could imagine creating some data and intentionally getting a bad model that should provide wrong interpretation regarding the data generative model.
* for the second point, we could introduce briefly the principle of the causality principle with an example and give an example of what it takes to make such causality deduction.

@dsleo @ogrisel might want to add some thoughts regarding this proposal.

Hints from issue discussion:
Thanks for opening the issue ! 

For the first point, we could give an example with synthetic data to show how controlling over different sets of variables leads to different coefficients which don't recover the "causal ones" _even_ if the learnt model has very good fit. I can try to craft an example.

For the second point, we could add to the [common pitfall in interpretation of coefficients of linear model](https://scikit-learn.org/stable/auto_examples/inspection/plot_linear_model_coefficient_interpretation.html) a first warning on the non causal-interpretation of the coefficients. Something along the following lines could be added at the end of the first paragraph:
_Coefficients in linear models have no causal interpretation. They merely give indication on how the model relies on the features not of their importance in the system they seek to describe or predict._
We could also add a word of warning in the [PDP/ICE doc](https://scikit-learn.org/stable/auto_examples/inspection/plot_partial_dependence.html#sphx-glr-auto-examples-inspection-plot-partial-dependence-py).

For the last point, what do you have in mind by "the causality principle" ? How do you see the scope of a dedicated tutorial/documentation on causality ?

Commit message: Fix pep8 violations

Patch to review:
diff --git a/examples/inspection/plot_causal_interpretation.py b/examples/inspection/plot_causal_interpretation.py
new file mode 100644
index 0000000000000..49be8a661b907
--- /dev/null
+++ b/examples/inspection/plot_causal_interpretation.py
@@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
+
+clf_without = LinearRegression()
+clf_without.fit(np.stack([experiences, parent_hourly_wages, college_degrees],
+                axis=1), hourly_wages)
+
+
+print('College degree coefficient with ability control: {}'.format(
+    clf_with.coef_[2]))
+
+print('College degree coefficient without ability control: {}'.format(
+    clf_without.coef_[2]))
+
+# %%
+# Lessons learned
+# ---------------
+#
+# Machine learning models are not reliable for inferring causation causal effect.
+# While we showed this with a linear model, OVB can strike any type of models.
+#
+# Whenever interpreting a coefficient or a change in predictions
+# brought about by a change in one of the features, it is important to think of
+# potentially unobsered variables and whether they could be correlated with both
+# the feature in question and the target variable.
+# Such variables are called `Confounding Variables
+# <https://en.wikipedia.org/wiki/Confounding>`_.
+# To avoid these, researchers usually conduct experiments in which
+# the treatment variable (e.g. college degree) is randomized. When an experiment
+# is prohibitively expensive or unethical, researchers
+# can somethimes use other causal inference techniques such
+# as `Instrumental Variables
+# <https://en.wikipedia.org/wiki/Instrumental_variables_estimation>`_ (IV)
+# estimations.
diff --git a/examples/inspection/plot_linear_model_coefficient_interpretation.py b/examples/inspection/plot_linear_model_coefficient_interpretation.py
index 459b180f00e36..f9b831ab6d8e4 100644
--- a/examples/inspection/plot_linear_model_coefficient_interpretation.py
+++ b/examples/inspection/plot_linear_model_coefficient_interpretation.py
@@ -653,7 +653,49 @@
 # information on top of the other features. Additionnaly, this selection
 # is unstable for correlated features, and should be interpreted with
 # caution.
+
+# %%
+# Causal interpretation
+# ---------------------
+#
+# Policy makers might want to know the effect of education on wage to
+# assess whether or not a a certain policy designed to entice people to
+# pursue more education would make economic sense.
+# While machine Learning models are great for measuring statistical
+# associations, they are generally unable to infer causal effects.
+#
+# It might be tempting to look at the coefficient of education on wage
+# from our last model (or any model for that matter) and conclude that it
+# captures the true effect of a change in the standardized education variable
+# on wages.
+#
+# Unfortunately there are likely unobserved confounding variables that either
+# inflate or deflate that coefficient. A confounding variable is a variable that
+# causes both EDUCATION and WAGE. One example of such variable is ability.
+# Presumably, more able people are more likely to pursue education while at the
+# same time being more likely to earn a higher hourly wage at any level of
+# education. In this case, ability induces a positive `Omitted Variable Bias
+# <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+# on the EDUCATION coefficient, thereby exaggerating the effect of education
+# on wages.
 #
+# See the :ref:`sphx_glr_auto_examples_inspection_plot_causal_interpretation.py`
+# for a simulated case of ability OVB.
+
+# %%
+# Warning: data and model quality
+# -------------------------------
+#
+# Keep in mind that the outcome y and features X are the product
+# of a Data Generating Process that is hidden from us. Machine
+# learning models are trained to approximate the unobserved
+# mathematical function that links X to y from sample data. As a
+# result, any interpretation made about a model may not necessarily
+# generalize to the true Data Generating Process. This is especially
+# true when the model is of bad quality or when the sample data is
+# not representative of the population.
+
+# %%
 # Lessons learned
 # ---------------
 #
@@ -669,3 +711,7 @@
 #   coefficients could significantly vary from one another.
 # * Inspecting coefficients across the folds of a cross-validation loop
 #   gives an idea of their stability.
+# * Coefficient are unlikely to have any causal meaning. They tend
+#   to be biased by unobserved confounders.
+# * Inspection tools may not necessarily provide insights on the true
+#   Data Generating Process.


Code review comments:
diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
You can remove the underscore in the naming. The plot will be nicer then.

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
I think it would be better to create a variable `X` (a dataframe) an `y`. In the second case, we could easily do `X.drop(columns="ability")`

@user1:
It would be a bit more readable than the columns stacking.

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
+
+clf_without = LinearRegression()
+clf_without.fit(np.stack([experiences, parent_hourly_wages, college_degrees],
+                axis=1), hourly_wages)
+
+
+print('College degree coefficient with ability control: {}'.format(
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
From what I recall, we also wanted to show that the model fitting was quite reasonable indeed.
So we might want to add a train/test split and just get the score on the test.

@user1:
Here you can use an f-string instead of `format`.

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
+
+clf_without = LinearRegression()
+clf_without.fit(np.stack([experiences, parent_hourly_wages, college_degrees],
+                axis=1), hourly_wages)
+
+
+print('College degree coefficient with ability control: {}'.format(
+    clf_with.coef_[2]))
+
+print('College degree coefficient without ability control: {}'.format(
+    clf_without.coef_[2]))
+
+# %%
+# Lessons learned
+# ---------------
+#
+# Machine learning models are not reliable for inferring causation causal effect.
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
maybe 

```suggestion
# Machine learning models are not designed for estimation of causal effects.
```

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
```suggestion
n_samples = 10_000
```

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
+
+clf_without = LinearRegression()
+clf_without.fit(np.stack([experiences, parent_hourly_wages, college_degrees],
+                axis=1), hourly_wages)
+
+
+print('College degree coefficient with ability control: {}'.format(
+    clf_with.coef_[2]))
+
+print('College degree coefficient without ability control: {}'.format(
+    clf_without.coef_[2]))
+
+# %%
+# Lessons learned
+# ---------------
+#
+# Machine learning models are not reliable for inferring causation causal effect.
+# While we showed this with a linear model, OVB can strike any type of models.
+#
+# Whenever interpreting a coefficient or a change in predictions
+# brought about by a change in one of the features, it is important to think of
+# potentially unobsered variables and whether they could be correlated with both
+# the feature in question and the target variable.
+# Such variables are called `Confounding Variables
+# <https://en.wikipedia.org/wiki/Confounding>`_.
+# To avoid these, researchers usually conduct experiments in which
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
Small clarification as we can't really avoid these but merely get around. 
```suggestion
# In order to still estimate causal effect in the presence of confounding
# variables, researchers usually resort to experiments in which
```

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
+clf_with.fit(np.stack([experiences, parent_hourly_wages, college_degrees,
+             abilities], axis=1), hourly_wages)
+
+clf_without = LinearRegression()
+clf_without.fit(np.stack([experiences, parent_hourly_wages, college_degrees],
+                axis=1), hourly_wages)
+
+
+print('College degree coefficient with ability control: {}'.format(
+    clf_with.coef_[2]))
+
+print('College degree coefficient without ability control: {}'.format(
+    clf_without.coef_[2]))
+
+# %%
+# Lessons learned
+# ---------------
+#
+# Machine learning models are not reliable for inferring causation causal effect.
+# While we showed this with a linear model, OVB can strike any type of models.
+#
+# Whenever interpreting a coefficient or a change in predictions
+# brought about by a change in one of the features, it is important to think of
+# potentially unobsered variables and whether they could be correlated with both
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
typo 

```suggestion
# potentially unobserved variables and whether they could be correlated with both
```

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
The best here is to make a run of black on the example for the coding style

```suggestion
hourly_wages = (
    0.2 * experiences + parent_hourly_wages + 2
    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N
)
```

diff_hunk: @@ -653,7 +653,49 @@
 # information on top of the other features. Additionnaly, this selection
 # is unstable for correlated features, and should be interpreted with
 # caution.
+
+# %%
+# Causal interpretation
+# ---------------------
+#
+# Policy makers might want to know the effect of education on wage to
+# assess whether or not a a certain policy designed to entice people to
+# pursue more education would make economic sense.
+# While machine Learning models are great for measuring statistical
+# associations, they are generally unable to infer causal effects.
+#
+# It might be tempting to look at the coefficient of education on wage
+# from our last model (or any model for that matter) and conclude that it
+# captures the true effect of a change in the standardized education variable
+# on wages.
+#
+# Unfortunately there are likely unobserved confounding variables that either
+# inflate or deflate that coefficient. A confounding variable is a variable that
+# causes both EDUCATION and WAGE. One example of such variable is ability.
+# Presumably, more able people are more likely to pursue education while at the
+# same time being more likely to earn a higher hourly wage at any level of
+# education. In this case, ability induces a positive `Omitted Variable Bias
+# <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+# on the EDUCATION coefficient, thereby exaggerating the effect of education
+# on wages.
 #
+# See the :ref:`sphx_glr_auto_examples_inspection_plot_causal_interpretation.py`
+# for a simulated case of ability OVB.
+
+# %%
+# Warning: data and model quality
+# -------------------------------
+#
+# Keep in mind that the outcome y and features X are the product
+# of a Data Generating Process that is hidden from us. Machine
path: examples/inspection/plot_linear_model_coefficient_interpretation.py
text: @user1:
```suggestion
# of a data generating process that is hidden from us. Machine
```

diff_hunk: @@ -0,0 +1,126 @@
+"""
+================================================
+Causal Interpretation of Machine Learning Models
+================================================
+
+Machine Learning models are great for measuring statistical associations.
+Unfortunately, unless we're willing to make strong assumptions about the data,
+those models are unable to infer causal effects.
+
+To illustrate this, we will simulate a situation in which we try to answer
+one of the most important questions in economics of education:
+what is the causal effect of earning a college degree on hourly wages?
+Although the answer to this question is crucial to policy makers,
+`Omitted-Variable Biases <https://en.wikipedia.org/wiki/Omitted-variable_bias>`_ (OVB)
+prevent us from identifying that causal effect.
+
+
+.. contents::
+   :local:
+   :depth: 1
+"""
+
+print(__doc__)
+
+import numpy as np
+from sklearn.linear_model import LinearRegression
+
+# %%
+# The dataset: simulated hourly wages
+# -----------------------------------
+#
+# The data generating process is laid out in the code below.
+# Work experience in years and a measure of ability are drawn
+# from a Normal distribution;
+# the hourly wage of one of the parents is drawn from Beta distribution.
+# We then create an indicator of college degree which is positivily
+# impacted by ability and parental hourly wage. Finally, we model
+# hourly wages as a linear function of all the previous variables and a
+# random component. Note that all variables have a positivie effect on
+# hourly wages.
+
+N = 10000
+experiences = np.random.normal(20, 10, size=N).astype(int)
+experiences[experiences < 0] = 0
+abilities = np.random.normal(0, 0.15, size=N)
+parent_hourly_wages = 50 * np.random.beta(2, 8, size=N)
+parent_hourly_wages[parent_hourly_wages < 0] = 0
+
+
+college_degrees = (9 * abilities + 0.02 * parent_hourly_wages
+                   + np.random.randn() > 0.7).astype(int)
+
+hourly_wages = 0.2 * experiences + parent_hourly_wages + 2 \
+    * college_degrees + 5 * abilities + np.random.normal(0, 1, size=N)
+
+hourly_wages[hourly_wages < 0] = 0
+
+# %%
+# Description of the simulated data
+# ---------------------------------
+#
+# The following plot shows the distribution of each variable,
+# and pairwise scatter plots. Key to our OVB story is the positive
+# relationship between ability and college degree.
+
+import pandas as pd
+import seaborn as sns
+
+df = pd.DataFrame({
+    'hourly_wage': hourly_wages,
+    'experience': experiences,
+    'parent_hourly_wage': parent_hourly_wages,
+    'college_degree': college_degrees,
+    'ability': abilities
+    })
+
+_ = sns.pairplot(df, diag_kind='kde')
+
+# %%
+# Predicting income with and without the ability feature
+# ------------------------------------------------------
+#
+# Let's now train two :class:`~sklearn.linear_model.LinearRegression`
+# models to predict our hourly wage. We include the ability feature
+# in the first model and show that our estimate of the college degree
+# coefficient is close to 2 which is the true causal effect from our
+# data generating process. In real life, intellectual ability is either
+# never observed or poorly measured (e.g. IQ score).
+# Researchers are forced to "omit" the ability feature from their models,
+# thereby inflating the estimate via a positive OVB.
+
+clf_with = LinearRegression()
path: examples/inspection/plot_causal_interpretation.py
text: @user1:
usually, we use `clf` for classifier. It might be better to use `regressor_with_ability` and `regressor_without_ability`