variance decomposition issues across several models

[pierpaolocreanza]

Hi there,

Thank you for putting together this package! I'm trying to estimate several two-way models for my application and I get puzzling results. Here's a summary by case.

1. Standard AKM, without collapsing at the spell level and using leave-out-observation methods for HE.

HE-corrected estimates do not sum to one. The more suspicious change relative to HO-correction is var(eps), but is not obvious that that's where the problem is.

These is my pytwoway setup for these results:

# Prepare input DataFrame for PyTwoWay
    input_df = pd.DataFrame({
        'i': panel['cluster_id'].values,
        'j': panel['hgroup_id_temp'].values,
        'y': panel['outcome'].values,
        't': panel['year'].values
    })
    
    # Set parameters
    clean_params = bpd.clean_params({
        'connectedness': 'leave_out_observation',
        'drop_single_stayers': True,
        'collapse_at_connectedness_measure': False,
        'copy': False
    })

    bdf = bpd.BipartiteDataFrame(input_df, track_id_changes=True)
    
    # Clean the data
    bdf = bdf.clean(clean_params)

    # Estimation params
    fe_params = tw.fe_params({
        'ncore': 12,
        'ndraw_trace_sigma_2': 200,
        'ndraw_trace_ho': 200,
        'he': True,
        'ndraw_trace_he': 200,
        'ndraw_lev_he': 1800,
        'attach_fe_estimates': True,
        'preconditioner': 'ichol'
    })

    # Run model and save results
    
    # Initialize and fit the FE estimator with controls
    fe_params['Q_var'] = [tw.Q.VarPsi(), tw.Q.VarAlpha()]  
    fe_params['Q_cov'] = [tw.Q.CovPsiAlpha()]
    
    fe_estimator = tw.FEEstimator(bdf, fe_params)
    fe_estimator.fit()

2. Standard AKM, collapsing at the spell level and using leave-out-spell methods for HE.

I run exactly the same code on the same data as above but change the clean_params to:

clean_params = bpd.clean_params({
        'connectedness': 'leave_out_spell',
        'drop_single_stayers': True,
        'collapse_at_connectedness_measure': True,
        'copy': False
    })

Yet, my results are even more puzzling. var(eps) now exceeds var(y) for all models---even FE!---, and by a large amount. var(alpha)_he is negative, which clearly cannot be.

3. BLM with no controls and without collapsing at the spell level (analogous to case 1)

My setup is very basic:

input_df = pd.DataFrame({
        'i': panel['cluster_id'],
        'j': panel['hgroup_id_temp'],
        'y': panel['patents'],
        't': panel['year']
    })

    nl = 3 # Number of worker types
    nk = 4 # Number of firm types
    blm_params = tw.blm_params({
        'nl': nl, 'nk': nk,
        'a1_mu': 0, 'a1_sig': 1.5, 'a2_mu': 0, 'a2_sig': 1.5,
        's1_low': 0.5, 's1_high': 1.5, 's2_low': 0.5, 's2_high': 1.5
    })
    cluster_params = bpd.cluster_params({
        'measures': bpd.measures.CDFs(),
        'grouping': bpd.grouping.KMeans(n_clusters=nk),
        'is_sorted': True,
        'copy': False
    })
    clean_params = bpd.clean_params({
        'drop_returns': 'returners',
        'copy': False
    })
    bdf = bpd.BipartiteDataFrame(input_df, track_id_changes=True)
    bdf = bdf.clean(clean_params).collapse(is_sorted=True)
    bdf = bdf.cluster(cluster_params)
    bdf = bdf.to_eventstudy(is_sorted=True, copy=False)

    movers = bdf.get_worker_m(is_sorted=True)
    jdata = bdf.loc[movers,:]
    sdata = bdf.loc[~movers,:]

    blm_fit = tw.BLMEstimator(blm_params)
    blm_fit.fit(jdata=jdata, sdata=sdata, n_init=100, n_best=5, ncore=12)
    blm_model = blm_fit.model

    var_decomp = tw.BLMVarianceDecomposition(blm_params)
    var_decomp.fit(jdata=jdata, sdata=sdata, blm_model=blm_model, n_samples=100, Q_var=[tw.Q.VarPsi(), tw.Q.VarAlpha()], ncore=12)

    results = {k: np.mean(v) for k, v in var_decomp.res.items()}
    results_df = pd.DataFrame.from_dict(results, orient='index', columns=['mean_value'])
    results_df.to_csv(temp_dir / "blm_variance_decomposition.csv")

Once again, var(eps) exceeds var(y):

What do you think could be going on? And do you have any suggestions for fixing it?
I haven't noticed a big difference by increasing the number of draws and iterations from the default to the level I have now. They're not as high as they can be---I kept them manageable as I was still tying out different approaches. Do you think more draws/iterations will help? How high would you go? For reference, in my non-collapsed connected set I have about 3 million observations for 265k workers and 30k firms.

Could it be a version compatibility issues with some of the package dependencies?

[adamoppenheimer]

Hi Pier,

Thank you for using PyTwoWay!

Sorry for not directly answering all your questions here, instead I'll reference some previous issues people have posted:

Please check out #48 for a discussion of negative bias-corrected variance estimates.

Please check out #52 to see how to deal with the variance decomposition on BLM with collapsed data.

Finally, regarding your main issue about var(eps), there is a fundamental issue with how to properly deal with the residual for collapsed data.

For background, the issue derives from which var(eps) you are interested in. You may be interested in var(eps) from the collapsed data, or from the uncollapsed data. This could explain why in your collapsed data, var(y) is less than var(eps), since the var(eps) is for the uncollapsed data, and collapsing the data will mechanically reduce var(y).

Once you have determined the var(eps) you want, if you are interested in the variance of the residual in the uncollapsed data, the current code on GitHub makes the assumption that the errors are uncorrelated within-spell. This is a strong assumption, as it scales the residual by the weight in the collapsed data. You can see this here when the biased sigma^2 (var(eps)) is computed. The same assumption is used for estimating the HE var(eps). Unfortunately I can't remember how the HO correction deals with this, but I don't think it requires this assumption.

In some new, unpublished code I made multiple options about how to deal with this, but I don't know when I will publish it.

Regardless of the current options, I believe the more reasonable approach to handle this is to estimate the total variance explained by the explanatory variables, and estimate var(eps) as the residual. This way you don't need to make an assumption about the within-spell correlation of residuals. If you are interested in the uncollapsed data results, then you should include weights during estimation on the collapsed data, which will ensure that the var(alpha), var(psi), etc. are scaled correctly for the uncollapsed data. Then compute var(y) for the uncollapsed data, and take the residual to back out var(eps). If instead you are interested in the results for the collapsed data, you should consider estimating without weights, and back out var(eps) using the residual from var(y) for the collapsed data.

Finally, it could be helpful to simulate data from a process with a similar structure to your real data, but where you know the underlying parameters. Then you can estimate using PyTwoWay and verify that the estimates are correct.

I hope this helps, but please let me know if anything is unclear and if you have any more questions!

Best,
Adam

[pierpaolocreanza]

Hi Adam,

Thank you very much for your prompt and helpful reply! I appreciate it. A few follow-ups:

(1) Replacing w1 and w2 with all ones did solve the variance decomposition inconsistency in BLM, thank you! Are there cases where one does want to weight? Or should this be the default?

(2) For the HE-correction in the uncollapsed panel not summing to 1, I've gone through some of the suggestions.
(a) I'm taking a while to make attrition plots, but I'll post an update if anything interesting emerges.
(b) Increasing the trace and leverage draws further did not change much.
(c) You also suggested in #48 to explicitly include var(alpha), which is already what I do. Would I be justified in deriving var(eps) as a residual, or your argument above serial correlation only applies to collapsed data? If var(eps) and var(alpha) [or any other term] seem incompatible with each other, is there an argument for 'believing' one over the other?

(3) More generally, since I believe serial correlation to be strong in my application, your advice to estimate everything at the collapsed (weighted) level, and then back out var(eps) as a residual seems great. Thank you! Is there a reference that justifies this or is it more of ad hoc solution? Sorry if this question is a bit outside the scope of the package, but I'm really new to these issues!

Best,
Pier

[adamoppenheimer]

Hi Pier,

(1) That's great that it's working now. From what I understand, there should never be a case where you want to include weights in the variance decomposition for BLM.

(2c) For uncollapsed data, I believe the estimator for var(eps) should be correct. The issue with collapsed data is about rescaling it to represent the variance in the original data, which relies on assumptions about serial correlation. For uncollapsed data, you don't rescale it.

Also, I came up with an alternative proposal to back out var(alpha) when the bias-correction indicates it is negative. You can include the Q matrix to estimate var(psi), var(alpha), cov(psi, alpha), and var(psi + alpha) (this file has all the Q matrices you can compute). So if var(alpha) is negative, and you believe it is unreliable, you can potentially back out an alternative estimate by taking var(psi + alpha) - var(psi) - 2cov(psi, alpha). Ideally, this should exactly equal var(alpha), but I am not sure that the bias corrections have the property that var(psi + alpha) = var(psi) + var(alpha) + 2cov(psi, alpha) exactly, the relationship might hold only in expectation. So it's possible this will lead to a non-negative estimate for var(alpha), but it will still be unbiased (maybe it would be better to take the average over the two estimates). There may be a legitimate argument behind this approach, since var(psi + alpha) and (var(psi) + 2cov(psi, alpha)) could potentially be less subject to limited mobility bias than var(psi) or var(alpha) independently (I am not sure if this is true though).

(3) There may be a reference on this, but I haven't seen any. However, I just took a quick glance at my derivations and I realize that the bias-corrections rely on estimates of the variance of the residuals, which will depend on serial correlation. So if you believe the errors are highly correlated within-spell, it would probably make sense to use the unweighted estimators on the collapsed data, since the average over the spell is already highly informative. I think weights would make the most sense if you believe that longer spells are more informative about fixed effects, but with highly correlated errors they won't be. Without weights, you can do the full analysis using the unweighted var(y) on the collapsed data.

Ideally you might want to estimate the FE on the collapsed data without weights (since the additional observations within a spell aren't informative) but the var(eps) on the uncollapsed data so that you can account for the high correlation within-spell. Unfortunately I don't have code to do this at the moment.

I hope this helps, and again please let me know if you have any more questions!

Best,
Adam