From 656fba3368821191937ae1cc0a00712617754194 Mon Sep 17 00:00:00 2001 From: Matthieu Stigler Date: Sat, 5 Dec 2020 15:29:36 -0800 Subject: [PATCH 01/87] lfe-Ex.Rout.save: Plain and blind update --- tests/Examples/lfe-Ex.Rout.save | 98 ++++++++++----------------------- 1 file changed, 29 insertions(+), 69 deletions(-) diff --git a/tests/Examples/lfe-Ex.Rout.save b/tests/Examples/lfe-Ex.Rout.save index b571d86..5c37589 100644 --- a/tests/Examples/lfe-Ex.Rout.save +++ b/tests/Examples/lfe-Ex.Rout.save @@ -59,9 +59,9 @@ Loading required package: Matrix > # find bias corrections > bccorr(est) corr v1 v2 cov d1 d2 - 0.008569979 0.933723478 1.297425684 0.009432572 0.082334457 0.061969275 + 0.008017737 0.934873578 1.301275975 0.008843272 0.081184357 0.058118985 d12 --0.004721002 +-0.004131701 > > > @@ -679,21 +679,21 @@ F-statistic(proj model): 252.8 on 2 and 499 DF, p-value: < 2.2e-16 > # find bias corrections, there's little bias in this example > print(yv <- fevcov(est, lhs='y')) id firm foo -id 1.031163938 0.006729247 -0.006303323 -firm 0.006729247 1.046116404 0.010550689 -foo -0.006303323 0.010550689 1.058185345 +id 1.031372636 0.006691582 -0.006322827 +firm 0.006691582 1.049086528 0.010535564 +foo -0.006322827 0.010535564 1.057014174 attr(,"bias") id firm foo -id 0.0307580151 -0.0002094436 -0.0001106874 -firm -0.0002094436 0.0224115429 -0.0001244024 -foo -0.0001106874 -0.0001244024 0.0132579840 +id 3.054932e-02 -0.0001717782 -9.118383e-05 +firm -1.717782e-04 0.0194414181 -1.092769e-04 +foo -9.118383e-05 -0.0001092769 1.442915e-02 > ## Here's how to compute the unbiased correlation matrix: > cm <- cov2cor(yv) > structure(cm,bias=NULL) id firm foo -id 1.000000000 0.006479069 -0.006034272 -firm 0.006479069 1.000000000 0.010027899 -foo -0.006034272 0.010027899 1.000000000 +id 1.000000000 0.006433027 -0.006055682 +firm 0.006433027 1.000000000 0.010004876 +foo -0.006055682 0.010004876 1.000000000 > > > @@ -718,11 +718,11 @@ foo -0.006034272 0.010027899 1.000000000 > #both bootstrap and computed se: > cbind(getfe(est,ef=efactory(est,'ref'),se=TRUE), fse=fixedse(est)) effect obs comp fe idx se fse -f.1 0.9569311 185 1 f 1 0.07807322 0.07535548 -f.2 1.9949066 223 1 f 2 0.06163607 0.06860181 -f.3 3.0326896 185 1 f 3 0.07650146 0.07532290 -f.4 4.0294033 196 1 f 4 0.07595171 0.07317074 -f.5 4.9413841 211 1 f 5 0.07077708 0.07054858 +f.1 0.9569311 185 1 f 1 0.07749880 0.07535548 +f.2 1.9949066 223 1 f 2 0.06038432 0.06860181 +f.3 3.0326896 185 1 f 3 0.07696251 0.07532290 +f.4 4.0294033 196 1 f 4 0.07527526 0.07317074 +f.5 4.9413841 211 1 f 5 0.07015206 0.07054858 > #compare with lm: > summary(lm(y ~x+f-1)) @@ -1114,24 +1114,24 @@ Levels: 1 2 3 [1] 0.6807816+0i > # mctrace is not really useful for small problems. > mctrace(fun,ncol(A),tol=0.05) -[1] 0.6590496 +[1] 0.6600197 attr(,"sd") -[1] 0.03281585 +[1] 0.03292458 attr(,"iterations") -[1] 11102 +[1] 11040 > # try a larger problem (3000x3000): > f1 <- factor(sample(1500,3000,replace=TRUE)) > f2 <- factor(sample(1500,3000,replace=TRUE)) > fl <- list(f1,f2) > mctrace(fl,tol=-5) -[1] 503.8349 +[1] 499.8301 attr(,"sd") -[1] 4.485407 +[1] 3.43084 attr(,"iterations") -[1] 28 +[1] 56 > # exact: > length(f1) - nlevels(f1) - nlevels(f2) + nlevels(compfactor(fl)) -[1] 503 +[1] 495 > > > @@ -1183,13 +1183,6 @@ F-statistic(proj model): 1.883 on 2 and 97 DF, p-value: 0.1577 > nlexpect(est, quote(c(x1=x1,x2=x2,var=c((x1-pt1)^2,(x2-pt2)^2)))) x1 x2 var.x1 var.x2 0.17450420 0.46517350 0.02611881 1.22460757 -> ## No test: -> # the covariance matrix: -> nlexpect(est, tcrossprod(as.matrix(c(x1-pt1,x2-pt2)))) - x1 x2 -x1 0.02609831 -0.122408 -x2 -0.12240804 1.224730 -> ## End(No test) > #Wald test of single variable > waldtest(est, ~x1)['p.F'] p.F @@ -1209,39 +1202,6 @@ x2 -0.12240804 1.224730 x1 a 0.8587288 b 0.5324731 -> ## No test: -> -> # Non-linear test: -> -> # A simple one, what's the probability that product x1*x2 is between 0 and |E(x1)|? -> nlexpect(est, x1*x2 > 0 & x1*x2 < abs(pt1), vectorize=TRUE, method='divonne') - x1 -0.3914772 -> # Then a more complicated one with the expected value of a polynomal in the coefficients -> f <- function(x) c(poly=x[['x1']]*(6*x[['x1']]-x[['x2']]^2)) -> # This is the linearized test: -> waldtest(est, f)['p.F'] - p.F -0.7432511 -> # In general, for a function f, the non-linear Wald test is something like -> # the following: -> # expected value of function -> Ef <- nlexpect(est, f, coefs=c('x1','x2')) -> # point value of function -> Pf <- f(c(pt1,pt2)) -> # similar to a Wald test, but non-linear: -> nlexpect(est, function(x) (f(x)-Ef)^2 > Pf^2, c('x1','x2'), vectorize=TRUE) - poly -0.6851808 -> # one-sided -> nlexpect(est, function(x) f(x)-Ef > abs(Pf), c('x1','x2'), vectorize=TRUE) - poly -0.267744 -> # other sided -> nlexpect(est, function(x) f(x)-Ef < -abs(Pf), c('x1','x2'), vectorize=TRUE) - poly -0.4193062 -> ## End(No test) > > > @@ -1284,14 +1244,14 @@ b 0.5324731 > # estimate the covariance matrix of the fixed effects > fevcov(est, alpha) id firm foo -id 0.93372348 0.01224916 -0.11164451 -firm 0.01224916 4.95014733 0.04605155 -foo -0.11164451 0.04605155 6.48715940 +id 0.93487358 0.01198024 -0.11085834 +firm 0.01198024 4.95403947 0.04691949 +foo -0.11085834 0.04691949 6.49011310 attr(,"bias") id firm foo -id 0.082334457 -0.004334716 -0.002472491 -firm -0.004334716 0.060817078 -0.002995544 -foo -0.002472491 -0.002995544 0.044123684 +id 0.081184357 -0.004065798 -0.003258662 +firm -0.004065798 0.056924931 -0.003863477 +foo -0.003258662 -0.003863477 0.041169982 > # estimate variances of the diagonal > varvars(est, alpha) [1] 0.008598521 0.043820907 0.057656021 From 046b7b706006690136508218888b280e2d519ef2 Mon Sep 17 00:00:00 2001 From: Matthieu Stigler Date: Sat, 5 Dec 2020 15:38:18 -0800 Subject: [PATCH 02/87] Update some urls from https to https (R devel on R-hub was complaining about it) --- R/condfstat.R | 4 ++-- R/lfe-package.R | 8 ++++---- man/cgsolve.Rd | 2 +- man/condfstat.Rd | 4 ++-- man/lfe-package.Rd | 12 ++++++------ 5 files changed, 15 insertions(+), 15 deletions(-) diff --git a/R/condfstat.R b/R/condfstat.R index a124de8..f8c0179 100644 --- a/R/condfstat.R +++ b/R/condfstat.R @@ -140,10 +140,10 @@ ivbootstrap <- function(z, x, y, quantiles=0.95, N=100L, cluster=NULL) { #' @references Sanderson, E. and F. Windmeijer (2014) \cite{A weak instrument #' F-test in linear IV models with multiple endogenous variables}, Journal of #' Econometrics, 2015. -#' \url{http://www.sciencedirect.com/science/article/pii/S0304407615001736} +#' \url{https://www.sciencedirect.com/science/article/pii/S0304407615001736} #' #' Stock, J.H. and M. Yogo (2004) \cite{Testing for weak instruments in linear -#' IV regression}, \url{http://ssrn.com/abstract=1734933} in +#' IV regression}, \url{https://ssrn.com/abstract=1734933} in #' \cite{Identification and inference for econometric models: Essays in honor #' of Thomas Rothenberg}, 2005. #' @examples diff --git a/R/lfe-package.R b/R/lfe-package.R index 1521fe5..eb3faa3 100644 --- a/R/lfe-package.R +++ b/R/lfe-package.R @@ -135,7 +135,7 @@ #' #' Croissant, Y. and G. Millo (2008) \cite{Panel Data Econometrics in R: The #' plm Package}, Journal of Statistical Software, 27(2). -#' \url{http://www.jstatsoft.org/v27/i02/} +#' \url{https://www.jstatsoft.org/v27/i02/} #' #' Gaure, S. (2013) \cite{OLS with Multiple High Dimensional Category #' Variables.} Computational Statistics and Data Analysis, 66:8--18, 2013 @@ -152,17 +152,17 @@ #' Guimaraes, P. and Portugal, P. (2010) \cite{A simple feasible procedure to #' fit models with high-dimensional fixed effects.} The Stata Journal, #' 10(4):629--649, 2010. -#' \url{http://www.stata-journal.com/article.html?article=st0212} +#' \url{https://www.stata-journal.com/article.html?article=st0212} #' #' Ouazad, A. (2008) \cite{A2REG: Stata module to estimate models with two #' fixed effects.} Statistical Software Components S456942, Boston College #' Department of Economics. -#' \url{http://ideas.repec.org/c/boc/bocode/s456942.html} +#' \url{https://ideas.repec.org/c/boc/bocode/s456942.html} #' #' Sanderson, E. and F. Windmeijer (2014) \cite{A weak instrument F-test in #' linear IV models with multiple endogenous variables}, Journal of #' Econometrics, 2015. -#' \url{http://www.sciencedirect.com/science/article/pii/S0304407615001736} +#' \url{https://www.sciencedirect.com/science/article/pii/S0304407615001736} #' @keywords regression models #' @examples #' diff --git a/man/cgsolve.Rd b/man/cgsolve.Rd index 5c7a5a2..7e52fd2 100644 --- a/man/cgsolve.Rd +++ b/man/cgsolve.Rd @@ -74,7 +74,7 @@ Thus, the test is neither deterministic nor perfect. \references{ Kaasschieter, E. (1988) \cite{A practical termination criterion for the conjugate gradient method}, BIT Numerical Mathematics, -28(2):308-322. \url{http://link.springer.com/article/10.1007\%2FBF01934094} +28(2):308-322. \url{https://link.springer.com/article/10.1007\%2FBF01934094} } \seealso{ \code{\link{kaczmarz}} diff --git a/man/condfstat.Rd b/man/condfstat.Rd index 4903eea..e40e838 100644 --- a/man/condfstat.Rd +++ b/man/condfstat.Rd @@ -93,10 +93,10 @@ condfstat(est, quantiles=c(0.05, 0.95)) Sanderson, E. and F. Windmeijer (2014) \cite{A weak instrument F-test in linear IV models with multiple endogenous variables}, Journal of Econometrics, 2015. -\url{http://www.sciencedirect.com/science/article/pii/S0304407615001736} +\url{https://www.sciencedirect.com/science/article/pii/S0304407615001736} Stock, J.H. and M. Yogo (2004) \cite{Testing for weak instruments in linear -IV regression}, \url{http://ssrn.com/abstract=1734933} in +IV regression}, \url{https://ssrn.com/abstract=1734933} in \cite{Identification and inference for econometric models: Essays in honor of Thomas Rothenberg}, 2005. } diff --git a/man/lfe-package.Rd b/man/lfe-package.Rd index d743227..a28abbf 100644 --- a/man/lfe-package.Rd +++ b/man/lfe-package.Rd @@ -150,16 +150,16 @@ J.R. Stat. Soc.(A) 171(3), 673--697. Cornelissen, T. (2008) \cite{The stata command felsdvreg to fit a linear model with two high-dimensional fixed effects.} Stata Journal, 8(2):170--189, 2008. -\url{http://econpapers.repec.org/RePEc:tsj:stataj:v:8:y:2008:i:2:p:170-189} +\url{https://econpapers.repec.org/RePEc:tsj:stataj:v:8:y:2008:i:2:p:170-189} Correia, S. (2014) \cite{REGHDFE: Stata module to perform linear or instrumental-variable regression absorbing any number of high-dimensional fixed effects}, Statistical Software Components, Boston College Department -of Economics. \url{http://econpapers.repec.org/RePEc:boc:bocode:s457874} +of Economics. \url{https://econpapers.repec.org/RePEc:boc:bocode:s457874} Croissant, Y. and G. Millo (2008) \cite{Panel Data Econometrics in R: The plm Package}, Journal of Statistical Software, 27(2). -\url{http://www.jstatsoft.org/v27/i02/} +\url{https://www.jstatsoft.org/v27/i02/} Gaure, S. (2013) \cite{OLS with Multiple High Dimensional Category Variables.} Computational Statistics and Data Analysis, 66:8--18, 2013 @@ -176,17 +176,17 @@ fixed-effects linear regression}, Stat 3(1):379-390, 2014. Guimaraes, P. and Portugal, P. (2010) \cite{A simple feasible procedure to fit models with high-dimensional fixed effects.} The Stata Journal, 10(4):629--649, 2010. -\url{http://www.stata-journal.com/article.html?article=st0212} +\url{https://www.stata-journal.com/article.html?article=st0212} Ouazad, A. (2008) \cite{A2REG: Stata module to estimate models with two fixed effects.} Statistical Software Components S456942, Boston College Department of Economics. -\url{http://ideas.repec.org/c/boc/bocode/s456942.html} +\url{https://ideas.repec.org/c/boc/bocode/s456942.html} Sanderson, E. and F. Windmeijer (2014) \cite{A weak instrument F-test in linear IV models with multiple endogenous variables}, Journal of Econometrics, 2015. -\url{http://www.sciencedirect.com/science/article/pii/S0304407615001736} +\url{https://www.sciencedirect.com/science/article/pii/S0304407615001736} } \keyword{models} \keyword{regression} From dbbcb4eaea749dae2dcd071e1857308dd8e071f6 Mon Sep 17 00:00:00 2001 From: Matthieu Stigler Date: Sun, 6 Dec 2020 18:53:34 -0800 Subject: [PATCH 03/87] Incoporate change made by CRAN maintianers --- R/startup.R | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/R/startup.R b/R/startup.R index b888d93..92421be 100644 --- a/R/startup.R +++ b/R/startup.R @@ -29,7 +29,8 @@ setoption <- function(...) { if(is.na(cr)) cr <- as.integer(Sys.getenv('OMP_NUM_THREADS')) if(is.na(cr)) cr <- as.integer(Sys.getenv('OMP_THREAD_LIMIT')) if(is.na(cr)) cr <- as.integer(Sys.getenv('NUMBER_OF_PROCESSORS')) - if(is.na(cr)) cr <- parallel::detectCores(all.tests=TRUE) + ## all.tests=TRUE is unsafe for batch use (addefd by ripley) + if(is.na(cr)) cr <- parallel::detectCores() #all.tests=TRUE) if(is.na(cr) || cr < 1) { cr <- 1 } From 59f952e8ae6e2576d19445d0b612ad80daad6477 Mon Sep 17 00:00:00 2001 From: Matthieu Stigler Date: Sun, 6 Dec 2020 20:27:31 -0800 Subject: [PATCH 04/87] fixes sprintf issue --- tests/degenerate.R | 3 ++- tests/degenerate.Rout.save | 6 ++++-- 2 files changed, 6 insertions(+), 3 deletions(-) diff --git a/tests/degenerate.R b/tests/degenerate.R index e59aed2..779f83c 100644 --- a/tests/degenerate.R +++ b/tests/degenerate.R @@ -31,6 +31,7 @@ ideff <- alpha[paste('id',id,sep='.'),'effect'] firmeff <- alpha[paste('firm',firm,sep='.'),'effect'] ## verify that id and firm coefficients are 1 -cat('accuracy:',sprintf('%.8e',coef(lm(y ~ x + x2 + x3 + x4 + x5 + x6 + ideff + firmeff-1))[7:8],'\n')) +co <- coef(lm(y ~ x + x2 + x3 + x4 + x5 + x6 + ideff + firmeff-1))[7:8] +cat('accuracy:', formatC(co, digits=8, format="f"),'\n') diff --git a/tests/degenerate.Rout.save b/tests/degenerate.Rout.save index 879b7f2..519ce89 100644 --- a/tests/degenerate.Rout.save +++ b/tests/degenerate.Rout.save @@ -79,8 +79,10 @@ F-statistic(proj model): 1.21e+03 on 6 and 1980 DF, p-value: <2e-16 > firmeff <- alpha[paste('firm',firm,sep='.'),'effect'] > > ## verify that id and firm coefficients are 1 -> cat('accuracy:',sprintf('%.8e',coef(lm(y ~ x + x2 + x3 + x4 + x5 + x6 + ideff + firmeff-1))[7:8],'\n')) -accuracy: 1.00000000e+00 1.00000000e+00> +> co <- coef(lm(y ~ x + x2 + x3 + x4 + x5 + x6 + ideff + firmeff-1))[7:8] +> cat('accuracy:', formatC(co, digits=8, format="f"),'\n') +accuracy: 1.00000000 1.00000000 +> > > > proc.time() From 1e6a3a79eee0941934b4795e644d410f3114847a Mon Sep 17 00:00:00 2001 From: Matthieu Stigler Date: Sun, 6 Dec 2020 20:28:49 -0800 Subject: [PATCH 05/87] remove examples/lfe-Ex --- tests/Examples/lfe-Ex.Rout.save | 1333 ------------------------------- 1 file changed, 1333 deletions(-) delete mode 100644 tests/Examples/lfe-Ex.Rout.save diff --git a/tests/Examples/lfe-Ex.Rout.save b/tests/Examples/lfe-Ex.Rout.save deleted file mode 100644 index 5c37589..0000000 --- a/tests/Examples/lfe-Ex.Rout.save +++ /dev/null @@ -1,1333 +0,0 @@ - -R Under development (unstable) (2019-12-10 r77548) -- "Unsuffered Consequences" -Copyright (C) 2019 The R Foundation for Statistical Computing -Platform: x86_64-pc-linux-gnu (64-bit) - -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. - - Natural language support but running in an English locale - -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. - -> pkgname <- "lfe" -> source(file.path(R.home("share"), "R", "examples-header.R")) -> options(warn = 1) -> library('lfe') -Loading required package: Matrix -> -> base::assign(".oldSearch", base::search(), pos = 'CheckExEnv') -> base::assign(".old_wd", base::getwd(), pos = 'CheckExEnv') -> cleanEx() -> nameEx("bccorr") -> ### * bccorr -> -> flush(stderr()); flush(stdout()) -> -> ### Name: bccorr -> ### Title: Compute limited mobility bias corrected correlation between -> ### fixed effects -> ### Aliases: bccorr -> -> ### ** Examples -> -> x <- rnorm(500) -> x2 <- rnorm(length(x)) -> -> ## create individual and firm -> id <- factor(sample(40,length(x),replace=TRUE)) -> firm <- factor(sample(30,length(x),replace=TRUE,prob=c(2,rep(1,29)))) -> foo <- factor(sample(20,length(x),replace=TRUE)) -> ## effects -> id.eff <- rnorm(nlevels(id)) -> firm.eff <- rnorm(nlevels(firm)) -> foo.eff <- rnorm(nlevels(foo)) -> ## left hand side -> y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] + foo.eff[foo] + rnorm(length(x)) -> -> # make a data frame -> fr <- data.frame(y,x,x2,id,firm,foo) -> ## estimate and print result -> est <- felm(y ~ x+x2|id+firm+foo, data=fr, keepX=TRUE) -> # find bias corrections -> bccorr(est) - corr v1 v2 cov d1 d2 - 0.008017737 0.934873578 1.301275975 0.008843272 0.081184357 0.058118985 - d12 --0.004131701 -> -> -> -> cleanEx() -> nameEx("btrap") -> ### * btrap -> -> flush(stderr()); flush(stdout()) -> -> ### Name: btrap -> ### Title: Bootstrap standard errors for the group fixed effects -> ### Aliases: btrap -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=2) -> ## create covariates -> x <- rnorm(3000) -> x2 <- rnorm(length(x)) -> -> ## create individual and firm -> id <- factor(sample(700,length(x),replace=TRUE)) -> firm <- factor(sample(300,length(x),replace=TRUE)) -> -> ## effects -> id.eff <- rlnorm(nlevels(id)) -> firm.eff <- rexp(nlevels(firm)) -> -> ## left hand side -> y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] + rnorm(length(x)) -> -> ## estimate and print result -> est <- felm(y ~ x+x2 | id + firm) -> summary(est) - -Call: - felm(formula = y ~ x + x2 | id + firm) - -Residuals: - Min 1Q Median 3Q Max --2.5692 -0.5356 0.0000 0.5265 3.0282 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 0.98931 0.02125 46.55 <2e-16 *** -x2 0.25710 0.02188 11.75 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 0.981 on 2006 degrees of freedom -Multiple R-squared(full model): 0.8987 Adjusted R-squared: 0.8485 -Multiple R-squared(proj model): 0.5385 Adjusted R-squared: 0.3101 -F-statistic(full model):17.92 on 993 and 2006 DF, p-value: < 2.2e-16 -F-statistic(proj model): 1170 on 2 and 2006 DF, p-value: < 2.2e-16 - - -> ## extract the group effects -> alpha <- getfe(est) -> head(alpha) - effect obs comp fe idx -id.1 0.7558736 6 1 id 1 -id.2 -0.7764543 1 1 id 2 -id.3 0.9832472 7 1 id 3 -id.4 2.7358165 4 1 id 4 -id.5 0.9211335 4 1 id 5 -id.6 3.8248845 2 1 id 6 -> ## bootstrap standard errors -> head(btrap(alpha,est)) - effect obs comp fe idx se -id.1 0.7558736 6 1 id 1 0.4565457 -id.2 -0.7764543 1 1 id 2 0.8725093 -id.3 0.9832472 7 1 id 3 0.3565344 -id.4 2.7358165 4 1 id 4 0.4586499 -id.5 0.9211335 4 1 id 5 0.4740637 -id.6 3.8248845 2 1 id 6 0.6129610 -> -> ## bootstrap some differences -> ef <- function(v,addnames) { -+ w <- c(v[2]-v[1],v[3]-v[2],v[3]-v[1]) -+ if(addnames) { -+ names(w) <-c('id2-id1','id3-id2','id3-id1') -+ attr(w,'extra') <- list(note=c('line1','line2','line3')) -+ } -+ w -+ } -> # check that it's estimable -> is.estimable(ef,est$fe) -[1] TRUE -> -> head(btrap(alpha,est,ef=ef)) - effect note se -id2-id1 -1.5323279 line1 0.8175766 -id3-id2 1.7597015 line2 0.8622066 -id3-id1 0.2273736 line3 0.4678413 -> options(oldopts) -> -> -> -> -> cleanEx() -> nameEx("cgsolve") -> ### * cgsolve -> -> flush(stderr()); flush(stdout()) -> -> ### Name: cgsolve -> ### Title: Solve a symmetric linear system with the conjugate gradient -> ### method -> ### Aliases: cgsolve -> -> ### ** Examples -> -> -> N <- 100000 -> # create some factors -> f1 <- factor(sample(34000,N,replace=TRUE)) -> f2 <- factor(sample(25000,N,replace=TRUE)) -> # a matrix of dummies, which probably is rank deficient -> B <- makeDmatrix(list(f1,f2)) -> dim(B) -[1] 100000 56753 -> # create a right hand side -> b <- as.matrix(B %*% rnorm(ncol(B))) -> # solve B' B x = B' b -> sol <- cgsolve(crossprod(B), crossprod(B, b), eps=-1e-2) -> #verify solution -> sqrt(sum((B %*% sol - b)^2)) -[1] 0.0007802951 -> -> -> -> -> cleanEx() -> nameEx("chainsubset") -> ### * chainsubset -> -> flush(stderr()); flush(stdout()) -> -> ### Name: chainsubset -> ### Title: Chain subset conditions -> ### Aliases: chainsubset -> -> ### ** Examples -> -> set.seed(48) -> N <- 10000 -> dat <- data.frame(y=rnorm(N), x=rnorm(N)) -> # It's not the same as and'ing the conditions: -> felm(y ~ x,data=dat,subset=chainsubset(x < mean(y), y < 2*mean(x))) -(Intercept) x - -2.058171 -0.002087 -> felm(y ~ x,data=dat,subset=chainsubset(y < 2*mean(x), x < mean(y))) -(Intercept) x - -0.74578 0.03168 -> felm(y ~ x,data=dat,subset=(x < mean(y)) & (y < 2*mean(x))) -(Intercept) x - -0.791573 0.007835 -> lm(y ~ x, data=dat, subset=chainsubset(x < mean(y), x > mean(y))) - -Call: -lm(formula = y ~ x, data = dat, subset = chainsubset(x < mean(y), - x > mean(y))) - -Coefficients: -(Intercept) x - 0.0586 -14.1016 - -> -> -> -> cleanEx() -> nameEx("compfactor") -> ### * compfactor -> -> flush(stderr()); flush(stdout()) -> -> ### Name: compfactor -> ### Title: Find the connected components -> ### Aliases: compfactor -> -> ### ** Examples -> -> -> ## create two factors -> f1 <- factor(sample(300,400,replace=TRUE)) -> f2 <- factor(sample(300,400,replace=TRUE)) -> -> ## find the components -> cf <- compfactor(list(f1=f1,f2=f2)) -> -> ## show the third largest component -> fr <- data.frame(f1,f2,cf) -> fr[cf==3,] - f1 f2 cf -7 277 111 3 -28 198 115 3 -52 141 65 3 -55 108 219 3 -204 108 65 3 -246 141 264 3 -278 108 115 3 -301 277 264 3 -325 141 185 3 -352 109 65 3 -396 277 63 3 -> -> -> -> -> cleanEx() -> nameEx("condfstat") -> ### * condfstat -> -> flush(stderr()); flush(stdout()) -> -> ### Name: condfstat -> ### Title: Compute conditional F statistic for weak instruments in an -> ### IV-estimation with multiple endogenous variables -> ### Aliases: condfstat -> -> ### ** Examples -> -> -> z1 <- rnorm(4000) -> z2 <- rnorm(length(z1)) -> u <- rnorm(length(z1)) -> # make x1, x2 correlated with errors u -> -> x1 <- z1 + z2 + 0.2*u + rnorm(length(z1)) -> x2 <- z1 + 0.94*z2 - 0.3*u + rnorm(length(z1)) -> y <- x1 + x2 + u -> est <- felm(y ~ 1 | 0 | (x1 | x2 ~ z1 + z2)) -> summary(est) - -Call: - felm(formula = y ~ 1 | 0 | (x1 | x2 ~ z1 + z2)) - -Residuals: - Min 1Q Median 3Q Max --4.3972 -0.7110 -0.0092 0.7191 3.5336 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -(Intercept) -0.01607 0.01962 -0.819 0.412825 -`x1(fit)` 0.90461 0.24097 3.754 0.000177 *** -`x2(fit)` 1.11122 0.25643 4.333 1.5e-05 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.047 on 3997 degrees of freedom -Multiple R-squared(full model): 0.8992 Adjusted R-squared: 0.8992 -Multiple R-squared(proj model): 0.8992 Adjusted R-squared: 0.8992 -F-statistic(full model):1.479e+04 on 2 and 3997 DF, p-value: < 2.2e-16 -F-statistic(proj model): 1.479e+04 on 2 and 3997 DF, p-value: < 2.2e-16 -F-statistic(endog. vars):1.479e+04 on 2 and 3997 DF, p-value: < 2.2e-16 - - -> ## Not run: -> ##D summary(est$stage1, lhs='x1') -> ##D summary(est$stage1, lhs='x2') -> ## End(Not run) -> -> # the joint significance of the instruments in both the first stages are ok: -> t(sapply(est$stage1$lhs, function(lh) waldtest(est$stage1, ~z1|z2, lhs=lh))) - p chi2 df1 p.F F df2 -x1 0 8416.522 2 0 4208.261 3997 -x2 0 6889.635 2 0 3444.817 3997 -> # everything above looks fine, t-tests for instruments, -> # as well as F-tests for excluded instruments in the 1st stages. -> # The conditional F-test reveals that the instruments are jointly weak -> # (it's close to being only one instrument, z1+z2, for both x1 and x2) -> condfstat(est, quantiles=c(0.05, 0.95)) - x1 x2 -iid F 8.023097 8.021403 -attr(,"df1") -[1] 1 -attr(,"quantiles") - 5% 95% -x1 0.4076302 1.238373 -x2 0.7415684 1.635039 -attr(,"quantiles")attr(,"q") -[1] 0.05 0.95 -attr(,"quantiles")attr(,"samples") -[1] 100 -> -> -> -> -> cleanEx() -> nameEx("demeanlist") -> ### * demeanlist -> -> flush(stderr()); flush(stdout()) -> -> ### Name: demeanlist -> ### Title: Centre vectors on multiple groups -> ### Aliases: demeanlist -> -> ### ** Examples -> -> oldopts <- options(lfe.threads=1) -> ## create a matrix -> mtx <- data.frame(matrix(rnorm(999),ncol=3)) -> # a list of factors -> rgb <- c('red','green','blue') -> fl <- replicate(4, factor(sample(rgb,nrow(mtx),replace=TRUE)), simplify=FALSE) -> names(fl) <- paste('g',seq_along(fl),sep='') -> # centre on all means -> mtx0 <- demeanlist(mtx,fl) -> head(data.frame(mtx0,fl)) - X1 X2 X3 g1 g2 g3 g4 -1 -0.64196076 1.4456118 -1.3025993 red red green blue -2 -0.06326532 -0.2158289 0.6045219 green green green red -3 -0.40641864 0.7223192 0.2344265 blue blue blue green -4 1.40753491 0.8459852 2.1454416 blue red red blue -5 0.09988421 -0.2905892 -1.0429431 blue blue red blue -6 -0.75980638 -0.1095204 1.1490805 blue green green green -> # verify that the group means for the columns are zero -> lapply(fl, function(f) apply(mtx0,2,tapply,f,mean)) -$g1 - X1 X2 X3 -blue 4.029092e-12 -4.113887e-13 -3.656474e-10 -green 3.004903e-12 6.395388e-14 3.751943e-10 -red -7.680146e-12 3.528239e-13 -5.681341e-11 - -$g2 - X1 X2 X3 -blue -4.816087e-18 6.232643e-18 2.922735e-17 -green -2.429030e-18 7.222454e-18 -2.173773e-17 -red 6.679396e-18 5.807147e-18 1.077004e-17 - -$g3 - X1 X2 X3 -blue -2.727461e-11 7.652167e-13 3.084120e-10 -green -3.646924e-11 -4.633995e-14 -3.716751e-10 -red 6.153743e-11 -5.825537e-13 1.427556e-10 - -$g4 - X1 X2 X3 -blue 2.109630e-13 -4.238365e-13 1.771802e-10 -green -1.727496e-13 1.232453e-13 1.812780e-11 -red -1.337501e-14 2.238489e-13 -1.548618e-10 - -> options(oldopts) -> -> -> -> cleanEx() -> nameEx("efactory") -> ### * efactory -> -> flush(stderr()); flush(stdout()) -> -> ### Name: efactory -> ### Title: Create estimable function -> ### Aliases: efactory -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=1) -> id <- factor(sample(5000,50000,replace=TRUE)) -> firm <- factor(sample(3000,50000,replace=TRUE)) -> fl <- list(id=id,firm=firm) -> obj <- list(fe=fl,cfactor=compfactor(fl)) -> ## the trivial least-norm transformtion, which by the way is non-estimable -> print(ef <- efactory(obj,'ln')) -function (v, addnames) -{ - if (addnames) { - names(v) <- nm - attr(v, "extra") <- list(obs = obs, comp = comp, fe = fef, - idx = idx) - } - v -} - - -> is.estimable(ef,fl) -Warning in is.estimable(ef, fl) : - non-estimable function, largest error 0.003 in coordinate 1050 ("id.1050") -[1] FALSE -> ## then the default -> print(ef <- efactory(obj,'ref')) -function (v, addnames) -{ - esum <- sum(v[extrarefs]) - df <- v[refsubs] - sub <- ifelse(is.na(df), 0, df) - df <- v[refsuba] - add <- ifelse(is.na(df), 0, df + esum) - v <- v - sub + add - if (addnames) { - names(v) <- nm - attr(v, "extra") <- list(obs = obs, comp = comp, fe = fef, - idx = idx) - } - v -} - - -attr(,"verified") -[1] TRUE -> is.estimable(ef,fl) -[1] TRUE -> # get the names of the coefficients, i.e. the nm-variable in the function -> head(evalq(nm,environment(ef))) -[1] "id.1" "id.2" "id.3" "id.4" "id.5" "id.6" -> options(oldopts) -> -> -> -> -> cleanEx() -> nameEx("felm") -> ### * felm -> -> flush(stderr()); flush(stdout()) -> -> ### Name: felm -> ### Title: Fit a linear model with multiple group fixed effects -> ### Aliases: felm -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=1) -> -> ## Simulate data -> -> # Covariates -> x <- rnorm(1000) -> x2 <- rnorm(length(x)) -> # Individuals and firms -> id <- factor(sample(20,length(x),replace=TRUE)) -> firm <- factor(sample(13,length(x),replace=TRUE)) -> # Effects for them -> id.eff <- rnorm(nlevels(id)) -> firm.eff <- rnorm(nlevels(firm)) -> # Left hand side -> u <- rnorm(length(x)) -> y <- x + 0.5*x2 + id.eff[id] + firm.eff[firm] + u -> -> ## Estimate the model and print the results -> est <- felm(y ~ x + x2 | id + firm) -> summary(est) - -Call: - felm(formula = y ~ x + x2 | id + firm) - -Residuals: - Min 1Q Median 3Q Max --3.4329 -0.6925 -0.0261 0.6946 2.7795 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 0.96790 0.03120 31.02 <2e-16 *** -x2 0.48603 0.03123 15.56 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.009 on 966 degrees of freedom -Multiple R-squared(full model): 0.7456 Adjusted R-squared: 0.7369 -Multiple R-squared(proj model): 0.5549 Adjusted R-squared: 0.5397 -F-statistic(full model):85.78 on 33 and 966 DF, p-value: < 2.2e-16 -F-statistic(proj model): 602.1 on 2 and 966 DF, p-value: < 2.2e-16 - - -> -> ## Not run: -> ##D # Compare with lm -> ##D summary(lm(y ~ x + x2 + id + firm-1)) -> ## End(Not run) -> -> ## Example with 'reverse causation' (IV regression) -> -> # Q and W are instrumented by x3 and the factor x4. -> x3 <- rnorm(length(x)) -> x4 <- sample(12,length(x),replace=TRUE) -> Q <- 0.3*x3 + x + 0.2*x2 + id.eff[id] + 0.3*log(x4) - 0.3*y + rnorm(length(x),sd=0.3) -> W <- 0.7*x3 - 2*x + 0.1*x2 - 0.7*id.eff[id] + 0.8*cos(x4) - 0.2*y+ rnorm(length(x),sd=0.6) -> # Add them to the outcome variable -> y <- y + Q + W -> -> ## Estimate the IV model and report robust SEs -> ivest <- felm(y ~ x + x2 | id + firm | (Q|W ~ x3 + factor(x4))) -> summary(ivest, robust=TRUE) - -Call: - felm(formula = y ~ x + x2 | id + firm | (Q | W ~ x3 + factor(x4))) - -Residuals: - Min 1Q Median 3Q Max --3.4927 -0.6913 -0.0241 0.6976 2.7567 - -Coefficients: - Estimate Robust s.e t value Pr(>|t|) -x 0.89435 0.16419 5.447 6.5e-08 *** -x2 0.48430 0.03160 15.326 < 2e-16 *** -`Q(fit)` 1.03746 0.10833 9.577 < 2e-16 *** -`W(fit)` 0.97835 0.04679 20.908 < 2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.017 on 964 degrees of freedom -Multiple R-squared(full model): 0.6976 Adjusted R-squared: 0.6867 -Multiple R-squared(proj model): 0.6043 Adjusted R-squared: 0.59 -F-statistic(full model, *iid*): 71.4 on 35 and 964 DF, p-value: < 2.2e-16 -F-statistic(proj model): 460.9 on 4 and 964 DF, p-value: < 2.2e-16 -F-statistic(endog. vars):681.6 on 2 and 964 DF, p-value: < 2.2e-16 - - -> condfstat(ivest) - Q W -iid F 40.07999 66.90775 -attr(,"df1") -[1] 11 -> -> ## Not run: -> ##D # Compare with the not instrumented fit: -> ##D summary(felm(y ~ x + x2 + Q + W | id + firm)) -> ## End(Not run) -> -> ## Example with multiway clustering -> -> # Create a large cluster group (500 clusters) and a small one (20 clusters) -> cl1 <- factor(sample(rep(1:500, length.out=length(x)))) -> cl2 <- factor(sample(rep(1:20, length.out=length(x)))) -> # Function for adding clustered noise to our outcome variable -> cl_noise <- function(cl) { -+ obs_per_cluster <- length(x)/nlevels(cl) -+ unlist(replicate(nlevels(cl), rnorm(obs_per_cluster, mean=rnorm(1), sd=runif(1)), simplify=FALSE)) -+ } -> # New outcome variable -> y_cl <- x + 0.5*x2 + id.eff[id] + firm.eff[firm] + cl_noise(cl1) + cl_noise(cl2) -> -> ## Estimate and print the model with cluster-robust SEs (default) -> est_cl <- felm(y_cl ~ x + x2 | id + firm | 0 | cl1 + cl2) -> summary(est_cl) - -Call: - felm(formula = y_cl ~ x + x2 | id + firm | 0 | cl1 + cl2) - -Residuals: - Min 1Q Median 3Q Max --4.8756 -0.9697 0.0232 1.0041 5.2598 - -Coefficients: - Estimate Cluster s.e. t value Pr(>|t|) -x 0.97330 0.05007 19.44 <2e-16 *** -x2 0.53359 0.03043 17.54 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.528 on 966 degrees of freedom -Multiple R-squared(full model): 0.5746 Adjusted R-squared: 0.5601 -Multiple R-squared(proj model): 0.3633 Adjusted R-squared: 0.3415 -F-statistic(full model, *iid*):39.54 on 33 and 966 DF, p-value: < 2.2e-16 -F-statistic(proj model): 252.8 on 2 and 499 DF, p-value: < 2.2e-16 - - -> -> ## Not run: -> ##D # Print ordinary standard errors: -> ##D summary(est_cl, robust = FALSE) -> ##D # Match cluster-robust SEs from Stata's reghdfe package: -> ##D summary(felm(y_cl ~ x + x2 | id + firm | 0 | cl1 + cl2, cmethod="reghdfe")) -> ## End(Not run) -> -> options(oldopts) -> -> -> -> -> cleanEx() -> nameEx("fevcov") -> ### * fevcov -> -> flush(stderr()); flush(stdout()) -> -> ### Name: fevcov -> ### Title: Compute limited mobility bias corrected covariance matrix -> ### between fixed effects -> ### Aliases: fevcov -> -> ### ** Examples -> -> -> x <- rnorm(5000) -> x2 <- rnorm(length(x)) -> -> ## create individual and firm -> id <- factor(sample(40,length(x),replace=TRUE)) -> firm <- factor(sample(30,length(x),replace=TRUE,prob=c(2,rep(1,29)))) -> foo <- factor(sample(20,length(x),replace=TRUE)) -> ## effects -> id.eff <- rnorm(nlevels(id)) -> firm.eff <- runif(nlevels(firm)) -> foo.eff <- rchisq(nlevels(foo),df=1) -> ## left hand side -> id.m <- id.eff[id] -> firm.m <- firm.eff[firm] -> foo.m <- foo.eff[foo] -> # normalize them -> id.m <- id.m/sd(id.m) -> firm.m <- firm.m/sd(firm.m) -> foo.m <- foo.m/sd(foo.m) -> y <- x + 0.25*x2 + id.m + firm.m + foo.m + rnorm(length(x),sd=2) -> z <- x + 0.5*x2 + 0.7*id.m + 0.5*firm.m + 0.3*foo.m + rnorm(length(x),sd=2) -> # make a data frame -> fr <- data.frame(y,z,x,x2,id,firm,foo) -> ## estimate and print result -> est <- felm(y|z ~ x+x2|id+firm+foo, data=fr, keepX=TRUE) -> # find bias corrections, there's little bias in this example -> print(yv <- fevcov(est, lhs='y')) - id firm foo -id 1.031372636 0.006691582 -0.006322827 -firm 0.006691582 1.049086528 0.010535564 -foo -0.006322827 0.010535564 1.057014174 -attr(,"bias") - id firm foo -id 3.054932e-02 -0.0001717782 -9.118383e-05 -firm -1.717782e-04 0.0194414181 -1.092769e-04 -foo -9.118383e-05 -0.0001092769 1.442915e-02 -> ## Here's how to compute the unbiased correlation matrix: -> cm <- cov2cor(yv) -> structure(cm,bias=NULL) - id firm foo -id 1.000000000 0.006433027 -0.006055682 -firm 0.006433027 1.000000000 0.010004876 -foo -0.006055682 0.010004876 1.000000000 -> -> -> -> -> cleanEx() -> nameEx("fixedse") -> ### * fixedse -> -> flush(stderr()); flush(stdout()) -> -> ### Name: fixedse -> ### Title: Compute standard errors for fixed effects -> ### Aliases: fixedse -> ### Keywords: internal -> -> ### ** Examples -> -> x <- rnorm(1000) -> f <- factor(sample(5,1000,replace=TRUE)) -> y <- x + (1:5)[f] + rnorm(1000) -> est <- felm(y ~ x | f, keepX=TRUE) -> #both bootstrap and computed se: -> cbind(getfe(est,ef=efactory(est,'ref'),se=TRUE), fse=fixedse(est)) - effect obs comp fe idx se fse -f.1 0.9569311 185 1 f 1 0.07749880 0.07535548 -f.2 1.9949066 223 1 f 2 0.06038432 0.06860181 -f.3 3.0326896 185 1 f 3 0.07696251 0.07532290 -f.4 4.0294033 196 1 f 4 0.07527526 0.07317074 -f.5 4.9413841 211 1 f 5 0.07015206 0.07054858 -> #compare with lm: -> summary(lm(y ~x+f-1)) - -Call: -lm(formula = y ~ x + f - 1) - -Residuals: - Min 1Q Median 3Q Max --3.5561 -0.6437 -0.0299 0.7166 2.9214 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 1.01522 0.03135 32.38 <2e-16 *** -f1 0.95693 0.07536 12.70 <2e-16 *** -f2 1.99491 0.06860 29.08 <2e-16 *** -f3 3.03269 0.07532 40.26 <2e-16 *** -f4 4.02940 0.07317 55.07 <2e-16 *** -f5 4.94138 0.07055 70.04 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.024 on 994 degrees of freedom -Multiple R-squared: 0.9216, Adjusted R-squared: 0.9211 -F-statistic: 1947 on 6 and 994 DF, p-value: < 2.2e-16 - -> -> -> -> cleanEx() -> nameEx("getfe") -> ### * getfe -> -> flush(stderr()); flush(stdout()) -> -> ### Name: getfe -> ### Title: Retrieve the group fixed effects -> ### Aliases: getfe -> ### Keywords: models regression -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=2) -> ## create covariates -> x <- rnorm(4000) -> x2 <- rnorm(length(x)) -> -> ## create individual and firm -> id <- factor(sample(500,length(x),replace=TRUE)) -> firm <- factor(sample(300,length(x),replace=TRUE)) -> -> ## effects -> id.eff <- rlnorm(nlevels(id)) -> firm.eff <- rexp(nlevels(firm)) -> -> ## left hand side -> y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] + rnorm(length(x)) -> -> ## estimate and print result -> est <- felm(y ~ x+x2 | id + firm) -> summary(est) - -Call: - felm(formula = y ~ x + x2 | id + firm) - -Residuals: - Min 1Q Median 3Q Max --2.94796 -0.59482 -0.01038 0.59008 3.15936 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 1.01409 0.01687 60.11 <2e-16 *** -x2 0.26766 0.01752 15.28 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 0.9892 on 3200 degrees of freedom -Multiple R-squared(full model): 0.9149 Adjusted R-squared: 0.8937 -Multiple R-squared(proj model): 0.5457 Adjusted R-squared: 0.4323 -F-statistic(full model):43.06 on 799 and 3200 DF, p-value: < 2.2e-16 -F-statistic(proj model): 1922 on 2 and 3200 DF, p-value: < 2.2e-16 - - -> ## extract the group effects -> alpha <- getfe(est,se=TRUE) -> -> ## find some estimable functions, with standard errors, we don't get -> ## names so we must precompute some numerical indices in ef -> idx <- match(c('id.5','id.6','firm.11','firm.12'),rownames(alpha)) -> alpha[idx,] - effect obs comp fe idx se -id.5 2.8685887 6 1 id 5 0.4564400 -id.6 2.1252963 7 1 id 6 0.4294043 -firm.11 0.2437394 18 1 firm 11 0.2809155 -firm.12 0.2925972 14 1 firm 12 0.3584084 -> ef <- function(v,addnames) { -+ w <- c(v[idx[[2]]]-v[idx[[1]]],v[idx[[4]]]+v[idx[[1]]], -+ v[idx[[4]]]-v[idx[[3]]]) -+ if(addnames) names(w) <-c('id6-id5','f12+id5','f12-f11') -+ w -+ } -> getfe(est,ef=ef,se=TRUE) - effect se -id6-id5 -0.74329244 0.5178066 -f12+id5 3.16118596 0.4907870 -f12-f11 0.04885782 0.3275784 -> options(oldopts) -> ## Not run: -> ##D summary(lm(y ~ x+x2+id+firm-1)) -> ## End(Not run) -> -> -> -> -> cleanEx() -> nameEx("is.estimable") -> ### * is.estimable -> -> flush(stderr()); flush(stdout()) -> -> ### Name: is.estimable -> ### Title: Verify estimability of function -> ### Aliases: is.estimable -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=1) -> ## create individual and firm -> id <- factor(sample(5000,50000,replace=TRUE)) -> firm <- factor(sample(3000,50000,replace=TRUE)) -> -> ## create some estimable functions. It's faster to -> ## use numerical indices in ef rather than strings, and the input v -> ## to ef has no names, we have to add them when requested -> ef <- function(v,addnames) { -+ w <- c(v[6]-v[5],v[7000]+v[5],v[7000]-v[6000]) -+ if(addnames) names(w) <-c('id6-id5','f2k+id5','f2k-f1k') -+ w -+ } -> is.estimable(ef,list(id=id,firm=firm)) -[1] TRUE -> -> ## Then make an error; in the last coordinate, sum two firms -> ef <- function(v,addnames) { -+ w <- c(v[6]-v[5],v[7000]+v[5],v[7000]+v[6000]) -+ if(addnames) names(w) <-c('id6-id5','f2k+id5','f2k-f1k') -+ w -+ } -> is.estimable(ef, list(id=id,firm=firm), keepdiff=TRUE) -Warning in is.estimable(ef, list(id = id, firm = firm), keepdiff = TRUE) : - non-estimable function, largest error 0.01 in coordinate 3 ("f2k-f1k") -[1] FALSE -attr(,"diff") - id6-id5 f2k+id5 f2k-f1k -2.797201e-11 1.212475e-11 1.436988e-02 -> options(oldopts) -> -> -> -> -> cleanEx() -> nameEx("kaczmarz") -> ### * kaczmarz -> -> flush(stderr()); flush(stdout()) -> -> ### Name: kaczmarz -> ### Title: Solve a linear system defined by factors -> ### Aliases: kaczmarz -> -> ### ** Examples -> -> -> ## create factors -> f1 <- factor(sample(24000,100000,replace=TRUE)) -> f2 <- factor(sample(20000,length(f1),replace=TRUE)) -> f3 <- factor(sample(10000,length(f1),replace=TRUE)) -> f4 <- factor(sample(8000,length(f1),replace=TRUE)) -> ## the matrix of dummies -> D <- makeDmatrix(list(f1,f2,f3,f4)) -> dim(D) -[1] 100000 61499 -> ## an x -> truex <- runif(ncol(D)) -> ## and the right hand side -> R <- as.vector(D %*% truex) -> ## solve it -> sol <- kaczmarz(list(f1,f2,f3,f4),R) -> ## verify that the solution solves the system Dx = R -> sqrt(sum((D %*% sol - R)^2)) -[1] 3.600146e-07 -> ## but the solution is not equal to the true x, because the system is -> ## underdetermined -> sqrt(sum((sol - truex)^2)) -[1] 52.45481 -> ## moreover, the solution from kaczmarz has smaller norm -> sqrt(sum(sol^2)) < sqrt(sum(truex^2)) -[1] TRUE -> -> -> -> -> cleanEx() -> nameEx("lfe-package") -> ### * lfe-package -> -> flush(stderr()); flush(stdout()) -> -> ### Name: lfe-package -> ### Title: Overview. Linear Group Fixed Effects -> ### Aliases: lfe-package lfe -> ### Keywords: models regression -> -> ### ** Examples -> -> -> oldopts <- options(lfe.threads=1) -> x <- rnorm(1000) -> x2 <- rnorm(length(x)) -> id <- factor(sample(10,length(x),replace=TRUE)) -> firm <- factor(sample(3,length(x),replace=TRUE,prob=c(2,1.5,1))) -> year <- factor(sample(10,length(x),replace=TRUE,prob=c(2,1.5,rep(1,8)))) -> id.eff <- rnorm(nlevels(id)) -> firm.eff <- rnorm(nlevels(firm)) -> year.eff <- rnorm(nlevels(year)) -> y <- x + 0.25*x2 + id.eff[id] + firm.eff[firm] + -+ year.eff[year] + rnorm(length(x)) -> est <- felm(y ~ x+x2 | id + firm + year) -> summary(est) - -Call: - felm(formula = y ~ x + x2 | id + firm + year) - -Residuals: - Min 1Q Median 3Q Max --3.4685 -0.6450 -0.0293 0.6631 2.8981 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 0.91611 0.03103 29.520 < 2e-16 *** -x2 0.24114 0.03095 7.792 1.68e-14 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.008 on 977 degrees of freedom -Multiple R-squared(full model): 0.7621 Adjusted R-squared: 0.7567 -Multiple R-squared(proj model): 0.4893 Adjusted R-squared: 0.4778 -F-statistic(full model):142.3 on 22 and 977 DF, p-value: < 2.2e-16 -F-statistic(proj model): 468.1 on 2 and 977 DF, p-value: < 2.2e-16 -*** Standard errors may be too high due to more than 2 groups and exactDOF=FALSE - - -> -> getfe(est,se=TRUE) - effect obs comp fe idx se -id.1 -0.25923494 88 1 id 1 0.13632132 -id.2 1.06062345 90 1 id 2 0.14546483 -id.3 0.06034638 90 1 id 3 0.13201034 -id.4 0.44271187 109 1 id 4 0.13256707 -id.5 0.51696693 98 1 id 5 0.13117504 -id.6 1.43021599 108 1 id 6 0.14482875 -id.7 0.23001980 100 1 id 7 0.13899292 -id.8 1.07498270 112 1 id 8 0.13169287 -id.9 -0.38411917 97 1 id 9 0.14582033 -id.10 1.62621598 108 1 id 10 0.12815771 -firm.1 0.00000000 415 1 firm 1 0.00000000 -firm.2 -0.15451360 359 1 firm 2 0.07219465 -firm.3 2.59917570 226 1 firm 3 0.08069538 -year.1 0.00000000 166 2 year 1 0.00000000 -year.2 -1.90629543 144 2 year 2 0.11552284 -year.3 -0.72206632 90 2 year 3 0.13077982 -year.4 -1.34267982 89 2 year 4 0.13507944 -year.5 -1.84370897 87 2 year 5 0.13840235 -year.6 0.39746394 81 2 year 6 0.13503649 -year.7 -0.62958841 82 2 year 7 0.15321151 -year.8 -0.33821454 72 2 year 8 0.15098197 -year.9 -1.77057443 90 2 year 9 0.12325007 -year.10 0.15294690 99 2 year 10 0.12299566 -> # compare with an ordinary lm -> summary(lm(y ~ x+x2+id+firm+year-1)) - -Call: -lm(formula = y ~ x + x2 + id + firm + year - 1) - -Residuals: - Min 1Q Median 3Q Max --3.4685 -0.6450 -0.0293 0.6631 2.8981 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -x 0.91611 0.03103 29.520 < 2e-16 *** -x2 0.24114 0.03095 7.792 1.68e-14 *** -id1 -0.25923 0.13589 -1.908 0.056720 . -id2 1.06062 0.13104 8.094 1.71e-15 *** -id3 0.06035 0.13485 0.448 0.654599 -id4 0.44271 0.12720 3.480 0.000523 *** -id5 0.51697 0.13192 3.919 9.52e-05 *** -id6 1.43022 0.12637 11.318 < 2e-16 *** -id7 0.23002 0.13224 1.739 0.082278 . -id8 1.07498 0.12360 8.697 < 2e-16 *** -id9 -0.38412 0.13391 -2.869 0.004213 ** -id10 1.62622 0.12621 12.885 < 2e-16 *** -firm2 -0.15451 0.07361 -2.099 0.036071 * -firm3 2.59918 0.08399 30.945 < 2e-16 *** -year2 -1.90630 0.11514 -16.557 < 2e-16 *** -year3 -0.72207 0.13338 -5.414 7.78e-08 *** -year4 -1.34268 0.13293 -10.101 < 2e-16 *** -year5 -1.84371 0.13476 -13.682 < 2e-16 *** -year6 0.39746 0.13803 2.880 0.004070 ** -year7 -0.62959 0.13669 -4.606 4.65e-06 *** -year8 -0.33821 0.14309 -2.364 0.018290 * -year9 -1.77057 0.13274 -13.339 < 2e-16 *** -year10 0.15295 0.12844 1.191 0.234006 ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.008 on 977 degrees of freedom -Multiple R-squared: 0.7679, Adjusted R-squared: 0.7625 -F-statistic: 140.6 on 23 and 977 DF, p-value: < 2.2e-16 - -> options(oldopts) -> -> -> -> -> cleanEx() -> nameEx("makeDmatrix") -> ### * makeDmatrix -> -> flush(stderr()); flush(stdout()) -> -> ### Name: makeDmatrix -> ### Title: Make sparse matrix of dummies from factor list -> ### Aliases: makeDmatrix -> -> ### ** Examples -> -> -> fl <- lapply(1:3, function(i) factor(sample(3,10,replace=TRUE))) -> fl -[[1]] - [1] 1 3 1 2 1 3 3 2 2 3 -Levels: 1 2 3 - -[[2]] - [1] 3 1 1 1 2 2 2 2 3 1 -Levels: 1 2 3 - -[[3]] - [1] 3 1 1 1 1 2 1 1 2 2 -Levels: 1 2 3 - -> makeDmatrix(fl, weights=seq(0.1,1,0.1)) -10 x 9 sparse Matrix of class "dgCMatrix" - f1.1 f1.2 f1.3 f2.1 f2.2 f2.3 f3.1 f3.2 f3.3 - [1,] 0.1 . . . . 0.1 . . 0.1 - [2,] . . 0.2 0.2 . . 0.2 . . - [3,] 0.3 . . 0.3 . . 0.3 . . - [4,] . 0.4 . 0.4 . . 0.4 . . - [5,] 0.5 . . . 0.5 . 0.5 . . - [6,] . . 0.6 . 0.6 . . 0.6 . - [7,] . . 0.7 . 0.7 . 0.7 . . - [8,] . 0.8 . . 0.8 . 0.8 . . - [9,] . 0.9 . . . 0.9 . 0.9 . -[10,] . . 1.0 1.0 . . . 1.0 . -> -> -> -> -> cleanEx() -> nameEx("mctrace") -> ### * mctrace -> -> flush(stderr()); flush(stdout()) -> -> ### Name: mctrace -> ### Title: Compute trace of a large matrix by sample means -> ### Aliases: mctrace -> -> ### ** Examples -> -> -> A <- matrix(rnorm(25),5) -> fun <- function(x) A %*% x -> sum(diag(A)) -[1] 0.6807816 -> sum(eigen(A,only.values=TRUE)$values) -[1] 0.6807816+0i -> # mctrace is not really useful for small problems. -> mctrace(fun,ncol(A),tol=0.05) -[1] 0.6600197 -attr(,"sd") -[1] 0.03292458 -attr(,"iterations") -[1] 11040 -> # try a larger problem (3000x3000): -> f1 <- factor(sample(1500,3000,replace=TRUE)) -> f2 <- factor(sample(1500,3000,replace=TRUE)) -> fl <- list(f1,f2) -> mctrace(fl,tol=-5) -[1] 499.8301 -attr(,"sd") -[1] 3.43084 -attr(,"iterations") -[1] 56 -> # exact: -> length(f1) - nlevels(f1) - nlevels(f2) + nlevels(compfactor(fl)) -[1] 495 -> -> -> -> -> cleanEx() -> nameEx("nlexpect") -> ### * nlexpect -> -> flush(stderr()); flush(stdout()) -> -> ### Name: nlexpect -> ### Title: Compute expectation of a function of the coefficients. -> ### Aliases: nlexpect -> -> ### ** Examples -> -> -> N <- 100 -> x1 <- rnorm(N) -> # make some correlation -> x2 <- 0.1*rnorm(N) + 0.1*x1 -> y <- 0.1*x1 + x2 + rnorm(N) -> summary(est <- felm(y ~ x1 + x2)) - -Call: - felm(formula = y ~ x1 + x2) - -Residuals: - Min 1Q Median 3Q Max --2.94359 -0.43645 0.00202 0.63692 2.63941 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -(Intercept) 0.02535 0.10519 0.241 0.810 -x1 0.17458 0.15997 1.091 0.278 -x2 0.46533 1.09479 0.425 0.672 - -Residual standard error: 1.043 on 97 degrees of freedom -Multiple R-squared(full model): 0.03737 Adjusted R-squared: 0.01752 -Multiple R-squared(proj model): 0.03737 Adjusted R-squared: 0.01752 -F-statistic(full model):1.883 on 2 and 97 DF, p-value: 0.1577 -F-statistic(proj model): 1.883 on 2 and 97 DF, p-value: 0.1577 - - -> pt1 <- coef(est)['x1'] -> pt2 <- coef(est)['x2'] -> # expected values of coefficients, should match the summary -> # and variance, i.e. square of standard errors in the summary -> nlexpect(est, quote(c(x1=x1,x2=x2,var=c((x1-pt1)^2,(x2-pt2)^2)))) - x1 x2 var.x1 var.x2 -0.17450420 0.46517350 0.02611881 1.22460757 -> #Wald test of single variable -> waldtest(est, ~x1)['p.F'] - p.F -0.2778482 -> # the same with nlexpect, i.e. probability for observing abs(x1)>abs(pt1) conditional -> # on E(x1) = 0. -> nlexpect(est, (x1-pt1)^2 > pt1^2, tol=1e-7, vectorize=TRUE) - x1 -0.2778482 -> # which is the same as -> 2*nlexpect(est, x1*sign(pt1) < 0) - x1 -0.2779329 -> -> # Here's a multivalued, vectorized example -> nlexpect(est, rbind(a=x1*x2 < pt1, b=x1*x2 > 0), vectorize=TRUE, method='divonne') - x1 -a 0.8587288 -b 0.5324731 -> -> -> -> -> cleanEx() -> nameEx("varvars") -> ### * varvars -> -> flush(stderr()); flush(stdout()) -> -> ### Name: varvars -> ### Title: Compute the variance of the fixed effect variance estimate -> ### Aliases: varvars -> -> ### ** Examples -> -> -> x <- rnorm(500) -> x2 <- rnorm(length(x)) -> -> ## create individual and firm -> id <- factor(sample(40,length(x),replace=TRUE)) -> firm <- factor(sample(30,length(x),replace=TRUE,prob=c(2,rep(1,29)))) -> foo <- factor(sample(20,length(x),replace=TRUE)) -> ## effects -> id.eff <- rnorm(nlevels(id)) -> firm.eff <- rnorm(nlevels(firm)) -> foo.eff <- rnorm(nlevels(foo)) -> ## left hand side -> id.m <- id.eff[id] -> firm.m <- 2*firm.eff[firm] -> foo.m <- 3*foo.eff[foo] -> y <- x + 0.25*x2 + id.m + firm.m + foo.m + rnorm(length(x)) -> -> # make a data frame -> fr <- data.frame(y,x,x2,id,firm,foo) -> ## estimate and print result -> est <- felm(y ~ x+x2|id+firm+foo, data=fr, keepX=TRUE) -> alpha <- getfe(est) -> # estimate the covariance matrix of the fixed effects -> fevcov(est, alpha) - id firm foo -id 0.93487358 0.01198024 -0.11085834 -firm 0.01198024 4.95403947 0.04691949 -foo -0.11085834 0.04691949 6.49011310 -attr(,"bias") - id firm foo -id 0.081184357 -0.004065798 -0.003258662 -firm -0.004065798 0.056924931 -0.003863477 -foo -0.003258662 -0.003863477 0.041169982 -> # estimate variances of the diagonal -> varvars(est, alpha) -[1] 0.008598521 0.043820907 0.057656021 -> -> -> -> -> cleanEx() -> nameEx("waldtest") -> ### * waldtest -> -> flush(stderr()); flush(stdout()) -> -> ### Name: waldtest -> ### Title: Compute Wald test for joint restrictions on coefficients -> ### Aliases: waldtest -> -> ### ** Examples -> -> -> x <- rnorm(10000) -> x2 <- rnorm(length(x)) -> y <- x - 0.2*x2 + rnorm(length(x)) -> #Also works for lm -> summary(est <- lm(y ~ x + x2 )) - -Call: -lm(formula = y ~ x + x2) - -Residuals: - Min 1Q Median 3Q Max --3.6477 -0.6843 0.0158 0.6744 3.6442 - -Coefficients: - Estimate Std. Error t value Pr(>|t|) -(Intercept) 0.007554 0.010078 0.75 0.454 -x 0.999265 0.009956 100.37 <2e-16 *** -x2 -0.199066 0.010172 -19.57 <2e-16 *** ---- -Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 - -Residual standard error: 1.008 on 9997 degrees of freedom -Multiple R-squared: 0.5108, Adjusted R-squared: 0.5107 -F-statistic: 5219 on 2 and 9997 DF, p-value: < 2.2e-16 - -> # We do not reject the true values -> waldtest(est, ~ x-1|x2+0.2|`(Intercept)`) - p chi2 df1 p.F F df2 - 0.9019639 0.5757362 3.0000000 0.9019611 0.1919121 9997.0000000 -attr(,"formula") -~x - 1 | x2 + 0.2 | `(Intercept)` - -> # The Delta-method coincides when the function is linear: -> waldtest(est, function(x) x - c(0, 1, -0.2)) - p chi2 df1 p.F F df2 - 0.9019639 0.5757362 3.0000000 0.9019611 0.1919121 9997.0000000 -attr(,"formula") -function (x) -x - c(0, 1, -0.2) - -> -> -> -> -> ### *