Repository: alteryx/evalml (python)

Issue:
Add support for additional estimators for multiseries problems
Currently, our model performance is relatively limited with just the VARMAX regressor. We can generate better results if we included additional non-time-series-native estimators for multiseries problems.

Pull request description:
Added support for additional estimators for multiseries datasets
Resolves #4386 

Patch:
diff --git a/docs/source/release_notes.rst b/docs/source/release_notes.rst
index 4dda45a7..994d705a 100644
--- a/docs/source/release_notes.rst
+++ b/docs/source/release_notes.rst
@@ -2,6 +2,7 @@ Release Notes
 -------------
 **Future Releases**
     * Enhancements
+        * Added support for additional estimators for multiseries datasets :pr:`4385`
     * Fixes
         * Fixed bug in `_downcast_nullable_y` causing woodwork initialization issues :pr:`4369`
         * Fixed multiseries prediction interval labels :pr:`4377`
diff --git a/evalml/pipelines/components/estimators/regressors/catboost_regressor.py b/evalml/pipelines/components/estimators/regressors/catboost_regressor.py
index 1a5945e3..a1d854f4 100644
--- a/evalml/pipelines/components/estimators/regressors/catboost_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/catboost_regressor.py
@@ -47,10 +47,12 @@ class CatBoostRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(
diff --git a/evalml/pipelines/components/estimators/regressors/decision_tree_regressor.py b/evalml/pipelines/components/estimators/regressors/decision_tree_regressor.py
index 5bd70668..df6af93e 100644
--- a/evalml/pipelines/components/estimators/regressors/decision_tree_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/decision_tree_regressor.py
@@ -55,10 +55,12 @@ class DecisionTreeRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(
diff --git a/evalml/pipelines/components/estimators/regressors/elasticnet_regressor.py b/evalml/pipelines/components/estimators/regressors/elasticnet_regressor.py
index 417fffc5..8cab2ad1 100644
--- a/evalml/pipelines/components/estimators/regressors/elasticnet_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/elasticnet_regressor.py
@@ -33,10 +33,12 @@ class ElasticNetRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(
diff --git a/evalml/pipelines/components/estimators/regressors/et_regressor.py b/evalml/pipelines/components/estimators/regressors/et_regressor.py
index c5991db8..81673715 100644
--- a/evalml/pipelines/components/estimators/regressors/et_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/et_regressor.py
@@ -56,10 +56,12 @@ class ExtraTreesRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(
diff --git a/evalml/pipelines/components/estimators/regressors/lightgbm_regressor.py b/evalml/pipelines/components/estimators/regressors/lightgbm_regressor.py
index 7070f09e..d371267d 100644
--- a/evalml/pipelines/components/estimators/regressors/lightgbm_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/lightgbm_regressor.py
@@ -68,7 +68,10 @@ class LightGBMRegressor(Estimator):
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
     ]
-    """[ProblemTypes.REGRESSION]"""
+    """[
+        ProblemTypes.REGRESSION,
+        ProblemTypes.TIME_SERIES_REGRESSION,
+    ]"""
 
     SEED_MIN = 0
     SEED_MAX = SEED_BOUNDS.max_bound
diff --git a/evalml/pipelines/components/estimators/regressors/linear_regressor.py b/evalml/pipelines/components/estimators/regressors/linear_regressor.py
index f255761d..15c95d46 100644
--- a/evalml/pipelines/components/estimators/regressors/linear_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/linear_regressor.py
@@ -28,10 +28,12 @@ class LinearRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(self, fit_intercept=True, n_jobs=-1, random_seed=0, **kwargs):
diff --git a/evalml/pipelines/components/estimators/regressors/rf_regressor.py b/evalml/pipelines/components/estimators/regressors/rf_regressor.py
index 3a2939ff..5ab13317 100644
--- a/evalml/pipelines/components/estimators/regressors/rf_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/rf_regressor.py
@@ -37,10 +37,12 @@ class RandomForestRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     def __init__(
diff --git a/evalml/pipelines/components/estimators/regressors/xgboost_regressor.py b/evalml/pipelines/components/estimators/regressors/xgboost_regressor.py
index 3a63a9f5..3f06e431 100644
--- a/evalml/pipelines/components/estimators/regressors/xgboost_regressor.py
+++ b/evalml/pipelines/components/estimators/regressors/xgboost_regressor.py
@@ -40,10 +40,12 @@ class XGBoostRegressor(Estimator):
     supported_problem_types = [
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]
     """[
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     ]"""
 
     # xgboost supports seeds from -2**31 to 2**31 - 1 inclusive. these limits ensure the random seed generated below
diff --git a/evalml/pipelines/components/transformers/preprocessing/drop_nan_rows_transformer.py b/evalml/pipelines/components/transformers/preprocessing/drop_nan_rows_transformer.py
index 3c52c164..b71d4eb9 100644
--- a/evalml/pipelines/components/transformers/preprocessing/drop_nan_rows_transformer.py
+++ b/evalml/pipelines/components/transformers/preprocessing/drop_nan_rows_transformer.py
@@ -1,4 +1,5 @@
 """Transformer to drop rows specified by row indices."""
+import pandas as pd
 from woodwork import init_series
 
 from evalml.pipelines.components.transformers import Transformer
@@ -43,12 +44,24 @@ class DropNaNRowsTransformer(Transformer):
         y_t = infer_feature_types(y) if y is not None else None
 
         X_t_schema = X_t.ww.schema
+        y_t_logical = None
+        y_t_semantic = None
         if y_t is not None:
-            y_t_logical = y_t.ww.logical_type
+            if isinstance(y_t, pd.DataFrame):
+                y_t_logical = y_t.ww.logical_types
+            else:
+                y_t_logical = y_t.ww.logical_type
             y_t_semantic = y_t.ww.semantic_tags
 
         X_t, y_t = drop_rows_with_nans(X_t, y_t)
         X_t.ww.init_with_full_schema(X_t_schema)
         if y_t is not None:
-            y_t = init_series(y_t, logical_type=y_t_logical, semantic_tags=y_t_semantic)
+            if isinstance(y_t, pd.DataFrame):
+                y_t.ww.init(logical_types=y_t_logical, semantic_tags=y_t_semantic)
+            else:
+                y_t = init_series(
+                    y_t,
+                    logical_type=y_t_logical,
+                    semantic_tags=y_t_semantic,
+                )
         return X_t, y_t
diff --git a/evalml/pipelines/components/transformers/preprocessing/time_series_featurizer.py b/evalml/pipelines/components/transformers/preprocessing/time_series_featurizer.py
index f8124710..bbd35cc1 100644
--- a/evalml/pipelines/components/transformers/preprocessing/time_series_featurizer.py
+++ b/evalml/pipelines/components/transformers/preprocessing/time_series_featurizer.py
@@ -127,16 +127,32 @@ class TimeSeriesFeaturizer(Transformer):
         if self.time_index is None:
             raise ValueError("time_index cannot be None!")
 
-        # For the multiseries case, where we only want the start delay lag for the baseline
-        if isinstance(y, pd.DataFrame):
-            self.statistically_significant_lags = [self.start_delay]
-        else:
-            self.statistically_significant_lags = self._find_significant_lags(
-                y,
-                conf_level=self.conf_level,
-                start_delay=self.start_delay,
-                max_delay=self.max_delay,
+        if y is None:
+            # Set lags to all possible lag values
+            self.statistically_significant_lags = np.arange(
+                self.start_delay,
+                self.start_delay + self.max_delay + 1,
             )
+        else:
+            # For the multiseries case, each series ID has individualized lag values
+            if isinstance(y, pd.Series) or isinstance(y, np.ndarray):
+                y = pd.DataFrame(y)
+
+            self.statistically_significant_lags = {}
+            for column in y.columns:
+                self.statistically_significant_lags[
+                    column
+                ] = self._find_significant_lags(
+                    y[column],
+                    conf_level=self.conf_level,
+                    start_delay=self.start_delay,
+                    max_delay=self.max_delay,
+                )
+                if len(y.columns) == 1:
+                    self.statistically_significant_lags = (
+                        self.statistically_significant_lags[column]
+                    )
+                    return self
         return self
 
     @staticmethod
@@ -160,31 +176,28 @@ class TimeSeriesFeaturizer(Transformer):
     @staticmethod
     def _find_significant_lags(y, conf_level, start_delay, max_delay):
         all_lags = np.arange(start_delay, start_delay + max_delay + 1)
-        if y is not None:
-            # Compute the acf and find its peaks
-            acf_values, ci_intervals = acf(
-                y,
-                nlags=len(y) - 1,
-                fft=True,
-                alpha=conf_level,
-            )
-            peaks, _ = find_peaks(acf_values)
-            # Significant lags are the union of:
-            # 1. the peaks (local maxima) that are significant
-            # 2. The significant lags among the first 10 lags.
-            # We then filter the list to be in the range [start_delay, start_delay + max_delay]
-            index = np.arange(len(acf_values))
-            significant = np.logical_or(ci_intervals[:, 0] > 0, ci_intervals[:, 1] < 0)
-            first_significant_10 = index[:10][significant[:10]]
-            significant_lags = (
-                set(index[significant]).intersection(peaks).union(first_significant_10)
-            )
-            # If no lags are significant get the first lag
-            significant_lags = sorted(significant_lags.intersection(all_lags)) or [
-                start_delay,
-            ]
-        else:
-            significant_lags = all_lags
+        # Compute the acf and find its peaks
+        acf_values, ci_intervals = acf(
+            y,
+            nlags=len(y) - 1,
+            fft=True,
+            alpha=conf_level,
+        )
+        peaks, _ = find_peaks(acf_values)
+        # Significant lags are the union of:
+        # 1. the peaks (local maxima) that are significant
+        # 2. The significant lags among the first 10 lags.
+        # We then filter the list to be in the range [start_delay, start_delay + max_delay]
+        index = np.arange(len(acf_values))
+        significant = np.logical_or(ci_intervals[:, 0] > 0, ci_intervals[:, 1] < 0)
+        first_significant_10 = index[:10][significant[:10]]
+        significant_lags = (
+            set(index[significant]).intersection(peaks).union(first_significant_10)
+        )
+        # If no lags are significant get the first lag
+        significant_lags = sorted(significant_lags.intersection(all_lags)) or [
+            start_delay,
+        ]
         return significant_lags
 
     def _compute_rolling_transforms(self, X, y, original_features):
@@ -234,7 +247,25 @@ class TimeSeriesFeaturizer(Transformer):
             col = data[col_name]
             if categorical_columns and col_name in categorical_columns:
                 col = X_categorical[col_name]
-            for t in self.statistically_significant_lags:
+            # Lags are stored in a dict for multiseries problems
+            # Returns the lags corresponding to the series ID value
+            if isinstance(self.statistically_significant_lags, dict):
+                from evalml.pipelines.utils import MULTISERIES_SEPARATOR_SYMBOL
+
+                col_series_id = (
+                    MULTISERIES_SEPARATOR_SYMBOL
+                    + col_name.split(MULTISERIES_SEPARATOR_SYMBOL)[-1]
+                )
+                for (
+                    series_id_target_name,
+                    lag_list,
+                ) in self.statistically_significant_lags.items():
+                    if series_id_target_name.endswith(col_series_id):
+                        lags = lag_list
+                        break
+            else:
+                lags = self.statistically_significant_lags
+            for t in lags:
                 lagged_features[self.df_colname_prefix.format(col_name, t)] = col.shift(
                     t,
                 )
diff --git a/evalml/pipelines/multiseries_regression_pipeline.py b/evalml/pipelines/multiseries_regression_pipeline.py
index 948ce040..df9bce66 100644
--- a/evalml/pipelines/multiseries_regression_pipeline.py
+++ b/evalml/pipelines/multiseries_regression_pipeline.py
@@ -83,6 +83,7 @@ class MultiseriesRegressionPipeline(TimeSeriesRegressionPipeline):
 
         self.component_graph.fit(X_unstacked, y_unstacked)
         self.input_feature_names = self.component_graph.input_feature_names
+        self.series_id_target_names = y_unstacked.columns
 
     def predict_in_sample(
         self,
@@ -144,7 +145,7 @@ class MultiseriesRegressionPipeline(TimeSeriesRegressionPipeline):
         ]
         y_overlapping_features = [
             feature
-            for feature in y_train_unstacked.columns
+            for feature in self.series_id_target_names
             if feature in y_unstacked.columns
         ]
         y_unstacked = y_unstacked[y_overlapping_features]
@@ -154,7 +155,6 @@ class MultiseriesRegressionPipeline(TimeSeriesRegressionPipeline):
         y_train_unstacked = infer_feature_types(y_train_unstacked)
         X_unstacked = infer_feature_types(X_unstacked)
         y_unstacked = infer_feature_types(y_unstacked)
-
         unstacked_predictions = super().predict_in_sample(
             X_unstacked,
             y_unstacked,
@@ -163,16 +163,50 @@ class MultiseriesRegressionPipeline(TimeSeriesRegressionPipeline):
             objective,
             calculating_residuals,
         )
+        unstacked_predictions = unstacked_predictions[
+            [
+                series_id_target
+                for series_id_target in self.series_id_target_names
+                if series_id_target in unstacked_predictions.columns
+            ]
+        ]
+
+        # Add `time_index` column to index for generating stacked datetime column in `stack_data()`
+        unstacked_predictions.index = X_unstacked[self.time_index]
         stacked_predictions = stack_data(
             unstacked_predictions,
-            include_series_id=include_series_id,
+            include_series_id=True,
             series_id_name=self.series_id,
         )
+        # Move datetime index into separate date column to use when merging later
+        stacked_predictions = stacked_predictions.reset_index(drop=False)
+
+        sp_dtypes = {
+            self.time_index: X[self.time_index].dtype,
+            self.series_id: X[self.series_id].dtype,
+            self.input_target_name: y.dtype,
+        }
+        stacked_predictions = stacked_predictions.astype(sp_dtypes)
+
+        # Order prediction based on input (date, series_id)
+        output_cols = (
+            [self.series_id, self.input_target_name]
+            if include_series_id
+            else [self.input_target_name]
+        )
+        stacked_predictions = pd.merge(
+            X,
+            stacked_predictions,
+            on=[self.time_index, self.series_id],
+        )[output_cols]
 
         # Index will start at the unstacked index, so we need to reset it to the original index
         stacked_predictions.index = X.index
-        stacked_predictions = infer_feature_types(stacked_predictions)
-        return stacked_predictions
+
+        if not include_series_id:
+            return infer_feature_types(stacked_predictions[self.input_target_name])
+        else:
+            return infer_feature_types(stacked_predictions)
 
     def get_forecast_period(self, X):
         """Generates all possible forecasting time points based on latest data point in X.
diff --git a/evalml/pipelines/utils.py b/evalml/pipelines/utils.py
index 26ffb946..9b73e40c 100644
--- a/evalml/pipelines/utils.py
+++ b/evalml/pipelines/utils.py
@@ -132,6 +132,7 @@ def _get_datetime(X, y, problem_type, estimator_class, sampler_name=None):
     if add_datetime_featurizer and estimator_class.model_family not in [
         ModelFamily.ARIMA,
         ModelFamily.PROPHET,
+        ModelFamily.VARMAX,
     ]:
         components.append(DateTimeFeaturizer)
     return components
@@ -298,13 +299,7 @@ def _get_preprocessing_components(
     Returns:
         list[Transformer]: A list of applicable preprocessing components to use with the estimator.
     """
-    if is_multiseries(problem_type):
-        if include_decomposer:
-            components_functions = [_get_decomposer]
-        else:
-            return []
-
-    elif is_time_series(problem_type):
+    if is_time_series(problem_type):
         components_functions = [
             _get_label_encoder,
             _get_drop_all_null,
@@ -1508,7 +1503,7 @@ def stack_X(X, series_id_name, time_index, starting_index=None, series_id_values
         time_index (str): The name of the time index column.
         starting_index (int): The starting index to use for the stacked DataFrame. If None, the starting index
             will match that of the input data. Defaults to None.
-        series_id_values (set, list): The unique values of a series ID, used to generate the index. If None, values will
+        series_id_values (list): The unique values of a series ID, used to generate the index. If None, values will
             be generated from X column values. Required if X only has time index values and no exogenous values.
             Defaults to None.
 
@@ -1516,14 +1511,21 @@ def stack_X(X, series_id_name, time_index, starting_index=None, series_id_values
         pd.DataFrame: The restacked features.
     """
     original_columns = set()
-    series_ids = series_id_values or set()
-    if series_id_values is None:
+    if series_id_values is not None:
+        series_ids = series_id_values
+    else:
+        # Using list to maintain order (vs. a set)
+        series_ids = list()
         for col in X.columns:
             if col == time_index:
                 continue
             separated_name = col.split(MULTISERIES_SEPARATOR_SYMBOL)
             original_columns.add(MULTISERIES_SEPARATOR_SYMBOL.join(separated_name[:-1]))
-            series_ids.add(separated_name[-1])
+            series_ids.append(separated_name[-1])
+        # Remove duplicates while maintaining insertion order
+        # Need order to match series ID labels correctly when restacking columns
+        seen = set()
+        series_ids = [val for val in series_ids if not (val in seen or seen.add(val))]
 
     if len(series_ids) == 0:
         raise ValueError(
@@ -1542,7 +1544,7 @@ def stack_X(X, series_id_name, time_index, starting_index=None, series_id_values
         restacked_X = pd.DataFrame(
             {
                 time_index: time_index_col,
-                series_id_name: sorted(list(series_ids)) * len(X),
+                series_id_name: list(series_ids) * len(X),
             },
             index=stacked_index,
         )
diff --git a/evalml/preprocessing/utils.py b/evalml/preprocessing/utils.py
index dc17e75e..013dc71f 100644
--- a/evalml/preprocessing/utils.py
+++ b/evalml/preprocessing/utils.py
@@ -77,9 +77,9 @@ def split_multiseries_data(X, y, series_id, time_index, **kwargs):
         X_unstacked, y_unstacked, problem_type="time series regression", **kwargs
     )
 
-    # Get unique series value from X if there is only the time_index column
+    # Get unique series values (as a list to maintain order) from X if there is only the time_index column
     # Otherwise, this information is generated in `stack_X` from the column values
-    series_id_values = set(X[series_id]) if len(X_unstacked.columns) == 1 else None
+    series_id_values = X[series_id].unique() if len(X_unstacked.columns) == 1 else None
 
     X_train = stack_X(
         X_train_unstacked,


Test patch:
diff --git a/evalml/tests/automl_tests/test_default_algorithm.py b/evalml/tests/automl_tests/test_default_algorithm.py
index 31b8a166..826b5253 100644
--- a/evalml/tests/automl_tests/test_default_algorithm.py
+++ b/evalml/tests/automl_tests/test_default_algorithm.py
@@ -670,7 +670,7 @@ def test_default_algorithm_multiseries_time_series(
     )
 
     first_batch = algo.next_batch()
-    assert len(first_batch) == 2
+    assert len(first_batch) == 8
     pipeline = first_batch[0]
     assert pipeline.model_family == ModelFamily.VARMAX
     assert pipeline.parameters["pipeline"] == search_parameters["pipeline"]
@@ -679,8 +679,8 @@ def test_default_algorithm_multiseries_time_series(
 
     long_explore = algo.next_batch()
     long_estimators = set([pipeline.estimator.name for pipeline in long_explore])
-    assert len(long_explore) == 100
-    assert len(long_estimators) == 1
+    assert len(long_explore) == 300
+    assert len(long_estimators) == 3
 
 
 @pytest.mark.parametrize(
diff --git a/evalml/tests/automl_tests/test_iterative_algorithm.py b/evalml/tests/automl_tests/test_iterative_algorithm.py
index 3030c099..4fc4f0f5 100644
--- a/evalml/tests/automl_tests/test_iterative_algorithm.py
+++ b/evalml/tests/automl_tests/test_iterative_algorithm.py
@@ -18,7 +18,6 @@ from evalml.pipelines.components import (
     DateTimeFeaturizer,
     EmailFeaturizer,
     NaturalLanguageFeaturizer,
-    STLDecomposer,
     TimeSeriesFeaturizer,
     URLFeaturizer,
 )
@@ -98,7 +97,6 @@ def test_iterative_algorithm_init(
     assert algo.batch_number == 0
     assert algo.default_max_batches == 1
     estimators = get_estimators(problem_type)
-    decomposer = [STLDecomposer] if is_regression(problem_type) else []
     assert len(algo.allowed_pipelines) == len(
         [
             make_pipeline(
@@ -106,11 +104,15 @@ def test_iterative_algorithm_init(
                 y,
                 estimator,
                 problem_type,
+                include_decomposer=include_decomposer,
                 parameters=search_parameters,
             )
             for estimator in estimators
-        ]
-        + decomposer,
+            # Generate both decomposer and non-decomposer pipelines when problem type is multiseries time series reg.
+            for include_decomposer in (
+                [True, False] if is_regression(problem_type) else [False]
+            )
+        ],
     )
 
 
diff --git a/evalml/tests/component_tests/test_decision_tree_regressor.py b/evalml/tests/component_tests/test_decision_tree_regressor.py
index 882c8813..ba14cb87 100644
--- a/evalml/tests/component_tests/test_decision_tree_regressor.py
+++ b/evalml/tests/component_tests/test_decision_tree_regressor.py
@@ -14,6 +14,7 @@ def test_problem_types():
     assert set(DecisionTreeRegressor.supported_problem_types) == {
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     }
 
 
diff --git a/evalml/tests/component_tests/test_drop_nan_rows_transformer.py b/evalml/tests/component_tests/test_drop_nan_rows_transformer.py
index 20214a3c..ad3d6a77 100644
--- a/evalml/tests/component_tests/test_drop_nan_rows_transformer.py
+++ b/evalml/tests/component_tests/test_drop_nan_rows_transformer.py
@@ -27,8 +27,9 @@ def test_drop_rows_transformer():
     assert_frame_equal(fit_transformed_X, X_expected)
 
 
+@pytest.mark.parametrize("y_is_df", [True, False])
 @pytest.mark.parametrize("null_value", [pd.NA, np.NaN])
-def test_drop_rows_transformer_retain_ww_schema(null_value):
+def test_drop_rows_transformer_retain_ww_schema(null_value, y_is_df):
     # Expecting float because of np.NaN values
     X = pd.DataFrame(
         {"a column": [null_value, 2, 3, 4], "another col": ["a", null_value, "c", "d"]},
@@ -46,20 +47,47 @@ def test_drop_rows_transformer_retain_ww_schema(null_value):
     )
     X_expected_schema = X.ww.schema
 
-    y = pd.Series([3, 2, 1, null_value])
-    y = init_series(y, logical_type="IntegerNullable", semantic_tags="y_custom_tag")
+    if y_is_df:
+        y = pd.DataFrame(
+            {"series_a": [3, 2, 1, null_value], "series_b": [1, null_value, 3, 4]},
+        )
+        y.ww.init()
+        y.ww.set_types(
+            logical_types={
+                "series_a": "IntegerNullable",
+                "series_b": "IntegerNullable",
+            },
+            semantic_tags={"series_a": "custom_tag_a", "series_b": "custom_tag_b"},
+        )
 
-    y_expected = pd.Series([1], index=[2])
-    y_expected = init_series(
-        y_expected,
-        logical_type="IntegerNullable",
-        semantic_tags="y_custom_tag",
-    )
+        y_expected = pd.DataFrame({"series_a": [1], "series_b": [3]}, index=[2])
+        y_expected.ww.init()
+        y_expected.ww.set_types(
+            logical_types={
+                "series_a": "IntegerNullable",
+                "series_b": "IntegerNullable",
+            },
+            semantic_tags={"series_a": "custom_tag_a", "series_b": "custom_tag_b"},
+        )
+    else:
+        y = pd.Series([3, 2, 1, null_value])
+        y = init_series(y, logical_type="IntegerNullable", semantic_tags="y_custom_tag")
+
+        y_expected = pd.Series([1], index=[2])
+        y_expected = init_series(
+            y_expected,
+            logical_type="IntegerNullable",
+            semantic_tags="y_custom_tag",
+        )
     y_expected_schema = y.ww.schema
 
     drop_rows_transformer = DropNaNRowsTransformer()
     transformed_X, transformed_y = drop_rows_transformer.fit_transform(X, y)
     assert_frame_equal(transformed_X, X_expected)
-    assert_series_equal(transformed_y, y_expected)
     assert _schema_is_equal(transformed_X.ww.schema, X_expected_schema)
+
+    if y_is_df:
+        assert_frame_equal(transformed_y, y_expected)
+    else:
+        assert_series_equal(transformed_y, y_expected)
     assert transformed_y.ww.schema == y_expected_schema
diff --git a/evalml/tests/component_tests/test_en_regressor.py b/evalml/tests/component_tests/test_en_regressor.py
index 064ce573..e4eb4188 100644
--- a/evalml/tests/component_tests/test_en_regressor.py
+++ b/evalml/tests/component_tests/test_en_regressor.py
@@ -24,6 +24,7 @@ def test_problem_types():
     assert set(ElasticNetRegressor.supported_problem_types) == {
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     }
 
 
diff --git a/evalml/tests/component_tests/test_et_regressor.py b/evalml/tests/component_tests/test_et_regressor.py
index c570dd7c..aa6a4af4 100644
--- a/evalml/tests/component_tests/test_et_regressor.py
+++ b/evalml/tests/component_tests/test_et_regressor.py
@@ -14,6 +14,7 @@ def test_problem_types():
     assert set(ExtraTreesRegressor.supported_problem_types) == {
         ProblemTypes.REGRESSION,
         ProblemTypes.TIME_SERIES_REGRESSION,
+        ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
     }
 
 
diff --git a/evalml/tests/component_tests/test_time_series_featurizer.py b/evalml/tests/component_tests/test_time_series_featurizer.py
index 703132b4..8f58906c 100644
--- a/evalml/tests/component_tests/test_time_series_featurizer.py
+++ b/evalml/tests/component_tests/test_time_series_featurizer.py
@@ -990,7 +990,7 @@ def test_featurizer_y_dataframe(multiseries_ts_data_unstacked):
     featurizer = TimeSeriesFeaturizer(time_index="date", gap=1, forecast_horizon=5)
     featurizer.fit(X, y)
 
-    assert featurizer.statistically_significant_lags == [6]
+    assert featurizer.statistically_significant_lags == {col: [6] for col in y.columns}
 
     expected_y_cols = [
         f"target{MULTISERIES_SEPARATOR_SYMBOL}{i}_delay_6" for i in range(y.shape[1])
diff --git a/evalml/tests/conftest.py b/evalml/tests/conftest.py
index 80036d07..553ee9fe 100644
--- a/evalml/tests/conftest.py
+++ b/evalml/tests/conftest.py
@@ -987,6 +987,7 @@ def X_y_based_on_pipeline_or_problem_type(X_y_binary, X_y_multi, X_y_regression)
             ProblemTypes.TIME_SERIES_BINARY: "binary",
             ProblemTypes.TIME_SERIES_MULTICLASS: "multiclass",
             ProblemTypes.TIME_SERIES_REGRESSION: "regression",
+            ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION: "regression",
         }
         pipeline_classes = {
             BinaryClassificationPipeline: "binary",
diff --git a/evalml/tests/pipeline_tests/regression_pipeline_tests/test_multiseries_regression_pipeline.py b/evalml/tests/pipeline_tests/regression_pipeline_tests/test_multiseries_regression_pipeline.py
index e777c539..fb5be37b 100644
--- a/evalml/tests/pipeline_tests/regression_pipeline_tests/test_multiseries_regression_pipeline.py
+++ b/evalml/tests/pipeline_tests/regression_pipeline_tests/test_multiseries_regression_pipeline.py
@@ -120,7 +120,7 @@ def test_multiseries_pipeline_predict_in_sample(
         range(55, 65),
         index=range(90, 100),
         name="target",
-        dtype="float64",
+        dtype="int64",
     )
     if include_series_id:
         expected = pd.concat([X_holdout["series_id"], expected], axis=1)
@@ -147,7 +147,6 @@ def test_multiseries_pipeline_predict_in_sample_series_out_of_order(
 
     # Reorder rows but keep ordered by date
     # Store ordered series ID values to compare to output later
-    X_holdout_series_id = X_holdout["series_id"]
     X_index = X_holdout.index
     X_holdout = X_holdout.sample(frac=1).sort_values(by="date")
     y_holdout = y_holdout.reindex(X_holdout.index)
@@ -165,17 +164,39 @@ def test_multiseries_pipeline_predict_in_sample_series_out_of_order(
         y_train=y_train,
         include_series_id=include_series_id,
     )
+
     expected = pd.Series(
         range(55, 65),
         index=range(90, 100),
         name="target",
-        dtype="float64",
+        dtype="int64",
+    )
+    expected = pd.concat(
+        [
+            X_holdout["date"],
+            pd.Series(
+                [0, 1, 2, 3, 4] * 2,
+                name="series_id",
+                dtype=int,
+                index=range(90, 100),
+            ),
+            expected,
+        ],
+        axis=1,
     )
+    expected = pd.merge(
+        infer_feature_types(X_holdout[["date", "series_id"]]),
+        expected,
+        on=["date", "series_id"],
+    )
+    expected = expected.drop("date", axis=1)
+    expected.index = range(90, 100)
+
     if include_series_id:
-        expected = pd.concat([X_holdout_series_id, expected], axis=1)
         expected = infer_feature_types(expected)
         pd.testing.assert_frame_equal(y_pred, expected)
     else:
+        expected = expected["target"]
         pd.testing.assert_series_equal(y_pred, expected)
 
 
@@ -209,7 +230,7 @@ def test_multiseries_pipeline_predict(
             range(55, 65),
             index=range(90, 100),
             name="target",
-            dtype="float64",
+            dtype="int64",
         )
     # Only the first predicted value is present in the delayed features
     else:
@@ -217,7 +238,7 @@ def test_multiseries_pipeline_predict(
             [85, 86, 87, 88, 89, 0, 0, 0, 0, 0],
             index=range(90, 100),
             name="target",
-            dtype="float64",
+            dtype="int64",
         )
     pd.testing.assert_series_equal(y_pred, expected)
 
diff --git a/evalml/tests/pipeline_tests/test_pipeline_utils.py b/evalml/tests/pipeline_tests/test_pipeline_utils.py
index 5a9d4b16..44641ef4 100644
--- a/evalml/tests/pipeline_tests/test_pipeline_utils.py
+++ b/evalml/tests/pipeline_tests/test_pipeline_utils.py
@@ -146,7 +146,7 @@ def test_make_pipeline(
             datetime = (
                 [DateTimeFeaturizer]
                 if estimator_class.model_family
-                not in [ModelFamily.ARIMA, ModelFamily.PROPHET]
+                not in [ModelFamily.ARIMA, ModelFamily.PROPHET, ModelFamily.VARMAX]
                 and "dates" in column_names
                 else []
             )
@@ -170,25 +170,22 @@ def test_make_pipeline(
             )
 
             if is_time_series(problem_type):
-                if is_multiseries(problem_type):
-                    expected_components = dfs + decomposer + [estimator_class]
-                else:
-                    expected_components = (
-                        dfs
-                        + label_encoder
-                        + email_featurizer
-                        + url_featurizer
-                        + drop_null
-                        + natural_language_featurizer
-                        + imputer
-                        + delayed_features
-                        + decomposer
-                        + datetime
-                        + ohe
-                        + drop_nan_rows_transformer
-                        + standard_scaler
-                        + [estimator_class]
-                    )
+                expected_components = (
+                    dfs
+                    + label_encoder
+                    + email_featurizer
+                    + url_featurizer
+                    + drop_null
+                    + natural_language_featurizer
+                    + imputer
+                    + delayed_features
+                    + decomposer
+                    + datetime
+                    + ohe
+                    + drop_nan_rows_transformer
+                    + standard_scaler
+                    + [estimator_class]
+                )
             else:
                 expected_components = (
                     dfs
@@ -624,7 +621,7 @@ def test_get_estimators():
                 problem_type=ProblemTypes.MULTISERIES_TIME_SERIES_REGRESSION,
             ),
         )
-        == 1
+        == 5
     )
 
     assert len(get_estimators(problem_type=ProblemTypes.BINARY, model_families=[])) == 0
@@ -1474,13 +1471,11 @@ def test_stack_data_noop():
     pd.testing.assert_series_equal(stack_data(series_y), series_y)
 
 
-@pytest.mark.parametrize("series_id_values_type", [set, list])
 @pytest.mark.parametrize("no_features", [True, False])
 @pytest.mark.parametrize("starting_index", [None, 1, 132])
 def test_stack_X(
     starting_index,
     no_features,
-    series_id_values_type,
     multiseries_ts_data_stacked,
     multiseries_ts_data_unstacked,
 ):
@@ -1491,7 +1486,7 @@ def test_stack_X(
         X_expected.index = X_expected.index + starting_index
 
     if no_features:
-        series_id_values = series_id_values_type(str(i) for i in range(0, 5))
+        series_id_values = list(str(i) for i in range(0, 5))
         X = pd.DataFrame(X["date"])
         X_expected = X_expected[["date", "series_id"]]