Different results from random effects plm (R) and xtreg (Stata)

Question

I'm trying to convert the analysis from an old research paper from Stata to R. However, I ran into a problem that I have been unable to fix.

When I compare the coefficients from the two programs, they are not the same, even though the input is the same. I found this thread which describes a similar problem but for fixed effects: Difference between fixed effects models in R (plm) and Stata (xtreg)

However, the answer there gets a lot smaller difference than I get, accounting just for the difference in how plm and xtreg handle year effects.

For example, using the V-dem v9 Country-Year Full+Others dataset https://www.v-dem.net/en/data/data-version-9/, I ran this:

    library(plm)
Vdemv9 <- readRDS("./Country_Year_V-dem_Full+others_R_v9/V-Dem-CY-Full+Others-v9.rds")

model2 <- plm(v2x_polyarchy ~ v2elembaut+v2elrgstry,
                data = Vdemv9,
                model = "random",
                index = c("country_id","year"))
    summary(model2)

## Results:
Coefficients:
             Estimate Std. Error z-value  Pr(>|z|)    
(Intercept) 0.3735057  0.0059080  63.221 < 2.2e-16 ***
v2elembaut  0.1105280  0.0020646  53.534 < 2.2e-16 ***
v2elrgstry  0.0600031  0.0023033  26.051 < 2.2e-16 ***

Stata gives me the following results:

xtset country_id year, yearly
xtreg v2x_polyarchy v2elembaut v2elrgstry

## Results
v2x_polyar~y |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
  v2elembaut |   .1105945   .0020701    53.43   0.000     .1065372    .1146518
  v2elrgstry |   .0601527   .0023079    26.06   0.000     .0556292    .0646761
       _cons |   .3733406   .0062298    59.93   0.000     .3611304    .3855508

Am I doing something wrong? If not, is this something I need to worry about? The difference is small at only .0000665 for the coefficient for v2elembaut, but I would have expected it to not be there at all.

Answer 1

If your data are unbalanced, you will get different results from plm and Stata for the Swamy-Arora methods for random effects. For balanced data, you should get identical estimates.

This difference is due to a slight variation in the implementation: see plm's vignette https://cran.r-project.org/web/packages/plm/vignettes/B_plmFunction.html, section "Unbalanced Panels" and especially the paper by Cottrell (2017) mentioned there with more details (and Stata's documentation for xtreg, but it is discussed in Cottrell (2017)).

Best to my knowledge, plm and EViews are the only major implementations that use the formulae as given by Baltagi/Chang (1994) for the unbalanced case of Swamy-Arora per default. Stata and gretl use a slightly different implementation, also when using Stata's sa option to specifially request Swamy-Arora. However, it is possible to replicate Stata's results (using the sa option) with plm (see plm`s vignette).