Handle np.nan / missing values in SplineTransformer

I think it would be quite natural to add an option to `SplineTransformer` to accept inputs with missing values as follows:

- `handle_missing="error"` the default (keep current behavior)
- `handle_missing="zero"`/"constant": encode missing values by setting all output features for that input column to 0 (or some other constant, see discussion below),
- `handle_missing="indicator"`: append an extra binary feature as missingness indicator and encode missing values as 0 on the remaining output features.

Note that `handle_missing="indicator"` would be different and statistically more meaningful than  `SimpleImputer(strategy="mean", add_indicator=True)` with `SplineTransformer` and furthermore would make for leaner ML pipelines (better UX).

I am not sure if we need to add the `handle_missing="zero"` option. It would break the property the sum of output values of a given `SplineTransformer` encoding always sum to 1 while `handle_missing="indicator"` would preserve this property (in addition to make missingness more explicit to the downstream model in case missingness is informative one way or another).

If we want to preserve the sum to 1 property while not adding an explicit missingness indicator feature, maybe we could instead provide `handle_missing="constant"` (not sure about the name) that would encode missing values as `1 / n_outputs` to preserve the "sum to 1" property. Not entirely sure if this would result in a more interesting prior than the zero encoding.

COMMENT:
Here is an example.

Let's consider a classification problem where the input feature has some completely at random missing values and the target is only statistically associated with the missingness pattern but is fundamentally independent on the underlying actual values.

```python
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd

from sklearn.preprocessing import SplineTransformer
from sklearn.linear_model import LogisticRegression
from sklearn.impute import SimpleImputer
from sklearn.pipeline import make_pipeline
from sklearn.inspection import permutation_importance

rng = np.random.default_rng(0)

X = rng.uniform(size=(100_000, 1))

# add some missing values at random, independently on the original non
# missing values:
missing_mask = rng.uniform(size=X.shape) > 0.99
X[missing_mask] = np.nan
X = pd.DataFrame(X, columns=["x"])

# target is only based on the presence or absence of missing values, not on the
# actual value.
y = missing_mask.any(axis=1)
```

By making the spline transformer able to encode missingness by itself, it's quite easy to fit downstream model with builtin feature selection (e.g. a linear classifier with an elastic net penalty) that would naturally only pickup that missing indicator.

```python
n_knots = 15

lr = LogisticRegression(
    C=1e-2,
    penalty="elasticnet",
    solver="saga",
    l1_ratio=0.5,
    max_iter=1000,
    random_state=0,
)

spline_indicator = make_pipeline(
    SplineTransformer(n_knots=n_knots, handle_missing="indicator"),
    lr,
).fit(X, y)


coef = spline_indicator[-1].coef_.ravel()
plt.bar(np.arange(len(coef)), coef)
```
![image](https://github.com/scikit-learn/scikit-learn/assets/89061/f7143fba-1b01-4afb-ae34-fe2fc2e43edb)

Note that only the last coef (the one for the missingness indicator) is non-zero.

However, if you use the `SimpleImputer` instead you have two problem:

- first, you need a more complex pipeline to avoid spline encoding the missingness indicator feature generated by the previous step. I agree that this point could probably be fixed if we had a better pipeline API.
- second, and this is the main motivation for my issue, the imputation by the mean would not be neutral: it would collapse missingness information with the spline feature(s) that happen to be non-zero for the empirical mean of the non-missing values of that column. As a result a downstream model with builtin feature selection would be confused and lead to a suboptimal model:

```python
from sklearn.compose import ColumnTransformer

impute_indicator_then_spline = make_pipeline(
    SimpleImputer(strategy="mean", add_indicator=True),
    ColumnTransformer(
        [
            ("spline", SplineTransformer(n_knots=n_knots), [0]),
        ],
        remainder="passthrough",
    ),
    lr,
).fit(X, y)


coef = impute_indicator_then_spline[-1].coef_.ravel()
plt.bar(np.arange(len(coef)), coef);
```
![image](https://github.com/scikit-learn/scikit-learn/assets/89061/b339f2ef-130e-41ea-b685-eb5ba4e175cb)

Here you can see that elastic net would naturally pickup two features:
- the last one (the missingness indicator of the first stage of the pipeline),
- the spline feature that best encode the mean of the non-missing values (because of imputationby the mean in the first step of the pipeline).

Permutation importance confirm the problem:

<details>

- for the first model, only the expected feature is important:

```python
X_transformed = spline_indicator[:-1].transform(X)

result = permutation_importance(
    spline_indicator[-1],
    X_transformed,
    y,
    n_repeats=10,
    random_state=42,
    n_jobs=2,
    scoring="neg_log_loss",
)
sorted_importances_idx = result.importances_mean.argsort()
importances = pd.DataFrame(
    result.importances[sorted_importances_idx].T,
    columns=np.arange(len(coef))[sorted_importances_idx],
)
ax = importances.plot.box(vert=False, whis=10)
ax.set_title("Permutation Importances")
ax.axvline(x=0, color="k", linestyle="--")
ax.set_xlabel("Increase in log-loss")
ax.figure.tight_layout()
```
![image](https://github.com/scikit-learn/scikit-learn/assets/89061/75018af4-ca47-493f-9a5c-5c06b7c381cc)


- the second model also attributes some importance to the spline of the mean:

```python
X_transformed = impute_indicator_then_spline[:-1].transform(X)

result = permutation_importance(
    impute_indicator_then_spline[-1],
    X_transformed,
    y,
    n_repeats=10,
    random_state=42,
    n_jobs=2,
    scoring="neg_log_loss",
)
sorted_importances_idx = result.importances_mean.argsort()
importances = pd.DataFrame(
    result.importances[sorted_importances_idx].T,
    columns=np.arange(len(coef))[sorted_importances_idx],
)
ax = importances.plot.box(vert=False, whis=10)
ax.set_title("Permutation Importances")
ax.axvline(x=0, color="k", linestyle="--")
ax.set_xlabel("Increase in log-loss")
ax.figure.tight_layout()
```
![image](https://github.com/scikit-learn/scikit-learn/assets/89061/6a88f688-89c9-40e3-a396-fc7810b89e2c)

</details>

So in conclusion, even if we have a better pipeline API, it would not be possible to build the best model by composing `SimpleImputer(strategy="mean", add_indicator=True)` with our current `SplineTransformer`, hence my expresssed need for #28043.

Labels: Enhancement