Several issues

[simon-lowe]

Hi,

First I wanted to thank you for the great package. It is super useful! I've been using it for two projects, so I have a bunch of questions/issues and I didn't want to open a separate issue for each, but if you prefer that I can also do that.

Small issues:

• On my personal computer, on windows, I have not been able to install pytwoway, even after a couple of hours of trying because of the dependency quadprog. This seems to be a known issue, and I tried the fixes, but havent been able to make it work. In particular, I think quadprog now made a version with a prebuilt wheel and is recommending to packages which use quadprog to put that as a dependency.
• Wanted to mention 2 warnings which seem to have no effect on the commands running:

1. I get the warning: "future version 'df.iloc[:,i] = newvals' will attempt to set the values in place instead of always setting in a new array". Just wanted to give a heads-up in case that will change something.
2. When I run a model with fit(), it works, but I do get the following warning: "sys:1: ResourceWarning: unclosed socket" (and then a bit more detail)

New feature maybe?:

• Is it currently possible to add a third fixed-effect (an origin firm fixed-effect) like in It Ain't Where You're From, It's Where You're At: Hiring Origins, Firm Heterogeneity, and Wages (Sabrina Di Addario, Patrick Kline, Raffaele Saggio, and Mikkel Sølvsten) Journal of Econometrics, 232 (April 2023), pp. 340-374.. In particular, does simply adding it as a control work (the connectedness might not be correct then?)?

(Apologies in what follows I don't have reproducible examples, because I'm using confidential admin data)
Bigger issues:

• If I use summary on the cleaned data, and in the next line run FEEstimator, the variances of y are different (see examples below). This is also the case with the simulated datasets in the examples. This only happens with collapse_at_connectedness_measure: True. Is that because different weights are used?
• Negative implied alphas or the terms don't add up to the variance. I am using a binary variable as a dependent variable (like in Employers and Unemployment Insurance Take-up (Lachowska, Sorkin and Woodbury)). Below are 3 different versions of the summary and estimation results:

1. with collapse_at_connectedness_measure: True
   
   
2. with collapse_at_connectedness_measure: False
   
   
3. with collapse_at_connectedness_measure: False and controls
   
   

Sorry about the very long post and I hope that I didn't do something stupid!

Thanks a lot for your help!

[adamoppenheimer]

Hi Simon,

Thanks for using the package and for reaching out!

Small issues:

• I will have to talk with Professor Lamadon about updating quadprog. The reason pytwoway requires a specific version is that he requires a Conda version of quadprog, and since they don't officially release that he needs to compile it himself (I can't since I don't have a Windows computer). If you are open to publishing an updated version of quadprog on Conda then I can update the dependency, but I will warn you that Professor Lamadon has told me it is a difficult process.
• Thank you for letting me know about the warnings. I am going to open an issue for the first one, the second one I've been getting for a very long time. I've looked it up online, I think it's very difficult to diagnose and it hasn't caused any problems for me yet.

New Feature:

• There is already an issue open for this, but thank you for the paper, I've been looking for a good resource on how to estimate this. So far I have written some sample code to estimate the connected set (not leave-one-out) that seems to work well in testing, but I've only tried it on simulated data. The code is based on the papers by Weeks & Williams (1964) and Godolphin & Godolphin (2001). This code is a simple implementation and depends on sparsity, so if you are including a categorical control with a small number of categories, the code will not work (since it constructs the XX' matrix (not X'X), so if there are only a small number of categories many of the matrix elements will become non-zero). When I have time I'll look at the paper you linked, and also try to optimize the code I already wrote. I think this is an important feature to have.
• If you try out my code and have any feedback, I would greatly appreciate it!
• While running the code I pasted can be helpful, in my experience the FE estimator will raise an error when estimating leverages if one of the leverages equals one (which happens when the set isn't leave-one-out connected). So if you are using a categorical control and the HE correction still works, it is very likely that the set is leave-one-out connected.

Bigger Issues:

• Summary gives an unweighted variance whereas the FEEstimator gives the weighted variance. When you collapse, the spells receive a weight equal to their length so that's why you are getting different values.
• Another user had an issue with negative estimates, but with the CRE estimator. The correspondence was over email, so I'll paste my response:

---

Professor Lamadon said you could bound the CRE variance below at 0, but the fact that it's negative points to other issues. I personally don't think bounding it below seems like the right decision, since it is unlikely the true variance is actually 0.

He said what you should do is to run attrition plots to investigate the estimators' behaviors. I just updated the attrition notebook (also available on the github under docs/source/notebooks), so you should be able to just fill it in with the data you are investigating to generate the attrition plots (I also updated the Monte Carlo, FE, and CRE notebooks, so you may want to check out the new versions of those on the documentation site as well).

---

• The values don't sum to the variance since they are approximations, if you want the estimates to be better then you should increase the number of trace and leverage draws.

---

Please let me know if this helps address your questions, and please let me know whatever new questions you may have!

Best,
Adam

---

Here's the code to compute the connected components:

'''
Connected components for arbitrarily many categorical control variables.
'''
import numpy as np
import pandas as pd
import bipartitepandas as bpd
from bipartitepandas.util import to_list
from scipy.sparse import csc_matrix, hstack
from scipy.sparse.linalg import svds, eigsh
from igraph import Graph

def SPxD(sp, diag):
    '''
    Faster product of a sparse matrix and a diagonal, i.e. take sp @ diag.

    Arguments:
        sp (CSC Matrix): sparse matrix
        diag (NumPy Array or float): diagonal entries or multiplicative factor

    Returns:
        (CSC Matrix): product of sparse matrix and diagonal
    '''
    if isinstance(diag, (float, int)):
        # If multiplicative factor
        if diag == 1:
            return sp
        return diag * sp

    # If diagonal entries
    return sp.multiply(diag)

def n_connected_components(bdf, categorical_controls, method='approximate', n_categories=None, null_space_dim=1, warn=2, rng=None):
    '''
    Compute the n-way connected components of a dataframe.

    Arguments:
        bdf (BipartiteDataFrame): (collapsed) long format data
        categorical_controls (list): categorical controls
        method (str): if 'approximate', use the method of Weeks & Williams (1964) - this guarantees connectedness, but does not guarantee that different subsets are disconnected; if 'exact', use the method of Godolphin & Godolphin (2001) - this guarantees both connectedness and that different subsets are disconnected
        n_categories (dict of ints or None): link each categorical variable to the number of categories it has, for computational efficiency, using 'i' for worker ids and 'j' for firm ids; None is equivalent to {}
        null_space_dim (int): used when `method='exact'`; dimension of the null space, it is recommended to slowly increment this (by increasing the value) until the singular value warning disappears; then once the function outputs the correct number of 0 singular values, set `null_space_dim` to the correct value and set `warn=0`
        warn (int): used when `method='exact'`; if 2, raise an error if all the singular values are 0; if 1 or 2, print how many singular values are 0; if 0, do nothing
        rng (np.random.Generator or None): used when `method='exact'`; NumPy random number generator; None is equivalent to np.random.default_rng(None)

    Returns:
        (list of lists): each entry of the primary list gives a list of the observations included in a particular connected component
    '''
    if rng is None:
        rng = np.random.default_rng(None)

    categorical_controls = to_list(categorical_controls)

    ## Set n_categories ##
    if n_categories is None:
        n_categories = {}
    # Make sure not to overwrite original
    n_categories = n_categories.copy()
    if 'i' not in n_categories.keys():
        n_categories['i'] = bdf.n_workers()
    if 'j' not in n_categories.keys():
        n_categories['j'] = bdf.n_firms()
    for cat_control in categorical_controls:
        if cat_control not in n_categories.keys():
            n_categories[cat_control] = bdf.n_unique_ids(cat_control)

    ## Indices for each category ##
    cat_indices = {}
    cum_cols = 0

    ## General variables ##
    n_obs = len(bdf)
    ones = np.ones(n_obs, dtype=int)
    idx = np.arange(n_obs)

    ## Create sparse matrices ##
    X1 = csc_matrix((ones, (idx, bdf.loc[:, 'j'].to_numpy())), shape=(n_obs, n_categories['j']))
    X2 = csc_matrix((ones, (idx, bdf.loc[:, 'i'].to_numpy())), shape=(n_obs, n_categories['i']))

    # Update indices
    cat_indices['j'] = (cum_cols, cum_cols + n_categories['j'])
    cum_cols += n_categories['j']

    for cat_control in categorical_controls:
        # Control columns
        Xc = csc_matrix((ones, (idx, bdf.loc[:, cat_control].to_numpy())), shape=(n_obs, n_categories[cat_control]))
        X1 = hstack([X1, Xc])

        # Update indices
        cat_indices[cat_control] = (cum_cols, cum_cols + n_categories[cat_control])
        cum_cols += n_categories[cat_control]
    del Xc

    if method == 'approximate':
        ### Weeks & Williams (1964) ###
        X = hstack([X2, X1])
        del X1, X2

        # How many groups in common are required to consider observations adjacent
        adjacency_thres = len(n_categories.keys()) - 1

        ## Adjacency matrix ##
        adjacency_matrix = ((X @ X.T) == adjacency_thres).astype(int, copy=False)

        ## Construct graph ##
        G = Graph.Adjacency(adjacency_matrix)

        ## Find connected components ##
        cc_list = sorted(G.components(mode='weak'), reverse=True, key=len)

        return cc_list

    ### Godolphin & Godolphin (2001) ###
    bdf = bdf.copy(deep=True)
    ## Information Matrix ##
    N = X1.T @ X2
    C = X1.T @ X1 - SPxD(N, 1 / np.bincount(bdf.loc[:, 'i'].to_numpy())) @ N.T
    del X1, X2, N

    ## Singular Values ##
    # svecs, svals, _ = svds(C, k=null_space_dim, which='SM', return_singular_vectors='u', random_state=rng)
    evals, evecs = eigsh(C, k=null_space_dim, sigma=1e-10, which='LM')
    svals, svecs = evals, evecs

    n_zero_svals = np.sum(np.isclose(svals, 0))
    if n_zero_svals == len(n_categories.keys()) - 1:
        # Matrix is full rank
        return [np.arange(n_obs)]
    if (warn == 2) and n_zero_svals == null_space_dim:
        raise ValueError('All singular values are 0.')
    elif warn in [1, 2]:
        print(f'{n_zero_svals} singular values are 0.')

    ## Null Space ##
    kernel_basis = svecs[:, np.nonzero(np.isclose(svals, 0))[0]]

    # # Normalize kernel basis
    # kernel_basis -= kernel_basis.min(axis=0)
    # kernel_basis = MxD(kernel_basis, 1 / kernel_basis.max(axis=0))

    # Round (don't want to consider basis vectors different if they're extremely close)
    kernel_basis = np.round(kernel_basis, 6)

    ## Connectedness ##
    fully_identified = []
    for cat_control in ['j'] + categorical_controls:
        # Connected sets for each control variable are determined by unique rows of the kernel basis
        # Source for sparse unique rows: https://stackoverflow.com/a/52891452/17333120
        control_idx = cat_indices[cat_control]
        kernel_subset = kernel_basis[control_idx[0]: control_idx[1], :]
        unique_basis_vectors, control_concomps = np.unique(kernel_subset, axis=0, return_inverse=True)
        if (len(unique_basis_vectors) == 1) and (cat_control != 'j'):
            # If entire category is identified (ignore j)
            fully_identified.append(cat_control)
        else:
            # Create new category based on connected components
            bdf.loc[:, cat_control] = control_concomps[bdf.loc[:, cat_control].to_numpy()]

            # Update n_categories
            n_categories[cat_control] = len(unique_basis_vectors)
    if len(fully_identified) > 0:
        ## If some columns are fully identified, ignore them ##
        # Make sure not to overwrite original
        categorical_controls = categorical_controls.copy()
        for col in fully_identified:
            # Delete fully identified columns
            del categorical_controls[categorical_controls.index(col)]
            del n_categories[col]

    if len(categorical_controls) == 0:
        print('All components can be found by directly computing largest set of connected workers + firms.')
        return None

    return n_connected_components(bdf=bdf, categorical_controls=categorical_controls, method='approximate', n_categories=n_categories)

##############
## Simulate ##
##############
rng = np.random.default_rng(1234)
sim_params = bpd.sim_params(
    {
        'n_workers': 500,
        'firm_size': 10,
        'p_move': 0.03
    }
)
clean_params = bpd.clean_params(
    {
        'connectedness': 'connected',
        'copy': True
    }
)
sim_data = bpd.SimBipartite(sim_params).simulate(rng)

n_control1 = 15
n_control2 = 8
sim_data['worker_control'] = rng.integers(n_control1, size=len(sim_data['i'].unique()))[sim_data['i']]
sim_data['firm_control'] = rng.integers(n_control2, size=len(sim_data['j'].unique()))[sim_data['j']]

bdf = bpd.BipartiteDataFrame(
    sim_data,
    custom_categorical_dict={'worker_control': True, 'firm_control': True},
    custom_dtype_dict={'worker_control': 'categorical', 'firm_control': 'categorical'},
    custom_how_collapse_dict={'worker_control': 'first', 'firm_control': 'first'}
).clean()

categorical_controls = ['worker_control', 'firm_control']
n_categories = {'i': bdf.n_workers(), 'j': bdf.n_firms(), 'worker_control': n_control1, 'firm_control': n_control2}

cc_list1 = n_connected_components(bdf, categorical_controls, method='approximate', n_categories=n_categories, null_space_dim=1, warn=2, rng=None)

cc_list2 = n_connected_components(bdf, categorical_controls, method='exact', n_categories=n_categories, null_space_dim=30, warn=2, rng=None)

###################
## Connected set ##
###################
bdf2 = bdf.loc[cc_list2[0], :].clean()

categorical_controls2 = ['worker_control', 'firm_control']
n_categories2 = {'i': bdf2.n_workers(), 'j': bdf2.n_firms(), 'worker_control': bdf2.n_unique_ids('worker_control'), 'firm_control': bdf2.n_unique_ids('firm_control')}

cc_list2_1 = n_connected_components(bdf2, categorical_controls2, method='approximate', n_categories=n_categories2, null_space_dim=1, warn=2, rng=None)

cc_list2_2 = n_connected_components(bdf2, categorical_controls2, method='exact', n_categories=n_categories2, null_space_dim=7, warn=1, rng=None)

[simon-lowe]

Hi Adam,

Thanks a lot for the fast answer! I will take a look at the things you suggested and get back to you.

Cheers,
Simon

[adamoppenheimer]

Hi Simon,

In case you want to check out the papers I used for the connectedness, here are the links (I had to use another paper from Godolphin to find their algorithm, since the 2001 paper isn't available online):
Weeks & Williams (1964)
Godolphin (2004)

If you choose to read them and have any insights please let me know! I was a bit lost in the Godolphin paper since they use terminology I had never seen before, so I may be misinterpreting what their algorithm does. But like I said before, when I experimented with some sample data it does seem to find connected components that are larger than those found by the Weeks & Williams algorithm.

Best,
Adam

Edit: I also took a look at the DKSS paper. Their algorithm is very clever, but from what I can tell what it is doing is computing the largest connected component that is connected to the categories they normalize. I don't think this is a sensible approach in general, since it is largely based on chance whether you happen to normalize the correct categories. Because of this, unless someone informs me that my intuition on this is incorrect, I don't think I'm going to implement their algorithm. However, one part of their algorithm is useful for checking that the data is full rank, and I plan to implement it into pytwoway - I was looking for something like this! Thanks for sharing the paper!

[adamoppenheimer]

Hi Simon,

I was thinking about this a bit and for the negative alphas, you could always add the direct computation of var(alpha) by adding tw.Q.VarAlpha() to your Q matrix for FE without controls, maybe if you compute it this way you'll get something positive.

Also, I happened to stumble upon this paper today, if you check out the bottom of page 7/top of page 8, they discuss negative variance estimates and an alternative consistent estimator where the variances will always be non-negative (but you lose unbiasedness, and you'd need to think about how to move the variance from alpha into the other terms).

Best,
Adam

[adamoppenheimer]

Hi Simon,

Just following up on this.

Best,
Adam