added group_gaussian family

Author: mathias-von-ottenbreit

API_REFERENCE.md

@@ -14,7 +14,7 @@ The learning rate. Must be greater than zero and not more than one. The higher t
 Used to randomly split training observations into training and validation if ***validation_set_indexes*** is not specified when fitting.
 
 #### family (default = "gaussian")
-Determines the loss function used. Allowed values are "gaussian", "binomial", "poisson", "gamma" and "tweedie". This is used together with ***link_function***. 
+Determines the loss function used. Allowed values are "gaussian", "binomial", "poisson", "gamma", "tweedie" and "group_gaussian". This is used together with ***link_function***. When ***family*** is "group_gaussian" then the "group" argument in the ***fit*** method must be provided. In the latter case APLR will try to minimize group MSE when training the model.
 
 #### link_function (default = "identity")
 Determines how the linear predictor is transformed to predictions. Allowed values are "identity", "logit" and "log". For an ordinary regression model use ***family*** "gaussian" and ***link_function*** "identity". For logistic regression use ***family*** "binomial" and ***link_function*** "logit". For a multiplicative model use the "log" ***link_function***. The "log" ***link_function*** often works best with a "poisson", "gamma" or "tweedie" ***family***, depending on the data. The ***family*** "poisson", "gamma" or "tweedie" should only be used with the "log" ***link_function***. Inappropriate combinations of ***family*** and ***link_function*** may result in a warning message when fitting the model and/or a poor model fit. Please note that values other than "identity" typically require a significantly higher ***m*** (or ***v***) in order to converge.
@@ -55,7 +55,7 @@ Specifies the variance power for the "tweedie" ***family***.
 #### validation_tuning_metric (default = "default")
 Specifies which metric to use for validating the model and tuning ***m***. Available options are "default" (using the same methodology as when calculating the training error), "mse", "mae", "negative_gini" and "rankability". The default is often a choice that fits well with respect to the ***family*** chosen. However, if you want to use ***family*** or ***tweedie_power*** as tuning parameters then the default is not suitable. "rankability" uses a methodology similar to the one described in https://towardsdatascience.com/how-to-calculate-roc-auc-score-for-regression-models-c0be4fdf76bb except that the metric is inverted and can be weighted by sample weights.
 
-## 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]=[])
+## 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))
 
 ***This method fits the model to data.***
 
@@ -82,6 +82,9 @@ An optional list of integers specifying the indexes of predictors (columns) in *
 #### 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.
 
+#### 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).
+
 
 ## Method: predict(X:npt.ArrayLike, cap_predictions_to_minmax_in_training:bool=True)
 

aplr/aplr.py

@@ -48,9 +48,9 @@ def __set_params_cpp(self):
         self.APLRRegressor.tweedie_power=self.tweedie_power
         self.APLRRegressor.validation_tuning_metric=self.validation_tuning_metric
 
-    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]=[]):
+    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)):
         self.__set_params_cpp()
-        self.APLRRegressor.fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints)
+        self.APLRRegressor.fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group)
 
     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

@@ -46,15 +46,20 @@ class APLRRegressor
     std::vector<size_t> predictor_indexes;
     std::vector<size_t> prioritized_predictors_indexes;
     std::vector<int> monotonic_constraints;
+    VectorXi group_train;
+    VectorXi group_validation;
+    std::set<int> unique_groups_train;
+    std::set<int> unique_groups_validation;
 
     //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 std::vector<int> &monotonic_constraints, const VectorXi &group);
     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 define_training_and_validation_sets(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight, const std::vector<size_t> &validation_set_indexes);
+    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);
     bool check_if_base_term_has_only_one_unique_value(size_t base_term);
     void add_term_to_terms_eligible_current(Term &term);
@@ -146,7 +151,8 @@ class APLRRegressor
     APLRRegressor(const APLRRegressor &other);
     ~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 std::vector<size_t> &prioritized_predictors_indexes={}, const std::vector<int> &monotonic_constraints={},
+        const VectorXi &group=VectorXi(0));
     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);
@@ -207,13 +213,13 @@ 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 std::vector<int> &monotonic_constraints, const VectorXi &group)
 {
     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);
-    define_training_and_validation_sets(X,y,sample_weight,validation_set_indexes);
+    validate_input_to_fit(X,y,sample_weight,X_names,validation_set_indexes,prioritized_predictors_indexes,monotonic_constraints,group);
+    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);
     execute_boosting_steps();
@@ -239,7 +245,9 @@ void APLRRegressor::throw_error_if_family_does_not_exist()
     else if(family=="gamma")
         family_exists=true;
     else if(family=="tweedie")
-        family_exists=true;        
+        family_exists=true;
+    else if(family=="group_gaussian")
+        family_exists=true;
     if(!family_exists)
         throw std::runtime_error("Family "+family+" is not available in APLR.");   
 }
@@ -268,7 +276,7 @@ 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 std::vector<size_t> &prioritized_predictors_indexes, const std::vector<int> &monotonic_constraints, const VectorXi &group)
 {
     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.");
@@ -281,6 +289,8 @@ void APLRRegressor::validate_input_to_fit(const MatrixXd &X,const VectorXd &y,co
     throw_error_if_monotonic_constraints_has_invalid_indexes(X, monotonic_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()};
+    if(group_is_of_incorrect_size) throw std::runtime_error("When family is group_gaussian then y and group must have the same number of rows.");
 }
 
 void APLRRegressor::throw_error_if_validation_set_indexes_has_invalid_indexes(const VectorXd &y, const std::vector<size_t> &validation_set_indexes)
@@ -381,7 +391,8 @@ void APLRRegressor::throw_error_if_sample_weight_contains_invalid_values(const V
     }
 }
 
-void APLRRegressor::define_training_and_validation_sets(const MatrixXd &X,const VectorXd &y,const VectorXd &sample_weight, const std::vector<size_t> &validation_set_indexes)
+void APLRRegressor::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)
 {
     size_t y_size{static_cast<size_t>(y.size())};
     std::vector<size_t> train_indexes;
@@ -440,6 +451,16 @@ void APLRRegressor::define_training_and_validation_sets(const MatrixXd &X,const
             sample_weight_train[i]=sample_weight[train_indexes[i]];
         }
     }
+    bool groups_are_provided{group.size()>0};
+    if(groups_are_provided)
+    {
+        group_train.resize(train_indexes.size());
+        for (size_t i = 0; i < train_indexes.size(); ++i)
+        {
+            group_train[i]=group[train_indexes[i]];
+        }
+        unique_groups_train = get_unique_integers(group_train);
+    }
     //Populating test matrices
     for (size_t i = 0; i < validation_indexes.size(); ++i)
     {
@@ -454,6 +475,15 @@ void APLRRegressor::define_training_and_validation_sets(const MatrixXd &X,const
             sample_weight_validation[i]=sample_weight[validation_indexes[i]];
         }
     }
+    if(groups_are_provided)
+    {
+        group_validation.resize(validation_indexes.size());
+        for (size_t i = 0; i < validation_indexes.size(); ++i)
+        {
+            group_validation[i]=group[validation_indexes[i]];
+        }
+        unique_groups_validation = get_unique_integers(group_validation);
+    }
 }
 
 void APLRRegressor::scale_training_observations_if_using_log_link_function()
@@ -561,6 +591,21 @@ VectorXd APLRRegressor::calculate_neg_gradient_current()
         output=(y_train.array() - predictions_current.array()) / predictions_current.array() / predictions_current.array();
     else if(family=="tweedie")
         output=(y_train.array()-predictions_current.array()).array() * predictions_current.array().pow(-tweedie_power);
+    else if(family=="group_gaussian")
+    {
+        GroupData group_residuals_and_count{calculate_group_errors_and_count(y_train,predictions_current,group_train,unique_groups_train)};
+
+        for(int unique_group_value:unique_groups_train)
+        {
+            group_residuals_and_count.error[unique_group_value] /= group_residuals_and_count.count[unique_group_value];
+        }
+
+        output = VectorXd(y_train.rows());
+        for (Eigen::Index i = 0; i < y_train.size(); ++i)
+        {
+            output[i] = group_residuals_and_count.error[group_train[i]];
+        }
+    }
 
     if(link_function!="identity")
         output=output.array()*differentiate_predictions().array();
@@ -1004,7 +1049,7 @@ void APLRRegressor::calculate_and_validate_validation_error(size_t boosting_step
 void APLRRegressor::calculate_validation_error(size_t boosting_step, const VectorXd &predictions)
 {
     if(validation_tuning_metric=="default")
-        validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,family,tweedie_power),sample_weight_validation);
+        validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,family,tweedie_power,group_validation,unique_groups_validation),sample_weight_validation);
     else if(validation_tuning_metric=="mse")
         validation_error_steps[boosting_step]=calculate_mean_error(calculate_errors(y_validation,predictions,sample_weight_validation,FAMILY_GAUSSIAN),sample_weight_validation);
     else if(validation_tuning_metric=="mae")
@@ -1274,6 +1319,10 @@ void APLRRegressor::cleanup_after_fit()
     predictor_indexes.clear();
     prioritized_predictors_indexes.clear();
     monotonic_constraints.clear();
+    group_train.resize(0);
+    group_validation.resize(0);
+    unique_groups_train.clear();
+    unique_groups_validation.clear();
 }
 
 VectorXd APLRRegressor::predict(const MatrixXd &X, bool cap_predictions_to_minmax_in_training)

cpp/functions.h

@@ -1,6 +1,5 @@
 #pragma once
 #include <limits>
-#include "../dependencies/eigen-master/Eigen/Dense"
 #include <numeric> //std::iota
 #include <algorithm> //std::sort, std::stable_sort
 #include <vector>
@@ -9,6 +8,9 @@
 #include <thread>
 #include <future>
 #include <random>
+#include <set>
+#include <map>
+#include "../dependencies/eigen-master/Eigen/Dense"
 #include "constants.h"
 
 using namespace Eigen;
@@ -37,6 +39,12 @@ static bool is_approximately_zero(TReal a, TReal tolerance = std::numeric_limits
     return false;
 }
 
+std::set<int> get_unique_integers(const VectorXi &int_vector)
+{
+    std::set<int> unique_integers{int_vector.begin(),int_vector.end()};
+    return unique_integers;
+}
+
 double set_error_to_infinity_if_invalid(double error)
 {
     bool error_is_invalid{!std::isfinite(error)};
@@ -77,7 +85,49 @@ VectorXd calculate_tweedie_errors(const VectorXd &y,const VectorXd &predicted,do
     return errors;
 }
 
-VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),const std::string &family="gaussian",double tweedie_power=1.5)
+struct GroupData
+{
+    std::map<int,double> error;
+    std::map<int,size_t> count;
+};
+
+GroupData calculate_group_errors_and_count(const VectorXd &y,const VectorXd &predicted,const VectorXi &group, const std::set<int> &unique_groups)
+{
+    GroupData group_data;
+    for(int unique_group_value:unique_groups)
+    {
+        group_data.error[unique_group_value]=0.0;
+        group_data.count[unique_group_value]=0;
+    }
+    for (Eigen::Index i = 0; i < group.size(); ++i)
+    {
+        group_data.error[group[i]] += y[i] - predicted[i];
+        group_data.count[group[i]] += 1;
+    } 
+    return group_data;
+}
+
+VectorXd calculate_group_gaussian_errors(const VectorXd &y,const VectorXd &predicted,const VectorXi &group, const std::set<int> &unique_groups)
+{
+    GroupData group_residuals_and_count{calculate_group_errors_and_count(y,predicted,group,unique_groups)};
+
+    for(int unique_group_value:unique_groups)
+    {
+        group_residuals_and_count.error[unique_group_value] *= group_residuals_and_count.error[unique_group_value];
+        group_residuals_and_count.error[unique_group_value] /= group_residuals_and_count.count[unique_group_value];
+    }
+
+    VectorXd errors(y.rows());
+    for (Eigen::Index i = 0; i < y.size(); ++i)
+    {
+        errors[i] = group_residuals_and_count.error[group[i]];
+    }
+    
+    return errors;
+}
+
+
+VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const VectorXd &sample_weight=VectorXd(0),const std::string &family="gaussian",double tweedie_power=1.5, const VectorXi &group=VectorXi(0), const std::set<int> &unique_groups={})
 {   
     VectorXd errors;
     if(family=="gaussian")
@@ -90,7 +140,9 @@ VectorXd calculate_errors(const VectorXd &y,const VectorXd &predicted,const Vect
         errors=calculate_gamma_errors(y,predicted);
     else if(family=="tweedie")
         errors=calculate_tweedie_errors(y,predicted,tweedie_power);
-    
+    else if(family=="group_gaussian")
+        errors=calculate_group_gaussian_errors(y,predicted,group,unique_groups);
+
     if(sample_weight.size()>0)
         errors=errors.array()*sample_weight.array();
 

cpp/pythonbinding.cpp

@@ -22,7 +22,8 @@ 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::call_guard<py::scoped_ostream_redirect,py::scoped_estream_redirect>())
+            py::arg("monotonic_constraints")=std::vector<int>(),py::arg("group")=VectorXi(0),
+            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"))
         .def("calculate_local_feature_importance",&APLRRegressor::calculate_local_feature_importance,py::arg("X"))