Random Effects Results in R (plm) Cannot Replicate Results in Stata

Question

I have been working on migrating a current project from Stata to R, where I have encountered difficulties with differing results of random effects regressions.

I have panel data from an experiment where the treatment dummy is perfectly correlated with the group indicator because it is time-invariant. This means that a fixed effects regression of the outcome variable on the treatment dummy is not possible – however, a random-effects regression should be, since it only partially time-demeans the data. I am willing to assume that the treatment dummy and other covariates are not correlated with the group-specific error.

In Stata, this worked without a problem. The random-effects regression of the continuous outcome variable on the treatment dummy gives a result that makes sense, and the fixed effects regression omits the treatment dummy, exactly as expected.

However, in R, using the plm package, it did not work. I have received the error message "empty model." Curiously, this is not the case if the model does not include the treatment-dummy but other variables as regressors that are not perfectly correlated with the group indicator. In this case, plm's default method "swar" gives the same results as Stata.

I have tried to use other methods that are supplied by plm, and only the "walhus" method does work. In the case of a regression with the treatment dummy as a covariate, this gives the same result on the coefficients as Stata. However, it gives different results for models without the treatment dummy. These differences are not huge but considerable.

So in conclusion, I am able to replicate Stata's results in R, but with different methods where Stata uses only one. I have not found an explanation for that behavior in the Stata Documentation or in the plm paper in the Journal of Statistical Software. The plm paper gives sources for the different methods for RE (that supposedly differ in their estimation of theta) but does not explain the differences itself. The original sources for "swar" and "walhus" are Econometrica papers from the late 60s / early 70s. Quite frankly, I was not able to find a solution in these either. I have also found this question on Stackexchange, but I believe that this is a different issue.

Any help or ideas would be much appreciated! This has already taken an immense ammount of time and I find it to be really troubling.

P.S. I cannot share the original data, but I have created a dataset with similar properties with which these problems can be replicated. I have put it into a dropbox, as .Rdata and .dta.

The "original" Stata code:

xtset GroupID Round


xtreg outcome Treatment, re
------------------------------------------------------------------------------
     outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   Treatment |   36.93656    5.97516     6.18   0.000     25.22546    48.64766
       _cons |   51.16955   4.225076    12.11   0.000     42.88855    59.45055
-------------+----------------------------------------------------------------


xtreg outcome X1, re
------------------------------------------------------------------------------
     outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
          X1 |  -.0278302   .1193763    -0.23   0.816    -.2618033     .206143
       _cons |   70.84536   6.953707    10.19   0.000     57.21635    84.47438
-------------+----------------------------------------------------------------

The corresponding R-code:

library(plm)
testdata <- pdata.frame(testdata, index=c("GroupID","Round"))


Model1 <- plm(outcome ~ Treatment, data = testdata, model="random", random.method="swar") 
summary(Model1) # This doesn’t work
Error in plm.fit(data, model = models[1], effect = effect) : empty model


Model2 <- plm(outcome ~ Treatment, data = testdata, model="random", random.method="walhus") 
summary(Model2) # This gives the same results as Stata
            Estimate Std. Error z-value  Pr(>|z|)    
(Intercept)  51.1695     4.2251 12.1109 < 2.2e-16 ***
Treatment    36.9366     5.9752  6.1817 6.342e-10 ***


Model3 <- plm(outcome ~ X1, data = testdata, model="random", random.method="swar")
summary(Model3) # This gives the same results as Stata
            Estimate Std. Error z-value Pr(>|z|)    
(Intercept) 70.84536    6.95371 10.1881   <2e-16 ***
X1          -0.02783    0.11938 -0.2331   0.8157    


Model4 <- plm(outcome ~ X1, data = testdata, model="random", random.method="walhus")
summary(Model4) # This gives slightly different results than Stata
             Estimate Std. Error z-value Pr(>|z|)    
(Intercept) 70.682277   7.003460 10.0925   <2e-16 ***
X1          -0.024072   0.119074 -0.2022   0.8398    

EDIT: I have tried something else and found that plm's default method "swar" does also work for a model that includes both the time-invariant treatment-dummy and a time-varying continuous covariate:

Model1.2 <- plm(outcome ~ Treatment + X1, data = testdata, model="random", random.method="swar")
summary(Model1.2) # This somehow works
             Estimate Std. Error z-value  Pr(>|z|)    
(Intercept) 14.906599  11.284649  1.3210    0.1865    
Treatment   36.835123   6.075290  6.0631 1.335e-09 ***
X1          -0.012018   0.108785 -0.1105    0.9120  

This gives the same results on the coefficients (but not the intercept) as Stata:

xtreg outcome Treatment X1, re
------------------------------------------------------------------------------
     outcome |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   Treatment |   36.83512    6.07529     6.06   0.000     24.92777    48.74247
          X1 |   -.012018   .1087849    -0.11   0.912    -.2252326    .2011965
       _cons |   51.74172   6.697543     7.73   0.000     38.61478    64.86866
-------------+----------------------------------------------------------------

Answer 1

While the question seems as a software question at first, there is some statistics behind it (and, thus, I deem this to be on-topic for xvalidated):

The random effect estimator as per Swamy-Arora uses the variation of the associated within model and the associated between model. For a plm-based exposition see one of the package's vignettes https://cran.rstudio.com/web/packages/plm/vignettes/B_plmFunction.html, section "Unbalanced panels" (but this is not specific to unbalances panels). Any good text book about panel models will cover this, e.g., Wooldridge or Baltagi. Other random effect estimators like Wallace-Hussain use slightly other "base models" (but Amemiya's estimator uses the within model twice), see Baltagi's text book for an overview.

Now, looking at the software implementation for plm if model = "swar": The function estimates a within model first. This fails (correctly) for the specific example you have as there is no within variation of the only covariate (Treatment, as you observed correctly). The function then does not continue to estimate the between model. Stata does continue for these data (and also gretl) and gives an output. Thus, the model you want to estimate is equivalent to the between specification. The between model can be estimated by:

plm(outcome ~ Treatment, data = testdata, model = "between") 

# Coefficients:
#              Estimate Std. Error t-value  Pr(>|t|)    
# (Intercept)  51.1695     3.7313 13.7135 5.722e-11 ***
# Treatment    36.9366     5.2769  6.9997 1.555e-06 ***

-> You get the same estimates for the coefficient as Stata gives (the difference for the standard errors is due to some adjustment either specific to RE models or to Stata, I suppose. Also, for RE models z values are usually applied as the finite distribution is usually not known).