Clustered standard errors are completely different in R than in STATA

Question

I'm trying to reproduce a study in R. Here are its core elements:

• study wants to measure the effect of a transit strike on highway
  delay

• independent variables:

  strike: binary

  dateresidual: difference from the start of strike (negative for pre-strike, positive for during strike)

  dateresidual*strike

  4 binary variables for 4/5 weekdays

• each observation is measured by one of the thousands of road sensors (sensorid) for a particular hour of the day. So each vds would have multiple observations per day.

• The STATA code ran this with cluster(sensorid) and absorb(sensorid), meaning the standard errors are clustered at the sensor level and sensor id is the fixed effect.

• The regression has a weight for highway length/total flow

  areg delay strike dateresidual datestrike mon tue wed thu [aw=weight], cluster(sensorid) absorb(sensorid)

• I have tried to run this in r using plm. I was able to get the exact same estimates:

plm1<-plm(delay~strike+dateresidual+datestrike+mon+tue+wed+thu,mydata,model="within",index=c("sensorid"))

However, when I tried to run the clustered standard errors at sensor id, the standard errors are way off from the stata results and the effects are no longer significant.

coeftest(plm1,vcovHC)

Could you tell me what I should tweak in coeftest to represent what the code in STATA does? Or should I use a different package?

Thank you very much!

Answer 1

First, Stata uses a finite sample correction that R does not use when clustering. Second, areg is designed for datasets with many groups, but not a number that grows with the sample size. One example is states in the US. Your plm is much more like xtreg, fe. Although the point estimates produced by areg and xtreg, fe are the same, the estimated VCEs differ with clustering because the commands make different assumptions about whether the number of groups/sensors increases with the sample size. When you cluster with xtreg, fe, the asymptotics relies on the number of groups going to infinity. So this is not an apples to apples comparison.

In your setting, xtreg, fe seems more suitable since many sensors could be added. If you have to replicate areg's output, you can use felm.

Here is an econometrically stupid example demonstrating these claims. Here I am using Roger Newson's rsource to run R from within Stata, but it is not strictly necessary:

. rsource, terminator(END_OF_R)
Assumed R program path: "/usr/local/bin/R"

Beginning of R output
> suppressPackageStartupMessages({
+         require(plm)
+         require(lmtest)
+         library(foreign)
+         library(lfe)
+         data(Cigar)
+ })
> xtregfe <- plm(sales ~ price, model = 'within', data = Cigar)
> G <- length(unique(Cigar$state))
> c <- G/(G - 1)
> coeftest(xtregfe,c * vcovHC(xtregfe, type = "HC1", cluster = "group"))

t test of coefficients:

      Estimate Std. Error t value  Pr(>|t|)    
price -0.20984    0.03575 -5.8697 5.503e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> areg<-felm(sales ~ price | state | 0 | state, Cigar)
> coeftest(areg)

t test of coefficients:

       Estimate Std. Error t value  Pr(>|t|)    
price -0.209840   0.036348 -5.7731 9.672e-09 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

> write.dta(Cigar,"~/Desktop/Cigar.dta")
> 
End of R output

. 
. use "~/Desktop/Cigar.dta", clear
(Written by R.              )

. xtset state year
       panel variable:  state (strongly balanced)
        time variable:  year, 63 to 92
                delta:  1 unit

. xtreg sales price, fe cluster(state)

Fixed-effects (within) regression               Number of obs     =      1,380
Group variable: state                           Number of groups  =         46

R-sq:                                           Obs per group:
     within  = 0.2559                                         min =         30
     between = 0.1007                                         avg =       30.0
     overall = 0.0969                                         max =         30

                                                F(1,45)           =      34.45
corr(u_i, Xb)  = 0.0329                         Prob > F          =     0.0000

                                 (Std. Err. adjusted for 46 clusters in state)
------------------------------------------------------------------------------
             |               Robust
       sales |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       price |  -.2098402   .0357498    -5.87   0.000     -.281844   -.1378364
       _cons |   138.3669   2.456008    56.34   0.000     133.4202    143.3135
-------------+----------------------------------------------------------------
     sigma_u |  25.678164
     sigma_e |  15.174773
         rho |  .74116131   (fraction of variance due to u_i)
------------------------------------------------------------------------------

. areg sales price, absorb(state) cluster(state)

Linear regression, absorbing indicators         Number of obs     =      1,380
Absorbed variable: state                        No. of categories =         46
                                                F(   1,     45)   =      33.33
                                                Prob > F          =     0.0000
                                                R-squared         =     0.7682
                                                Adj R-squared     =     0.7602
                                                Root MSE          =    15.1748

                                 (Std. Err. adjusted for 46 clusters in state)
------------------------------------------------------------------------------
             |               Robust
       sales |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
       price |  -.2098402   .0363482    -5.77   0.000    -.2830493   -.1366312
       _cons |   138.3669   2.497119    55.41   0.000     133.3374    143.3963
------------------------------------------------------------------------------

As you can see, areg/felm give you a price coefficient of -0.20984 with a clustered standard error of 0.03635. The panel fixed effect approaches both give you -0.20984, but with a smaller CSE of 0.03575.