Allow boosting of estimators in scikit-learn pipelines using residuals

### Describe the workflow you want to enable

I want to be able to use multiple estimators in one pipeline.  E.g.

```python
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.datasets import make_regression

X, y = make_regression(n_samples=100, n_features=10, noise=0.1)

pipeline = Pipeline([
    ('linear', LinearRegression()),
    ('gbm', GradientBoostingRegressor())
])

pipeline.fit(X, y)
```

This currently doesn't work: 
```python
TypeError: All intermediate steps should be transformers and implement fit and transform or be the string 'passthrough' 'LinearRegression()' (type <class 'sklearn.linear_model._base.LinearRegression'>) doesn't
```

This is for good reason: naively passing predict from the linear model to the fit for the gbm will lead to overfitting!

### Describe your proposed solution

Allow boosting estimators in scikit learn pipeline by subtracting the predict output from y at `fit` time and then adding the predict outputs at transform time. This enhancement would allow multiple regression estimators to be used consecutively, where each subsequent estimator fits on the residuals of the previous ones.

In my example above, at fit time the linear step proceeds as normal.  Then the gbm step detects 2 chained estimators, and subtracts the prediction from the linear step from y, and then fits to the residual.

At transform time, the predictions from the linear model would be added to the predictions from the gbm model.

To simplify things, I think you still wouldn't be allowed to connect an estimator to a transformer.

### Describe alternatives you've considered, if relevant

Vecstack is a good alternative, but its slow.  You have to cross-validate each estimator, which gets expensive if you chain more than 2 in a row.

```python
from sklearn.pipeline import Pipeline
from sklearn.linear_model import LinearRegression, Ridge
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.datasets import make_regression
from vecstack import StackingTransformer

X, y = make_regression(n_samples=100, n_features=10, noise=0.1)

pipeline = Pipeline([
    ('linear', StackingTransformer([('lr', LinearRegression())], regression=True, shuffle=True)),
    ('gbm', GradientBoostingRegressor())
])

pipeline.fit(X, y)
y_pred = pipeline.predict(X)
```

### Additional context

I love being able to use stacking inside sklearn pipelines.  It would be very cool to be able to use boosting too.  I haven't been able to figure out how to make a boosting class that adjusts y as it moves through the pipeline, which is why I think this feature would need to be implemented in the core library.

COMMENT:
> @a-wozniakowski That works if I want to only boost one estimator, and use the same preprocessing for each one. It does not work if I want to boost _different_ estimators, with varied preprocessing.
> 
> For example, I would like to have a pipeline that starts with a linear regression model that uses one-hot encoding. I would then like to git to add a Decision Tree that uses ordinal encoding, fit on the residuals of the linear model.

Is the following implementation close to the example, namely:
```
from sklearn.base import BaseEstimator
from sklearn.datasets import make_regression
from sklearn.ensemble import GradientBoostingRegressor
from sklearn.linear_model import LinearRegression
from sklearn.metrics import mean_squared_error
from sklearn.pipeline import FeatureUnion, Pipeline
from sklearn.preprocessing import StandardScaler, QuantileTransformer

X, y = make_regression(n_samples=100, n_features=10, noise=10.0, random_state=42)

class Regressor(BaseEstimator):
    def __init__(self, n_features):
        self.n_features = n_features
        self.linear_reg = LinearRegression()
        self.boosting_reg = GradientBoostingRegressor(
            n_estimators=1,
            init="zero",
            random_state=42
        )

    def fit(self, X, y):
        self.linear_reg.fit(X=X[:, :self.n_features], y=y)
        y_pred = self.linear_reg.predict(X=X[:, :self.n_features])
        residual = y - y_pred
        self.boosting_reg.fit(X=X[:, self.n_features:], y=residual)
        return self

    def predict(self, X):
        y_pred_linear_reg = self.linear_reg.predict(X=X[:, :self.n_features])
        residual = self.boosting_reg.predict(X=X[:, self.n_features:])
        y_pred = y_pred_linear_reg + residual
        return y_pred_linear_reg, y_pred


union = FeatureUnion(
    transformer_list=[("ss", StandardScaler()),
                      ("qt", QuantileTransformer(n_quantiles=2))]
)

X_union = union.fit_transform(X=X)

n_features = 10
pipe = Pipeline(
    steps=[("reg", Regressor(n_features=n_features))]
).fit(X=X_union, y=y)
y_pred_linear_reg, y_pred = pipe.predict(X=X_union)
print(mean_squared_error(y_true=y, y_pred=y_pred_linear_reg))
print(mean_squared_error(y_true=y, y_pred=y_pred))

reg = GradientBoostingRegressor(
    n_estimators=65,
    init=None,
    random_state=42
).fit(X=X_union, y=y)
y_pred = reg.predict(X=X_union)
print(mean_squared_error(y_true=y, y_pred=y_pred))
```
where the print statements result in:
```
92.6124323169397
89.77572408088261
89.90872773258083
```

Note that the preprocessing was done with `StandardScaler` and `QuantileTransformer`, but those choices were arbitrary. Also, I put 65 decision trees in the latter as that resulted in a roughly similar training error as the twicing method, which was implemented in `Regressor`.

Labels: New Feature, Needs Decision