Pooling levels of categorical variables for regression trees

I have a data set I would like to do a regression analysis for. There are many features of both categorical and continuous types. One of the categorical features has many (>75) levels so this is an issue. I have reason to believe that some of the levels are essentially the same. I intend to use decision trees with some ensemble method (ie Bagging or boosting).

I would like to try to pool/cluster the levels of the problematic feature to improve performance. I realize that theoretically if the ensemble/number of leaves is large enough this is not necessary but I am already having computational issues.

Is there a standard method to combine levels which perform the same?

-------------Edit--------------

I think I found a method which would work. The idea would be to use use the proximity matrix. This is essentially the N_Obs by N_Obs matrix for the fraction of out of bag trees where the observations where in the same terminal node. We can then aggregate this into a level by level matrix where the elements are the average of the fraction in the proximity matrix. We would then pool all levels together when they past a threshold and see if this improves RMSE. It is likely best to take a step-wise iterative approach to find the optimal pooling but I might just take the threshold as the average of the diagonal. This should give a reasonable threshold because it represents how often each level is in the same terminal node as itself. Comments welcome, I will report back on results.

• Can you pool based on domain knowlege? This is a better first step than a supervised machine learning approach, as it cannot overfit or leak. – Matthew Drury Mar 30 '17 at 0:07
• We could guess at it a bit but it would be hard. – Keith Mar 30 '17 at 16:23
• I'm not advocating guessing, but using your understanding of the business or scientific problem to group factors conceptually is a fine thing to do, and is safer than relying on your data to do the work. – Matthew Drury Mar 30 '17 at 18:26
• Have a look at stats.stackexchange.com/questions/227125/… and links therein. Maybe some similar ideas can be used with trees? – kjetil b halvorsen May 17 '17 at 11:02
• Possible duplicate of Principled way of collapsing categorical variables with many categories – shadowtalker Aug 22 '17 at 15:31

I have implemented my solution to this. I wrote two functions:

prox_matrix(df, target, features, cluster_dimension,trees = 10)

Parameters

• df: Input dataframe
• target: Dependant variable you are trying to predict with the random forrest
• features: List of independent variables
• cluster_dimension: Dimension you would like to cluster/pool to add to your list of features
• trees: the number of trees to use in your Random Forest

Returns

• D: DataFrame of the proximity matrix for the cluster_dimension

Code Below

def prox_matrix(df, target, features, cluster_dimension,trees = 10):
#https://www.stat.berkeley.edu/~breiman/RandomForests/cc_home.htm#prox

from sklearn.ensemble import RandomForestRegressor
import numpy as np
import pandas as pd

#initialize datframe for independant variables
independant = pd.DataFrame()

#Handle Categoricals: This should really be added to RandomForestRegressor
for column,data_type in df[features].dtypes.iteritems():
try:
independant[column] = pd.to_numeric(df[column],downcast = 'integer')
except ValueError:
contains_nulls = df[column].isnull().values.any()
dummies = pd.get_dummies(df[column],prefix=column,dummy_na=contains_nulls,drop_first=True)
independant[dummies.columns] = dummies

if len(independant.index) != len(df.index):
raise Exception('independant variables not stored properly')

#train Model
clf = RandomForestRegressor(n_estimators=trees, n_jobs=-1)
clf.fit(independant, df[target])

#Final leaf for each tree
leaves = clf.apply(independant)
#value in cluster dimension
labels = df[cluster_dimension].values

numerator_matrix = {}
for i,value_i in enumerate(labels):
for j,value_j in enumerate(labels):
if i >= j:
numerator_matrix[(value_i,value_j)] = numerator_matrix.get((value_i,value_j), 0) + np.count_nonzero(leaves[i]==leaves[j])
numerator_matrix[(value_j,value_i)] = numerator_matrix[(value_i,value_j)]

#normalize by the total number of possible matchnig leaves
prox_matrix = {key: 1.0 - float(x)/(trees*np.count_nonzero(labels==key[0])*np.count_nonzero(labels==key[1])) for key, x in numerator_matrix.iteritems()}

#make sorted dataframe
levels = np.unique(labels)
D = pd.DataFrame(data=[[ prox_matrix[(i,j)] for i in levels] for j in levels],index=levels,columns=levels)

return D


kMedoids(D, k, tmax=100)

Parameters

• D: Proximity/distance matrix
• k: Number of clusters
• tmax: Maximum number of iterations to check for convergence of clustering

Returns

• M: List of mediods
• C: Dictionary mapping the clustered levels to each mediod
• S: Silhouette of each cluster for evaluation of performance

Code Below

def kMedoids(D, k, tmax=100):
#https://www.researchgate.net/publication/272351873_NumPy_SciPy_Recipes_for_Data_Science_k-Medoids_Clustering
import numpy as np
import pandas as pd

# determine dimensions of distance matrix D
m, n = D.shape

if m != n:
raise Exception('matrix not symmetric')

if sum(D.columns.values != D.index.values):
raise Exception('rows and columns do not match')

if k > n:
raise Exception('too many medoids')

#Some distance matricies will not have a 0 diagonal
Dtemp =D.copy()
np.fill_diagonal(Dtemp.values,0)

# randomly initialize an array of k medoid indices
M = list(Dtemp.sample(k).index.values)

# initialize a dictionary to represent clusters
Cnew = {}

for t in xrange(tmax):
# determine mapping to clusters
J = Dtemp.loc[M].idxmin(axis='index')
#Fill dictionary with cluster members
C = {kappa: J[J==kappa].index.values for kappa in J.unique()}
# update cluster medoids
Cnew = {Dtemp.loc[C[kappa],C[kappa]].mean().idxmin() : C[kappa] for kappa in C.keys()}
#Update mediod list
M = Cnew.keys()

# check for convergence (ie same clusters)
if set(C.keys()) == set(Cnew.keys()):
if not sum(set(C[kappa]) != set(Cnew[kappa]) for kappa in C.keys()): break
else:
print('did not converge')

#Calculate silhouette
S = {}
for kappa_same in Cnew.keys():
a = Dtemp.loc[Cnew[kappa_same],Cnew[kappa_same]].mean().mean()
b = np.min([Dtemp.loc[Cnew[kappa_other],Cnew[kappa_same]].mean().mean() for kappa_other in Cnew.keys() if kappa_other!=kappa_same])
S[kappa_same] = (b - a) / max(a, b)

# return results
return M, Cnew, S


Notes:

1. There are links to theory documentation in the code
2. I used all records not strictly the OOB records. Follow up here
3. The prox_matrix() method is very slow. I have done a few things to speed it up but most of the cost comes from the double loop. Updates welcome.
4. The diagonal of the proximity matrix need not be zeros. I force this in the KMedoids method so that I get convergence.
• would you have a vignette or blogpost that demonstrates the application and performance of this code? – Rahul Jul 18 '17 at 2:09
• Sorry. I don't blog. It is pretty straight forward. 1) Make a proximity matrix based on the feature you want to cluster 2) Cluster using K-medoids . All the concepts are relatively well documented defined elsewhere. It is just the code that was needed. – Keith Jul 18 '17 at 4:45