While replicating the IV model in R using the *fixest* library ([[Statistical modeling of IV]]) and in Python using the *linearmodels* library ([[Replication of IV in Python]]), we had a discrepancy in the way the standard errors are corrected for within-state correlation (clustered). Let's start with the simplest model, without any clustering first (so the standard errors are assumed to be [[Independent and identically distributed random variables|i.i.d.]]): ```r df <- readr::read_csv('cpi-policy-lobbyist-as-iv_v202.csv') fixest::etable(fixest::feols(consumer_price_index ~ esworker_health_bill | as.factor(state) + as.factor(date), data = df), vcov = "iid") fixest::feols(cons.. Dependent Var.: consumer_price_index esworker_health_bill 10.68*** (1.801) Fixed-Effects: -------------------- as.factor(state) Yes as.factor(date) Yes ____________________ ____________________ S.E. type IID Observations 302 R2 0.98019 Within R2 0.14473 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` In *R/fixest*, we then apply the cluster correction to the standard errors. Actually, all we need to do is remove the last argument, because the standard errors are clustered by the group variable "state" by default. The following is the equivalent of replacing the last argument value with "cluster" (`vcov = "cluster"`). ```r fixest::etable(fixest::feols(consumer_price_index ~ esworker_health_bill | as.factor(state) + as.factor(date), data = df)) fixest::feols(cons.. Dependent Var.: consumer_price_index esworker_health_bill 10.68* (3.897) Fixed-Effects: -------------------- as.factor(state) Yes as.factor(date) Yes ____________________ ____________________ S.E.: Clustered by: state) Observations 302 R2 0.98019 Within R2 0.14473 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` Ideally, with the assumptions discussed in [[Data and conceptual model]], this is our (correctly) corrected standard error, which takes into account the within-state correlations of the consumer price index. To have a referee for this comparison, we have also [[Naïve model in Stata|replicated the model in Stata]] using the `xtreg` function. Why do this? Because Stata's `vce(cluster clustvar)` has a well-documented behavior as discussed in [[Clustered standard errors]] and performs the intended correction we want to apply here. Therefore, we know that this is the correct standard error unless Stata screwed up its code for the function. But then we also know that Stata didn't screw up its code, because the results are identical to the *R/fixest* results. Now let's turn to Python and take a look at the naïve [[Independent and identically distributed random variables|i.i.d.]] model first: ```python import pandas as pd from linearmodels import PanelOLS # Load the data df = pd.read_csv('cpi-policy-lobbyist-as-iv_v202.csv') # Convert date to datetime df['date'] = pd.to_datetime(df['date']) # Set state and date as index df = df.set_index(['state', 'date']) # Create the panel regression model model = PanelOLS.from_formula('consumer_price_index ~ esworker_health_bill + EntityEffects + TimeEffects', data=df) # Fit the model results = model.fit() # Display the results summary print(results) Parameter Estimates ================================================================= Parameter Std. Err. T-stat P-value ----------------------------------------------------------------- esworker_health_bill 10.683 1.8007 5.9328 0.0000 ================================================================= F-test for Poolability: 91.695 P-value: 0.0000 Distribution: F(92,208) Included effects: Entity, Time ``` So far, so good. The standard error is the same as in the previous (R) output. But we know that the errors are not [[Independent and identically distributed random variables|i.i.d.]], so we need to do something about it. Let's cluster the standard errors using the group (called entity in *Python/linearmodels*): state. ```python # Fit the model with clustered standard errors by state results = model.fit(cov_type = 'clustered', cluster_entity = True) Parameter Estimates ================================================================= Parameter Std. Err. T-stat P-value ----------------------------------------------------------------- esworker_health_bill 10.683 3.7121 2.8779 0.0044 ================================================================= F-test for Poolability: 91.695 P-value: 0.0000 Distribution: F(92,208) Included effects: Entity, Time ``` The effect size is again consistent, as expected, but the standard error (3.71) is lower than what we believe to be the true corrected standard error, taking into account the correlation within states (3.90). This is a larger difference than might otherwise be attributed to computational differences, such as the use of different matrix inversion algorithms. **So what's going on here?** Well, the culprit is a behind-the-scenes computational correction called small sample adjustment. Small sample adjustment originates from the inherent trade-off of using a sample for population-level inference. The larger the sample, the smaller the impact of the correction on the standard error, because the adjustment is a function of the sample size and the number of parameters estimated. While this part is fairly straightforward, the adjustment itself can be applied in many different ways. This is where *R/fixest* and *Python/linearmodels* differ. So how do we know? Let's remove the small sample correction from both estimations. First, in *R/fixest*, the function to manage the small sample correction is `ssc()`. We need to make two changes here: 1. `adj = False` 2. `cluster.adj = False` The first argument to the function, `adj`, is a binary decision whether to apply the correction (in the form of `(n - 1) / (n - K)` with a sample size of `n` and the number of estimated parameters `K`). There is another argument we need to use here because we are clustering the standard errors. It's `cluster.adj` and it adds another term to the functional form of the correction, `G/(G-1)`, where `G` is the number of unique elements of the cluster variable. Note that `cluster.adj` must be set to `False` in addition to `adj`, because setting `adj` to `False` removes `(n - 1) / (n - K)` from the variance-covariance matrix correction, but keeps `G/(G-1)`. In other words, the first argument multiplies the variance-covariance matrix by "the `n` correction", and the second argument multiplicatively adds "the `G` correction": $(n - 1) / (n - K) * G/(G-1)$ Now let's go back to the original estimate and remove both small sample corrections: ```r fixest::etable(fixest::feols(consumer_price_index ~ esworker_health_bill | as.factor(state) + as.factor(date), data = df), vcov = "cluster", ssc = fixest::ssc(adj = FALSE, cluster.adj = FALSE)) fixest::feols(cons.. Dependent Var.: consumer_price_index esworker_health_bill 10.68* (3.081) Fixed-Effects: -------------------- as.factor(state) Yes as.factor(date) Yes ____________________ ____________________ S.E.: Clustered by: state) Observations 302 R2 0.98019 Within R2 0.14473 --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 ``` So with clustering, but without the small sample correction, the standard error in the R/fixest output is 3.081. Now let's turn to Python/linearmodels. As in R/fixest, removing the small sample correction requires multiple arguments. We need to change the values of the following three arguments in the `fit()` function (by default they are all `True`): 1. Set `debiased = False` to remove the small sample correction, `n / (n - K)`, from the estimation. 2. Set `auto_df = False` to disable the automatic adjustment of the degrees of freedom for fixed effects (`n - K - G`). 3. Set `count_effects = False` to disable explicit counting of fixed effects in the degrees of freedom adjustment. This argument is used only when `auto_df=False`, and also results in the `n - K - G` correction when set to True (default). The difference between using `auto_df` and `count_effects` is a slight difference in how `G` is counted. Now let's apply the changes in Python: ```python # Fit the model with clustered standard errors by state but without the small sample correction results = model.fit(cov_type = 'clustered', cluster_entity = True, debiased=False, auto_df=False, count_effects=False) Parameter Estimates ================================================================= Parameter Std. Err. T-stat P-value ----------------------------------------------------------------- esworker_health_bill 10.683 3.0807 3.4678 0.0005 ================================================================= F-test for Poolability: 91.695 P-value: 0.0000 Distribution: F(92,208) Included effects: Entity, Time ``` Voilà! The standard error is exactly the same as in *R/fixest*. The mystery is solved.