I'm trying to do as much as I can with R instead of Stata, but Stata has positive network externalities in the economics profession and it's hard to move away if you have coauthors using it. Anyway, one of the most common regressions I have to run is a fixed effects regression with clustered standard errors. Suppose that Y is your dependent variable, X is an explanatory variable and F is a categorical variable that defines your fixed effects. You also want to cluster your standard errors according to F. If F has a lot of categories, demeaning instead of including the fixed effects explicitly has to be used so as not to have an excessive number of regressors. This would be the code in Stata:

xtset F

xtreg Y X, fe cluster(F) robust dfadj

The "dfadj" option is something I discovered recently that is needed to get the correct standard errors (it stands for "degrees of freedom adjustment" and at least for me, it was quite hidden in the manuals).

To produce the same results in R, you have to use the felm function of the lfe package. This is the code:

model <- felm(Y ~ X | F | 0 | F)

summary(model)

In the felm function arguments are separated by "|". The first argument is the equation to be estimated, the next one is the categorical variable that defines the fixed effects to demean the variables. The third one, in this case "0", could be used to introduce the instruments in instrumental variable estimation, and the last one defines the clustering of the standard errors. This R function produces exactly the same results as the Stata code above.