Repository: scikit-learn/scikit-learn
Pull request: ENH kind accept a list of str in plot_partial_dependence
closes #19305 

This PR proposes that:

* `kind` supports a list of str of the same length as `features` to specify different PD plots for each axis.
* When a `kind="both"` and a 2-way partial dependence is requested in feature, a warning is raised instead of not making any plot. The warning will stipulate to pass a list to silence it. 

Problem statement:
Limitation of ICE when required for 1D and 2D interactions
The following code requested to display ICE will fail:

```python
features = ["EDUCATION", "AGE", ("EDUCATION", "AGE")]
fig, ax = plt.subplots(figsize=(15, 5))
_ = plot_partial_dependence(
    model.named_steps['histgradientboostingregressor'],
    X_train_preprocessed, features, feature_names=feature_names,
    n_jobs=-1, grid_resolution=20, ax=ax, kind="both",
)
```

It is due to the fact that we are requesting a 2D interaction (`EDUCATION`, `AGE`) for which we cannot return ICE.
I think that we should support this use-case by only returning the PDP plot, potentially with a warning.
It is a better user experience than having to call the function twice.

@NicolasHug @thomasjpfan Do you think it would be a sensible thing to do?

Hints from issue discussion:
I would be okay with no warning and document that the 2d interactions do not have ICE.
I would warn (to educate users that 2-way ICEs are impossible to render in a 2d plot) and also allow for `kind` to be a list of strings. But no strong opinion. Either way, I agree we should document this.
take
> I would warn (to educate users that 2-way ICEs are impossible to render in a 2d plot) and also allow for kind to be a list of strings. 

I think this the best way. One can silence the warning by passing a list to `kind`.

Commit message: TST

Patch to review:
diff --git a/doc/whats_new/v1.1.rst b/doc/whats_new/v1.1.rst
index f755aeba20030..0599141300de7 100644
--- a/doc/whats_new/v1.1.rst
+++ b/doc/whats_new/v1.1.rst
@@ -113,10 +113,10 @@ Changelog
 
 - |Fix| :class:`decomposition.FastICA` now validates input parameters in `fit` instead of `__init__`.
   :pr:`21432` by :user:`Hannah Bohle <hhnnhh>` and :user:`Maren Westermann <marenwestermann>`.
-  
+
 - |Fix| :class:`decomposition.FactorAnalysis` now validates input parameters
   in `fit` instead of `__init__`.
-  :pr:`21713` by :user:`Haya <HayaAlmutairi>` and 
+  :pr:`21713` by :user:`Haya <HayaAlmutairi>` and
   :user:`Krum Arnaudov <krumeto>`.
 
 - |Fix| :class:`decomposition.KernelPCA` now validates input parameters in
@@ -178,6 +178,16 @@ Changelog
   multilabel classification.
   :pr:`19689` by :user:`Guillaume Lemaitre <glemaitre>`.
 
+:mod:`sklearn.inspection`
+.........................
+
+- |Enhancement| In
+  :meth:`~sklearn.inspection.PartialDependenceDisplay.from_estimator` and
+  :meth:`~sklearn.inspection.PartialDependenceDisplay.from_predictions`, allow
+  `kind` to accept a list of strings to specify for each feature interaction
+  which type of plot a user wants to draw.
+  :pr:`19438` by :user:`Guillaume Lemaitre <glemaitre>`.
+
 :mod:`sklearn.linear_model`
 ...........................
 
diff --git a/examples/inspection/plot_partial_dependence.py b/examples/inspection/plot_partial_dependence.py
index 4d22077fdbc09..e281f1e73c88e 100644
--- a/examples/inspection/plot_partial_dependence.py
+++ b/examples/inspection/plot_partial_dependence.py
@@ -108,24 +108,18 @@
 # We will plot the partial dependence, both individual (ICE) and averaged one
 # (PDP). We limit to only 50 ICE curves to not overcrowd the plot.
 
-import matplotlib.pyplot as plt
-from sklearn.inspection import partial_dependence
 from sklearn.inspection import PartialDependenceDisplay
 
+common_params = {"subsample": 50, "n_jobs": 3, "grid_resolution": 20, "random_state": 0}
+
 print("Computing partial dependence plots...")
 tic = time()
-features = ["MedInc", "AveOccup", "HouseAge", "AveRooms"]
 display = PartialDependenceDisplay.from_estimator(
     est,
     X_train,
-    features,
+    features=["MedInc", "AveOccup", "HouseAge", "AveRooms"],
     kind="both",
-    subsample=50,
-    n_jobs=3,
-    grid_resolution=20,
-    random_state=0,
-    ice_lines_kw={"color": "tab:blue", "alpha": 0.2, "linewidth": 0.5},
-    pd_line_kw={"color": "tab:orange", "linestyle": "--"},
+    **common_params,
 )
 print(f"done in {time() - tic:.3f}s")
 display.figure_.suptitle(
@@ -168,14 +162,9 @@
 display = PartialDependenceDisplay.from_estimator(
     est,
     X_train,
-    features,
+    features=["MedInc", "AveOccup", "HouseAge", "AveRooms"],
     kind="both",
-    subsample=50,
-    n_jobs=3,
-    grid_resolution=20,
-    random_state=0,
-    ice_lines_kw={"color": "tab:blue", "alpha": 0.2, "linewidth": 0.5},
-    pd_line_kw={"color": "tab:orange", "linestyle": "--"},
+    **common_params,
 )
 print(f"done in {time() - tic:.3f}s")
 display.figure_.suptitle(
@@ -223,20 +212,23 @@
 # Another consideration is linked to the performance to compute the PDPs. With
 # the tree-based algorithm, when only PDPs are requested, they can be computed
 # on an efficient way using the `'recursion'` method.
+import matplotlib.pyplot as plt
 
-features = ["AveOccup", "HouseAge", ("AveOccup", "HouseAge")]
 print("Computing partial dependence plots...")
 tic = time()
 _, ax = plt.subplots(ncols=3, figsize=(9, 4))
+
+# Note that we could have called the method `from_estimator` three times and
+# provide one feature, one kind of plot, and one axis for each call.
 display = PartialDependenceDisplay.from_estimator(
     est,
     X_train,
-    features,
-    kind="average",
-    n_jobs=3,
-    grid_resolution=20,
+    features=["AveOccup", "HouseAge", ("AveOccup", "HouseAge")],
+    kind=["both", "both", "average"],
     ax=ax,
+    **common_params,
 )
+
 print(f"done in {time() - tic:.3f}s")
 display.figure_.suptitle(
     "Partial dependence of house value on non-location features\n"
@@ -257,9 +249,9 @@
 #
 # Let's make the same partial dependence plot for the 2 features interaction,
 # this time in 3 dimensions.
-
 import numpy as np
 from mpl_toolkits.mplot3d import Axes3D
+from sklearn.inspection import partial_dependence
 
 fig = plt.figure()
 
diff --git a/sklearn/inspection/_plot/partial_dependence.py b/sklearn/inspection/_plot/partial_dependence.py
index e20205e86e51e..33a6340fa6752 100644
--- a/sklearn/inspection/_plot/partial_dependence.py
+++ b/sklearn/inspection/_plot/partial_dependence.py
@@ -232,18 +232,37 @@ def plot_partial_dependence(
 
         .. versionadded:: 0.22
 
-    kind : {'average', 'individual', 'both'}, default='average'
+    kind : {'average', 'individual', 'both'} or list of such str, \
+            default='average'
         Whether to plot the partial dependence averaged across all the samples
         in the dataset or one line per sample or both.
 
         - ``kind='average'`` results in the traditional PD plot;
-        - ``kind='individual'`` results in the ICE plot.
+        - ``kind='individual'`` results in the ICE plot;
+        - ``kind='both'`` results in plotting both the ICE and PD on the same
+          plot.
 
-       Note that the fast ``method='recursion'`` option is only available for
-       ``kind='average'``. Plotting individual dependencies requires using the
-       slower ``method='brute'`` option.
+        A list of such strings can be provided to specify `kind` on a per-plot
+        basis. The length of the list should be the same as the number of
+        interaction requested in `features`.
+
+        .. note::
+           ICE ('individual' or 'both') makes it impossible rendering 2-ways
+           interactions plot. As a result, an error will be raised and a
+           regular 'average' plot will be made instead.
+
+        .. note::
+           The fast ``method='recursion'`` option is only available for
+           ``kind='average'``. Plotting individual dependencies requires using
+           the slower ``method='brute'`` option.
 
         .. versionadded:: 0.24
+           Add `kind` parameter with `'average'`, `'individual'`, and `'both'`
+           options.
+
+        .. versionadded:: 1.1
+           Add the possibility to pass a list of string specifying `kind`
+           for each plot.
 
     subsample : float, int or None, default=1000
         Sampling for ICE curves when `kind` is 'individual' or 'both'.
@@ -371,6 +390,15 @@ def _plot_partial_dependence(
     if len(set(feature_names)) != len(feature_names):
         raise ValueError("feature_names should not contain duplicates.")
 
+    # expand kind to always be a list of str
+    kind_ = [kind] * len(features) if isinstance(kind, str) else kind
+    if len(kind_) != len(features):
+        raise ValueError(
+            "When `kind` is provided as a list of strings, it should contain "
+            f"as many elements as `features`. `kind` contains {len(kind_)} "
+            f"element(s) and `features` contains {len(features)} element(s)."
+        )
+
     def convert_feature(fx):
         if isinstance(fx, str):
             try:
@@ -380,8 +408,8 @@ def convert_feature(fx):
         return int(fx)
 
     # convert features into a seq of int tuples
-    tmp_features = []
-    for fxs in features:
+    tmp_features, ice_for_two_way_pd = [], []
+    for kind_plot, fxs in zip(kind_, features):
         if isinstance(fxs, (numbers.Integral, str)):
             fxs = (fxs,)
         try:
@@ -396,13 +424,26 @@ def convert_feature(fx):
                 "Each entry in features must be either an int, "
                 "a string, or an iterable of size at most 2."
             )
-        if kind != "average" and np.size(fxs) > 1:
-            raise ValueError(
-                "It is not possible to display individual effects for more "
-                f"than one feature at a time. Got: features={features}."
-            )
+        # store the information if 2-way PD was requested with ICE to later
+        # raise a warning
+        ice_for_two_way_pd.append(kind_plot != "average" and np.size(fxs) > 1)
+
         tmp_features.append(fxs)
 
+    if any(ice_for_two_way_pd):
+        # raise an error an be specific regarding the parameter values
+        # when 1- and 2-way PD were requested
+        kind_ = [
+            "average" if forcing_average else kind_plot
+            for forcing_average, kind_plot in zip(ice_for_two_way_pd, kind_)
+        ]
+        raise ValueError(
+            "You requested an ICE plot with 2-way feature interactions. "
+            "This is impossible to render. You need to explicitely request "
+            "a PD plot ('average') for the 2-way feature interactions. "
+            f"Therefore for features={features!r}, you need to provide "
+            f"kind={kind_!r}"
+        )
     features = tmp_features
 
     # Early exit if the axes does not have the correct number of axes
@@ -442,9 +483,9 @@ def convert_feature(fx):
             method=method,
             grid_resolution=grid_resolution,
             percentiles=percentiles,
-            kind=kind,
+            kind=kind_plot,
         )
-        for fxs in features
+        for kind_plot, fxs in zip(kind_, features)
     )
 
     # For multioutput regression, we can only check the validity of target
@@ -455,7 +496,7 @@ def convert_feature(fx):
     pd_result = pd_results[0]  # checking the first result is enough
     n_tasks = (
         pd_result.average.shape[0]
-        if kind == "average"
+        if kind_[0] == "average"
         else pd_result.individual.shape[0]
     )
     if is_regressor(estimator) and n_tasks > 1:
@@ -467,9 +508,9 @@ def convert_feature(fx):
 
     # get global min and max average predictions of PD grouped by plot type
     pdp_lim = {}
-    for pdp in pd_results:
+    for idx, pdp in enumerate(pd_results):
         values = pdp["values"]
-        preds = pdp.average if kind == "average" else pdp.individual
+        preds = pdp.average if kind_[idx] == "average" else pdp.individual
         min_pd = preds[target_idx].min()
         max_pd = preds[target_idx].max()
         n_fx = len(values)
@@ -555,18 +596,37 @@ class PartialDependenceDisplay:
     deciles : dict
         Deciles for feature indices in ``features``.
 
-    kind : {'average', 'individual', 'both'}, default='average'
+    kind : {'average', 'individual', 'both'} or list of such str, \
+            default='average'
         Whether to plot the partial dependence averaged across all the samples
         in the dataset or one line per sample or both.
 
         - ``kind='average'`` results in the traditional PD plot;
-        - ``kind='individual'`` results in the ICE plot.
+        - ``kind='individual'`` results in the ICE plot;
+        - ``kind='both'`` results in plotting both the ICE and PD on the same
+          plot.
+
+        A list of such strings can be provided to specify `kind` on a per-plot
+        basis. The length of the list should be the same as the number of
+        interaction requested in `features`.
 
-       Note that the fast ``method='recursion'`` option is only available for
-       ``kind='average'``. Plotting individual dependencies requires using the
-       slower ``method='brute'`` option.
+        .. note::
+           ICE ('individual' or 'both') makes it impossible rendering 2-ways
+           interactions plot. As a result, an error will be raised and a
+           regular 'average' plot will be made instead.
+
+        .. note::
+           The fast ``method='recursion'`` option is only available for
+           ``kind='average'``. Plotting individual dependencies requires using
+           the slower ``method='brute'`` option.
 
         .. versionadded:: 0.24
+           Add `kind` parameter with `'average'`, `'individual'`, and `'both'`
+           options.
+
+        .. versionadded:: 1.1
+           Add the possibility to pass a list of string specifying `kind`
+           for each plot.
 
     subsample : float, int or None, default=1000
         Sampling for ICE curves when `kind` is 'individual' or 'both'.
@@ -1026,6 +1086,7 @@ def _plot_average_dependence(
 
     def _plot_one_way_partial_dependence(
         self,
+        kind,
         preds,
         avg_preds,
         feature_values,
@@ -1042,6 +1103,8 @@ def _plot_one_way_partial_dependence(
 
         Parameters
         ----------
+        kind : str
+            The kind of partial plot to draw.
         preds : ndarray of shape \
                 (n_instances, n_grid_points) or None
             The predictions computed for all points of `feature_values` for a
@@ -1071,7 +1134,7 @@ def _plot_one_way_partial_dependence(
         """
         from matplotlib import transforms  # noqa
 
-        if self.kind in ("individual", "both"):
+        if kind in ("individual", "both"):
             self._plot_ice_lines(
                 preds[self.target_idx],
                 feature_values,
@@ -1082,9 +1145,9 @@ def _plot_one_way_partial_dependence(
                 ice_lines_kw,
             )
 
-        if self.kind in ("average", "both"):
+        if kind in ("average", "both"):
             # the average is stored as the last line
-            if self.kind == "average":
+            if kind == "average":
                 pd_line_idx = pd_plot_idx
             else:
                 pd_line_idx = pd_plot_idx * n_lines + n_ice_lines
@@ -1119,7 +1182,7 @@ def _plot_one_way_partial_dependence(
         else:
             ax.set_yticklabels([])
 
-        if pd_line_kw.get("label", None) and self.kind != "individual":
+        if pd_line_kw.get("label", None) and kind != "individual":
             ax.legend()
 
     def _plot_two_way_partial_dependence(
@@ -1263,50 +1326,57 @@ def plot(
         import matplotlib.pyplot as plt  # noqa
         from matplotlib.gridspec import GridSpecFromSubplotSpec  # noqa
 
+        if isinstance(self.kind, str):
+            kind = [self.kind] * len(self.features)
+        else:
+            kind = self.kind
+
+        if len(kind) != len(self.features):
+            raise ValueError(
+                "When `kind` is provided as a list of strings, it should "
+                "contain as many elements as `features`. `kind` contains "
+                f"{len(kind)} element(s) and `features` contains "
+                f"{len(self.features)} element(s)."
+            )
+
+        valid_kinds = {"average", "individual", "both"}
+        if any([k not in valid_kinds for k in kind]):
+            raise ValueError(
+                f"Values provided to `kind` must be one of: {valid_kinds!r} or a list"
+                f" of such values. Currently, kind={self.kind}"
+            )
+
         if line_kw is None:
             line_kw = {}
         if ice_lines_kw is None:
             ice_lines_kw = {}
         if pd_line_kw is None:
             pd_line_kw = {}
-        if contour_kw is None:
-            contour_kw = {}
 
         if ax is None:
             _, ax = plt.subplots()
 
+        if contour_kw is None:
+            contour_kw = {}
         default_contour_kws = {"alpha": 0.75}
         contour_kw = {**default_contour_kws, **contour_kw}
 
-        default_line_kws = {
-            "color": "C0",
-            "label": "average" if self.kind == "both" else None,
-        }
-        if self.kind in ("individual", "both"):
-            default_ice_lines_kws = {"alpha": 0.3, "linewidth": 0.5}
-        else:
-            default_ice_lines_kws = {}
-
-        ice_lines_kw = {
-            **default_line_kws,
-            **line_kw,
-            **default_ice_lines_kws,
-            **ice_lines_kw,
-        }
-        del ice_lines_kw["label"]
-
-        pd_line_kw = {**default_line_kws, **line_kw, **pd_line_kw}
-
         n_features = len(self.features)
-        if self.kind in ("individual", "both"):
-            n_ice_lines = self._get_sample_count(len(self.pd_results[0].individual[0]))
-            if self.kind == "individual":
-                n_lines = n_ice_lines
-            else:
-                n_lines = n_ice_lines + 1
-        else:
+        is_average_plot = [kind_plot == "average" for kind_plot in kind]
+        if all(is_average_plot):
+            # only average plots are requested
             n_ice_lines = 0
             n_lines = 1
+        else:
+            # we need to determine the number of ICE samples computed
+            ice_plot_idx = is_average_plot.index(False)
+            n_ice_lines = self._get_sample_count(
+                len(self.pd_results[ice_plot_idx].individual[0])
+            )
+            if any([kind_plot == "both" for kind_plot in kind]):
+                n_lines = n_ice_lines + 1  # account for the average line
+            else:
+                n_lines = n_ice_lines
 
         if isinstance(ax, plt.Axes):
             # If ax was set off, it has most likely been set to off
@@ -1326,7 +1396,7 @@ def plot(
             n_rows = int(np.ceil(n_features / float(n_cols)))
 
             self.axes_ = np.empty((n_rows, n_cols), dtype=object)
-            if self.kind == "average":
+            if all(is_average_plot):
                 self.lines_ = np.empty((n_rows, n_cols), dtype=object)
             else:
                 self.lines_ = np.empty((n_rows, n_cols, n_lines), dtype=object)
@@ -1355,7 +1425,7 @@ def plot(
             self.bounding_ax_ = None
             self.figure_ = ax.ravel()[0].figure
             self.axes_ = ax
-            if self.kind == "average":
+            if all(is_average_plot):
                 self.lines_ = np.empty_like(ax, dtype=object)
             else:
                 self.lines_ = np.empty(ax.shape + (n_lines,), dtype=object)
@@ -1368,22 +1438,62 @@ def plot(
         self.deciles_vlines_ = np.empty_like(self.axes_, dtype=object)
         self.deciles_hlines_ = np.empty_like(self.axes_, dtype=object)
 
-        for pd_plot_idx, (axi, feature_idx, pd_result) in enumerate(
-            zip(self.axes_.ravel(), self.features, self.pd_results)
+        for pd_plot_idx, (axi, feature_idx, pd_result, kind_plot) in enumerate(
+            zip(self.axes_.ravel(), self.features, self.pd_results, kind)
         ):
             avg_preds = None
             preds = None
             feature_values = pd_result["values"]
-            if self.kind == "individual":
+            if kind_plot == "individual":
                 preds = pd_result.individual
-            elif self.kind == "average":
+            elif kind_plot == "average":
                 avg_preds = pd_result.average
-            else:  # kind='both'
+            else:  # kind_plot == 'both'
                 avg_preds = pd_result.average
                 preds = pd_result.individual
 
             if len(feature_values) == 1:
+                # define the line-style for the current plot
+                default_line_kws = {
+                    "color": "C0",
+                    "label": "average" if kind_plot == "both" else None,
+                }
+                if kind_plot in "individual":
+                    default_ice_lines_kws = {"alpha": 0.3, "linewidth": 0.5}
+                    default_pd_lines_kws = {}
+                elif kind_plot in "both":
+                    # by default, we need to distinguish the average line from
+                    # the individual lines via color and line style
+                    default_ice_lines_kws = {
+                        "alpha": 0.3,
+                        "linewidth": 0.5,
+                        "color": "tab:blue",
+                    }
+                    default_pd_lines_kws = {
+                        "color": "tab:orange",
+                        "linestyle": "--",
+                    }
+                else:
+                    default_ice_lines_kws = {}
+                    default_pd_lines_kws = {}
+
+                ice_lines_kw = {
+                    **default_line_kws,
+                    **default_ice_lines_kws,
+                    **line_kw,
+                    **ice_lines_kw,
+                }
+                del ice_lines_kw["label"]
+
+                pd_line_kw = {
+                    **default_line_kws,
+                    **default_pd_lines_kws,
+                    **line_kw,
+                    **pd_line_kw,
+                }
+
                 self._plot_one_way_partial_dependence(
+                    kind_plot,
                     preds,
                     avg_preds,
                     feature_values[0],
diff --git a/sklearn/inspection/_plot/tests/test_plot_partial_dependence.py b/sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
index ed5e4388eef48..d14657856be45 100644
--- a/sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
+++ b/sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
@@ -159,6 +159,8 @@ def test_plot_partial_dependence_kind(
     assert disp.contours_[0, 1] is None
     assert disp.contours_[0, 2] is None
 
+    assert disp.kind == kind
+
 
 @pytest.mark.filterwarnings("ignore:A Bunch will be returned")
 @pytest.mark.parametrize(
@@ -592,13 +594,8 @@ def test_plot_partial_dependence_dataframe(
         ),
         (
             dummy_classification_data,
-            {"features": [(1, 2)], "kind": "individual"},
-            "It is not possible to display individual effects for more than one",
-        ),
-        (
-            dummy_classification_data,
-            {"features": [(1, 2)], "kind": "both"},
-            "It is not possible to display individual effects for more than one",
+            {"features": [1, 2], "kind": ["both"]},
+            "When `kind` is provided as a list of strings, it should contain",
         ),
         (
             dummy_classification_data,
@@ -610,6 +607,16 @@ def test_plot_partial_dependence_dataframe(
             {"features": [1], "subsample": 1.2},
             r"When a floating-point, subsample=1.2 should be in the \(0, 1\) range",
         ),
+        (
+            dummy_classification_data,
+            {"features": [1], "kind": "foo"},
+            "Values provided to `kind` must be one of",
+        ),
+        (
+            dummy_classification_data,
+            {"features": [0, 1], "kind": ["foo", "individual"]},
+            "Values provided to `kind` must be one of",
+        ),
     ],
 )
 def test_plot_partial_dependence_error(
@@ -729,6 +736,71 @@ def test_partial_dependence_overwrite_labels(
             assert legend_text[0].get_text() == label
 
 
+def test_partial_dependence_kind_list(
+    plot_partial_dependence,
+    pyplot,
+    clf_diabetes,
+    diabetes,
+):
+    """Check that we can provide a list of strings to kind parameter."""
+    matplotlib = pytest.importorskip("matplotlib")
+
+    disp = plot_partial_dependence(
+        clf_diabetes,
+        diabetes.data,
+        features=[0, 2, (1, 2)],
+        grid_resolution=20,
+        kind=["both", "both", "average"],
+    )
+
+    for idx in [0, 1]:
+        assert all(
+            [
+                isinstance(line, matplotlib.lines.Line2D)
+                for line in disp.lines_[0, idx].ravel()
+            ]
+        )
+        assert disp.contours_[0, idx] is None
+
+    assert disp.contours_[0, 2] is not None
+    assert all([line is None for line in disp.lines_[0, 2].ravel()])
+
+
+@pytest.mark.parametrize(
+    "features, kind",
+    [
+        ([0, 2, (1, 2)], "individual"),
+        ([0, 2, (1, 2)], "both"),
+        ([(0, 1), (0, 2), (1, 2)], "individual"),
+        ([(0, 1), (0, 2), (1, 2)], "both"),
+        ([0, 2, (1, 2)], ["individual", "individual", "individual"]),
+        ([0, 2, (1, 2)], ["both", "both", "both"]),
+    ],
+)
+def test_partial_dependence_kind_error(
+    plot_partial_dependence,
+    pyplot,
+    clf_diabetes,
+    diabetes,
+    features,
+    kind,
+):
+    """Check that we raise an informative error when 2-way PD is requested
+    together with 1-way PD/ICE"""
+    warn_msg = (
+        "You requested an ICE plot with 2-way feature interactions. "
+        "This is impossible to render. You need to explicitely request "
+    )
+    with pytest.raises(ValueError, match=warn_msg):
+        plot_partial_dependence(
+            clf_diabetes,
+            diabetes.data,
+            features=features,
+            grid_resolution=20,
+            kind=kind,
+        )
+
+
 @pytest.mark.filterwarnings("ignore:A Bunch will be returned")
 @pytest.mark.parametrize(
     "line_kw, pd_line_kw, ice_lines_kw, expected_colors",
@@ -773,7 +845,7 @@ def test_plot_partial_dependence_lines_kw(
     if pd_line_kw is not None and "linestyle" in pd_line_kw:
         assert line.get_linestyle() == pd_line_kw["linestyle"]
     else:
-        assert line.get_linestyle() == "-"
+        assert line.get_linestyle() == "--"
 
     line = disp.lines_[0, 0, 0]
     assert line.get_color() == expected_colors[1]
@@ -781,3 +853,34 @@ def test_plot_partial_dependence_lines_kw(
         assert line.get_linestyle() == ice_lines_kw["linestyle"]
     else:
         assert line.get_linestyle() == "-"
+
+
+def test_partial_dependence_display_wrong_len_kind(
+    pyplot,
+    clf_diabetes,
+    diabetes,
+):
+    """Check that we raised an error when the `kind` is a list and does not have the
+    right length.
+
+    This bug cannot be trigger apart of creating a `PartialDependenceDisplay` without
+    calling `.from_estimator`.
+    """
+    disp = PartialDependenceDisplay.from_estimator(
+        clf_diabetes,
+        diabetes.data,
+        features=[0, 2],
+        grid_resolution=20,
+        kind="average",
+    )
+
+    # alter `kind` to be a list that does not have the same length than the
+    # length of `features`
+    disp.kind = ["average"]
+    err_msg = (
+        r"When `kind` is provided as a list of strings, it should contain as many"
+        r" elements as `features`. `kind` contains 1 element\(s\) and `features`"
+        r" contains 2 element\(s\)."
+    )
+    with pytest.raises(ValueError, match=err_msg):
+        disp.plot()


Code review comments:
diff_hunk: @@ -232,18 +232,37 @@ def plot_partial_dependence(
 
         .. versionadded:: 0.22
 
-    kind : {'average', 'individual', 'both'}, default='average'
+    kind : {'average', 'individual', 'both'} or list of such str, \
+            default='average'
         Whether to plot the partial dependence averaged across all the samples
         in the dataset or one line per sample or both.
 
         - ``kind='average'`` results in the traditional PD plot;
-        - ``kind='individual'`` results in the ICE plot.
+        - ``kind='individual'`` results in the ICE plot;
+        - ``kind='both'`` results in plotting both the ICE and PD on the same
+          plot.
 
-       Note that the fast ``method='recursion'`` option is only available for
-       ``kind='average'``. Plotting individual dependencies requires using the
-       slower ``method='brute'`` option.
+        A list of such strings can be provided to specify `kind` on a per-plot
+        basis. The length of the list should be the same as the number of
+        interaction requested in `features`.
+
+        .. note::
+           ICE ('individual' or 'both') makes it impossible rendering 2-ways
+           interactions plot. As a result, an error will be raised and a
+           regular 'average' plot will be made instead.
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
If an error is raised (as opposed to a warning as previously), the "average" plot will no longer be used automatically anymore. This phrasing needs to be adjusted:

```suggestion
           ICE ('individual' or 'both') is not a valid option for 2-ways
           interactions plot. As a result, an error will be raised.
           2-ways interaction plots should always be configured to
           use the 'average' kind instead.
```

diff_hunk: @@ -555,18 +596,37 @@ class PartialDependenceDisplay:
     deciles : dict
         Deciles for feature indices in ``features``.
 
-    kind : {'average', 'individual', 'both'}, default='average'
+    kind : {'average', 'individual', 'both'} or list of such str, \
+            default='average'
         Whether to plot the partial dependence averaged across all the samples
         in the dataset or one line per sample or both.
 
         - ``kind='average'`` results in the traditional PD plot;
-        - ``kind='individual'`` results in the ICE plot.
+        - ``kind='individual'`` results in the ICE plot;
+        - ``kind='both'`` results in plotting both the ICE and PD on the same
+          plot.
+
+        A list of such strings can be provided to specify `kind` on a per-plot
+        basis. The length of the list should be the same as the number of
+        interaction requested in `features`.
 
-       Note that the fast ``method='recursion'`` option is only available for
-       ``kind='average'``. Plotting individual dependencies requires using the
-       slower ``method='brute'`` option.
+        .. note::
+           ICE ('individual' or 'both') makes it impossible rendering 2-ways
+           interactions plot. As a result, an error will be raised and a
+           regular 'average' plot will be made instead.
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
```suggestion
           ICE ('individual' or 'both') is not a valid option for 2-ways
           interactions plot. As a result, an error will be raised.
           2-ways interaction plots should always be configured to
           use the 'average' kind instead.
```

diff_hunk: @@ -729,6 +736,71 @@ def test_partial_dependence_overwrite_labels(
             assert legend_text[0].get_text() == label
 
 
+def test_partial_dependence_kind_list(
+    plot_partial_dependence,
+    pyplot,
+    clf_diabetes,
+    diabetes,
+):
+    """Check that we can provide a list of strings to kind parameter."""
+    matplotlib = pytest.importorskip("matplotlib")
+
+    disp = plot_partial_dependence(
+        clf_diabetes,
+        diabetes.data,
+        features=[0, 2, (1, 2)],
+        grid_resolution=20,
+        kind=["both", "both", "average"],
+    )
+
+    for idx in [0, 1]:
+        assert all(
+            [
+                isinstance(line, matplotlib.lines.Line2D)
+                for line in disp.lines_[0, idx].ravel()
+            ]
+        )
+        assert disp.contours_[0, idx] is None
+
+    assert disp.contours_[0, 2] is not None
+    assert all([line is None for line in disp.lines_[0, 2].ravel()])
+
+
+@pytest.mark.parametrize(
+    "features, kind",
+    [
+        ([0, 2, (1, 2)], "individual"),
+        ([0, 2, (1, 2)], "both"),
+        ([(0, 1), (0, 2), (1, 2)], "individual"),
+        ([(0, 1), (0, 2), (1, 2)], "both"),
+        ([0, 2, (1, 2)], ["individual", "individual", "individual"]),
+        ([0, 2, (1, 2)], ["both", "both", "both"]),
+    ],
+)
+def test_partial_dependence_kind_error(
+    plot_partial_dependence,
+    pyplot,
+    clf_diabetes,
+    diabetes,
+    features,
+    kind,
+):
+    """Check that we raise an informative error when 2-way PD is requested
+    together with 1-way PD/ICE"""
+    warn_msg = (
+        "You requested an ICE plot with 2-way feature interactions. "
+        "This is impossible to render. You need to explicitely request "
path: sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
text: @user1:
This needs would need to be adapted accordingly.

```suggestion
        "ICE plot cannot be rendered for 2-way feature interactions. 2-way "
        "feature interactions mandates PD plots using the 'average' kind"
```

I am not sure if regexp escaping is needed and if the suggestion above is black formatted though.

diff_hunk: @@ -178,6 +178,16 @@ Changelog
   multilabel classification.
   :pr:`19689` by :user:`Guillaume Lemaitre <glemaitre>`.
 
+:mod:`sklearn.inspection`
+.........................
+
+- |Enhancement| In
+  :meth:`~sklearn.inspection.PartialDependenceDisplay.from_estimator` and
+  :meth:`~sklearn.inspection.PartialDependenceDisplay.from_predictions`, allow
+  `kind` to accept a list of strings to specify for each feature interaction
+  which type of plot a user wants to draw.
path: doc/whats_new/v1.1.rst
text: @user1:
```suggestion
  `kind` to accept a list of strings to specify  which type of
   plot to draw for each feature interaction.
```

diff_hunk: @@ -1263,50 +1326,57 @@ def plot(
         import matplotlib.pyplot as plt  # noqa
         from matplotlib.gridspec import GridSpecFromSubplotSpec  # noqa
 
+        if isinstance(self.kind, str):
+            kind = [self.kind] * len(self.features)
+        else:
+            kind = self.kind
+
+        if len(kind) != len(self.features):
+            raise ValueError(
+                "When `kind` is provided as a list of strings, it should "
+                "contain as many elements as `features`. `kind` contains "
+                f"{len(kind)} element(s) and `features` contains "
+                f"{len(self.features)} element(s)."
+            )
+
+        valid_kinds = {"average", "individual", "both"}
+        if any([k not in valid_kinds for k in kind]):
+            raise ValueError(
+                f"Values provided to `kind` must be one of: {valid_kinds!r} or a list"
+                f" of such values. Currently, kind={self.kind}"
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
```suggestion
                f" of such values. Currently, kind={self.kind!r}"
```

diff_hunk: @@ -396,13 +424,26 @@ def convert_feature(fx):
                 "Each entry in features must be either an int, "
                 "a string, or an iterable of size at most 2."
             )
-        if kind != "average" and np.size(fxs) > 1:
-            raise ValueError(
-                "It is not possible to display individual effects for more "
-                f"than one feature at a time. Got: features={features}."
-            )
+        # store the information if 2-way PD was requested with ICE to later
+        # raise a warning
+        ice_for_two_way_pd.append(kind_plot != "average" and np.size(fxs) > 1)
+
         tmp_features.append(fxs)
 
+    if any(ice_for_two_way_pd):
+        # raise an error an be specific regarding the parameter values
+        # when 1- and 2-way PD were requested
+        kind_ = [
+            "average" if forcing_average else kind_plot
+            for forcing_average, kind_plot in zip(ice_for_two_way_pd, kind_)
+        ]
+        raise ValueError(
+            "You requested an ICE plot with 2-way feature interactions. "
+            "This is impossible to render. You need to explicitely request "
+            "a PD plot ('average') for the 2-way feature interactions. "
+            f"Therefore for features={features!r}, you need to provide "
+            f"kind={kind_!r}"
path: sklearn/inspection/_plot/partial_dependence.py
text: @user1:
```suggestion
            "ICE plot cannot be rendered for 2-way feature interactions. "
            "2-way feature interactions mandates PD plots using the "
            "'average' kind: "
            f"features={features!r} should be configured to use "
            f"kind={kind_!r} explicitly."
```

diff_hunk: @@ -773,11 +845,42 @@ def test_plot_partial_dependence_lines_kw(
     if pd_line_kw is not None and "linestyle" in pd_line_kw:
         assert line.get_linestyle() == pd_line_kw["linestyle"]
     else:
-        assert line.get_linestyle() == "-"
+        assert line.get_linestyle() == "--"
 
     line = disp.lines_[0, 0, 0]
     assert line.get_color() == expected_colors[1]
     if ice_lines_kw is not None and "linestyle" in ice_lines_kw:
         assert line.get_linestyle() == ice_lines_kw["linestyle"]
     else:
         assert line.get_linestyle() == "-"
+
+
+def test_partial_dependence_display_wrong_len_kind(
+    pyplot,
+    clf_diabetes,
+    diabetes,
+):
+    """Check that we raised an error when the `kind` is a list and does not have the
+    right length.
+
+    This bug cannot be trigger apart of creating a `PartialDependenceDisplay` without
+    calling `.from_estimator`.
+    """
+    disp = PartialDependenceDisplay.from_estimator(
+        clf_diabetes,
+        diabetes.data,
+        features=[0, 2],
+        grid_resolution=20,
+        kind="average",
path: sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
text: @user1:
Why not pass a list with the wrong length directly instead of changing it later?

```suggestion
        kind=["average"],  # len(kind) != len(features)
```

diff_hunk: @@ -773,11 +845,42 @@ def test_plot_partial_dependence_lines_kw(
     if pd_line_kw is not None and "linestyle" in pd_line_kw:
         assert line.get_linestyle() == pd_line_kw["linestyle"]
     else:
-        assert line.get_linestyle() == "-"
+        assert line.get_linestyle() == "--"
 
     line = disp.lines_[0, 0, 0]
     assert line.get_color() == expected_colors[1]
     if ice_lines_kw is not None and "linestyle" in ice_lines_kw:
         assert line.get_linestyle() == ice_lines_kw["linestyle"]
     else:
         assert line.get_linestyle() == "-"
+
+
+def test_partial_dependence_display_wrong_len_kind(
+    pyplot,
+    clf_diabetes,
+    diabetes,
+):
+    """Check that we raised an error when the `kind` is a list and does not have the
+    right length.
+
+    This bug cannot be trigger apart of creating a `PartialDependenceDisplay` without
+    calling `.from_estimator`.
+    """
+    disp = PartialDependenceDisplay.from_estimator(
+        clf_diabetes,
+        diabetes.data,
+        features=[0, 2],
+        grid_resolution=20,
+        kind="average",
+    )
+
+    # alter `kind` to be a list that does not have the same length than the
+    # length of `features`
path: sklearn/inspection/_plot/tests/test_plot_partial_dependence.py
text: @user1:
```suggestion
    # alter `kind` to be a list with a length different from length of `features`
```