The HRT for mixed data types (Python implementation)


In my last post I showed how the holdout random test (HRT) could be used to obtain valid p-values for any machine learning model by sampling from the conditional distribution of the design matrix. Like the permutation-type approaches used to assess variable importance for decision trees, this method sees whether a measure of performance accuracy declines when a column of the data has its values shuffled. However these ad-hoc permutation approaches lack statistical rigor and will not obtain a valid inferential assessment, even asymptotically, as the non-permuted columns of the data are not conditioned on. For example, if two features are correlated with the data, but only one has a statistical relationship with the response, then naive permutation approaches will often find the correlated noise column to be significant simply by it riding on the statistical coattails of the true variable. The HRT avoids this issue by fully conditioning on the data.

One simple way of learning the conditional distribution of the design matrix is to assume a multivariate Gaussian distribution but simply estimating the precision matrix. However when the columns of the data are not Gaussian or not continuous then this learned distribution will prove a poor estimate of the conditional relationship of the data. The goal is this post is two-fold. First, show how to fit a marginal regression model to each column of the data (regularized Gaussian and Binomial regressions are used). Second, a python implementation will be used to complement the R code used previously. While this post will use an un-tuned random forest classifier, any machine learning model can be used for the training set of the data.

# import the necessary modules
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.datasets import make_classification
import seaborn as sns

Split a dataset into a tranining and a test folder

In the code blocks below we load a real and synthetic dataset to highlight the HRT at the bottom of the script. Note that the column types of each data need to be defined in the cn_type variable.

Option 1: South African Heart Dataset

link_data = ""
dat_sah = pd.read_csv(link_data)
# Extract the binary response and then drop
y_sah = dat_sah['chd']
# one-hot encode famhist
dat_sah['famhist'] = pd.get_dummies(dat_sah['famhist'])['Present']
# Convert the X matrix to a numpy array
X_sah = np.array(dat_sah)
cn_type_sah = np.where(dat_sah.columns=='famhist','binomial','gaussian')
# Do a train/test split
idx = np.arange(len(y_sah))
idx_test = np.where((idx % 5) == 0)[0]
idx_train = np.where((idx % 5) != 0)[0]
X_train_sah = X_sah[idx_train]
X_test_sah = X_sah[idx_test]
y_train_sah = y_sah[idx_train]
y_test_sah = y_sah[idx_test]

Option 2: Non-linear decision boundary dataset

# ---- Random circle data ---- #
n_circ = 1000
X_circ = np.random.randn(n_circ,5)
X_circ = X_circ + np.random.randn(n_circ,1)
y_circ = np.where(np.apply_along_axis(arr=X_circ[:,0:2],axis=1,func1d= lambda x: np.sqrt(np.sum(x**2)) ) > 1.2,1,0)

cn_type_circ = np.repeat('gaussian',X_circ.shape[1])

idx = np.arange(n_circ)
idx_test = np.where((idx % 5) == 0)[0]
idx_train = np.where((idx % 5) != 0)[0]

X_train_circ = X_circ[idx_train]
X_test_circ = X_circ[idx_test]
y_train_circ = y_circ[idx_train]
y_test_circ = y_circ[idx_test]



Function support

The code block below provides a wrapper to implement the HRT algorithm for a binary outcome using a single training and test split. See my previous post for generalizations of this method for cross-validation. The function also requires a cn_type argument to specify whether the column is continuous or Bernoulli. The glm_l2 function implements an L2-regularized generalized regression model for Gaussian and Binomial data using an iteratively re-weighted least squares method. This can generalized for elastic-net regularization as well as different generalized linear model classes. The dgp_fun function takes a model with with glm_l2 and will generate a new vector of the data conditional on the rest of the design matrix.


def hrt_bin_fun(X_train,y_train,X_test,y_test,cn_type):
    # Sigmoid function
    def sigmoid(x):
        return( 1/(1+np.exp(-x)) )
    # Sigmoid weightin
    def sigmoid_w(x):
        return( sigmoid(x)*(1-sigmoid(x)) )

    def glm_l2(resp,x,standardize,family='binomial',lam=0,add_int=True,tol=1e-4,max_iter=100):
        y = np.array(resp.copy())
        X = x.copy()
        n = X.shape[0]

        # Make sure all the response values are zeros or ones
        check1 = (~np.all(np.isin(np.array(resp),[0,1]))) & (family=='binomial')
        if check1:
            print('Error! Response variable is not all binary'); #return()
        # Make sure the family type is correct
        check2 = ~pd.Series(family).isin(['gaussian','binomial'])[0]
        if check2:
            print('Error! Family must be either gaussian or binoimal')

        # Normalize if requested
        if standardize:
            mu_X = X.mean(axis=0).reshape(1,X.shape[1])
            std_X = X.std(axis=0).reshape(1,X.shape[1])
            mu_X = np.repeat(0,p).reshape(1,X.shape[1])
            std_X = np.repeat(1,p).reshape(1,X.shape[1])

        X = (X - mu_X)/std_X

        # Add intercept
        if add_int:
            X = np.append(X,np.repeat(1,n).reshape(n,1),axis=1)

        # Calculate dimensions
        y = y.reshape(n,1)

        p = X.shape[1]
        # l2-regularization
        Lambda = n * np.diag(np.repeat(lam,p))

        bhat = np.repeat(0,X.shape[1])

        if family=='binomial':
            bb = np.log(np.mean(y)/(1-np.mean(y)))
            bb = np.mean(y)

        if add_int:
            bhat = np.append(bhat[1:p],bb).reshape(p,1)

        if family=='binomial':
            ii = 0
            diff = 1
            while( (ii < max_iter) & (diff > tol) ):
                ii += 1
                # Predicted probabilities
                eta =
                phat = sigmoid(eta)
                res = y - phat
                what = phat*(1-phat)
                # Adjusted response
                z = eta + res/what
                # Weighted-least squares
                bhat_new = np.linalg.inv( * what).T,X) + Lambda), * what).T, z) )
                diff = np.mean((bhat_new - bhat)**2)
                bhat = bhat_new.copy()
            sig2 = 0

            bhat = np.linalg.inv(,X) + Lambda ),, y) )
            # Calculate the standard error of the residuals
            res = y -,bhat)
            sig2 = np.sum(res**2) / (n - (p - add_int))

        # Separate the intercept
        if add_int:
            b0 = bhat[p-1][0]
            bhat2 = bhat[0:(p-1)].copy() / std_X.T # Extract intercept
            b0 = b0 - np.sum(bhat2 * mu_X.T)
            bhat2 = bhat.copy() / std_X.T
            b0 = 0

        # Create a dictionary to store the results
        ret_dict = {'b0':b0, 'bvec':bhat2, 'family':family, 'sig2':sig2, 'n':n}
        return ret_dict

    # mdl=mdl_lst[4].copy(); x = tmp_X.copy() 
    # Function to generate data from a fitted model
    def dgp_fun(mdl,x):
        tmp_n = mdl['n']
        tmp_family = mdl['family']
        tmp_sig2 = mdl['sig2']
        tmp_b0 = mdl['b0']
        tmp_bvec = mdl['bvec']
        # Fitted value
        fitted = np.squeeze(, tmp_bvec) + tmp_b0)

        if tmp_family=='gaussian':
            # Generate some noise
            noise = np.random.randn(tmp_n)*np.sqrt(tmp_sig2) + tmp_b0
            y_ret = fitted + noise
            y_ret = np.random.binomial(n=1,p=sigmoid(fitted),size=tmp_n)
        # Return

    # Logistic loss function
    def loss_binomial(y,yhat):
        ll = -1*np.mean(y*np.log(yhat) + (1-y)*np.log(1-yhat))
    # Loop through and fit a model to each column
    mdl_lst = []
    for cc in np.arange(len(cn_type)):
        tmp_y = X_test[:,cc]
        tmp_X = np.delete(X_test, cc, 1)
        tmp_family = cn_type[cc]
    # Fit random forest
    clf = RandomForestClassifier(n_estimators=100, max_depth=2, random_state=0), y_train)
    # Baseline predicted probabilities and logistic loss
    phat_baseline = clf.predict_proba(X_test)[:,1]    
    loss_baseline = loss_binomial(y=y_test,yhat=phat_baseline)

    pval_lst = []
    nsim = 250
    for cc in np.arange(len(cn_type)):
        print('Variable %i of %i' % (cc+1, len(cn_type)))
        mdl_cc = mdl_lst[cc]
        X_test_not_cc = np.delete(X_test, cc, 1)
        X_test_cc = X_test.copy()
        loss_lst = []
        for ii in range(nsim):
            xx_draw_test = dgp_fun(mdl=mdl_cc,x=X_test_not_cc)
            X_test_cc[:,cc] = xx_draw_test
            phat_ii = clf.predict_proba(X_test_cc)[:,1]
            loss_ii = loss_binomial(y=y_test,yhat=phat_ii)
        pval_cc = np.mean(np.array(loss_lst) <= loss_baseline)

    # Return p-values

Get the p-values for the different datasets

Now that the hrt_bin_fun has been defined, we can perform inference on the columns of the two datasets created above. The results below show that the sbp, tobacco, ldl, adiposity, and age are statistically significant features for the South African Heart Dataset. As expected, the first two variables, var1, and var2 from the non-linear decision boundary dataset are important as these are the two variables which define the decision boundary with the rest of the variables being noise variables.

pval_circ = hrt_bin_fun(X_train=X_train_circ,y_train=y_train_circ,X_test=X_test_circ,y_test=y_test_circ,cn_type=cn_type_circ)
pval_sah = hrt_bin_fun(X_train=X_train_sah,y_train=y_train_sah,X_test=X_test_sah,y_test=y_test_sah,cn_type=cn_type_sah)
pd.concat([pd.DataFrame({'vars':dat_sah.columns, 'pval':pval_sah, 'dataset':'SAH'}),
           pd.DataFrame({'vars':['var'+str(x) for x in np.arange(5)+1],'pval':pval_circ,'dataset':'NLP'})])
vars pval dataset
sbp 0.000 SAH
tobacco 0.000 SAH
ldl 0.000 SAH
adiposity 0.000 SAH
famhist 0.112 SAH
typea 1.000 SAH
obesity 1.000 SAH
alcohol 0.568 SAH
age 0.000 SAH
var1 0.000 NLP
var2 0.000 NLP
var3 0.412 NLP
var4 0.304 NLP
var5 0.192 NLP
Written on April 23, 2019