We’ll start with a probabilistic regression example on the diabetes dataset from sklearn:
from ngboost import NGBRegressor
from sklearn.datasets import load_diabetes
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error
X, Y = load_diabetes(return_X_y=True)
X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2)
ngb = NGBRegressor().fit(X_train, Y_train)
Y_preds = ngb.predict(X_test)
Y_dists = ngb.pred_dist(X_test)
# test Mean Squared Error
test_MSE = mean_squared_error(Y_preds, Y_test)
print('Test MSE', test_MSE)
# test Negative Log Likelihood
test_NLL = -Y_dists.logpdf(Y_test).mean()
print('Test NLL', test_NLL)[iter 0] loss=5.7687 val_loss=0.0000 scale=1.0000 norm=66.8065
[iter 100] loss=5.4122 val_loss=0.0000 scale=2.0000 norm=92.0139
[iter 200] loss=5.1742 val_loss=0.0000 scale=1.0000 norm=36.9730
[iter 300] loss=5.0397 val_loss=0.0000 scale=1.0000 norm=33.5210
[iter 400] loss=4.9400 val_loss=0.0000 scale=2.0000 norm=62.7756
Test MSE 3508.531873714694
Test NLL 6.035871970356794
Getting the estimated distributional parameters at a set of points is easy. This returns the predicted mean and standard deviation of the first five observations in the test set:
Y_dists[0:5].params{'loc': array([192.20800911, 105.3969346 , 125.9068854 , 90.40939382,
74.68291579]),
'scale': array([37.28866 , 26.1406384 , 28.36160015, 34.12044 , 24.49221427])}Distributions¶
NGBoost can be used with a variety of distributions, broken down into those for regression (support on an infinite set) and those for classification (support on a finite set).
Regression Distributions¶
| Distribution | Parameters | Implemented Scores | Reference |
|---|---|---|---|
Normal | loc, scale | LogScore, CRPScore | scipy.stats normal |
LogNormal | s, scale | LogScore, CRPScore | scipy.stats lognormal |
Exponential | scale | LogScore, CRPScore | scipy.stats exponential |
Regression distributions can be used through the NGBRegressor() constructor by passing the appropriate class as the Dist argument. Normal is the default.
from ngboost.distns import Exponential, Normal
X, Y = load_diabetes(return_X_y=True)
X_reg_train, X_reg_test, Y_reg_train, Y_reg_test = train_test_split(X, Y, test_size=0.2)
ngb_norm = NGBRegressor(Dist=Normal, verbose=False).fit(X_reg_train, Y_reg_train)
ngb_exp = NGBRegressor(Dist=Exponential, verbose=False).fit(X_reg_train, Y_reg_train)There are two prediction methods for NGBRegressor objects: predict(), which returns point predictions as one would expect from a standard regressor, and pred_dist(), which returns a distribution object representing the conditional distribution of at the points in the test set.
ngb_norm.predict(X_reg_test)[0:5]array([ 85.68847007, 244.23240372, 165.14133044, 173.26063598,
248.59937868])ngb_exp.predict(X_reg_test)[0:5]array([ 80.23765636, 236.49217244, 166.17076131, 139.26094569,
261.08260799])ngb_exp.pred_dist(X_reg_test)[0:5].params{'scale': array([ 80.23765636, 236.49217244, 166.17076131, 139.26094569,
261.08260799])}Survival Regression¶
NGBoost supports analyses of right-censored data. Any distribution that can be used for regression in NGBoost can also be used for survival analysis in theory, but this requires the implementation of the right-censored version of the appropriate score. At the moment, LogNormal and Exponential have these scores implemented. To do survival analysis, use NGBSurvival and pass both the time-to-event (or censoring) and event indicator vectors to fit():
import numpy as np
from ngboost import NGBSurvival
from ngboost.distns import LogNormal
X, Y = load_diabetes(return_X_y=True)
X_surv_train, X_surv_test, Y_surv_train, Y_surv_test = train_test_split(X, Y, test_size=0.2)
# introduce administrative censoring to simulate survival data
T_surv_train = np.minimum(Y_surv_train, 250) # time of an event or censoring
E_surv_train = Y_surv_train <= 250 # 1 if T[i] is the time of an event, 0 if it's a time of censoring
ngb = NGBSurvival(Dist=LogNormal).fit(X_surv_train, T_surv_train, E_surv_train)[iter 0] loss=5.0123 val_loss=0.0000 scale=2.0000 norm=1.5198
[iter 100] loss=4.6344 val_loss=0.0000 scale=1.0000 norm=0.5617
[iter 200] loss=4.4640 val_loss=0.0000 scale=1.0000 norm=0.5131
[iter 300] loss=4.3313 val_loss=0.0000 scale=2.0000 norm=0.9611
[iter 400] loss=4.2191 val_loss=0.0000 scale=1.0000 norm=0.4611
The scores currently implemented assume that the censoring is independent of survival, conditional on the observed predictors.
Classification Distributions¶
| Distribution | Parameters | Implemented Scores | Reference |
|---|---|---|---|
k_categorical(K) | p0, p1... p{K-1} | LogScore | Categorical distribution on Wikipedia |
Bernoulli | p | LogScore | Bernoulli distribution on Wikipedia |
Classification distributions can be used through the NGBClassifier() constructor by passing the appropriate class as the Dist argument. Bernoulli is the default and is equivalent to k_categorical(2).
from ngboost import NGBClassifier
from ngboost.distns import k_categorical, Bernoulli
from sklearn.datasets import load_breast_cancer
X, y = load_breast_cancer(return_X_y=True)
y[0:15] = 2 # artificially make this a 3-class problem instead of a 2-class problem
X_cls_train, X_cls_test, Y_cls_train, Y_cls_test = train_test_split(X, y, test_size=0.2)
ngb_cat = NGBClassifier(Dist=k_categorical(3), verbose=False) # tell ngboost that there are 3 possible outcomes
_ = ngb_cat.fit(X_cls_train, Y_cls_train) # Y should have only 3 values: {0,1,2}When using NGBoost for classification, the outcome vector Y must consist only of integers from 0 to K-1, where K is the total number of classes. This is consistent with the classification standards in sklearn.
NGBClassifier objects have three prediction methods: predict() returns the most likely class, predict_proba() returns the class probabilities, and pred_dist() returns the distribution object.
ngb_cat.predict(X_cls_test)[0:5]array([0, 1, 1, 1, 0])ngb_cat.predict_proba(X_cls_test)[0:5]array([[9.48163772e-01, 6.77474734e-03, 4.50614804e-02],
[3.28800280e-01, 6.48155593e-01, 2.30441272e-02],
[3.59619610e-03, 9.96151763e-01, 2.52041150e-04],
[3.59615491e-03, 9.96140351e-01, 2.63494289e-04],
[9.91620820e-01, 7.08525332e-03, 1.29392628e-03]])ngb_cat.pred_dist(X_cls_test)[0:5].params{'p0': array([0.94816377, 0.32880028, 0.0035962 , 0.00359615, 0.99162082]),
'p1': array([0.00677475, 0.64815559, 0.99615176, 0.99614035, 0.00708525]),
'p2': array([0.04506148, 0.02304413, 0.00025204, 0.00026349, 0.00129393])}Scores¶
NGBoost supports the log score (LogScore, also known as negative log-likelihood) and CRPS (CRPScore), although each score may not be implemented for each distribution. The score is specified by the Score argument in the constructor.
from ngboost.scores import LogScore, CRPScore
NGBRegressor(Dist=Exponential, Score=CRPScore, verbose=False).fit(X_reg_train, Y_reg_train)
NGBClassifier(Dist=k_categorical(3), Score=LogScore, verbose=False).fit(X_cls_train, Y_cls_train)Base Learners¶
NGBoost can be used with any sklearn regressor as the base learner, specified with the Base argument. The default is a depth-3 regression tree.
from sklearn.tree import DecisionTreeRegressor
learner = DecisionTreeRegressor(criterion='friedman_mse', max_depth=5)
NGBSurvival(Dist=Exponential, Score=CRPScore, Base=learner, verbose=False).fit(X_surv_train, T_surv_train, E_surv_train)Other Arguments¶
The learning rate, number of estimators, minibatch fraction, and column subsampling are also easily adjusted:
ngb = NGBRegressor(n_estimators=100, learning_rate=0.01,
minibatch_frac=0.5, col_sample=0.5)
ngb.fit(X_reg_train, Y_reg_train)[iter 0] loss=5.7068 val_loss=0.0000 scale=1.0000 norm=60.8766
Sample weights (for training) are set using the sample_weight argument to fit.
ngb = NGBRegressor(n_estimators=100, learning_rate=0.01,
minibatch_frac=0.5, col_sample=0.5)
weights = np.random.random(Y_reg_train.shape)
ngb.fit(X_reg_train, Y_reg_train, sample_weight=weights)[iter 0] loss=5.7231 val_loss=0.0000 scale=1.0000 norm=65.9016