<jupyter_start><jupyter_text>Non-linear Panel Models--- Note:- This notebook uses the "Crime from North Carolina" dataset from notebook 2- `statsmodels`'s conditional logit implementation fails to converge with the Panel101 "toy dataset"- The "toy dataset" does not have sufficient variation within groups. It also suffers from high multicollinearity between variables within groups (see Section 1.a below) Overview (and pre-analysis plan)In this notebook, we will:- Fit multiple panel models to non-linear data (i.e. data with binary outcomes)- Models include:    - Pooled logit    - Conditional logit (country fixed effects)    - Conditional logit (year fixed effects)    - LPM    - LPM (country fixed effects)    - LPM (year fixed effects)    - LPM (two-way fixed effects)<jupyter_code>import pandas as pd

import statsmodels.api as sm

from linearmodels import PanelOLS

from statsmodels.discrete.conditional_models import ConditionalLogit
from statsmodels.stats.outliers_influence import variance_inflation_factor<jupyter_output><empty_output><jupyter_text>1.a Data prep<jupyter_code># Load data (from http://www.princeton.edu/~otorres/LogitR101.pdf )

data = pd.read_stata('http://dss.princeton.edu/training/Panel101.dta')
data
# Seperate exogenous columns, endogenous columns, 
# country column, and year column

exog = data[['x1', 'x2', 'x3']]
endog = data['y_bin'].astype(int)
countries = data['country']
years = data['year']
# The "toy dataset" has very few observations per group

data.groupby(by='country').size()
# Get variance within each country for y_bin, x1, x2, and x3

within_country_var = data[['y_bin', 'x1', 'x2', 'x3', 'country']].groupby(by='country').var()
within_country_var
# Check multicollinearity within groups using variance inflation factor

def get_vif(data):
    """Returns vif of each column in `data` as a np.array.
    """
    vif = []
    exog = data.iloc[:,:-1].to_numpy()
    for i in range(len(data.columns)-1):
        vif.append(variance_inflation_factor(exog, i))
    vif = (pd.DataFrame(vif)
             .set_index(pd.Index(['x1', 'x2', 'x3'],
                                 name='exog_var')))
    return vif


vif_country_var = (
    data[['x1', 'x2', 'x3', 'country']]
        .groupby(by='country')
        .apply(lambda x: get_vif(x))
        .pipe(pd.DataFrame)
        .rename(columns={0:'vif'})
        .explode(column='vif')
        .reset_index(level=1)
        .pivot(columns='exog_var')
)
vif_country_var<jupyter_output><empty_output><jupyter_text>Note:- As we see above, countries C, G, E, and A and have high multicollinearity between variables x1 and x2- I've tried multiple hyperaparameters (intial values, numerical methods) for `statsmodels`'s conditional logit implementation- `statsmodels`'s conditional logit fails to converge across all my attempts- This example demonstrates that numerical issues **can differ across implementations** (e.g. R's clogit has no issues with this "toy dataset")- Let's use a real dataset instead: "Crime in North Carolina"  1.b Crime Dataset<jupyter_code>data = pd.read_csv('data/crime.csv', index_col=0)
data
# Let's create a indicator variable for crime rate (1 if greater than 0.025)

endog = data['crmrte'].apply(lambda x: 1 if x > 0.025 else 0)
# Balanced outcomes

endog.value_counts()
# Same exogenous variables as in notebook 2's model

exog = data[['density', 'taxpc', 'wcon', 'pctmin']]
# Set index (entity, time) for PanelOLS's API

data.loc[:, 'outcome'] = endog
panel_data = data.set_index(['county', 'year'])<jupyter_output><empty_output><jupyter_text>2. Non-linear panel models  2.1 Pooled logit<jupyter_code># Pooled logit

logit_mod = sm.Logit(endog=endog, exog=exog)
logit_res = logit_mod.fit()
logit_res.summary()<jupyter_output>Optimization terminated successfully.
         Current function value: 0.459943
         Iterations 8<jupyter_text>2.2 Conditional logit (county fixed effects)<jupyter_code># Conditional logit (country FE)

clogit_country_mod = ConditionalLogit(endog=endog,
                                      exog=exog,
                                      groups=data['county'])
clogit_country_res = clogit_country_mod.fit(method='nm')
clogit_country_res.summary()<jupyter_output>/home/topher-lo/repos/PyEC333/ec333_env/lib/python3.8/site-packages/statsmodels/discrete/conditional_models.py:80: UserWarning: Dropped 73 groups and 511 observations for having no within-group variance
  warnings.warn(msg)<jupyter_text>Note:- Try changing the numerical method to 'bfgs' (which is default method for `statsmodels`'s conditional logit), it will fail to converge  2.3 Conditional logit (year fixed effects)<jupyter_code># Conditional logit (year FE)

clogit_year_mod = ConditionalLogit(endog=endog,
                                   exog=exog,
                                   groups=data['year'])
clogit_year_res = clogit_year_mod.fit(method='nm', maxiter=500)
clogit_year_res.summary()<jupyter_output>/home/topher-lo/repos/PyEC333/ec333_env/lib/python3.8/site-packages/statsmodels/discrete/conditional_models.py:283: RuntimeWarning: overflow encountered in double_scalars
  v = f(t - 1, k) + f(t - 1, k - 1) * exb[t - 1]<jupyter_text>2.4 LPM<jupyter_code># Pooled LPM

lpm_mod = sm.OLS(endog=endog,
                 exog=exog)
lpm_res = lpm_mod.fit()
lpm_res.summary()<jupyter_output><empty_output><jupyter_text>2.5 LPM (county fixed effects)<jupyter_code># Entity fixed effects

entity_fe_formula = "outcome ~ density + taxpc + wcon + pctmin + EntityEffects"
mod_entity_fe = PanelOLS.from_formula(entity_fe_formula, data=panel_data, drop_absorbed=True)
res_entity_fe = mod_entity_fe.fit()
res_entity_fe.summary<jupyter_output>/home/topher-lo/repos/PyEC333/ec333_env/lib/python3.8/site-packages/linearmodels/panel/model.py:1743: AbsorbingEffectWarning: 
Variables have been fully absorbed and have removed from the regression:

pctmin

  warnings.warn(<jupyter_text>2.6 LPM (year fixed effects)<jupyter_code># Time period fixed effects

year_fe_formula = "outcome ~ density + taxpc + wcon + pctmin + TimeEffects"
mod_year_fe = PanelOLS.from_formula(year_fe_formula, data=panel_data, drop_absorbed=True)
res_year_fe = mod_year_fe.fit()
res_year_fe.summary<jupyter_output><empty_output><jupyter_text>2.7 LPM (county and year fixed effects)<jupyter_code># Entity and time period fixed effects

twoways_fe_formula = "outcome ~ density + taxpc + wcon + pctmin + TimeEffects + EntityEffects"
mod_twoways_fe = PanelOLS.from_formula(twoways_fe_formula, data=panel_data, drop_absorbed=True)
res_twoways_fe = mod_twoways_fe.fit()
res_twoways_fe.summary<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.11032082324455206, 'fraction_numerical': 0.01740314769975787, 'mean_word_length': 4.619897959183674, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 38, 'https://': 0, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 5, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '10', 'dir': 'jupyter-structured-clean-dedup', 'id': '354976', 'license': 'MIT', 'n_tokens_mistral': 2330, 'n_tokens_neox': 2209, 'n_words': 553, 'path': 'notebooks/8_nonlinear_panel.ipynb', 'repo_name': 'topher-lo/PyEC333', 'size': '6608'}