R Under development (unstable) (2024-02-22 r85974 ucrt) -- "Unsuffered Consequences"
Copyright (C) 2024 The R Foundation for Statistical Computing
Platform: x86_64-w64-mingw32/x64

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(lfe)
Loading required package: Matrix
> options(lfe.threads=2,digits=5,warn=1)
> set.seed(65318)
> x <- rnorm(500)
> x2 <- rnorm(length(x))
> 
> ## create individual and firm
> id <- factor(sample(10,length(x),replace=TRUE))
> firm <- factor(sample(6,length(x),replace=TRUE,prob=c(2,rep(1,5))))
> 
> ## effects
> id.eff <- rnorm(nlevels(id))
> firm.eff <- rnorm(nlevels(firm))
> 
> ## left hand side
> y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] + rnorm(length(x))
> y2 <- 0.1*x + 0.2*x2 + id.eff[id] + firm.eff[firm] + rnorm(length(x))
> 
> ## estimate
> smry_est <- summary(est <- felm(y ~ x+x2 | id + firm))
> getfe(est)
           effect obs comp   fe idx
id.1   -1.6363769  43    1   id   1
id.2    0.0029098  51    1   id   2
id.3   -0.7658186  50    1   id   3
id.4   -0.8531414  51    1   id   4
id.5    2.0125831  35    1   id   5
id.6    1.1118543  56    1   id   6
id.7    1.3748350  53    1   id   7
id.8   -0.8797889  47    1   id   8
id.9   -2.5330891  55    1   id   9
id.10  -0.3791573  59    1   id  10
firm.1  0.0000000 134    1 firm   1
firm.2 -0.7685150  77    1 firm   2
firm.3  0.7307192  81    1 firm   3
firm.4  0.7251559  70    1 firm   4
firm.5 -1.5010210  61    1 firm   5
firm.6  2.2462290  77    1 firm   6
> smry_est_lm <- summary(est_lm <- lm(y ~ x + x2 + id + firm))
> smry_est_lm

Call:
lm(formula = y ~ x + x2 + id + firm)

Residuals:
   Min     1Q Median     3Q    Max 
-3.360 -0.665 -0.030  0.764  2.618 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept)  -1.6364     0.1730   -9.46  < 2e-16 ***
x             0.9692     0.0488   19.88  < 2e-16 ***
x2            0.3346     0.0484    6.91  1.6e-11 ***
id2           1.6393     0.2148    7.63  1.2e-13 ***
id3           0.8706     0.2162    4.03  6.6e-05 ***
id4           0.7832     0.2158    3.63  0.00031 ***
id5           3.6490     0.2361   15.46  < 2e-16 ***
id6           2.7482     0.2105   13.05  < 2e-16 ***
id7           3.0112     0.2139   14.08  < 2e-16 ***
id8           0.7566     0.2193    3.45  0.00061 ***
id9          -0.8967     0.2111   -4.25  2.6e-05 ***
id10          1.2572     0.2087    6.02  3.4e-09 ***
firm2        -0.7685     0.1495   -5.14  4.0e-07 ***
firm3         0.7307     0.1479    4.94  1.1e-06 ***
firm4         0.7252     0.1540    4.71  3.3e-06 ***
firm5        -1.5010     0.1604   -9.36  < 2e-16 ***
firm6         2.2462     0.1500   14.97  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 1.03 on 483 degrees of freedom
Multiple R-squared:  0.787,	Adjusted R-squared:  0.78 
F-statistic:  112 on 16 and 483 DF,  p-value: <2e-16

> 
> all.equal(coef(smry_est), coef(smry_est_lm)[c("x", "x2"),])
[1] TRUE
> 
> ## vcov, confint
> all.equal(vcov(est), vcov(est_lm)[c("x", "x2"), c("x", "x2")])
[1] TRUE
> all.equal(confint(est), confint(est_lm, parm = c("x", "x2")))
[1] TRUE
> all.equal(stats::confint.default(est), stats::confint.default(est_lm, parm = c("x", "x2")))
[1] TRUE
> 
> ## should differ:
> isTRUE(all.equal(confint(est, type ="iid"), confint(est, type ="robust")))
[1] FALSE
> 
> ## multi response
> est_felm_multi <- felm(y + y2 ~ x+x2 | id + firm)
> est_lm_multi <- lm(cbind(y, y2) ~ x+x2 + id + firm)
> 
> all.equal(coef(est_felm_multi), coef(est_lm_multi)[c("x", "x2"),])
[1] TRUE
> all.equal(coef(summary(est_felm_multi, lhs = "y")),
+           coef(summary(est_lm_multi))[[1]][c("x", "x2"),])
[1] TRUE
> all.equal(coef(summary(est_felm_multi, lhs = "y2")),
+           coef(summary(est_lm_multi))[[2]][c("x", "x2"),])
[1] TRUE
> 
> all.equal(confint(est_felm_multi),
+           confint(est_lm_multi)[c("y:x", "y:x2", "y2:x", "y2:x2"),])
[1] TRUE
> 
> 
> proc.time()
   user  system elapsed 
   0.95    0.06    1.01