diff --git a/changelog.d/207.fixed.md b/changelog.d/207.fixed.md new file mode 100644 index 00000000..7dd1efda --- /dev/null +++ b/changelog.d/207.fixed.md @@ -0,0 +1,2 @@ +Per-variable QRF models now derive distinct seeds, so variables imputed together no longer share one random quantile per row and come out comonotonic. `QRF` also accepts a `seed` argument. +Derived seeds stay within the supported uint32 range and are used consistently during numeric and classification tuning and target-specific subsampling. diff --git a/microimpute/models/qrf.py b/microimpute/models/qrf.py index 8edc5d08..9e5560d3 100644 --- a/microimpute/models/qrf.py +++ b/microimpute/models/qrf.py @@ -10,7 +10,7 @@ from quantile_forest import RandomForestQuantileRegressor from sklearn.ensemble import RandomForestClassifier -from microimpute.config import VALIDATE_CONFIG +from microimpute.config import RANDOM_STATE, VALIDATE_CONFIG from microimpute.models.imputer import Imputer, ImputerResults try: @@ -607,6 +607,7 @@ def __init__( batch_size: Optional[int] = None, cleanup_interval: int = 10, max_train_samples: Optional[int] = None, + seed: Optional[int] = RANDOM_STATE, ) -> None: """Initialize the QRF model. @@ -618,8 +619,19 @@ def __init__( max_train_samples: If set, subsample X_train to at most this many rows before fitting. Reduces memory and training time while preserving sequential covariance structure. + seed: Base random seed. Each imputed variable is given a distinct + seed derived from it, so variables imputed together draw + independently. Pass None for non-reproducible draws. """ - super().__init__(log_level=log_level) + if seed is not None and ( + isinstance(seed, bool) + or not isinstance(seed, (int, np.integer)) + or not 0 <= seed < 2**32 + ): + raise ValueError( + f"seed must be an integer from 0 to 2**32 - 1, or None, got {seed!r}" + ) + super().__init__(log_level=log_level, seed=seed) self.models = {} self.log_level = log_level self.memory_efficient = memory_efficient @@ -695,20 +707,45 @@ def _encode_imputed_variable( return data + def _seed_for_variable(self, variable: str) -> Optional[int]: + """Derive a distinct seed for one imputed variable. + + Each per-variable model builds its own generator from the seed it is + given. Handing every variable the same seed makes them draw the same + random quantiles in the same row order, so variables imputed together + come out comonotonic regardless of their dependence in the donor. The + offset is shared with target-specific subsampling and wraps within + sklearn's uint32 seed range. + """ + if self.seed is None: + return None + try: + variable_offset = (self.imputed_variables or []).index(variable) + except ValueError as error: + # Falling back to offset 0 would hand this variable the same draws + # as the first target, which is the comonotonicity this method + # exists to prevent - and it would do so silently. + raise ValueError( + f"Cannot derive a seed for {variable!r}: it is not among the " + f"imputed variables {list(self.imputed_variables or [])}." + ) from error + return (int(self.seed) + variable_offset) % 2**32 + def _create_model_for_variable(self, variable: str, **kwargs) -> Any: """Create the appropriate model (classifier or regressor) based on variable type.""" categorical_targets = getattr(self, "categorical_targets", {}) boolean_targets = getattr(self, "boolean_targets", {}) + seed = self._seed_for_variable(variable) if variable in categorical_targets: # Use classifier for categorical targets - return _RandomForestClassifierModel(seed=self.seed, logger=self.logger) + return _RandomForestClassifierModel(seed=seed, logger=self.logger) elif variable in boolean_targets: # Use classifier for boolean targets - return _RandomForestClassifierModel(seed=self.seed, logger=self.logger) + return _RandomForestClassifierModel(seed=seed, logger=self.logger) else: # Use QRF for numeric targets - return _QRFModel(seed=self.seed, logger=self.logger) + return _QRFModel(seed=seed, logger=self.logger) def _fit_model( self, @@ -800,12 +837,7 @@ def _target_fit_data( self.max_train_samples is not None and len(target_train) > self.max_train_samples ): - try: - variable_offset = (self.imputed_variables or []).index(variable) - except ValueError: - variable_offset = 0 - seed = None if self.seed is None else self.seed + variable_offset - rng = np.random.default_rng(seed) + rng = np.random.default_rng(self._seed_for_variable(variable)) sel = rng.choice( len(target_train), size=self.max_train_samples, replace=False ) @@ -1465,7 +1497,9 @@ def objective(trial: optuna.Trial) -> float: y_val = X_val_fold[var] # Create and fit QRF model with trial parameters - model = _QRFModel(seed=self.seed, logger=self.logger) + model = _QRFModel( + seed=self._seed_for_variable(var), logger=self.logger + ) model.fit( X_train_augmented[encoded_predictors], X_train_fold[var], @@ -1623,7 +1657,7 @@ def objective(trial: optuna.Trial) -> float: # Create and fit RFC model with trial parameters model = _RandomForestClassifierModel( - seed=self.seed, logger=self.logger + seed=self._seed_for_variable(var), logger=self.logger ) # Determine variable type and fit appropriately diff --git a/tests/test_models/test_qrf.py b/tests/test_models/test_qrf.py index e3099f75..50d3f31b 100644 --- a/tests/test_models/test_qrf.py +++ b/tests/test_models/test_qrf.py @@ -1649,3 +1649,166 @@ def test_qrf_fit_predict_with_max_train_samples() -> None: assert result.shape == (n_test, 1) assert not result.isna().any().any() + + +def test_per_variable_models_draw_independently() -> None: + """Variables imputed together must not share one random quantile per row. + + Every per-variable model builds its own generator from the seed it is + given. When they all receive the same seed they draw the same quantiles in + the same row order, so the imputed variables come out comonotonic whatever + their dependence in the donor. + """ + rng = np.random.default_rng(7) + n = 1500 + predictors = pd.DataFrame( + {"inc": rng.normal(30, 8, n), "age": rng.normal(45, 12, n)} + ) + targets = ["savings", "property_wealth", "corporate_wealth"] + train = predictors.copy() + for variable in targets: + # Shared signal, independent shocks: the conditional dependence is nil. + train[variable] = 0.5 * predictors["inc"] + rng.normal(0, 10, n) + + test = pd.DataFrame({"inc": rng.normal(30, 8, 600), "age": rng.normal(45, 12, 600)}) + + model = QRF(log_level="WARNING") + imputations = model.fit(train, ["inc", "age"], targets).predict(test) + + ranks = imputations[targets].rank().corr() + off_diagonal = [ + abs(ranks.loc[a, b]) for i, a in enumerate(targets) for b in targets[i + 1 :] + ] + assert max(off_diagonal) < 0.4, ( + "imputed variables are comonotonic; per-variable models are sharing " + f"a seed (rank correlations {off_diagonal})" + ) + + +def test_seed_is_configurable_and_reproducible() -> None: + """QRF should accept a seed, and the same seed should reproduce draws.""" + rng = np.random.default_rng(3) + n = 400 + train = pd.DataFrame({"x": rng.normal(size=n)}) + train["y"] = train["x"] + rng.normal(0, 1, n) + test = pd.DataFrame({"x": rng.normal(size=120)}) + + first = QRF(log_level="WARNING", seed=1234).fit(train, ["x"], ["y"]).predict(test) + second = QRF(log_level="WARNING", seed=1234).fit(train, ["x"], ["y"]).predict(test) + + np.testing.assert_allclose(first["y"], second["y"]) + + +@pytest.mark.parametrize("target_type", ["numeric", "boolean", "categorical"]) +def test_qrf_max_seed_supports_multiple_targets(target_type: str) -> None: + """Every valid sklearn seed must support reproducible multi-target fits.""" + rng = np.random.default_rng(71) + data = pd.DataFrame( + {name: rng.normal(size=120) for name in ["x", "first", "second"]} + ) + if target_type == "boolean": + data[["first", "second"]] = data[["first", "second"]] > 0 + elif target_type == "categorical": + for name in ["first", "second"]: + data[name] = np.where(data[name] > 0, "yes", "no") + + predictions = [] + for _ in range(2): + fitted = QRF(seed=2**32 - 1).fit( + data, ["x"], ["first", "second"], n_estimators=12 + ) + seeds = [model.seed for model in fitted.models.values()] + assert len(set(seeds)) == 2 + assert all(0 <= seed < 2**32 for seed in seeds) + predictions.append(fitted.predict(data[["x"]].iloc[:20])) + pd.testing.assert_frame_equal(*predictions) + assert not predictions[0].isna().any().any() + + +@pytest.mark.parametrize("seed", [None, 0, 42]) +def test_qrf_child_seeds_preserve_existing_seed_values(seed) -> None: + """Ordinary seeds and entropy-based draws retain their established meaning.""" + model = QRF(seed=seed) + model.imputed_variables = ["first", "second"] + expected = [None, None] if seed is None else [seed, seed + 1] + assert [ + model._seed_for_variable(name) for name in model.imputed_variables + ] == expected + + +@pytest.mark.parametrize("seed", [-1, 2**32, 1.5]) +def test_qrf_invalid_base_seed_is_not_normalized(seed) -> None: + """An invalid seed fails at construction, not part-way through a fit.""" + with pytest.raises(ValueError, match="seed must be an integer"): + QRF(seed=seed) + + +def test_qrf_seed_for_unknown_variable_raises() -> None: + """Falling back to offset 0 would silently recreate comonotonicity.""" + model = QRF(seed=100) + model.imputed_variables = ["a", "b"] + assert [model._seed_for_variable(v) for v in ["a", "b"]] == [100, 101] + with pytest.raises(ValueError, match="not among the imputed variables"): + model._seed_for_variable("z") + + +@pytest.mark.parametrize("target_type", ["numeric", "boolean"]) +def test_qrf_tuning_uses_distinct_target_seeds(monkeypatch, target_type: str) -> None: + """Real Optuna folds must fit the same independent streams as final models.""" + from microimpute.models.qrf import _RandomForestClassifierModel + + rng = np.random.default_rng(73) + data = pd.DataFrame( + {name: rng.normal(size=100) for name in ["x", "first", "second"]} + ) + model = QRF(seed=123) + model.imputed_variables = ["first", "second"] + if target_type == "numeric": + internal_model = _QRFModel + tune = model._tune_qrf_hyperparameters + else: + data[["first", "second"]] = data[["first", "second"]] > 0 + model.boolean_targets = {"first": {}, "second": {}} + internal_model = _RandomForestClassifierModel + tune = model._tune_rfc_hyperparameters + + fitted_seeds = [] + original_fit = internal_model.fit + + def record_fit(self, X, y, **kwargs): + fitted_seeds.append((y.name, self.seed)) + return original_fit(self, X, y, **kwargs) + + monkeypatch.setattr(internal_model, "fit", record_fit) + tune(data, ["x"], ["first", "second"], n_cv_folds=2, n_trials=1) + assert fitted_seeds == [("first", 123), ("second", 124)] * 2 + + +def test_qrf_target_subsampling_uses_bounded_child_seed(monkeypatch) -> None: + """Filtered training rows and their model use one consistent child seed.""" + rng = np.random.default_rng(9) + data = pd.DataFrame( + {name: rng.normal(size=120) for name in ["x", "first", "second"]} + ) + fitted_indices = {} + original_fit = _QRFModel.fit + + def record_fit(self, X, y, **kwargs): + fitted_indices[y.name] = X.index.to_numpy() + return original_fit(self, X, y, **kwargs) + + monkeypatch.setattr(_QRFModel, "fit", record_fit) + QRF(seed=2**32 - 1, max_train_samples=50).fit( + data, + ["x"], + ["first", "second"], + target_filters={ + name: np.ones(len(data), dtype=bool) for name in ["first", "second"] + }, + n_estimators=12, + ) + # The second stream wraps to zero at the uint32 boundary. + expected_indices = np.random.default_rng(0).choice( + len(data), size=50, replace=False + ) + np.testing.assert_array_equal(fitted_indices["second"], expected_indices)