# Econometrics.jl vs. R's plm

**URL:** <https://discourse.julialang.org/t/econometrics-jl-vs-rs-plm/68345>\
**Category:** Statistics\
**Tags:** econometrics\
**Created:** [September 17, 2021, 4:52pm UTC](https://discourse.julialang.org/t/econometrics-jl-vs-rs-plm/68345 "2021-09-17T16:52:56Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![mthelm85](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mthelm85/32/224164_2.png) [@mthelm85](https://discourse.julialang.org/u/mthelm85)\
**Post date:** [September 17, 2021, 4:52pm UTC](https://discourse.julialang.org/t/econometrics-jl-vs-rs-plm/68345/1 "2021-09-17T16:52:56Z")

</div>

I’m working through the [Econometrics With R online book](https://www.econometrics-with-r.org/), specifically the example in [10.3 on fixed effects regression](https://www.econometrics-with-r.org/10-3-fixed-effects-regression.html), but I’m trying to follow along with [Econometrics.jl](https://nosferican.github.io/Econometrics.jl/stable/estimators/#The-Basics-1).

Here’s the full code to reproduce the example:

```julia
using DataFramesMeta
using Econometrics
using RCall

df = rcopy(R"
library('AER')
data('Fatalities')
")

df = @chain @rget(Fatalities) begin
    @transform(fatal_rate = :fatal ./ :pop .* 10_000)
end

model = fit(
    EconometricModel,
    @formula(fatal_rate ~ beertax + absorb(state)),
    df
)

# Outputs:

Continuous Response Model
Number of observations: 336
Null Loglikelihood: -287.50
Loglikelihood: 105.99
R-squared: 0.9050
Wald: 12.19 ∼ F(1, 287) ⟹ Pr > F = 0.0006
Formula: fatal_rate ~ 1 + beertax + absorb(state)
Variance Covariance Estimator: OIM
──────────────────────────────────────────────────────────────────────────
                 PE SE t-value Pr > |t| 2.50% 97.50%
──────────────────────────────────────────────────────────────────────────
(Intercept) 2.37707 0.0969699 24.5135 <1e-71 2.18621 2.56794
beertax -0.655874 0.18785 -3.49148 0.0006 -1.02561 -0.286135
──────────────────────────────────────────────────────────────────────────

```

The book shows the output from R’s `plm` package and there are some differences:

```julia
#> Estimate Std. Error t value Pr(>|t|)  
#> beertax -0.65587 0.28880 -2.271 0.02388 *

```

Note that:

1. There is no intercept
2. The beertax coefficient is the same but the standard error and subsequent t and p values are different

Can anyone explain the difference? FixedEffectModels.jl produces an output that is much closer to the `plm` output:

```julia
using FixedEffectModels

reg(df, @formula(fatal_rate ~ beertax + fe(state)), Vcov.cluster(:state))

                          Fixed Effect Model
=======================================================================
Number of obs: 336 Degrees of freedom: 2
R2: 0.905 R2 Adjusted: 0.905
F-Stat: 5.05015 p-value: 0.029
R2 within: 0.041 Iterations: 1
=======================================================================
fatal_rate | Estimate Std.Error t value Pr(>|t|) Lower 95% Upper 95%
-----------------------------------------------------------------------
beertax | -0.655874 0.291856 -2.24725 0.025 -1.22998 -0.0817668
=======================================================================

```

---

<div class="post-metadata">

**Author:** ![junder873](https://avatars.discourse-cdn.com/v4/letter/j/e95f7d/32.png) [@junder873](https://discourse.julialang.org/u/junder873)\
**Post date:** [September 17, 2021, 8:05pm UTC](https://discourse.julialang.org/t/econometrics-jl-vs-rs-plm/68345/2 "2021-09-17T20:05:34Z")

</div>

First, in terms of the intercept. There is no consistency between packages on what exactly the intercept means in fixed effect models. Since using fixed effects is equivalent to adding an indicator variable for each category, the most common method is to somehow average across these different variables. Therefore, depending on who you ask, the intercept does or does not have meaning.

For the difference in standard errors, it looks like there is a difference in the assumptions made. Economectrics.jl seems to be just using normal standard errors, these are typically the smallest. `plm`, at least according to the book chapter you cited, is using robust standard errors, which allow for heteroskedasticity in the model (basically, the variance in the error term is not consistent). FixedEffectModels is using clustered standard errors, which allow for errors to be correlated across observations.

Try rerunning the FixedEffectModels as:

```julia
reg(df, @formula(fatal_rate ~ beertax + fe(state)), Vcov.robust())

```

Which should produce the same result.

---

<div class="post-metadata">

**Author:** ![Nosferican](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/nosferican/32/9275_2.png) [@Nosferican](https://discourse.julialang.org/u/Nosferican)\
**Post date:** [September 18, 2021, 3:24pm UTC](https://discourse.julialang.org/t/econometrics-jl-vs-rs-plm/68345/3 "2021-09-18T15:24:21Z")

</div>

Just adding a couple notes to @junder873 answer. The intercept when using the within estimator can be arbitrary and without a consistent interpretation depending on how it is computed. Econometrics.jl approach follows the principal that when possible, the intercept reported should match the one you would get if you had run the dummy variables model. In that sense, the intercept keeps the same value and interpretation it would have regardless of whether you estimate the model with all the features or `absorb` them.

In terms of the variance covariance estimators, it is now known that using cluster-robust vce is a bad default since it really depends on the assumed sampling design and whether we believe the treatment assignment is correlated with the sampling. In other words, just because we have “groups”, it doesn’t mean we have clusters (i.e., groups with different marginal effects) and that the sampling design causes clusters to have non-representative representations. If we do believe that is the case, you can choose a cluster robust variance covariance estimator that’s more appropriate. Looking at the current version, I believe the default is the observed information matrix (assumes White noise). Then you have your heteroscedasticity consistent options HC0, HC1, HC2, HC3, HC4. I can’t recall why I removed the cluster robust one… At some point, I think maybe I opted to incorporate [https://github.com/gragusa/CovarianceMatrices.jl](https://github.com/gragusa/CovarianceMatrices.jl) and never did for Arrellano/Rogers CRVE… I probably should have added those back as the extended HC1 with some residual degree adjustment heuristics.