Sample subset characteristics

Question

Edit: It seems that I have made the question too general, so I will provide a specific example of the type of problem I am trying to solve.

I have a database that contains every item that is sold at a grocery store, along with each defining feature of the items (i.e, price, country of origin, food category, producer...). I also have a database with customer purchases, so for each customer it lists n items that were bought.

For each customer I want to understand "as best as I can" the underlying reasoning for why they chose that group of n items in a quantitative manner.

A core caveat is that this is not being asked from an academic or theoretical viewpoint. This is purely practical

Original question:

When drawing a random sample from data, it is typically tested to check if the sample is properly representative of the total population.

Assume a scenario where a subset exists within a population, you know that it was not selected at random and that the individual points where chosen due to some sort of criteria.

If the sample was not chosen at random then it must have some distinguishable features and bias when compared to the population.

Are there any specific quantitative methods used for decomposing the differences in between a subset and the population, outside of just plotting distributions, one vs the other.

Also are there any Python packages or tools for this?

In plain terms:

I have a basket with a thousand items and I know the features/characteristics of each item. Someone comes and picks 10 items based on some preferences/characteristics/bias. I now want to understand the underlying reasoning for why they chose that group of 10 items in quantitative manner"

Answer 1

You can do clustering and then select the subsets so you are sure that your subset has similar characteristics of main dataset and other subset.

For the purpose of train-test split, I usually split main data into different clusters, and then split each cluster to 80-20 for training-test sets using sklearn train_test_split(... stratify=y_clus).

You can use my code; however, it's not always returning the best results and I may need to check different random_state values to find the best model.

In the first step, you need to encode your categorical variables and scale the numerical ones.

from sklearn import decomposition, datasets, model_selection, preprocessing, metrics
from sklearn.preprocessing import StandardScaler, OneHotEncoder, MinMaxScaler, LabelEncoder
from sklearn.pipeline import Pipeline
from sklearn.compose import ColumnTransformer

categorical_features = ['gender', 'marital','province','agegroup','isdirector']

categorical_transformer = Pipeline(steps=[
    ('onehot', OneHotEncoder(handle_unknown='ignore'))])


numeric_features = [col for col in df2.columns[1:-1] if col not in categorical_features]
#numeric_features=[el for el in numeric_features if el!='age']

numeric_transformer = Pipeline(steps=
    ('scaler', StandardScaler()) 
    ])



preprocessor = ColumnTransformer(
    transformers=[
        ('num', numeric_transformer, numeric_features),
        ('cat', categorical_transformer, categorical_features)])


y_encoder = LabelEncoder()
y = y_encoder.fit_transform(df2['sales'])

X = df2[numeric_features + categorical_features]

and the second step is to call the dataset_builder().

_, y_train, _, y_test, _, y_val, X_train_sc, X_test_sc, X_val_sc = dataset_builder(X,y, do_clustering=True, 
              singleclass=singcls,dataset_type='TVT', random_state=rnd_data)

The skipped variables ( _ ) are X_train, X_test, X_val for the unscaled (original) X.

BUT HOW IT WORKS????

The code use following function to do the clustering. I modified the code found on SciPy Hierarchical Clustering and Dendrogram Tutorial

# hierarchical/agglomerative
from scipy.cluster.hierarchy import dendrogram, linkage, fcluster
import numpy as np
import warnings

def classclustering(X_sc,y=None, Z=None, nclusters=0, method='ward', metric='euclidean', maxdepth_show = 20,show_charts=True):
    """
    Z: linkage matrix
    method: The linkage algorithm to use. Please check <scipy.cluster.hierarchy.linkage>
            single, complete, weighted,centroid, median, ward

            Methods ‘centroid’, ‘median’ and ‘ward’ are correctly defined only if Euclidean pairwise metric is used.

    metric: Pairwise distances between observations in n-dimensional space. Please check  <scipy.spatial.distance.pdist>
            euclidean, minkowski, cityblock, seuclidean (standardized Euclidean), cosine, correlation, 
            hamming, jaccard, chebyshev, canberra, braycurtis, mahalanobis, yule, matching, dice, kulsinski, 
            rogerstanimoto, russellrao, sokalmichener, sokalsneath, wminkowski
    """
    def performclustering(X_sc, Z=None, nclusters=0, method='ward', metric='euclidean', maxdepth_show = 20):
        linked=Z
        if linked is not None: # use previous linkage for custom number of clusters.
            if nclusters<2:
                raise Exception("nclus must be greater than 1 when linkage matrix (Z) has been used!")

            clus=fcluster(linked, nclusters, criterion='maxclust')       

        else:

            # faster calculation by showing only the first 20 clusters, p=20

            linked = linkage(X_sc, method, metric)

            labelList = range(1, 11)

            if show_charts:
                plt.figure(figsize=(10, 7))
                dendrogram(linked,
                            orientation='top',
                            #labels=labelList,
                            distance_sort='descending',
                            truncate_mode='lastp',  # show only the last p merged clusters
                            p=maxdepth_show,        # show only the last p merged clusters
                            show_leaf_counts=True,  # otherwise numbers in brackets are counts
                            leaf_rotation=90.,
                            leaf_font_size=12.,
                            show_contracted=True    # to get a distribution impression in truncated branches         
                          )

                plt.show()


            # Elbow Method
            # calculating the best number of clusters. It's 4 or 6 for only numberical data, and 3 or 9 for all data

            last = linked[-20:, 2]
            last_rev = last[::-1]
            idxs = np.arange(1, len(last) + 1)

            acceleration = np.diff(last, 2)  # 2nd derivative of the distances
            acceleration_rev = acceleration[::-1]
            k = acceleration_rev.argmax() + 2  # if idx 0 is the max of this we want 2 clusters

            if show_charts:
                plt.plot(idxs, last_rev)
                plt.xticks(np.arange(min(idxs), max(idxs)+1, 2.0))
                plt.xlabel("Number of clusters")
                plt.plot(idxs[:-2] + 1, acceleration_rev)

                plt.show()


            if nclusters>0:
                print("\033[1;31;47m Warning....\n    ncluster has been set. Optimal number of clusters (%s) has been disabled!\n"%k+'\033[0m')
            else:
                nclusters=k

            if show_charts:
                print ("clusters:", nclusters)
            clus=fcluster(linked, nclusters, criterion='maxclust')

        return clus, linked, nclusters


    if y is None:  # single-class clustering
        if type(Z)==list:
            raise Exception("Multi-class clustering is not working with predefined Linkage Matrix (Z)!")
        else:
            clus,linked, nclus = performclustering(X_sc, Z, nclusters, method, metric, maxdepth_show)

    else:   # perform multi-class clustering
        if Z is not None:
            raise Exception("Multi-class clustering is not working with predefined Linkage Matrix (Z)!")
        else:
            y_classes = set(y)

            #clus_y=[]
            linked=[]
            if show_charts:
                print("===========================")

            clus= np.zeros(X_sc.shape[0],dtype=int)
            tmpclus_old=[0]
            nclus=0
            for cl in y_classes:

                if show_charts:
                    print("Cluster analysis for class: %s"%cl)

                mask = y==cl   # indices
                tmpclus, tmplinked, tmp_nclus = performclustering(X_sc[mask,:], Z, nclusters, method, metric, maxdepth_show)
                nclus += tmp_nclus
                #clus_y.append(tmpclus)
                linked.append(tmplinked)
                clus[mask]=tmpclus+max(tmpclus_old)
                tmpclus_old = tmpclus

                if show_charts:
                    print("===========================")

    return clus,linked, nclus

To use the function, you just need to feed it with scaled data if you have categorical variables. The function can do clustering based on X only, or doing clustering for each calsses in y (clustering for YES, NO, ... separately).

scaler = preprocessor.fit(X)
X_sc = scaler.transform(X)

# single-class clustering
clus,Z,nclus= classclustering(X_sc,show_charts=True)

# multi-class clustering
#clus,Z, nclus = classclustering(X_sc, y, show_charts=True)

The output would be something like this:

and number of clusters is the peak in orange line:

Now, if you are going to split your data into training-test (dataset_type='TT') or training-validation-test sets (dataset_type='TVT'), use following function:

import imblearn.over_sampling as OverSampler

X_labels = ''
categorical_features_onehot = ''


def dataset_builder(X,y, do_clustering=True, singleclass=True, dataset_type='TVT', random_state=2):

    X_train, X_val, X_test, y_train, y_val, y_test = [],[],[],[],[],[]
    dataset_type=dataset_type.lower()

    if dataset_type not in ['tt','tvt']:
        raise Exception("Unknown dataset_type!")

    if not do_clustering:
        if dataset_type=='tt':
            X_train, y_train, X_test, y_test, X_val,y_val = train_test_builder(X, y, validation_size=0, test_size=0.2, 
                                                                               random_state=random_state)
        else:
            X_train, y_train, X_test, y_test, X_val,y_val = train_test_builder(X, y, validation_size=0.15, test_size=0.15, 
                                                                               random_state=random_state)

    else:

        scaler = preprocessor.fit(X)
        X_sc = scaler.transform(X)

        if singleclass:
            # single-class clustering
            clus,Z,nclus= classclustering(X_sc,show_charts=False)
        else:
            # multi-class clustering
            clus,Z, nclus = classclustering(X_sc, y, show_charts=False)


        if dataset_type=='tt':
            for cl in set(clus):
                mask = clus==cl
                X_clus = X[mask]
                y_clus = y[mask]

                X_train_clus, y_train_clus, X_test_clus, y_test_clus, _, _ = train_test_builder(X_clus, y_clus, 
                                                                                    validation_size=0, test_size=0.2, 
                                                                                    random_state=random_state)

                X_train.append(X_train_clus)
                X_test.append(X_test_clus)
                y_train.append(y_train_clus)
                y_test.append(y_test_clus)


            # method 1.2, fastest
            X_train = np.concatenate(X_train,axis=0)
            X_test = np.concatenate(X_test,axis=0)
            y_train = np.concatenate(y_train,axis=0)
            y_test = np.concatenate(y_test,axis=0)

            # convert to dataframe
            X_train = pd.DataFrame(X_train,columns=X.columns)
            X_test = pd.DataFrame(X_test,columns=X.columns)
        else:
            for cl in set(clus):
                mask = clus==cl
                X_clus = X[mask]
                y_clus = y[mask]

                X_train_clus, y_train_clus, X_test_clus, y_test_clus, X_val_clus, y_val_clus = train_test_builder(X_clus, y_clus, 
                                                                      validation_size=0.15, test_size=0.15, random_state=random_state)


                X_train.append(X_train_clus)
                X_val.append(X_val_clus)
                X_test.append(X_test_clus)
                y_train.append(y_train_clus)
                y_val.append(y_val_clus)
                y_test.append(y_test_clus)


            global xt,xv,xtt
            xt,xv,xtt = X_train,X_val,X_test
            # method 1.2, fastest
            X_train = np.concatenate(X_train,axis=0)
            X_val = np.concatenate(X_val,axis=0)
            X_test = np.concatenate(X_test,axis=0)
            y_train = np.concatenate(y_train,axis=0)
            y_val = np.concatenate(y_val,axis=0)
            y_test = np.concatenate(y_test,axis=0)


            # convert to dataframe
            X_train = pd.DataFrame(X_train,columns=X.columns)
            X_val = pd.DataFrame(X_val,columns=X.columns)
            X_test = pd.DataFrame(X_test,columns=X.columns)

    # preprocessing based on X_train:
    scaler = preprocessor.fit(X_train)

    X_train_sc, X_test_sc, X_val_sc = [],[],[]

    X_train_sc = scaler.transform(X_train)
    X_test_sc = scaler.transform(X_test)
    if len(X_val)>0:
        X_val_sc = scaler.transform(X_val)

    # dummy categorical vars name created by preprocessor
    ohe=scaler.named_transformers_['cat']
    ohe=ohe.named_steps['onehot']
    global categorical_features_onehot
    categorical_features_onehot = ohe.get_feature_names(categorical_features)

    global X_labels
    X_labels = numeric_features+list(categorical_features_onehot)

    return X_train, y_train, X_test, y_test, X_val, y_val, X_train_sc, X_test_sc, X_val_sc

My code uses some global variales such as preprocessor, categorical_features_onehot (the label of dummy variables)

Answer 2

In the generality stated, your problem is completely hopeless. You need to make strong assumptions or have good prior knowledge both on the underlying population distribution and the possible selection algorithm. However, even under these restrictive conditions deciding about 10 out of 1000 will often be hard or impossible.

Hopeless in general

To see why this cannot work in general, just do a thought experiment. If someone presents you a pick of ten items out of 1000, they may have been picked uniform at random or picked according to the very deterministic and biased algorithm "always pick exactly those ten". Since the resulting sample is the same in both cases, you have no way to distinguish between the two algorithms.

Using prior knowledge

The thing you are looking for are tests for randomness or (more likely) tests for specific distributions. Look up "randomness tests" or "normality tests". These tests work typically by defining a test statistic, which is simply a specific function of your sample. Under strong assumptions on the features of the underlying population (such as normality) you can derive a distribution of the test statistic. Using this distribution you then derive an interval such that, for example, 90% of all possible subsets of 10 are within this interval. If the sample you want to test is outside this interval, you are "90% sure" it was not picked randomly according to the distribution. This is still very hard in general, since you need to come up with a test statistic which puts the "bad" samples at the ends of its distribution. It obviously helps a lot to come up with a good test statistic and the interval you need, if you have a prior idea about the bias and its strength as well.

Other ideas

As described above you have no chance with just a single trial (i.e. picking 10 items once). But maybe you can convince your adversary to do repeat experiments. An ad-hoc approach which may work, is to provide repeated subsamples to choose from. Give your adversary a randomly selected subsample of size 500, say. Let him choose 10 items and record the items. Repeat this a large number of times (maybe up to a few million times) and keep a tally of how often each item was picked. If the small samples are chosen uniform at random, every item of your large sample will have been picked on average the same number of times. If there is preference for certain items, these will be picked more often. With this approach you can easily reveal the "always pick those 10 items" algorithm. With other selection strategies your mileage may vary.

Answer 3

1st: The hope to understand the "underlying reasoning" of customers from just counting co-occurrences in baskets and completely disregarding the customers themselves is futile. There is a whole field of market research that carry out experiments, where intended customers are spending quite some time participating. These experiments are done for a reason.

2nd: random sampling is not the way to go, as groups are likely too small in a grocery store. Noise will be huge.

3rd: I'm pretty sure that you are missing important data, when referring to item features. Likely you do not have data on location or co-localization/distance etc. of the items in the store, do you?

Nonetheless there are ways to quantify relations between item purchases.

There is a field called "market basket analysis" in statistics, which does exactly that (see e.g. https://en.wikipedia.org/wiki/Affinity_analysis).

The procedure you would like to employ is called "association rule learning" (https://en.wikipedia.org/wiki/Association_rule_learning)

The mainly interesting terms for you will likely be "support", "confidence" and "lift":

"Support" quantifies the frequency of purchasing a set of particular items together divided by the number of all transactions

"Confidence" is the probability of purchasing item 2 conditional on having item 1 in the basket

"Lift" denotes the ratio of support to expected proportion of co-occurrence under the assumption that no deviation from random distribution exists.

As for practical advice on python packages: I rarely use python, but maybe the MLxtend package can by of use here.

And https://pbpython.com/market-basket-analysis.html could be a starting point to get into this topic.