0
$\begingroup$

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!

|cite|improve this question
$\endgroup$
  • 2
    $\begingroup$ Take a look here for one potential difference. Also, note that the spelling is Stata since it is not an acronym. $\endgroup$ – Dimitriy V. Masterov Jul 25 '18 at 23:11
1
$\begingroup$

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.

|cite|improve this answer
$\endgroup$
  • $\begingroup$ the areg code is from the original study, only the R code is mine tho. I don't know the exact reason why they chose areg. Actually the SE is still very off in R. For example in STATA, the st.error for strike is 0.038 but in R its 0.778. Would the difference in areg and xtreg create such a big difference? $\endgroup$ – Annie Fannie Jul 26 '18 at 0:43
  • $\begingroup$ I don't have your data or even complete code, so I cannot really help. The code above manages to replicate output to five digits. $\endgroup$ – Dimitriy V. Masterov Jul 26 '18 at 0:57
  • $\begingroup$ Here's the original study with the data and the code. The areg is on line 294 aeaweb.org/articles?id=10.1257/aer.104.9.2763 Thank you so much for helping me!! $\endgroup$ – Annie Fannie Jul 26 '18 at 1:10
  • $\begingroup$ I think there must be something wrong with how I applied it in coeftest, because the SE is so much bigger... $\endgroup$ – Annie Fannie Jul 26 '18 at 1:16
  • $\begingroup$ I don't have access to that journal, but maybe you can add the code they use and what your complete R code to the original post. $\endgroup$ – Dimitriy V. Masterov Jul 26 '18 at 1:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.