int_constraints

Author: mathias-von-ottenbreit

API_REFERENCE.md

@@ -59,7 +59,7 @@ Specifies which metric to use for validating the model and tuning ***m***. Avail
 Specifies the quantile to use when ***family*** is "quantile".
 
 
-## Method: fit(X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0))
+## Method: fit(X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0), interaction_constraints:List[int]=[])
 
 ***This method fits the model to data.***
 
@@ -84,11 +84,14 @@ An optional list of integers specifying the indexes of observations to be used f
 An optional list of integers specifying the indexes of predictors (columns) in ***X*** that should be prioritized. Terms of the prioritized predictors will enter the model as long as they reduce the training error and do not contain too few effective observations. They will also be updated more often.
 
 #### monotonic_constraints
-An optional list of integers specifying monotonic constraints on model terms. For example, if there are three predictors in ***X***, then monotonic_constraints = [1,0,-1] means that 1) the first predictor in ***X*** cannot be used in interaction terms and all terms using the first predictor in ***X*** as a main effect must have positive regression coefficients, 2) there are no monotonic constraints on terms using the second predictor in ***X***, and 3) the third predictor in ***X*** cannot be used in interaction terms and all terms using the third predictor in ***X*** as a main effect must have negative regression coefficients.
+An optional list of integers specifying monotonic constraints on model terms. For example, if there are three predictors in ***X***, then monotonic_constraints = [1,0,-1] means that 1) the first predictor in ***X*** cannot be used in interaction terms as a secondary effect and all terms using the first predictor in ***X*** as a main effect must have positive regression coefficients, 2) there are no monotonic constraints on terms using the second predictor in ***X***, and 3) the third predictor in ***X*** cannot be used in interaction terms as a secondary effect and all terms using the third predictor in ***X*** as a main effect must have negative regression coefficients.
 
 #### group
 A numpy vector of integers that is used when ***family*** is "group_gaussian". For example, ***group*** may represent year (could be useful in a time series model).
 
+#### interaction_constraints
+An optional list of integers specifying interaction constraints on model terms. For example, if there are three predictors in ***X***, then interaction_constraints = [1,0,2] means that 1) the first predictor in ***X*** cannot be used in interaction terms as a secondary effect, 2) there are no interaction constraints on terms using the second predictor in ***X***, and 3) the third predictor in ***X*** cannot be used in any interaction terms.
+
 
 ## Method: predict(X:npt.ArrayLike, cap_predictions_to_minmax_in_training:bool=True)
 

aplr/aplr.py

@@ -50,9 +50,9 @@ def __set_params_cpp(self):
         self.APLRRegressor.validation_tuning_metric=self.validation_tuning_metric
         self.APLRRegressor.quantile=self.quantile
 
-    def fit(self, X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[],group:npt.ArrayLike = np.empty(0)):
+    def fit(self, X:npt.ArrayLike, y:npt.ArrayLike, sample_weight:npt.ArrayLike = np.empty(0), X_names:List[str]=[], validation_set_indexes:List[int]=[], prioritized_predictors_indexes:List[int]=[], monotonic_constraints:List[int]=[], group:npt.ArrayLike = np.empty(0), interaction_constraints:List[int]=[]):
         self.__set_params_cpp()
-        self.APLRRegressor.fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group)
+        self.APLRRegressor.fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group,interaction_constraints)
 
     def predict(self, X:npt.ArrayLike, cap_predictions_to_minmax_in_training:bool=True)->npt.ArrayLike:
         return self.APLRRegressor.predict(X, cap_predictions_to_minmax_in_training)

cpp/APLRRegressor.h

@@ -50,17 +50,20 @@ class APLRRegressor
     VectorXi group_validation;
     std::set<int> unique_groups_train;
     std::set<int> unique_groups_validation;
+    std::vector<int> interaction_constraints;
 
     //Methods
     void validate_input_to_fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight,const std::vector<std::string> &X_names, 
         const std::vector<size_t> &validation_set_indexes, const std::vector<size_t> &prioritized_predictors_indexes,
-        const std::vector<int> &monotonic_constraints, const VectorXi &group);
+        const std::vector<int> &monotonic_constraints, const VectorXi &group, const std::vector<int> &interaction_constraints);
     void throw_error_if_validation_set_indexes_has_invalid_indexes(const VectorXd &y, const std::vector<size_t> &validation_set_indexes);
     void throw_error_if_prioritized_predictors_indexes_has_invalid_indexes(const MatrixXd &X, const std::vector<size_t> &prioritized_predictors_indexes);
     void throw_error_if_monotonic_constraints_has_invalid_indexes(const MatrixXd &X, const std::vector<int> &monotonic_constraints);
+    void throw_error_if_interaction_constraints_has_invalid_indexes(const MatrixXd &X, const std::vector<int> &interaction_constraints);
     void define_training_and_validation_sets(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight, 
         const std::vector<size_t> &validation_set_indexes, const VectorXi &group);
-    void initialize(const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints);
+    void initialize(const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints,
+                    const std::vector<int> &interaction_constraints);
     bool check_if_base_term_has_only_one_unique_value(size_t base_term);
     void add_term_to_terms_eligible_current(Term &term);
     VectorXd calculate_neg_gradient_current(const VectorXd &sample_weight_train);
@@ -153,7 +156,7 @@ class APLRRegressor
     ~APLRRegressor();
     void fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight=VectorXd(0),const std::vector<std::string> &X_names={},const std::vector<size_t> &validation_set_indexes={},
         const std::vector<size_t> &prioritized_predictors_indexes={}, const std::vector<int> &monotonic_constraints={},
-        const VectorXi &group=VectorXi(0));
+        const VectorXi &group=VectorXi(0), const std::vector<int> &interaction_constraints={});
     VectorXd predict(const MatrixXd &X, bool cap_predictions_to_minmax_in_training=true);
     void set_term_names(const std::vector<std::string> &X_names);
     MatrixXd calculate_local_feature_importance(const MatrixXd &X);
@@ -214,15 +217,16 @@ APLRRegressor::~APLRRegressor()
 //invalidating validation_ratio. The rest of indexes are used to train. 
 void APLRRegressor::fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight,const std::vector<std::string> &X_names,
     const std::vector<size_t> &validation_set_indexes,const std::vector<size_t> &prioritized_predictors_indexes, 
-    const std::vector<int> &monotonic_constraints, const VectorXi &group)
+    const std::vector<int> &monotonic_constraints, const VectorXi &group, const std::vector<int> &interaction_constraints)
 {
     throw_error_if_family_does_not_exist();
     throw_error_if_link_function_does_not_exist();
     throw_error_if_tweedie_power_is_invalid();
-    validate_input_to_fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group);
+    validate_input_to_fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group,
+                        interaction_constraints);
     define_training_and_validation_sets(X,y,sample_weight,validation_set_indexes,group);
     scale_training_observations_if_using_log_link_function();
-    initialize(prioritized_predictors_indexes, monotonic_constraints);
+    initialize(prioritized_predictors_indexes, monotonic_constraints, interaction_constraints);
     execute_boosting_steps();
     update_coefficients_for_all_steps();
     print_final_summary();
@@ -281,7 +285,8 @@ void APLRRegressor::throw_error_if_tweedie_power_is_invalid()
 
 void APLRRegressor::validate_input_to_fit(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight,
     const std::vector<std::string> &X_names, const std::vector<size_t> &validation_set_indexes, 
-    const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints, const VectorXi &group)
+    const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints, const VectorXi &group,
+    const std::vector<int> &interaction_constraints)
 {
     if(X.rows()!=y.size()) throw std::runtime_error("X and y must have the same number of rows.");
     if(X.rows()<2) throw std::runtime_error("X and y cannot have less than two rows.");
@@ -292,6 +297,7 @@ void APLRRegressor::validate_input_to_fit(const MatrixXd &X,const VectorXd &y,co
     throw_error_if_validation_set_indexes_has_invalid_indexes(y, validation_set_indexes);
     throw_error_if_prioritized_predictors_indexes_has_invalid_indexes(X, prioritized_predictors_indexes);
     throw_error_if_monotonic_constraints_has_invalid_indexes(X, monotonic_constraints);
+    throw_error_if_interaction_constraints_has_invalid_indexes(X, interaction_constraints);
     throw_error_if_response_contains_invalid_values(y);
     throw_error_if_sample_weight_contains_invalid_values(y, sample_weight);
     bool group_is_of_incorrect_size{family=="group_gaussian" && group.rows()!=y.rows()};
@@ -329,6 +335,24 @@ void APLRRegressor::throw_error_if_monotonic_constraints_has_invalid_indexes(con
         throw std::runtime_error("monotonic_constraints must either be empty or a vector with one integer for each column in X.");
 }
 
+void APLRRegressor::throw_error_if_interaction_constraints_has_invalid_indexes(const MatrixXd &X, const std::vector<int> &interaction_constraints)
+{
+    bool dimension_error{interaction_constraints.size()>0 && interaction_constraints.size() != X.cols()};
+    if(dimension_error)
+        throw std::runtime_error("interaction_constraints must either be empty or a vector with one integer for each column in X.");
+    bool value_error{false};
+    for(auto &constraint:interaction_constraints)    
+    {
+        if(constraint<0 || constraint>2)
+        {
+            value_error=true;
+            break;
+        }
+    }
+    if(value_error)
+        throw std::runtime_error("The value of each element in interaction_constraints must be 0 (no constraint), 1 (the predictor can only appear as a main effect in an interaction), or 2 (the predictor cannot appear in any interaction).");
+}
+
 void APLRRegressor::throw_error_if_response_contains_invalid_values(const VectorXd &y)
 {
     if(link_function=="logit" || family=="binomial")
@@ -507,7 +531,8 @@ void APLRRegressor::scale_training_observations_if_using_log_link_function()
     }
 }
 
-void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints)
+void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints,
+                                const std::vector<int> &interaction_constraints)
 {
     number_of_base_terms=static_cast<size_t>(X_train.cols());
 
@@ -547,6 +572,16 @@ void APLRRegressor::initialize(const std::vector<size_t> &prioritized_predictors
         }
     }
 
+    this->interaction_constraints=interaction_constraints;
+    bool interaction_constraints_provided{interaction_constraints.size()>0};
+    if(interaction_constraints_provided)
+    {
+        for(auto &term_eligible_current:terms_eligible_current)
+        {
+            term_eligible_current.set_interaction_constraint(interaction_constraints[term_eligible_current.base_term]);
+        }
+    }
+
     linear_predictor_current=VectorXd::Constant(y_train.size(),intercept);
     linear_predictor_null_model=linear_predictor_current;
     linear_predictor_current_validation=VectorXd::Constant(y_validation.size(),intercept);
@@ -811,6 +846,7 @@ void APLRRegressor::determine_interactions_to_consider(const std::vector<size_t>
     interactions_to_consider=std::vector<Term>();
     interactions_to_consider.reserve(static_cast<size_t>(X_train.cols())*terms.size());
     bool monotonic_constraints_provided{monotonic_constraints.size()>0};
+    bool interaction_constraints_provided{interaction_constraints.size()>0};
 
     VectorXi indexes_for_terms_to_consider_as_interaction_partners{find_indexes_for_terms_to_consider_as_interaction_partners()};
     for (auto &model_term_index:indexes_for_terms_to_consider_as_interaction_partners)
@@ -823,11 +859,19 @@ void APLRRegressor::determine_interactions_to_consider(const std::vector<size_t>
                 Term interaction{Term(new_term_index)};
                 if(monotonic_constraints_provided)
                     interaction.set_monotonic_constraint(monotonic_constraints[new_term_index]);
+                if(interaction_constraints_provided)
+                {
+                    interaction.set_interaction_constraint(interaction_constraints[new_term_index]);
+                    bool new_term_cannot_be_an_interaction{interaction.get_interaction_constraint()==2};
+                    if(new_term_cannot_be_an_interaction) 
+                        continue;
+                }                    
                 Term model_term_without_given_terms{terms[model_term_index]};
                 model_term_without_given_terms.given_terms.clear();
                 model_term_without_given_terms.cleanup_when_this_term_was_added_as_a_given_term();
                 Term model_term_with_added_given_term{terms[model_term_index]};
-                bool model_term_without_given_terms_can_be_a_given_term{model_term_without_given_terms.get_monotonic_constraint()==0};
+                bool model_term_without_given_terms_can_be_a_given_term{model_term_without_given_terms.get_monotonic_constraint()==0
+                    && model_term_without_given_terms.get_interaction_constraint()==0};
                 if(model_term_without_given_terms_can_be_a_given_term)
                     model_term_with_added_given_term.given_terms.push_back(model_term_without_given_terms);
                 add_necessary_given_terms_to_interaction(interaction, model_term_with_added_given_term);
@@ -1345,6 +1389,7 @@ void APLRRegressor::cleanup_after_fit()
     group_validation.resize(0);
     unique_groups_train.clear();
     unique_groups_validation.clear();
+    interaction_constraints.clear();
 }
 
 VectorXd APLRRegressor::predict(const MatrixXd &X, bool cap_predictions_to_minmax_in_training)

cpp/pythonbinding.cpp

@@ -24,7 +24,7 @@ PYBIND11_MODULE(aplr_cpp, m) {
             )
         .def("fit", &APLRRegressor::fit,py::arg("X"),py::arg("y"),py::arg("sample_weight")=VectorXd(0),py::arg("X_names")=std::vector<std::string>(),
             py::arg("validation_set_indexes")=std::vector<size_t>(),py::arg("prioritized_predictors_indexes")=std::vector<size_t>(),
-            py::arg("monotonic_constraints")=std::vector<int>(),py::arg("group")=VectorXi(0),
+            py::arg("monotonic_constraints")=std::vector<int>(),py::arg("group")=VectorXi(0), py::arg("interaction_constraints")=std::vector<int>(),
             py::call_guard<py::scoped_ostream_redirect,py::scoped_estream_redirect>())
         .def("predict", &APLRRegressor::predict,py::arg("X"),py::arg("bool cap_predictions_to_minmax_in_training")=true)
         .def("set_term_names", &APLRRegressor::set_term_names,py::arg("X_names"))

cpp/term.h

@@ -38,6 +38,7 @@ class Term
     VectorXd negative_gradient_discretized;
     std::vector<size_t> observations_in_bins;
     int monotonic_constraint;
+    int interaction_constraint;
 
     void calculate_error_where_given_terms_are_zero(const VectorXd &negative_gradient, const VectorXd &sample_weight);
     void initialize_parameters_in_estimate_split_point(size_t bins,double v,size_t min_observations_in_split);
@@ -88,6 +89,8 @@ class Term
     bool get_can_be_used_as_a_given_term();
     void set_monotonic_constraint(int constraint);
     int get_monotonic_constraint();
+    void set_interaction_constraint(int constraint);
+    int get_interaction_constraint();
 
     friend bool operator== (const Term &p1, const Term &p2);
     friend class APLRRegressor;
@@ -96,14 +99,15 @@ class Term
 Term::Term(size_t base_term,const std::vector<Term> &given_terms,double split_point,bool direction_right,double coefficient):
 name{""},base_term{base_term},given_terms{given_terms},split_point{split_point},direction_right{direction_right},coefficient{coefficient},
 split_point_search_errors_sum{std::numeric_limits<double>::infinity()},ineligible_boosting_steps{0},can_be_used_as_a_given_term{false},
-monotonic_constraint{0}
+monotonic_constraint{0},interaction_constraint{0}
 {
 }
 
 Term::Term(const Term &other):
 name{other.name},base_term{other.base_term},given_terms{other.given_terms},split_point{other.split_point},direction_right{other.direction_right},
 coefficient{other.coefficient},coefficient_steps{other.coefficient_steps},split_point_search_errors_sum{other.split_point_search_errors_sum},
-ineligible_boosting_steps{0},can_be_used_as_a_given_term{other.can_be_used_as_a_given_term},monotonic_constraint{other.monotonic_constraint}
+ineligible_boosting_steps{0},can_be_used_as_a_given_term{other.can_be_used_as_a_given_term},monotonic_constraint{other.monotonic_constraint},
+interaction_constraint{other.interaction_constraint}
 {
 }
 
@@ -679,6 +683,16 @@ int Term::get_monotonic_constraint()
     return monotonic_constraint;
 }
 
+void Term::set_interaction_constraint(int constraint)
+{
+    interaction_constraint = constraint;
+}
+
+int Term::get_interaction_constraint()
+{
+    return interaction_constraint;
+}
+
 
 std::vector<std::vector<size_t>> distribute_terms_indexes_to_cores(std::vector<size_t> &term_indexes,size_t n_jobs)
 {