diff --git a/feature_engine/_base_transformers/base_numerical.py b/feature_engine/_base_transformers/base_numerical.py
index 10b24e99a..fed663213 100644
--- a/feature_engine/_base_transformers/base_numerical.py
+++ b/feature_engine/_base_transformers/base_numerical.py
@@ -28,18 +28,18 @@ class BaseNumericalTransformer(
     variable transformers, discretisers, math combination.
     """
 
-    def fit(self, X: pd.DataFrame) -> pd.DataFrame:
+    def _fit_setup(self, X: pd.DataFrame):
         """
         Checks that input is a dataframe, finds numerical variables, or alternatively
-        checks that variables entered by the user are of type numerical.
+        checks that variables entered by the user are of type numerical, and checks
+        for NA and Inf. Does not assign any trailing-underscore attribute, so that
+        subclasses can defer attribute assignment until the rest of their fit logic
+        has completed successfully.
 
         Parameters
         ----------
         X : Pandas DataFrame
 
-        y : Pandas Series, np.array. Default = None
-            Parameter is necessary for compatibility with sklearn Pipeline.
-
         Raises
         ------
         TypeError
@@ -53,6 +53,9 @@ def fit(self, X: pd.DataFrame) -> pd.DataFrame:
         -------
         X : Pandas DataFrame
             The same dataframe entered as parameter
+
+        variables_ : List
+            The variables that were found or checked.
         """
 
         # check input dataframe
@@ -60,23 +63,24 @@ def fit(self, X: pd.DataFrame) -> pd.DataFrame:
 
         # find or check for numerical variables
         if self.variables is None:
-            self.variables_ = find_numerical_variables(
-                X, return_empty=self.return_empty
-            )
+            variables_ = find_numerical_variables(X, return_empty=self.return_empty)
         else:
-            self.variables_ = check_numerical_variables(X, self.variables)
+            variables_ = check_numerical_variables(X, self.variables)
 
         # check if dataset contains na or inf
-        _check_contains_na(X, self.variables_)
-        _check_contains_inf(X, self.variables_)
+        _check_contains_na(X, variables_)
+        _check_contains_inf(X, variables_)
 
-        # save input features
-        self.feature_names_in_ = X.columns.tolist()
+        return X, variables_
+
+    def _get_feature_names_in(self, X):
+        """Get the names and number of features in the train set (the dataframe
+        used during fit)."""
 
-        # save train set shape
+        self.feature_names_in_ = X.columns.tolist()
         self.n_features_in_ = X.shape[1]
 
-        return X
+        return self
 
     def _check_transform_input_and_state(self, X: pd.DataFrame) -> pd.DataFrame:
         """
diff --git a/feature_engine/_base_transformers/mixins.py b/feature_engine/_base_transformers/mixins.py
index 4d4b7d254..9207873be 100644
--- a/feature_engine/_base_transformers/mixins.py
+++ b/feature_engine/_base_transformers/mixins.py
@@ -1,4 +1,4 @@
-from typing import Dict, List, Union
+from typing import Dict, List, Tuple, Union
 
 import pandas as pd
 from numpy import ndarray
@@ -46,10 +46,14 @@ def transform_x_y(self, X: pd.DataFrame, y: pd.Series):
 
 
 class FitFromDictMixin:
-    def _fit_from_dict(self, X: pd.DataFrame, user_dict_: Dict) -> pd.DataFrame:
+    def _fit_from_dict(
+        self, X: pd.DataFrame, user_dict_: Dict
+    ) -> Tuple[pd.DataFrame, List[Union[str, int]]]:
         """
         Checks that input is a dataframe, checks that variables in the dictionary
-        entered by the user are of type numerical.
+        entered by the user are of type numerical. Does not assign any
+        trailing-underscore attribute, so that subclasses can defer attribute
+        assignment until the rest of their fit logic has completed successfully.
 
         Parameters
         ----------
@@ -71,25 +75,22 @@ def _fit_from_dict(self, X: pd.DataFrame, user_dict_: Dict) -> pd.DataFrame:
         -------
         X : Pandas DataFrame
             The same dataframe entered as parameter
+
+        variables_ : List
+            The variables in the dictionary.
         """
         # check input dataframe
         X = check_X(X)
 
         # find or check for numerical variables
         variables = list(user_dict_.keys())
-        self.variables_ = check_numerical_variables(X, variables)
+        variables_ = check_numerical_variables(X, variables)
 
         # check if dataset contains na or inf
-        _check_contains_na(X, self.variables_)
-        _check_contains_inf(X, self.variables_)
-
-        # save input features
-        self.feature_names_in_ = X.columns.tolist()
-
-        # save train set shape
-        self.n_features_in_ = X.shape[1]
+        _check_contains_na(X, variables_)
+        _check_contains_inf(X, variables_)
 
-        return X
+        return X, variables_
 
 
 class GetFeatureNamesOutMixin:
diff --git a/feature_engine/creation/cyclical_features.py b/feature_engine/creation/cyclical_features.py
index ea201e6ac..24018b0cd 100644
--- a/feature_engine/creation/cyclical_features.py
+++ b/feature_engine/creation/cyclical_features.py
@@ -155,11 +155,15 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
             It is not needed in this transformer. You can pass y or None.
         """
         if self.max_values is None:
-            X = super().fit(X)
-            self.max_values_ = X[self.variables_].max().to_dict()
+            X, variables_ = self._fit_setup(X)
+            max_values_ = X[variables_].max().to_dict()
         else:
-            super()._fit_from_dict(X, self.max_values)
-            self.max_values_ = self.max_values
+            X, variables_ = super()._fit_from_dict(X, self.max_values)
+            max_values_ = self.max_values
+
+        self.variables_ = variables_
+        self.max_values_ = max_values_
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/feature_engine/creation/geo_features.py b/feature_engine/creation/geo_features.py
index 78b977c47..bb2698d07 100644
--- a/feature_engine/creation/geo_features.py
+++ b/feature_engine/creation/geo_features.py
@@ -234,8 +234,8 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         # check input dataframe
         X = check_X(X)
 
-        # Store coordinate variables
-        self.variables_: List[Union[str, int]] = [
+        # Coordinate variables
+        variables: List[Union[str, int]] = [
             self.lat1,
             self.lon1,
             self.lat2,
@@ -243,17 +243,17 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         ]
 
         # Check all coordinate columns exist
-        missing = set(self.variables_) - set(X.columns)
+        missing = set(variables) - set(X.columns)
         if missing:
             raise ValueError(
                 f"Coordinate columns {missing} are not present in the dataframe."
             )
 
         # Check coordinate columns are numerical
-        check_numerical_variables(X, self.variables_)
+        check_numerical_variables(X, variables)
 
         # Check for missing values
-        _check_contains_na(X, self.variables_)
+        _check_contains_na(X, variables)
 
         # Validate coordinate ranges if enabled
         if self.validate_ranges:
@@ -269,6 +269,9 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
                         f"Longitude values in '{lon_col}' must be between -180 and 180."
                     )
 
+        # save coordinate variables
+        self.variables_ = variables
+
         # save input features
         self.feature_names_in_ = X.columns.tolist()
 
diff --git a/feature_engine/discretisation/arbitrary.py b/feature_engine/discretisation/arbitrary.py
index 5ceae738d..5776cd71d 100644
--- a/feature_engine/discretisation/arbitrary.py
+++ b/feature_engine/discretisation/arbitrary.py
@@ -152,10 +152,12 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
             y is not needed in this transformer. You can pass y or None.
         """
         # check input dataframe
-        X = super()._fit_from_dict(X, self.binning_dict)
+        X, variables_ = super()._fit_from_dict(X, self.binning_dict)
 
+        self.variables_ = variables_
         # for consistency with the rest of the discretisers, we add this attribute
         self.binner_dict_ = self.binning_dict
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/feature_engine/discretisation/base_discretiser.py b/feature_engine/discretisation/base_discretiser.py
index 76302ea07..6c61d05d3 100644
--- a/feature_engine/discretisation/base_discretiser.py
+++ b/feature_engine/discretisation/base_discretiser.py
@@ -10,7 +10,8 @@ class BaseDiscretiser(BaseNumericalTransformer):
     """
     Shared set-up checks and methods across numerical discretisers.
 
-    Important: inherits fit() functionality and tags from BaseNumericalTransformer.
+    Important: inherits _fit_setup(), _get_feature_names_in() and tags from
+    BaseNumericalTransformer. Subclasses implement fit() themselves.
     """
 
     def __init__(
diff --git a/feature_engine/discretisation/decision_tree.py b/feature_engine/discretisation/decision_tree.py
index 53378e8ca..8af4b9b60 100644
--- a/feature_engine/discretisation/decision_tree.py
+++ b/feature_engine/discretisation/decision_tree.py
@@ -225,7 +225,7 @@ def __init__(
         self.random_state = random_state
         self.return_empty = return_empty
 
-    def fit(self, X: pd.DataFrame, y: pd.Series):  # type: ignore
+    def fit(self, X: pd.DataFrame, y: pd.Series):
         """
         Fit one decision tree per variable to discretise with cross-validation and
         grid-search for hyperparameters.
@@ -252,7 +252,7 @@ def fit(self, X: pd.DataFrame, y: pd.Series):  # type: ignore
             check_classification_targets(y)
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
         if self.param_grid:
             param_grid = self.param_grid
@@ -262,7 +262,7 @@ def fit(self, X: pd.DataFrame, y: pd.Series):  # type: ignore
         binner_dict_ = {}
         scores_dict_ = {}
 
-        for var in self.variables_:
+        for var in variables_:
 
             if self.regression:
                 model = DecisionTreeRegressor(random_state=self.random_state)
@@ -280,7 +280,7 @@ def fit(self, X: pd.DataFrame, y: pd.Series):  # type: ignore
             scores_dict_[var] = tree_model.score(X[var].to_frame(), y)
 
         if self.bin_output != "prediction":
-            for var in self.variables_:
+            for var in variables_:
                 clf = binner_dict_[var].best_estimator_
                 threshold = clf.tree_.threshold
                 feature = clf.tree_.feature
@@ -291,6 +291,9 @@ def fit(self, X: pd.DataFrame, y: pd.Series):  # type: ignore
 
         self.binner_dict_ = binner_dict_
         self.scores_dict_ = scores_dict_
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
+
         return self
 
     def transform(self, X: pd.DataFrame) -> pd.DataFrame:
diff --git a/feature_engine/discretisation/equal_frequency.py b/feature_engine/discretisation/equal_frequency.py
index 2471bfe01..a2137870f 100644
--- a/feature_engine/discretisation/equal_frequency.py
+++ b/feature_engine/discretisation/equal_frequency.py
@@ -170,17 +170,21 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
-        self.binner_dict_ = {}
+        binner_dict_ = {}
 
-        for var in self.variables_:
+        for var in variables_:
             tmp, bins = pd.qcut(x=X[var], q=self.q, retbins=True, duplicates="drop")
 
             # Prepend/Append infinities to accommodate outliers
             bins = list(bins)
             bins[0] = float("-inf")
             bins[len(bins) - 1] = float("inf")
-            self.binner_dict_[var] = bins
+            binner_dict_[var] = bins
+
+        self.binner_dict_ = binner_dict_
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/discretisation/equal_width.py b/feature_engine/discretisation/equal_width.py
index 366378fa3..bab5c5396 100644
--- a/feature_engine/discretisation/equal_width.py
+++ b/feature_engine/discretisation/equal_width.py
@@ -179,12 +179,12 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
         # fit
-        self.binner_dict_ = {}
+        binner_dict_ = {}
 
-        for var in self.variables_:
+        for var in variables_:
             tmp, bins = pd.cut(
                 x=X[var],
                 bins=self.bins,
@@ -197,6 +197,10 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
             bins = list(bins)
             bins[0] = float("-inf")
             bins[len(bins) - 1] = float("inf")
-            self.binner_dict_[var] = bins
+            binner_dict_[var] = bins
+
+        self.binner_dict_ = binner_dict_
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/discretisation/geometric_width.py b/feature_engine/discretisation/geometric_width.py
index 0bd913e85..6da74339c 100644
--- a/feature_engine/discretisation/geometric_width.py
+++ b/feature_engine/discretisation/geometric_width.py
@@ -174,12 +174,12 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
         # fit
-        self.binner_dict_ = {}
+        binner_dict_ = {}
 
-        for var in self.variables_:
+        for var in variables_:
             min_, max_ = X[var].min(), X[var].max()
             increment = np.power(max_ - min_, 1.0 / self.bins)
             bins = np.r_[
@@ -187,6 +187,10 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
             ]
             bins = np.sort(bins)
             bins = list(bins)
-            self.binner_dict_[var] = bins
+            binner_dict_[var] = bins
+
+        self.binner_dict_ = binner_dict_
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/arbitrary_imputer.py b/feature_engine/imputation/arbitrary_imputer.py
index 058ca9d1d..333a81918 100644
--- a/feature_engine/imputation/arbitrary_imputer.py
+++ b/feature_engine/imputation/arbitrary_imputer.py
@@ -163,17 +163,19 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         # find or check for numerical variables
         # create the imputer dictionary
         if self.imputer_dict:
-            self.variables_ = check_numerical_variables(
+            variables_ = check_numerical_variables(
                 X, list(self.imputer_dict.keys())
             )
-            self.imputer_dict_ = self.imputer_dict
+            imputer_dict_ = self.imputer_dict
         else:
             if self.variables is None:
-                self.variables_ = find_numerical_variables(X, self.return_empty)
+                variables_ = find_numerical_variables(X, self.return_empty)
             else:
-                self.variables_ = check_numerical_variables(X, self.variables)
-            self.imputer_dict_ = {var: self.arbitrary_number for var in self.variables_}
+                variables_ = check_numerical_variables(X, self.variables)
+            imputer_dict_ = {var: self.arbitrary_number for var in variables_}
 
+        self.variables_ = variables_
+        self.imputer_dict_ = imputer_dict_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/categorical.py b/feature_engine/imputation/categorical.py
index 1c989c437..2cd4a00a3 100644
--- a/feature_engine/imputation/categorical.py
+++ b/feature_engine/imputation/categorical.py
@@ -181,22 +181,22 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         # select variables to encode
         if self.ignore_format is True:
             if self.variables is None:
-                self.variables_ = find_all_variables(X, self.return_empty)
+                variables_ = find_all_variables(X, self.return_empty)
             else:
-                self.variables_ = check_all_variables(X, self.variables)
+                variables_ = check_all_variables(X, self.variables)
         else:
             if self.variables is None:
-                self.variables_ = find_categorical_variables(X, self.return_empty)
+                variables_ = find_categorical_variables(X, self.return_empty)
             else:
-                self.variables_ = check_categorical_variables(X, self.variables)
+                variables_ = check_categorical_variables(X, self.variables)
 
         if self.imputation_method == "missing":
-            self.imputer_dict_ = {var: self.fill_value for var in self.variables_}
+            imputer_dict_ = {var: self.fill_value for var in variables_}
 
         elif self.imputation_method == "frequent":
             # if imputing only 1 variable:
-            if len(self.variables_) == 1:
-                var = self.variables_[0]
+            if len(variables_) == 1:
+                var = variables_[0]
                 mode_vals = X[var].mode()
 
                 # Some variables may contain more than 1 mode:
@@ -205,13 +205,13 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
                         f"The variable {var} contains multiple frequent categories."
                     )
 
-                self.imputer_dict_ = {var: mode_vals[0]}
+                imputer_dict_ = {var: mode_vals[0]}
 
             # imputing multiple variables:
             else:
                 # Returns a dataframe with 1 row if there is one mode per
                 # variable, or more rows if there are more modes:
-                mode_vals = X[self.variables_].mode()
+                mode_vals = X[variables_].mode()
 
                 # Careful: some variables contain multiple modes
                 if len(mode_vals) > 1:
@@ -225,8 +225,10 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
                         f"categories."
                     )
 
-                self.imputer_dict_ = mode_vals.iloc[0].to_dict()
+                imputer_dict_ = mode_vals.iloc[0].to_dict()
 
+        self.variables_ = variables_
+        self.imputer_dict_ = imputer_dict_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/drop_missing_data.py b/feature_engine/imputation/drop_missing_data.py
index e92b22f23..5be7a19cb 100644
--- a/feature_engine/imputation/drop_missing_data.py
+++ b/feature_engine/imputation/drop_missing_data.py
@@ -163,16 +163,15 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # find variables for which indicator should be added
         if self.variables is None:
-            self.variables_ = find_all_variables(X, self.return_empty)
+            variables_ = find_all_variables(X, self.return_empty)
         else:
-            self.variables_ = check_all_variables(X, self.variables)
+            variables_ = check_all_variables(X, self.variables)
 
         # If user passes a threshold, then missing_only is ignored:
         if self.threshold is None and self.missing_only is True:
-            self.variables_ = [
-                var for var in self.variables_ if X[var].isnull().sum() > 0
-            ]
+            variables_ = [var for var in variables_ if X[var].isnull().sum() > 0]
 
+        self.variables_ = variables_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/end_tail.py b/feature_engine/imputation/end_tail.py
index cd895e971..e52500056 100644
--- a/feature_engine/imputation/end_tail.py
+++ b/feature_engine/imputation/end_tail.py
@@ -187,35 +187,37 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # find or check for numerical variables
         if self.variables is None:
-            self.variables_ = find_numerical_variables(X, self.return_empty)
+            variables_ = find_numerical_variables(X, self.return_empty)
         else:
-            self.variables_ = check_numerical_variables(X, self.variables)
+            variables_ = check_numerical_variables(X, self.variables)
 
         # estimate imputation values
         if self.imputation_method == "max":
-            self.imputer_dict_ = (X[self.variables_].max() * self.fold).to_dict()
+            imputer_dict_ = (X[variables_].max() * self.fold).to_dict()
 
         elif self.imputation_method == "gaussian":
             if self.tail == "right":
-                self.imputer_dict_ = (
-                    X[self.variables_].mean() + self.fold * X[self.variables_].std()
+                imputer_dict_ = (
+                    X[variables_].mean() + self.fold * X[variables_].std()
                 ).to_dict()
             elif self.tail == "left":
-                self.imputer_dict_ = (
-                    X[self.variables_].mean() - self.fold * X[self.variables_].std()
+                imputer_dict_ = (
+                    X[variables_].mean() - self.fold * X[variables_].std()
                 ).to_dict()
 
         elif self.imputation_method == "iqr":
-            IQR = X[self.variables_].quantile(0.75) - X[self.variables_].quantile(0.25)
+            IQR = X[variables_].quantile(0.75) - X[variables_].quantile(0.25)
             if self.tail == "right":
-                self.imputer_dict_ = (
-                    X[self.variables_].quantile(0.75) + (IQR * self.fold)
+                imputer_dict_ = (
+                    X[variables_].quantile(0.75) + (IQR * self.fold)
                 ).to_dict()
             elif self.tail == "left":
-                self.imputer_dict_ = (
-                    X[self.variables_].quantile(0.25) - (IQR * self.fold)
+                imputer_dict_ = (
+                    X[variables_].quantile(0.25) - (IQR * self.fold)
                 ).to_dict()
 
+        self.variables_ = variables_
+        self.imputer_dict_ = imputer_dict_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/mean_median.py b/feature_engine/imputation/mean_median.py
index 7812c1217..049b768e3 100644
--- a/feature_engine/imputation/mean_median.py
+++ b/feature_engine/imputation/mean_median.py
@@ -138,17 +138,19 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # find or check for numerical variables
         if self.variables is None:
-            self.variables_ = find_numerical_variables(X, self.return_empty)
+            variables_ = find_numerical_variables(X, self.return_empty)
         else:
-            self.variables_ = check_numerical_variables(X, self.variables)
+            variables_ = check_numerical_variables(X, self.variables)
 
         # find imputation parameters: mean or median
         if self.imputation_method == "mean":
-            self.imputer_dict_ = X[self.variables_].mean().to_dict()
+            imputer_dict_ = X[variables_].mean().to_dict()
 
         elif self.imputation_method == "median":
-            self.imputer_dict_ = X[self.variables_].median().to_dict()
+            imputer_dict_ = X[variables_].median().to_dict()
 
+        self.variables_ = variables_
+        self.imputer_dict_ = imputer_dict_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/missing_indicator.py b/feature_engine/imputation/missing_indicator.py
index dba8d3ef6..012cf7b23 100644
--- a/feature_engine/imputation/missing_indicator.py
+++ b/feature_engine/imputation/missing_indicator.py
@@ -141,15 +141,14 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # find variables for which indicator should be added
         if self.variables is None:
-            self.variables_ = find_all_variables(X, self.return_empty)
+            variables_ = find_all_variables(X, self.return_empty)
         else:
-            self.variables_ = check_all_variables(X, self.variables)
+            variables_ = check_all_variables(X, self.variables)
 
         if self.missing_only is True:
-            self.variables_ = [
-                var for var in self.variables_ if X[var].isnull().sum() > 0
-            ]
+            variables_ = [var for var in variables_ if X[var].isnull().sum() > 0]
 
+        self.variables_ = variables_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/imputation/random_sample.py b/feature_engine/imputation/random_sample.py
index 9366c1120..bc11e0dac 100644
--- a/feature_engine/imputation/random_sample.py
+++ b/feature_engine/imputation/random_sample.py
@@ -198,26 +198,29 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # find variables to impute
         if self.variables is None:
-            self.variables_ = find_all_variables(X, self.return_empty)
+            variables_ = find_all_variables(X, self.return_empty)
         else:
-            self.variables_ = check_all_variables(X, self.variables)
+            variables_ = check_all_variables(X, self.variables)
 
         # take a copy of the selected variables
-        self.X_ = X[self.variables_].copy()
+        X_ = X[variables_].copy()
 
         # check the variables assigned to the random state
         if self.seed == "observation":
-            self.random_state = _check_variables_input_value(self.random_state)
-            if isinstance(self.random_state, (int, str)):
-                self.random_state = [self.random_state]
-            if self.random_state and any(
-                var for var in self.random_state if var not in X.columns
+            random_state = _check_variables_input_value(self.random_state)
+            if isinstance(random_state, (int, str)):
+                random_state = [random_state]
+            if random_state and any(
+                var for var in random_state if var not in X.columns
             ):
                 raise ValueError(
                     "There are variables assigned as random state which are not part "
                     "of the training dataframe."
                 )
+            self.random_state = random_state
 
+        self.variables_ = variables_
+        self.X_ = X_
         self._get_feature_names_in(X)
 
         return self
diff --git a/feature_engine/scaling/mean_normalization.py b/feature_engine/scaling/mean_normalization.py
index c13949c3e..620865735 100644
--- a/feature_engine/scaling/mean_normalization.py
+++ b/feature_engine/scaling/mean_normalization.py
@@ -130,18 +130,24 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
-        self.mean_ = X[self.variables_].mean().to_dict()
-        self.range_ = (X[self.variables_].max() - X[self.variables_].min()).to_dict()
+        X, variables_ = self._fit_setup(X)
+
+        mean_ = X[variables_].mean().to_dict()
+        range_ = (X[variables_].max() - X[variables_].min()).to_dict()
 
         # check for constant columns
-        constant_columns = [col for col, value in self.range_.items() if value == 0]
+        constant_columns = [col for col, value in range_.items() if value == 0]
         if constant_columns:
             raise ValueError(
                 f"The following variable(s) are constant: {constant_columns}. "
                 "Division by zero is not allowed. Please remove constant columns."
             )
 
+        self.variables_ = variables_
+        self.mean_ = mean_
+        self.range_ = range_
+        self._get_feature_names_in(X)
+
         return self
 
     def transform(self, X: pd.DataFrame) -> pd.DataFrame:
diff --git a/feature_engine/transformation/arcsin.py b/feature_engine/transformation/arcsin.py
index 04fd8ac86..da4045fa4 100644
--- a/feature_engine/transformation/arcsin.py
+++ b/feature_engine/transformation/arcsin.py
@@ -133,15 +133,18 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
         # check if the variables are in the correct range
-        if ((X[self.variables_] < 0) | (X[self.variables_] > 1)).any().any():
+        if ((X[variables_] < 0) | (X[variables_] > 1)).any().any():
             raise ValueError(
                 "Some variables contain values outside the possible range 0-1. "
                 "Can't apply the arcsin transformation. "
             )
 
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
+
         return self
 
     def transform(self, X: pd.DataFrame) -> pd.DataFrame:
diff --git a/feature_engine/transformation/arcsinh.py b/feature_engine/transformation/arcsinh.py
index dfc9ae2cb..92ebf2c0a 100644
--- a/feature_engine/transformation/arcsinh.py
+++ b/feature_engine/transformation/arcsinh.py
@@ -172,7 +172,10 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe and find/check numerical variables
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
+
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/feature_engine/transformation/boxcox.py b/feature_engine/transformation/boxcox.py
index dc3b9e50f..52d5feffb 100644
--- a/feature_engine/transformation/boxcox.py
+++ b/feature_engine/transformation/boxcox.py
@@ -147,12 +147,16 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
-        self.lambda_dict_ = {}
+        lambda_dict_ = {}
 
-        for var in self.variables_:
-            _, self.lambda_dict_[var] = stats.boxcox(X[var])
+        for var in variables_:
+            _, lambda_dict_[var] = stats.boxcox(X[var])
+
+        self.variables_ = variables_
+        self.lambda_dict_ = lambda_dict_
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/feature_engine/transformation/log.py b/feature_engine/transformation/log.py
index 4e4024ac8..9cd12c307 100644
--- a/feature_engine/transformation/log.py
+++ b/feature_engine/transformation/log.py
@@ -172,31 +172,33 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
 
         # check input dataframe
         if isinstance(self.C, dict):
-            X = super()._fit_from_dict(X, self.C)
+            X, variables_ = super()._fit_from_dict(X, self.C)
         else:
-            X = super().fit(X)
+            X, variables_ = self._fit_setup(X)
 
-        self.C_ = self.C
+        C_ = self.C
 
         # calculate C to add to each variable
         if self.C == "auto":
             # we add 0 to positive variables
-            c_dict = {var: 0 for var in self.variables_ if X[var].min() > 0}
+            c_dict = {var: 0 for var in variables_ if X[var].min() > 0}
 
             # we add the minimum plus 1 to non-positive variables
-            non_positive_vars = [
-                var for var in self.variables_ if var not in c_dict.keys()
-            ]
+            non_positive_vars = [var for var in variables_ if var not in c_dict.keys()]
             c_dict.update(dict(X[non_positive_vars].min(axis=0).abs() + 1))
-            self.C_ = c_dict  # type:ignore
+            C_ = c_dict  # type:ignore
 
         # C=0 is the original LogTransformer contract: no constant is added,
         # so fail fast at fit time exactly as before this class supported C.
-        if self.C_ == 0 and (X[self.variables_] <= 0).any().any():
+        if C_ == 0 and (X[variables_] <= 0).any().any():
             raise ValueError(
                 "Some variables contain zero or negative values, can't apply log"
             )
 
+        self.variables_ = variables_
+        self.C_ = C_
+        self._get_feature_names_in(X)
+
         return self
 
     def transform(self, X: pd.DataFrame) -> pd.DataFrame:
diff --git a/feature_engine/transformation/power.py b/feature_engine/transformation/power.py
index c6257cf29..89aea9bf2 100644
--- a/feature_engine/transformation/power.py
+++ b/feature_engine/transformation/power.py
@@ -132,7 +132,10 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        super().fit(X)
+        X, variables_ = self._fit_setup(X)
+
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/feature_engine/transformation/reciprocal.py b/feature_engine/transformation/reciprocal.py
index 936fef68e..22678544c 100644
--- a/feature_engine/transformation/reciprocal.py
+++ b/feature_engine/transformation/reciprocal.py
@@ -124,15 +124,18 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
         # check if the variables contain the value 0
-        if (X[self.variables_] == 0).any().any():
+        if (X[variables_] == 0).any().any():
             raise ValueError(
                 "Some variables contain the value zero, can't apply reciprocal "
                 "transformation."
             )
 
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
+
         return self
 
     def transform(self, X: pd.DataFrame) -> pd.DataFrame:
diff --git a/feature_engine/transformation/yeojohnson.py b/feature_engine/transformation/yeojohnson.py
index f723ebe61..82fa53dac 100644
--- a/feature_engine/transformation/yeojohnson.py
+++ b/feature_engine/transformation/yeojohnson.py
@@ -140,12 +140,16 @@ def fit(self, X: pd.DataFrame, y: Optional[pd.Series] = None):
         """
 
         # check input dataframe
-        X = super().fit(X)
+        X, variables_ = self._fit_setup(X)
 
-        self.lambda_dict_ = {}
+        lambda_dict_ = {}
 
-        for var in self.variables_:
-            _, self.lambda_dict_[var] = stats.yeojohnson(X[var])
+        for var in variables_:
+            _, lambda_dict_[var] = stats.yeojohnson(X[var])
+
+        self.variables_ = variables_
+        self.lambda_dict_ = lambda_dict_
+        self._get_feature_names_in(X)
 
         return self
 
diff --git a/tests/estimator_checks/non_fitted_error_checks.py b/tests/estimator_checks/non_fitted_error_checks.py
index cfa2f64e2..774fce592 100644
--- a/tests/estimator_checks/non_fitted_error_checks.py
+++ b/tests/estimator_checks/non_fitted_error_checks.py
@@ -89,3 +89,40 @@ def _implements_inverse_transform(estimator):
         # implemented", which is all this check needs to establish.
         return True
     return True
+
+
+def check_raises_non_fitted_error_when_fit_fails(estimator, X, y=None):
+    """
+    Check that if fit() raises partway through (after some trailing-underscore
+    attributes may already have been computed, but before fit has completed),
+    the transformer does not end up looking fitted: transform() must still
+    raise NotFittedError.
+
+    This guards against attributes like `variables_` or `imputer_dict_` being
+    assigned to `self` as soon as they are computed, rather than only once the
+    rest of fit()'s logic has completed successfully. sklearn's own
+    check_estimator suite does not cover this: check_fit_check_is_fitted only
+    tests "never fitted" and "successfully fitted with well-behaved data",
+    never "fit raises on bad input" - that scenario is inherently specific to
+    what makes a given transformer's own fit logic fail.
+
+    Parameters
+    ----------
+    estimator: feature-engine transformer instance.
+
+    X: pandas DataFrame
+        Input designed to make estimator.fit() raise partway through.
+
+    y: pandas Series, default=None
+        Target, if the estimator's fit() requires one.
+    """
+    transformer = clone(estimator)
+
+    with pytest.raises((ValueError, TypeError, KeyError)):
+        if y is not None:
+            transformer.fit(X, y)
+        else:
+            transformer.fit(X)
+
+    with pytest.raises(NotFittedError):
+        transformer.transform(X)
diff --git a/tests/test_base_transformers/test_base_numerical_transformer.py b/tests/test_base_transformers/test_base_numerical_transformer.py
index 88006114c..4934e4b57 100644
--- a/tests/test_base_transformers/test_base_numerical_transformer.py
+++ b/tests/test_base_transformers/test_base_numerical_transformer.py
@@ -11,6 +11,12 @@ def __init__(self):
         self.variables = None
         self.return_empty = False
 
+    def fit(self, X):
+        X, variables_ = self._fit_setup(X)
+        self.variables_ = variables_
+        self._get_feature_names_in(X)
+        return X
+
     def transform(self, X):
         return self._check_transform_input_and_state(X)
 
diff --git a/tests/test_creation/test_check_estimator_creation.py b/tests/test_creation/test_check_estimator_creation.py
index e3c22caa1..23dec93c3 100644
--- a/tests/test_creation/test_check_estimator_creation.py
+++ b/tests/test_creation/test_check_estimator_creation.py
@@ -13,6 +13,9 @@
     RelativeFeatures,
 )
 from tests.estimator_checks.estimator_checks import check_feature_engine_estimator
+from tests.estimator_checks.non_fitted_error_checks import (
+    check_raises_non_fitted_error_when_fit_fails,
+)
 
 sklearn_version = parse_version(parse_version(sklearn.__version__).base_version)
 
@@ -97,3 +100,38 @@ def test_geo_distance_transformer_in_pipeline():
     Xtp = pipe.fit_transform(X.copy(), y)
 
     pd.testing.assert_frame_equal(Xtt, Xtp)
+
+
+@pytest.mark.parametrize(
+    "estimator",
+    [
+        CyclicalFeatures(),
+        MathFeatures(variables=["feature_1", "feature_2"], func=["sum", "mean"]),
+        RelativeFeatures(
+            variables=["feature_1"], reference=["feature_2"], func=["div"]
+        ),
+        DecisionTreeFeatures(regression=False),
+        GeoDistanceFeatures(lat1="lat1", lon1="lon1", lat2="lat2", lon2="lon2"),
+    ],
+)
+def test_raises_non_fitted_error_when_error_during_fit(estimator):
+    y = pd.Series([0, 1, 0, 1, 0])
+
+    if isinstance(estimator, GeoDistanceFeatures):
+        # non-numerical coordinate columns: fails after variables_ would have
+        # been set, at the "check coordinate columns are numerical" step.
+        X = pd.DataFrame(
+            {
+                "lat1": ["a", "b"],
+                "lon1": ["c", "d"],
+                "lat2": ["e", "f"],
+                "lon2": ["g", "h"],
+            }
+        )
+        y = pd.Series([0, 1])
+    else:
+        # the named/expected variables aren't in the df (or aren't numerical):
+        # fails at variable selection.
+        X = pd.DataFrame({"cat1": ["a", "b", "c", "a", "b"]})
+
+    check_raises_non_fitted_error_when_fit_fails(estimator, X, y)
diff --git a/tests/test_discretisation/test_check_estimator_discretisers.py b/tests/test_discretisation/test_check_estimator_discretisers.py
index 87e175eac..a1f78c1e0 100644
--- a/tests/test_discretisation/test_check_estimator_discretisers.py
+++ b/tests/test_discretisation/test_check_estimator_discretisers.py
@@ -14,6 +14,9 @@
     GeometricWidthDiscretiser,
 )
 from tests.estimator_checks.estimator_checks import check_feature_engine_estimator
+from tests.estimator_checks.non_fitted_error_checks import (
+    check_raises_non_fitted_error_when_fit_fails,
+)
 
 sklearn_version = parse_version(parse_version(sklearn.__version__).base_version)
 
@@ -63,3 +66,13 @@ def test_transformers_within_pipeline(transformer):
     Xtp = pipe.fit_transform(X, y)
 
     pd.testing.assert_frame_equal(Xtt, Xtp)
+
+
+@pytest.mark.parametrize("estimator", _estimators)
+def test_raises_non_fitted_error_when_error_during_fit(estimator):
+    # no numerical variables in the df: fails at variable selection, before any
+    # of binner_dict_/scores_dict_/variables_ would be computed.
+    X = pd.DataFrame({"cat1": ["a", "b", "c", "a", "b"]})
+    y = pd.Series([0, 1, 0, 1, 0])
+
+    check_raises_non_fitted_error_when_fit_fails(estimator, X, y)
diff --git a/tests/test_imputation/test_check_estimator_imputers.py b/tests/test_imputation/test_check_estimator_imputers.py
index 8b47b6147..6f9d0c4fc 100644
--- a/tests/test_imputation/test_check_estimator_imputers.py
+++ b/tests/test_imputation/test_check_estimator_imputers.py
@@ -15,6 +15,9 @@
     RandomSampleImputer,
 )
 from tests.estimator_checks.estimator_checks import check_feature_engine_estimator
+from tests.estimator_checks.non_fitted_error_checks import (
+    check_raises_non_fitted_error_when_fit_fails,
+)
 
 _estimators = [
     MeanImputer(),
@@ -74,3 +77,27 @@ def test_transformers_in_pipeline_with_set_output_pandas(transformer):
     Xtp = pipe.fit_transform(X, y)
 
     pd.testing.assert_frame_equal(Xtt, Xtp)
+
+
+@pytest.mark.parametrize("estimator", _estimators)
+def test_raises_non_fitted_error_when_error_during_fit(estimator):
+    if estimator.__class__.__name__ in ["MeanImputer", "EndTailImputer"]:
+        # no numerical variables in the df: fails at variable selection.
+        X = pd.DataFrame({"cat1": ["a", "b", "c", "a", "b"]})
+    elif estimator.__class__.__name__ == "ArbitraryImputer":
+        X = pd.DataFrame({"cat1": ["a", "b", "c", "a", "b"]})
+    elif estimator.__class__.__name__ == "CategoricalImputer":
+        # equally frequent categories: fails after variables_ would have been
+        # selected, inside the "frequent" imputation logic itself.
+        estimator = estimator.__class__(imputation_method="frequent")
+        X = pd.DataFrame({"cat1": ["a", "a", "b", "b"]})
+    elif estimator.__class__.__name__ == "RandomSampleImputer":
+        # invalid random_state: fails after variables_/X_ would have been set.
+        estimator = RandomSampleImputer(seed="observation", random_state="not_a_col")
+        X = pd.DataFrame({"num1": [1.0, 2.0, 3.0, 4.0, 5.0]})
+    else:
+        # AddMissingIndicator, DropMissingData: no reachable failure point
+        # once variables are selected, so fail at input validation instead.
+        X = pd.DataFrame()
+
+    check_raises_non_fitted_error_when_fit_fails(estimator, X)
diff --git a/tests/test_scaling/test_mean_normalization.py b/tests/test_scaling/test_mean_normalization.py
index 997b396db..807a8a9fc 100644
--- a/tests/test_scaling/test_mean_normalization.py
+++ b/tests/test_scaling/test_mean_normalization.py
@@ -6,6 +6,9 @@
 
 from feature_engine.scaling import MeanNormalisationScaler, MeanNormalizationScaler
 from tests.estimator_checks.fit_functionality_checks import check_return_empty
+from tests.estimator_checks.non_fitted_error_checks import (
+    check_raises_non_fitted_error_when_fit_fails,
+)
 
 DEPRECATION_WARNING = (
     "MeanNormalizationScaler was deprecated in favour of "
@@ -157,6 +160,21 @@ def test_constant_columns_error(transformer_class):
         transformer.fit(df)
 
 
+def test_raises_non_fitted_error_when_error_during_fit(transformer_class):
+    # constant column: fails after mean_/range_ would have been computed, at
+    # the "check for constant columns" step - real regression guard for the
+    # deferred trailing-underscore attribute assignment.
+    df = pd.DataFrame(
+        {
+            "var1": [1.0, 2.0, 3.0],
+            "var2": [4.0, 5.0, 3.0],
+            "var3": [7.0, 7.0, 7.0],
+        }
+    )
+    transformer = make_transformer(transformer_class)
+    check_raises_non_fitted_error_when_fit_fails(transformer, df)
+
+
 def test_check_return_empty(transformer_class):
     transformer = make_transformer(transformer_class)
     if transformer_class is MeanNormalizationScaler:
diff --git a/tests/test_transformation/test_check_estimator_transformers.py b/tests/test_transformation/test_check_estimator_transformers.py
index 8f482e10d..6aac49791 100644
--- a/tests/test_transformation/test_check_estimator_transformers.py
+++ b/tests/test_transformation/test_check_estimator_transformers.py
@@ -16,6 +16,9 @@
     YeoJohnsonTransformer,
 )
 from tests.estimator_checks.estimator_checks import check_feature_engine_estimator
+from tests.estimator_checks.non_fitted_error_checks import (
+    check_raises_non_fitted_error_when_fit_fails,
+)
 
 _estimators = [
     BoxCoxTransformer(),
@@ -95,3 +98,30 @@ def test_transformers_in_pipeline_with_set_output_pandas(transformer):
     Xtp = pipe.fit_transform(X, y)
 
     pd.testing.assert_frame_equal(Xtt, Xtp)
+
+
+@pytest.mark.parametrize("estimator", _estimators)
+def test_raises_non_fitted_error_when_error_during_fit(estimator):
+    name = estimator.__class__.__name__
+
+    if name == "BoxCoxTransformer":
+        # non-positive values: boxcox itself raises after variables_ selection.
+        X = pd.DataFrame({"num1": [-1.0, 2.0, 3.0, 4.0, 5.0]})
+    elif name == "LogTransformer":
+        # default C=0: zero/negative values raise after variables_ selection.
+        X = pd.DataFrame({"num1": [-1.0, 2.0, 3.0, 4.0, 5.0]})
+    elif name == "ArcsinTransformer":
+        # values outside 0-1: raises after variables_ selection.
+        X = pd.DataFrame({"num1": [1.1, 2.0, 3.0, 4.0, 5.0]})
+    elif name == "ReciprocalTransformer":
+        # zero values: raises after variables_ selection.
+        X = pd.DataFrame({"num1": [0.0, 2.0, 3.0, 4.0, 5.0]})
+    else:
+        # LogCpTransformer (C="auto" never fails on the values themselves),
+        # ArcSinhTransformer, PowerTransformer, YeoJohnsonTransformer: none of
+        # these validate values beyond being numerical, so there is no
+        # reachable failure point once variables_ would be selected. Fail at
+        # variable selection instead.
+        X = pd.DataFrame({"cat1": ["a", "b", "c", "a", "b"]})
+
+    check_raises_non_fitted_error_when_fit_fails(estimator, X)