Interfacing with scikit-learn - survivalpredict documentation

Survivalpredict’s SklearnSurvivalPipeline class allows us to interface directly with scikit-learn’s greater ecosystem. There is a whole ecosystem centered around hyperparameter tuning, feature selection, and scaling out compute for scikit-learn compatible code that we can interface with.

import warnings
warnings.filterwarnings('ignore')

import numpy as np
import itertools
import pandas as pd

from sklearn.model_selection import GridSearchCV , cross_val_score
from sklearn.compose import make_column_transformer, make_column_selector
from sklearn.preprocessing import StandardScaler, OneHotEncoder

from mlxtend.feature_selection import SequentialFeatureSelector
from mlxtend.plotting import plot_sequential_feature_selection as plot_sfs

from survivalpredict.estimators import  CoxProportionalHazard,KaplanMeierSurvivalEstimator
from survivalpredict.metrics import integrated_brier_score_administrative_sklearn_scorer 
from survivalpredict.pipeline import SklearnSurvivalPipeline,build_sklearn_pipeline_target,make_sklearn_survival_pipeline
from survivalpredict.strata_preprocessing import StrataBuilderEncoder, StrataColumnTransformer, make_strata_column_transformer

from SurvSet.data import SurvLoader
loader = SurvLoader()

Let’s start by loading in the data, using SurvSet. Notice that we are keeping the data as a dataframe. This allows us to retain feature names, enabling us to call out specific strata feature names when tuning for the best-performing strata.

df = loader.load_dataset("support2")["df"]
df = df.dropna()

times = df["time"].to_numpy().astype(np.int64)
events = df["event"].to_numpy().astype(np.bool_)
X_df = df[list(set(df.columns).difference(set(("pid", "event", "time"))))]

Next, we are going to determine the maximum point in time to observe for evaluating the model. As time goes forward, a smaller and smaller percentage of individuals will still be in the study, including the tail-end points in the integrated Brier score, which will artificially inflate the score, not highlighting improvements in performance in points in time that include most individuals.

max_time = np.percentile(times, 85).round().astype(np.int64)

The ‘build_sklearn_pipeline_target’ function allows us to take non-feature vectors, like times, events, times_start(for left censoring), and predetermined strata into a single vector that can be passed into scikit-learn classes. The output of ‘build_sklearn_pipeline_target’ will function as the ‘y’ vector for the SklearnSurvivalPipeline class.

Next, we will build a pipeline that will turn ‘fac_sex’, ‘fac_income’ columns into the strata, then preprocess the rest of the features, and finally train a Cox model.

y = build_sklearn_pipeline_target(times, events)

strata_transformers = make_strata_column_transformer(
    (StrataBuilderEncoder(), ["fac_sex", "fac_income"])
)

column_transformers = make_column_transformer(
    (StandardScaler(), make_column_selector(dtype_include=np.number)),
    (OneHotEncoder(), make_column_selector(dtype_include=[object, "string"])),
)

pipeline_cox = make_sklearn_survival_pipeline(
    strata_transformers,
    column_transformers,
    CoxProportionalHazard(),
    max_time=max_time,
)

pipeline_cox

We can now use sk’s native cross_val_score, but we will need to pass in our own scorer

cross_val_score(
    pipeline_cox,
    X_df,
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
    error_score='raise',
).mean()

np.float64(-155.14165982224978)

#we can also build a cox model without strata

cross_val_score(
    make_sklearn_survival_pipeline(
        column_transformers,
        CoxProportionalHazard(),
        max_time=max_time,
    ),
    X_df,
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
).mean()

np.float64(-155.57987707062597)

cross_val_score(
    make_sklearn_survival_pipeline(column_transformers,KaplanMeierSurvivalEstimator(),max_time=max_time),
    X_df,
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
).mean()

np.float64(-179.90976086220107)

It looks like our first Cox models are not bad, we are performing better than KaplanMeier. We can now start hyperparameters and strata tuning.

#getting combinations of all categorical data

cat_cols = X_df.select_dtypes(pd.CategoricalDtype()).columns

#in case we wana look at strata combinations
cat_col_combinations = list(
    itertools.chain.from_iterable(
        itertools.combinations(cat_cols, i) for i in range(1, 3)
    )
)

#, this makes gridsearch run faster
#cat_col_combinations = cat_cols 

#you can take a look at all pram names and prameters with pipeline_cox.get_params()

param_grid = {
    "stratacolumntransformer__stratabuilderencoder__columns": cat_col_combinations,
    "coxproportionalhazard__alpha": [0, 10, 25, 50, 75, 100],
}

#running the sklearn native GridSearchCV, you can replace GridSearchCV with your favorite hyperprameter tunning class
gs = GridSearchCV(
    pipeline_cox,
    param_grid,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
    n_jobs=5
)

#not all strata will be present in each cv split, this will cause errors/warnings; this should not be an issue, as sparse strata are not a good strata regardless
import warnings

with warnings.catch_warnings():
    warnings.simplefilter("ignore")

    gs.fit(X_df, y)

cv_results = pd.DataFrame(gs.cv_results_)

cv_results.sort_values('rank_test_score').head(5)[['params','mean_test_score','std_test_score','rank_test_score']]

	params	mean_test_score	std_test_score	rank_test_score
101	{'coxproportionalhazard__alpha': 10, 'strataco...	-119.853315	5.188538	1
68	{'coxproportionalhazard__alpha': 10, 'strataco...	-120.893513	6.171047	2
2	{'coxproportionalhazard__alpha': 0, 'stratacol...	-121.189043	5.495687	3
134	{'coxproportionalhazard__alpha': 25, 'strataco...	-121.211370	7.280129	4
233	{'coxproportionalhazard__alpha': 50, 'strataco...	-121.580471	8.243867	5

#the 'best' parameters.
cv_results[cv_results['rank_test_score'] == 1]['params'].values[0]

{'coxproportionalhazard__alpha': 10,
 'stratacolumntransformer__stratabuilderencoder__columns': ('fac_sfdm2',
  'fac_sex')}

strata_transformers = make_strata_column_transformer(
    (StrataBuilderEncoder(), ('fac_sfdm2', 'fac_sex'))
)

sp_pipeline_cox = make_sklearn_survival_pipeline(
    strata_transformers,
    column_transformers,
    CoxProportionalHazard(alpha=10),
    max_time=max_time,
)

cross_val_score(
    sp_pipeline_cox,
    X_df,
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
).mean()

np.float64(-119.8533147038668)

X_np = column_transformers.fit_transform(X_df)

pipeline_cox2 = make_sklearn_survival_pipeline(
    CoxProportionalHazard(),
    max_time=max_time,
)

sfs = SequentialFeatureSelector(pipeline_cox2, 
                          k_features= 1, 
                          forward=False,  
                          scoring=integrated_brier_score_administrative_sklearn_scorer, 
                          cv=5, 
                          n_jobs=5)

sfs.fit(X_np,y)

SequentialFeatureSelector(estimator=SklearnSurvivalPipeline(max_time=818,
                                                            steps=[('coxproportionalhazard',
                                                                    CoxProportionalHazard())]),
                          forward=False, k_features=(1, 1), n_jobs=5,
                          scoring=make_scorer(integrated_brier_score_administrative_sklearn_metric, greater_is_better=False, response_method='predict'))

In a Jupyter environment, please rerun this cell to show the HTML representation or trust the notebook.
On GitHub, the HTML representation is unable to render, please try loading this page with nbviewer.org.

fig = plot_sfs(sfs.get_metric_dict(), kind='std_err')

../../_images/a335c72dd4cbd9c3aca6139af319287696a4632ec4bd7557e9417667d660b099.png

step_wise_feature_selection_results = pd.DataFrame.from_records(sfs.subsets_).T

step_wise_feature_selection_results.sort_values('avg_score',ascending=False).head(5)

	feature_idx	cv_scores	avg_score	feature_names
7	(0, 3, 4, 15, 17, 18, 19)	[-143.56667782179701, -136.25594055650689, -15...	-144.112512	(0, 3, 4, 15, 17, 18, 19)
6	(0, 3, 4, 17, 18, 19)	[-143.83516942714164, -135.87268885869483, -15...	-144.263809	(0, 3, 4, 17, 18, 19)
8	(0, 1, 3, 4, 15, 17, 18, 19)	[-141.46966297859967, -138.35774779284924, -15...	-144.270329	(0, 1, 3, 4, 15, 17, 18, 19)
9	(0, 1, 3, 4, 15, 17, 18, 19, 23)	[-149.19738821435038, -135.50584628943156, -15...	-144.318642	(0, 1, 3, 4, 15, 17, 18, 19, 23)
10	(0, 1, 3, 4, 15, 16, 17, 18, 19, 23)	[-149.00292841678186, -134.64053897377698, -15...	-144.356467	(0, 1, 3, 4, 15, 16, 17, 18, 19, 23)

best_features_idxs = list(step_wise_feature_selection_results.iloc[7]['feature_idx'])                 

best_features_names = column_transformers.get_feature_names_out()[list(best_features_idxs)]
best_features_names

array(['standardscaler__num_hrt', 'standardscaler__num_surv2m',
       'standardscaler__num_scoma', 'standardscaler__num_sps',
       'standardscaler__num_adlp', 'standardscaler__num_adls',
       'standardscaler__num_num_co', 'standardscaler__num_surv6m'],
      dtype=object)

cross_val_score(
    pipeline_cox2,
    X_np[:,best_features_idxs],
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
).mean()

np.float64(-144.27032908030452)

best_features_col_names = [i.split('__')[1] for i in column_transformers.get_feature_names_out()[list(best_features_idxs)]]
strat_cols = ['fac_sfdm2', 'fac_sex']
cols = best_features_col_names + strat_cols

X_df2 = X_df[cols]


strata_transformers = make_strata_column_transformer(
    (StrataBuilderEncoder(), strat_cols)
)

column_transformers = make_column_transformer(
    (StandardScaler(), best_features_col_names ),
)

pipeline_cox3 = make_sklearn_survival_pipeline(
    strata_transformers,
    column_transformers,
    CoxProportionalHazard(),
    max_time=max_time,
)

Using both the best strata and best features, we are able to get the best performing model so far.

cross_val_score(
    pipeline_cox3,
    X_df2,
    y,
    scoring=integrated_brier_score_administrative_sklearn_scorer,
).mean()

np.float64(-112.16788210122328)

	transformers transformers: list of tuples List of (name, transformer, columns) tuples specifying the transformer objects to be applied to subsets of the data. name : str Like in Pipeline and FeatureUnion, this allows the transformer and its parameters to be set using ``set_params`` and searched in grid search. transformer : {'drop', 'passthrough'} or estimator Estimator must support :term:`fit` and :term:`transform`. Special-cased strings 'drop' and 'passthrough' are accepted as well, to indicate to drop the columns or to pass them through untransformed, respectively. columns : str, array-like of str, int, array-like of int, array-like of bool, slice or callable Indexes the data on its second axis. Integers are interpreted as positional columns, while strings can reference DataFrame columns by name. A scalar string or int should be used where ``transformer`` expects X to be a 1d array-like (vector), otherwise a 2d array will be passed to the transformer. A callable is passed the input data `X` and can return any of the above. To select multiple columns by name or dtype, you can use :obj:`make_column_selector`.	[('standardscaler', ...), ('onehotencoder', ...)]
	remainder remainder: {'drop', 'passthrough'} or estimator, default='drop' By default, only the specified columns in `transformers` are transformed and combined in the output, and the non-specified columns are dropped. (default of ``'drop'``). By specifying ``remainder='passthrough'``, all remaining columns that were not specified in `transformers`, but present in the data passed to `fit` will be automatically passed through. This subset of columns is concatenated with the output of the transformers. For dataframes, extra columns not seen during `fit` will be excluded from the output of `transform`. By setting ``remainder`` to be an estimator, the remaining non-specified columns will use the ``remainder`` estimator. The estimator must support :term:`fit` and :term:`transform`. Note that using this feature requires that the DataFrame columns input at :term:`fit` and :term:`transform` have identical order.	'drop'
	sparse_threshold sparse_threshold: float, default=0.3 If the output of the different transformers contains sparse matrices, these will be stacked as a sparse matrix if the overall density is lower than this value. Use ``sparse_threshold=0`` to always return dense. When the transformed output consists of all dense data, the stacked result will be dense, and this keyword will be ignored.	0.3
	n_jobs n_jobs: int, default=None Number of jobs to run in parallel. ``None`` means 1 unless in a :obj:`joblib.parallel_backend` context. ``-1`` means using all processors. See :term:`Glossary ` for more details.	None
	transformer_weights transformer_weights: dict, default=None Multiplicative weights for features per transformer. The output of the transformer is multiplied by these weights. Keys are transformer names, values the weights.	None
	verbose verbose: bool, default=False If True, the time elapsed while fitting each transformer will be printed as it is completed.	False
	verbose_feature_names_out verbose_feature_names_out: bool, str or Callable[[str, str], str], default=True - If True, :meth:`ColumnTransformer.get_feature_names_out` will prefix all feature names with the name of the transformer that generated that feature. It is equivalent to setting `verbose_feature_names_out="{transformer_name}__{feature_name}"`. - If False, :meth:`ColumnTransformer.get_feature_names_out` will not prefix any feature names and will error if feature names are not unique. - If ``Callable[[str, str], str]``, :meth:`ColumnTransformer.get_feature_names_out` will rename all the features using the name of the transformer. The first argument of the callable is the transformer name and the second argument is the feature name. The returned string will be the new feature name. - If ``str``, it must be a string ready for formatting. The given string will be formatted using two field names: ``transformer_name`` and ``feature_name``. e.g. ``"{feature_name}__{transformer_name}"``. See :meth:`str.format` method from the standard library for more info. .. versionadded:: 1.0 .. versionchanged:: 1.6 `verbose_feature_names_out` can be a callable or a string to be formatted.	True
	force_int_remainder_cols force_int_remainder_cols: bool, default=False This parameter has no effect. .. note:: If you do not access the list of columns for the remainder columns in the `transformers_` fitted attribute, you do not need to set this parameter. .. versionadded:: 1.5 .. versionchanged:: 1.7 The default value for `force_int_remainder_cols` will change from `True` to `False` in version 1.7. .. deprecated:: 1.7 `force_int_remainder_cols` is deprecated and will be removed in 1.9.	'deprecated'

	copy copy: bool, default=True If False, try to avoid a copy and do inplace scaling instead. This is not guaranteed to always work inplace; e.g. if the data is not a NumPy array or scipy.sparse CSR matrix, a copy may still be returned.	True
	with_mean with_mean: bool, default=True If True, center the data before scaling. This does not work (and will raise an exception) when attempted on sparse matrices, because centering them entails building a dense matrix which in common use cases is likely to be too large to fit in memory.	True
	with_std with_std: bool, default=True If True, scale the data to unit variance (or equivalently, unit standard deviation).	True

	categories categories: 'auto' or a list of array-like, default='auto' Categories (unique values) per feature: - 'auto' : Determine categories automatically from the training data. - list : ``categories[i]`` holds the categories expected in the ith column. The passed categories should not mix strings and numeric values within a single feature, and should be sorted in case of numeric values. The used categories can be found in the ``categories_`` attribute. .. versionadded:: 0.20	'auto'
	drop drop: {'first', 'if_binary'} or an array-like of shape (n_features,), default=None Specifies a methodology to use to drop one of the categories per feature. This is useful in situations where perfectly collinear features cause problems, such as when feeding the resulting data into an unregularized linear regression model. However, dropping one category breaks the symmetry of the original representation and can therefore induce a bias in downstream models, for instance for penalized linear classification or regression models. - None : retain all features (the default). - 'first' : drop the first category in each feature. If only one category is present, the feature will be dropped entirely. - 'if_binary' : drop the first category in each feature with two categories. Features with 1 or more than 2 categories are left intact. - array : ``drop[i]`` is the category in feature ``X[:, i]`` that should be dropped. When `max_categories` or `min_frequency` is configured to group infrequent categories, the dropping behavior is handled after the grouping. .. versionadded:: 0.21 The parameter `drop` was added in 0.21. .. versionchanged:: 0.23 The option `drop='if_binary'` was added in 0.23. .. versionchanged:: 1.1 Support for dropping infrequent categories.	None
	sparse_output sparse_output: bool, default=True When ``True``, it returns a :class:`scipy.sparse.csr_matrix`, i.e. a sparse matrix in "Compressed Sparse Row" (CSR) format. .. versionadded:: 1.2 `sparse` was renamed to `sparse_output`	True
	dtype dtype: number type, default=np.float64 Desired dtype of output.	<class 'numpy.float64'>
	handle_unknown handle_unknown: {'error', 'ignore', 'infrequent_if_exist', 'warn'}, default='error' Specifies the way unknown categories are handled during :meth:`transform`. - 'error' : Raise an error if an unknown category is present during transform. - 'ignore' : When an unknown category is encountered during transform, the resulting one-hot encoded columns for this feature will be all zeros. In the inverse transform, an unknown category will be denoted as None. - 'infrequent_if_exist' : When an unknown category is encountered during transform, the resulting one-hot encoded columns for this feature will map to the infrequent category if it exists. The infrequent category will be mapped to the last position in the encoding. During inverse transform, an unknown category will be mapped to the category denoted `'infrequent'` if it exists. If the `'infrequent'` category does not exist, then :meth:`transform` and :meth:`inverse_transform` will handle an unknown category as with `handle_unknown='ignore'`. Infrequent categories exist based on `min_frequency` and `max_categories`. Read more in the :ref:`User Guide `. - 'warn' : When an unknown category is encountered during transform a warning is issued, and the encoding then proceeds as described for `handle_unknown="infrequent_if_exist"`. .. versionchanged:: 1.1 `'infrequent_if_exist'` was added to automatically handle unknown categories and infrequent categories. .. versionadded:: 1.6 The option `"warn"` was added in 1.6.	'error'
	min_frequency min_frequency: int or float, default=None Specifies the minimum frequency below which a category will be considered infrequent. - If `int`, categories with a smaller cardinality will be considered infrequent. - If `float`, categories with a smaller cardinality than `min_frequency * n_samples` will be considered infrequent. .. versionadded:: 1.1 Read more in the :ref:`User Guide `.	None
	max_categories max_categories: int, default=None Specifies an upper limit to the number of output features for each input feature when considering infrequent categories. If there are infrequent categories, `max_categories` includes the category representing the infrequent categories along with the frequent categories. If `None`, there is no limit to the number of output features. .. versionadded:: 1.1 Read more in the :ref:`User Guide `.	None
	feature_name_combiner feature_name_combiner: "concat" or callable, default="concat" Callable with signature `def callable(input_feature, category)` that returns a string. This is used to create feature names to be returned by :meth:`get_feature_names_out`. `"concat"` concatenates encoded feature name and category with `feature + "_" + str(category)`.E.g. feature X with values 1, 6, 7 create feature names `X_1, X_6, X_7`. .. versionadded:: 1.3	'concat'

	alpha	0.0
	max_iter	100
	ties	'breslow'
	tol	1e-09

	estimator	SklearnSurviv...nalHazard())])
	k_features	(1, ...)
	forward	False
	floating	False
	verbose	0
	scoring	make_scorer(i...hod='predict')
	cv	5
	n_jobs	5
	pre_dispatch	'2*n_jobs'
	clone_estimator	True
	fixed_features	None
	feature_groups	None

	steps	[('stratacolumntransformer', ...), ('columntransformer', ...), ...]
	max_time	818
	memory	None

	strata_transformers	[('stratabuilderencoder', ...)]
	stratabuilderencoder__columns	['fac_sex', 'fac_income']

Interfacing with scikit-learn¶