diff --git a/NAMESPACE b/NAMESPACE index f8320765..738602f4 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -4,6 +4,7 @@ S3method(boottest,felm) S3method(boottest,fixest) S3method(boottest,ivreg) S3method(boottest,lm) +S3method(boottest,plm) S3method(confint,boottest) S3method(glance,boottest) S3method(glance,mboottest) diff --git a/R/arg_checks.R b/R/arg_checks.R index e82c5d5c..a63e7e8d 100644 --- a/R/arg_checks.R +++ b/R/arg_checks.R @@ -131,6 +131,30 @@ check_boottest_args_plus <- function( } } + + if(inherits(object, "plm")){ + + if (!is.null(fe)) { + if (fe %in% param) { + rlang::abort(paste("The function argument fe =", fe, "is included in the + hypothesis (via the `param` argument). This is not allowed. + Please set fe to another factor variable or NULL."), + call. = FALSE + ) + } + if (!(fe %in% names(index(object)))) { + rlang::abort(paste( + "The fixed effect to be projected out in the bootstrap,", + fe, "is not included as a dedicated fixed effect + in the estimated model." + )) + } + } + + + } + + if (((1 - sign_level) * (B + 1)) %% 1 != 0) { rlang::inform( paste("Note: The bootstrap usually performs best when the confidence", diff --git a/R/boottest_fixest.R b/R/boottest_fixest.R index 8fb2289c..742a7adc 100644 --- a/R/boottest_fixest.R +++ b/R/boottest_fixest.R @@ -485,6 +485,7 @@ boottest.fixest <- function(object, type = type, engine = engine ) + full_enumeration <- enumerate$full_enumeration B <- enumerate$B diff --git a/R/boottest_plm.r b/R/boottest_plm.r new file mode 100644 index 00000000..afd74831 --- /dev/null +++ b/R/boottest_plm.r @@ -0,0 +1,641 @@ +#' Fast wild cluster bootstrap inference for object of class plm +#' +#' `boottest.plm` is a S3 method that allows for fast wild cluster +#' bootstrap inference for objects of class plm by implementing +#' fast wild bootstrap algorithms as developed in Roodman et al., 2019 +#' and MacKinnon, Nielsen & Webb (2022). +#' +#' @param object An object of class plm and estimated via `plm::plm()`. +#' Non-linear models are not supported. +#' @param clustid A character vector or rhs formula containing the names of the +#' cluster variables. If `NULL`, a heteroskedasticity-robust (HC1) wild bootstrap is run. +#' @param param A character vector or rhs formula. The name of the regression +#' coefficient(s) for which the hypothesis is to be tested +#' @param B Integer. The number of bootstrap iterations. When the number of +#' clusters is low, increasing `B` adds little additional runtime. +#' @param bootcluster A character vector or rhs formula of length 1. Specifies +#' the bootstrap clustering variable or variables. If more +#' than one variable is specified, then bootstrapping is clustered by the +#' intersections of +#' clustering implied by the listed variables. To mimic the behavior of +#' Stata's boottest command, +#' the default is to cluster by the intersection of all the variables +#' specified via the `clustid` argument, +#' even though that is not necessarily recommended (see the paper by +#' Roodman et al. (2019) cited below, section 4.2). +#' Other options include `"min"`, where bootstrapping is clustered by +#' the cluster variable with the fewest clusters. +#' Further, the subcluster bootstrap (MacKinnon & Webb, 2018) is +#' supported - see the `vignette("fwildclusterboot", package = +#' "fwildclusterboot")` for details. +#' @param fe A character vector or rhs formula of length one which contains +#' the name of the fixed effect to be projected +#' out in the bootstrap. Note: if regression weights are used, `fe` +#' needs to be `NULL`. +#' @param sign_level A numeric between 0 and 1 which sets the significance level +#' of the inference procedure. E.g. `sign_level = 0.05` +#' returns 0.95% confidence intervals. By default, `sign_level = 0.05`. +#' @param conf_int A logical vector. If `TRUE`, boottest computes confidence +#' intervals by test inversion. If `FALSE`, only the p-value is returned. +#' @param engine Character scalar. Either `"R"`, `"R-lean"` or `"WildBootTests.jl"`. +#' Controls if `boottest()` should run via its native R implementation +#' or `WildBootTests.jl`. +#' `"R"` is the default and implements the cluster bootstrap +#' as in Roodman (2019). `"WildBootTests.jl"` executes the +#' wild cluster bootstrap via the WildBootTests.jl +#' package. For it to run, Julia and WildBootTests.jl need +#' to be installed. +#' The `"R-lean"` algorithm is a memory friendly, but less +#' performant rcpp-armadillo-based implementation of the wild +#' cluster bootstrap. +#' Note that if no cluster is provided, `boottest()` always +#' defaults to the `"R-lean"` algorithm. You can set the employed +#' algorithm globally by using the +#' `setBoottest_engine()` function. +#' @param bootstrap_type Determines which wild cluster bootstrap type should be +#' run. Options are `"fnw11"`, `"11"`, `"13"`, `"31"` and `"33"` for the wild cluster +#' bootstrap and `"11"` and `"31"` for the heteroskedastic bootstrap. +#' For more information, see the details section. `"fnw11"` is the default for +#' the cluster bootstrap, which runs a `"11"` type +#' wild cluster bootstrap via the algorithm outlined in "fast and wild" +#' (Roodman et al. (2019)). `"11"` is the default for the heteroskedastic +#' bootstrap. +#' @param R Hypothesis Vector giving linear combinations of coefficients. +#' Must be either `NULL` or a vector of the same length as `param`. If `NULL`, +#' a vector of ones of length param. +#' @param r A numeric. Shifts the null hypothesis +#' H0: param = r vs H1: param != r +#' @param beta0 Deprecated function argument. Replaced by function argument 'r'. +#' @param type character or function. The character string specifies the type +#' of boostrap to use: One of `"rademacher"`, `"mammen"`, `"norm"` +#' and `"webb"`. Alternatively, type can be a function(n) for drawing +#' wild bootstrap factors. `"rademacher"` by default. +#' For the Rademacher distribution, if the number of replications `B` +#' exceeds the number of possible draw ombinations, 2^(#number of clusters), +#' then `boottest()` will use each possible combination once (enumeration). +#' @param impose_null Logical. Controls if the null hypothesis is imposed on +#' the bootstrap dgp or not. Null imposed `(WCR)` by default. +#' If `FALSE`, the null is not imposed `(WCU)` +#' @param p_val_type Character vector of length 1. Type of p-value. +#' By default `"two-tailed"`. Other options include `"equal-tailed"`, +#' `">"` and `"<"`. +#' @param tol Numeric vector of length 1. The desired accuracy +#' (convergence tolerance) used in the root finding procedure to find +#' the confidence interval. +#' `1e-6` by default. +#' @param maxiter Integer. Maximum number of iterations used in the root +#' finding procedure to find the confidence interval. +#' `10L` by default. +#' @param nthreads The number of threads. Can be: a) an integer lower than, +#' or equal to, the maximum number of threads; b) 0: meaning +#' all available threads will be used; c) a number strictly +#' between 0 and 1 which represents the fraction of all threads +#' to use. The default is to use 1 core. +#' @param ssc An object of class `boot_ssc.type` obtained with the function +#' [fwildclusterboot::boot_ssc()]. Represents how the small sample +#' adjustments are computed. The defaults are `adj = TRUE, fixef.K = "none", +#' cluster.adj = "TRUE", cluster.df = "conventional"`. +#' You can find more details in the help file for `boot_ssc()`. +#' The function is purposefully designed to mimic plm's +#' [plm::vcovSSC()] function. +#' @param getauxweights Logical. Whether to save auxilliary weight matrix (v) +#' @param floattype Float64 by default. Other option: Float32. Should floating +#' point numbers in Julia be represented as 32 or 64 bit? Only relevant when +#' 'engine = "WildBootTests.jl"' +#' @param maxmatsize `NULL` by default = no limit. Else numeric scalar to set +#' the maximum size of auxilliary weight matrix (v), in gigabytes. Only +#' relevant when `engine = "WildBootTests.jl"` +#' @param bootstrapc Logical scalar, `FALSE` by default. `TRUE` to request +#' bootstrap-c instead of bootstrap-t. Only relevant when +#' `engine = "WildBootTests.jl"` +#' @param sampling `"dqrng"` or `"standard"`. If `"dqrng"`, the 'dqrng' package is +#' used for random number generation (when available). If `"standard"`, +#' functions from the 'stats' package are used when available. +#' This argument is mostly a convenience to control random number generation in +#' a wrapper package around `fwildclusterboot`, `wildrwolf`. +#' I recommend to use the fast' option. +#' @param ... Further arguments passed to or from other methods. + +#' @importFrom dreamerr check_arg validate_dots + +#' @return An object of class `boottest` +#' +#' \item{p_val}{The bootstrap p-value.} +#' \item{conf_int}{The bootstrap confidence interval.} +#' \item{param}{The tested parameter.} +#' \item{N}{Sample size. Might differ from the regression sample size if +#' the cluster variables contain NA values.} +#' \item{boot_iter}{Number of Bootstrap Iterations.} +#' \item{clustid}{Names of the cluster Variables.} +#' \item{N_G}{Dimension of the cluster variables as used in boottest.} +#' \item{sign_level}{Significance level used in boottest.} +#' \item{type}{Distribution of the bootstrap weights.} +#' \item{impose_null}{Whether the null was imposed on the bootstrap dgp or not.} +#' \item{R}{The vector "R" in the null hypothesis of interest Rbeta = r.} +#' \item{r}{The scalar "r" in the null hypothesis of interest Rbeta = r.} +#' \item{point_estimate}{R'beta. A scalar: the constraints vector times the +#' regression coefficients.} +#' \item{grid_vals}{All t-statistics calculated while calculating the +#' confidence interval.} +#' \item{p_grid_vals}{All p-values calculated while calculating the confidence +#' interval.} +#' \item{t_stat}{The 'original' regression test statistics.} +#' \item{t_boot}{All bootstrap t-statistics.} +#' \item{regression}{The regression object used in boottest.} +#' \item{call}{Function call of boottest.} +#' \item{engine}{The employed bootstrap algorithm.} +#' \item{nthreads}{The number of threads employed.} +#' +#' @export +#' @method boottest plm +#' +#' @section Setting Seeds: +#' To guarantee reproducibility, you need to +#' set a global random seed via +#' + `set.seed()` when using +#' 1) the lean algorithm (via `engine = "R-lean"`) including the +#' heteroskedastic wild bootstrap +#' 2) the wild cluster bootstrap via `engine = "R"` with Mammen weights or +#' 3) `engine = "WildBootTests.jl"` +#' + `dqrng::dqset.seed()` when using `engine = "R"` for Rademacher, Webb +#' or Normal weights +#' +#' @section Confidence Intervals: +#' `boottest` computes confidence intervals by inverting p-values. +#' In practice, the following procedure is used: +#' \itemize{ +#' \item Based on an initial guess for starting values, calculate p-values +#' for 26 equal spaced points between the starting values. +#' \item Out of the 26 calculated p-values, find the two pairs of values x +#' for which the corresponding p-values px cross the significance +#' sign_level sign_level. +#' \item Feed the two pairs of x into an numerical root finding procedure +#' and solve for the root. boottest currently relies on +#' `stats::uniroot` and sets an absolute tolerance of 1e-06 and +#' stops the procedure after 10 iterations. +#' } +#' @section Standard Errors: +#' `boottest` does not calculate standard errors. +#' @section Stata, Julia and Python Implementations: +#' The fast wild cluster bootstrap algorithms are further implemented in the +#' following software packages: +#' \itemize{ +#' \item Stata:[boottest](https://github.com/droodman/boottest) +#' \item Julia:[WildBootTests.jl](https://github.com/droodman/WildBootTests.jl) +#' \item Python:[wildboottest](https://github.com/s3alfisc/wildboottest) +#' } +#' @references Roodman et al., 2019, "Fast and wild: Bootstrap inference in +#' STATA using boottest", The STATA Journal. +#' () +#' @references MacKinnon, James G., Morten Ørregaard Nielsen, and +#' Matthew D. Webb. Fast and reliable jackknife and bootstrap +#' methods for cluster-robust inference. No. 1485. 2022. +#' @references Cameron, A. Colin, Jonah B. Gelbach, and Douglas L. Miller. +#' "Bootstrap-based improvements for inference with clustered errors." +#' The Review of Economics and Statistics 90.3 (2008): 414-427. +#' @references Cameron, A.Colin & Douglas L. Miller. +#' "A practitioner's guide to cluster-robust inference" +#' Journal of Human Resources (2015) \doi{doi:10.3368/jhr.50.2.317} +#' @references Davidson & MacKinnon. "Wild Bootstrap Tests for IV regression" +#' Journal of Economics and Business Statistics (2010) +#' \doi{https://doi.org/10.1198/jbes.2009.07221} +#' @references MacKinnon, James G., and Matthew D. Webb. +#' "The wild bootstrap for few (treated) clusters. +#' " The Econometrics Journal 21.2 (2018): 114-135. +#' @references MacKinnon, James G., and Matthew D. Webb. +#' "Cluster-robust inference: A guide to empirical practice" +#' Journal of Econometrics (2022) +#' \doi{https://doi.org/10.1016/j.jeconom.2022.04.001} +#' @references MacKinnon, James. "Wild cluster bootstrap confidence intervals." +#' L'Actualite economique 91.1-2 (2015): 11-33. +#' @references Webb, Matthew D. Reworking wild bootstrap based inference for +#' clustered errors. No. 1315. Queen's Economics Department Working Paper, +#' 2013. + +#' @examples +#' \dontrun{ +#' library("plm") +#' data("EmplUK", package="plm") +#' +#' E <- pdata.frame(EmplUK, index=c("firm","year"), drop.index=TRUE, row.names=TRUE) +#' head(E) +#' +#' plm_fit <- plm(inv~value+capital, data = Grunfeld, model = "within") +#' +#' boot1 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = "firm" +#' ) +#' +#' boot2 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = c("firm", "year") +#' ) +#' +#' boot3 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = c("firm", "year"), +#' fe = "firm" +#' ) +#' +#' boot4 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = c("firm", "year"), +#' fe = "firm", +#' sign_level = 0.2, +#' r = 2 +#' ) +#' # test treatment + ideology1 = 2 +#' boot5 <- boottest(plm_fit, +#' B = 9999, +#' clustid = c("firm", "year"), +#' param = c("value", "capital"), +#' R = c(1, 1), +#' r = 2 +#' ) +#' summary(boot1) +#' print(boot1) +#' plot(boot1) +#' nobs(boot1) +#' pval(boot1) +#' confint(boot1) +#' generics::tidy(boot1) +#' +#' # run different bootstrap types following MacKinnon, Nielsen & Webb (2022): +#' +#' # default: the fnw algorithm +#' boot_fnw11 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = "firm", +#' bootstrap_type = "fnw11" +#' ) +#' +#' # WCR 31 +#' boot_WCR31 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = "firm", +#' bootstrap_type = "31" +#' ) +#' +#' # WCU33 +#' boot_WCR31 <- boottest(plm_fit, +#' B = 9999, +#' param = "value", +#' clustid = "firm", +#' bootstrap_type = "33", +#' impose_null = FALSE +#' ) +#' +#' +#' } +#' +boottest.plm <- function(object, + param, + B, + clustid = NULL, + bootcluster = "max", + fe = NULL, + sign_level = 0.05, + conf_int = TRUE, + R = NULL, + r = 0, + beta0 = NULL, + type = "rademacher", + impose_null = TRUE, + bootstrap_type = "fnw11", + p_val_type = "two-tailed", + tol = 1e-6, + maxiter = 10L, + sampling = "dqrng", + nthreads = getBoottest_nthreads(), + ssc = boot_ssc( + adj = TRUE, + fixef.K = "none", + cluster.adj = TRUE, + cluster.df = "conventional" + ), + engine = getBoottest_engine(), + floattype = "Float64", + maxmatsize = FALSE, + bootstrapc = FALSE, + getauxweights = FALSE, + ...) { + call <- match.call() + + dreamerr::validate_dots(stop = TRUE) + + # Step 1: check arguments of feols call + check_arg(object, "MBT class(plm)") + check_arg(clustid, "NULL | character scalar | character vector | formula") + check_arg(param, "MBT scalar character | character vector | formula") + check_arg(B, "MBT scalar integer GT{99}") + check_arg(impose_null, "logical scalar") + check_arg(bootstrap_type, "charin(11, 13, 31, 33, fnw11)") + + check_arg(sign_level, "scalar numeric GT{0} LT{1}") + check_arg(type, "charin(rademacher, mammen, norm, gamma, webb)") + check_arg(p_val_type, "charin(two-tailed, equal-tailed,>, <)") + + check_arg(conf_int, "logical scalar") + check_arg(R, "NULL| scalar numeric | numeric vector") + check_arg(r, "numeric scalar | NULL") + check_arg(fe, "character scalar | NULL | formula") + check_arg(bootcluster, "character vector | formula") + check_arg(tol, "numeric scalar GT{0}") + check_arg(maxiter, "scalar integer GT{5}") + check_arg(boot_ssc, "class(ssc) | class(boot_ssc)") + check_arg(engine, "charin(R, R-lean, WildBootTests.jl)") + + check_arg(floattype, "charin(Float32, Float64)") + check_arg(maxmatsize, "scalar integer | NULL") + check_arg(bootstrapc, "scalar logical") + + check_arg(sampling, "charin(dqrng, standard)") + + # remind packages users to set a global seed + inform_seed( + frequency_id = "seed-reminder-boot-plm", + engine = engine + ) + + model <- object$args$model + effect <- object$args$effect + supported_models <- c("pooling", "within") + + + if(!(model %in% supported_models)){ + rlang::abort( + "The model type is not supported. boottest() only supports + models 'pooling' and 'within' for plm::plm().", + use_cli_format = TRUE + ) + } + + + # plm specific tests: check if fixed effects in index + + if(!is.null(fe)){ + + if(model == "pooling"){ + rlang::abort( + "The 'fe' argument is not supported for model = 'pooling'.", + use_cli_format = TRUE + ) + } else { + + index_names <- names(index(object)) + if(!(fe %in% index_names)){ + rlang::abort( + paste("The fe variable", fe, "is not part of the model index"), + use_cli_format = TRUE + ) + } else { + + if(effect == "individual"){ + if(fe != index_names[1]){ + rlang::abort( + paste( + "When you specify 'effect=individual', the fe argument + can only be equal to", index_names[1], ", but it is", fe, "."), + use_cli_format = TRUE + ) + } + } + + if(effect == "time"){ + if(fe != index_names[2]){ + rlang::abort( + paste( + "When you specify 'effect=time', the fe argument + can only be equal to", index_names[2], ", but it is", fe, "."), + use_cli_format = TRUE + ) + } + } + + } + + } + + } + + + + if(bootstrap_type != "fnw11"){ + if(engine == "R"){ + if(conf_int){ + rlang::inform( + "Confidence Intervals are currently only supported for + the R engine with 'bootstrap_type = 'fnw11' '.", + use_cli_format = TRUE, + .frequency = "regularly", + .frequency_id = "CI only for fnw algo." + ) + } + } + } + + if (!is.null(beta0)) { + rlang::abort( + c( + "The function argument 'beta0' is deprecated. + Please use the function argument 'r' instead, + by which it is replaced." + ), + use_cli_format = TRUE + ) + } + + # parse formula inputs to strings + + if (inherits(clustid, "formula")) { + clustid <- attr(terms(clustid), "term.labels") + } + + if (inherits(bootcluster, "formula")) { + bootcluster <- attr(terms(bootcluster), "term.labels") + } + + if (inherits(param, "formula")) { + param <- attr(terms(param), "term.labels") + } + + if (inherits(fe, "formula")) { + fe <- attr(terms(fe), "term.labels") + } + + # drop singleton fixed effects? How are these handled + # by plm? + + # -------------------------------------------- + + # check appropriateness of nthreads + nthreads <- check_set_nthreads(nthreads) + + if (is.null(clustid)) { + heteroskedastic <- TRUE + if (engine == "R") { + # heteroskedastic models should always be run through R-lean + engine <- "R-lean" + } + } else { + heteroskedastic <- FALSE + } + + check_bootstrap_types( + param = param, + bootstrap_type = bootstrap_type + ) + + R_long <- process_R( + R = R, + param = param + ) + + + if (engine != "WildBootTests.jl") { + r_algo_checks( + R = R_long, + p_val_type = p_val_type, + conf_int = conf_int, + B = B + ) + } + + check_boottest_args_plus( + object = object, + R = R_long, + param = param, + sign_level = sign_level, + B = B, + fe = fe + ) + + # preprocess the data: Y, X, weights, fixed_effect + preprocess <- preprocess2.plm( + object = object, + clustid = clustid, + R = R_long, + param = param, + bootcluster = bootcluster, + fe = fe, + engine = engine, + bootstrap_type = bootstrap_type + ) + + enumerate <- + check_set_full_enumeration( + preprocess = preprocess, + heteroskedastic = heteroskedastic, + B = B, + type = type, + engine = engine + ) + + full_enumeration <- enumerate$full_enumeration + B <- enumerate$B + + N <- preprocess$N + k <- preprocess$k + G <- + vapply(preprocess$clustid, function(x) { + length(unique(x)) + }, numeric(1)) + vcov_sign <- preprocess$vcov_sign + + small_sample_correction <- + get_ssc( + boot_ssc_object = ssc, + N = N, + k = k, + G = G, + vcov_sign = vcov_sign, + heteroskedastic = heteroskedastic + ) + + # clustermin, clusteradj + + + clustid_dims <- preprocess$clustid_dims + # R*beta; + point_estimate <- + as.vector(object$coefficients[param] %*% preprocess$R0[param]) + + boot_vcov <- boot_coef <- NULL + + res <- + run_bootstrap( + object = object, + engine = engine, + preprocess = preprocess, + bootstrap_type = bootstrap_type, + B = B, + point_estimate = point_estimate, + impose_null = impose_null, + r = r, + sign_level = sign_level, + param = param, + p_val_type = p_val_type, + nthreads = nthreads, + type = type, + full_enumeration = full_enumeration, + small_sample_correction = small_sample_correction, + conf_int = conf_int, + maxiter = maxiter, + tol = tol, + clustid = clustid, + fe = fe, + R_long = R_long, + heteroskedastic = heteroskedastic, + ssc = ssc, + floattype = floattype , + bootstrapc = bootstrapc , + getauxweights = getauxweights , + maxmatsize = maxmatsize, + sampling = sampling, + bootcluster = bootcluster + + ) + + # collect results + res_final <- list( + point_estimate = point_estimate, + p_val = res$p_val, + conf_int = res$conf_int, + p_grid_vals = res$p_grid_vals, + grid_vals = res$grid_vals, + t_stat = res$t_stat, + t_boot = res$t_boot, + # regression = res$object, + param = param, + N = preprocess$N, + boot_iter = B, + clustid = clustid, + # depvar = depvar, + N_G = preprocess$N_G, + sign_level = sign_level, + call = call, + type = type, + impose_null = impose_null, + R = R_long, + r = r, + engine = engine, + nthreads = nthreads, + boot_vcov = boot_vcov, + boot_coef = boot_coef + ) + + class(res_final) <- "boottest" + invisible(res_final) + +} + + + + diff --git a/R/model_matrix.R b/R/model_matrix.R index a0b0362c..2dcff825 100644 --- a/R/model_matrix.R +++ b/R/model_matrix.R @@ -10,6 +10,54 @@ model_matrix <- function(object, ...) { UseMethod("model_matrix") } +model_matrix.plm <- function(object, type, collin.rm = TRUE, ...) { + #' Enhanced model.matrix for objects of type plm + #' @method model_matrix plm + #' @param object An object of class plm + #' @param collin.rm Should collinear variables be dropped? + #' @param type 'rhs' for right-hand side variables, 'fixef' for fixed effects + #' @param ... Other arguments + #' @noRd + + dreamerr::check_arg(type, "charin(rhs, fixef)") + + if (type == "rhs") { + mm <- model.matrix(object) + if (collin.rm == TRUE) { + bn <- names(na.omit(coef(object))) + mm <- mm[, colnames(mm) %in% bn] + } + + } else if (type == "fixef") { + + mm <- index(object) + model <- object$args$model + effect <- object$args$effect + + + if(model == "within"){ + if(effect == "individual"){ + mm <- mm[,1, drop = FALSE] + } else if(effect == "time"){ + mm <- mm[,2, drop = FALSE] + } + } else if (model == "between"){ + mm <- mm[,2, drop = FALSE] + } else{ + stop("The plm object needs to be estimated via 'within' or between + for use with fixed effects. ") + } + + # make sure all fixed effect variables are factors + i <- seq_along(mm) + mm[, i] <- lapply(i, function(x) { + factor(mm[, x]) + }) + } + + mm +} + model_matrix.lm <- function(object, collin.rm = TRUE, ...) { #' Enhanced model.matrix for objects of type lm diff --git a/R/preprocess2.R b/R/preprocess2.R index fcc09c3e..8ae9305f 100644 --- a/R/preprocess2.R +++ b/R/preprocess2.R @@ -13,6 +13,197 @@ preprocess2 <- function(object, ...) { } +preprocess2.plm <- + function(object, + clustid, + R, + param, + fe, + engine, + bootcluster, + bootstrap_type) { + #' preprocess data for objects of type plm + #' + #' @param object an object of type fixest + #' @param clustid a character string containing the name(s) of the + #' cluster variables + #' @param R Hypothesis Vector giving linear combinations of coefficients. + #' @param fe character vector. name of the fixed effect to be projected + #' out in the bootstrap + #' @param param character vector. names of the parameter(s) to test + #' @param engine The bootstrap algorithm to run. Either "R" or + #' "WildBootTests.jl" + #' @param bootcluster a character string containing the name(s) of the + #' bootcluster variables. Alternatively, "min" or "max" + #' @param bootstrap_type Which bootstrap type should be run? + #' '11', '13', '31' or 33'. For more infos, see the paper + #' by MacKinnon, Nielsen & Webb (2022) + #' + #' @noRd + #' + #' @method preprocess2 fixest + + + call <- object$call + call_env <- object$call_env + fml <- formula(object) + #fml <- Formula::as.Formula(fml) + + fml_full <- formula(fml, collapse = TRUE) + + N <- nobs(object) + # lm and felm don't drop NAs due to multicollinearity, while fixest does + k <- length(na.omit(coef(object))) + + # plm specific checks + if(!(object$args$model %in% c("pooling", "within", "between"))){ + rlang::abort("boottest() only supports estimation of 'pooling', 'within' + and 'between' models via the plm::plm() function.", + use_cli_format = TRUE + ) + } + + is_iv <- FALSE + has_fe <- FALSE + + if(object$args$model %in% c("within", "between")){ + has_fe <- TRUE + } + + mf <- model.frame(object) + Y <- as.vector(mf[,1]) + + X <- as.matrix(mf[,2:ncol(mf)]) + + if(object$args$model == "pooling"){ + X <- cbind(1, X) + colnames(X) <- names(coef(object)) + } + + weights <- weights(object) + if (is.null(weights)) { + has_weights <- FALSE + weights <- rep(1, N) + } else { + has_weights <- TRUE + } + + if (has_weights) { + if (!is.null(fe)) { + rlang::abort( + "boottest() unfortunately currently does not support WLS and fixed + effects. Please set fe = NULL to run a bootstrap with WLS.", + use_cli_format = TRUE + ) + } + } + + fixed_effect <- NULL + k2 <- 0 + W <- n_fe <- NULL + + + if (has_fe) { + # if(!is.null(fe)){ + get_fe <- transform_fe( + object = object, + X = X, + Y = Y, + fe = fe, + N = N, + has_weights = has_weights, + engine = engine, + bootstrap_type = bootstrap_type, + clustid_char = clustid + ) + + + X <- get_fe$X + Y <- get_fe$Y + fixed_effect <- get_fe$fixed_effect + W <- get_fe$W + n_fe <- get_fe$n_fe + k2 <- get_fe$k2 + + } + + + # get cluster variable + if (!is.null(clustid)) { + clustid_list <- get_cluster( + object = object, + clustid_char = clustid, + N = N, + bootcluster = bootcluster, + call_env = call_env + ) + + vcov_sign <- clustid_list$vcov_sign + clustid <- clustid_list$clustid + clustid_dims <- ncol(clustid) + N_G <- clustid_list$N_G + cluster_names <- clustid_list$cluster_names + + cluster_bootcluster <- clustid_list$cluster_bootcluster + bootcluster <- clustid_list$bootcluster + all_c <- clustid_list$all_c + } else { + vcov_sign <- + clustid_dims <- + clustid <- bootcluster <- N_G <- cluster_names <- NULL + cluster_bootcluster <- bootcluster <- all_c <- NULL + } + + + # iv prep + instruments <- X_exog <- X_endog <- NULL + if (!is.matrix(R)) { + R0 <- rep(0, length(colnames(X))) + R0[match(param, colnames(X))] <- R + names(R0) <- colnames(X) + } else { + q <- nrow(R) + p <- ncol(R) + R0 <- matrix(0, q, ncol(X)) + R0[, 1:p] <- R + } + + res <- list( + Y = Y, + X = X, + weights = weights, + fixed_effect = fixed_effect, + W = W, + n_fe = n_fe, + N = N, + k = k, + k2 = k2, + clustid = clustid, + vcov_sign = vcov_sign, + clustid_dims = clustid_dims, + N_G = N_G, + cluster_bootcluster = cluster_bootcluster, + bootcluster = bootcluster, + N_G_bootcluster = length(unique(bootcluster[[1]])), + R0 = R0, + # model_frame = model_frame, + X_exog = X_exog, + X_endog = X_endog, + instruments = instruments, + has_fe = has_fe, + all_c = all_c + ) + + if (is_iv) { + class(res) <- c("preprocess", "iv") + } else { + class(res) <- c("preprocess", "ols") + } + + res + } + + preprocess2.fixest <- function(object, @@ -21,7 +212,7 @@ preprocess2.fixest <- param, fe, engine, - bootcluster, + bootcluster, bootstrap_type) { #' preprocess data for objects of type fixest #' @@ -36,7 +227,7 @@ preprocess2.fixest <- #' "WildBootTests.jl" #' @param bootcluster a character string containing the name(s) of the #' bootcluster variables. Alternatively, "min" or "max" - #' @param bootstrap_type Which bootstrap type should be run? + #' @param bootstrap_type Which bootstrap type should be run? #' '11', '13', '31' or 33'. For more infos, see the paper #' by MacKinnon, Nielsen & Webb (2022) #' @@ -61,11 +252,11 @@ preprocess2.fixest <- rlang::abort( "boottest() only supports OLS estimation via fixest::feols() - it does not support non-linear models computed via e.g. fixest::fepois() - or fixest::feglm.", + or fixest::feglm.", use_cli_format = TRUE ) } - + # if (!is.null(object$is_sunab)) { # if(object$is_sunab == TRUE){ # rlang::abort( @@ -120,7 +311,7 @@ preprocess2.fixest <- if (!is.null(fe)) { rlang::abort( "boottest() unfortunately currently does not support WLS and fixed - effects. Please set fe = NULL to run a bootstrap with WLS.", + effects. Please set fe = NULL to run a bootstrap with WLS.", use_cli_format = TRUE ) } @@ -140,8 +331,8 @@ preprocess2.fixest <- fe = fe, N = N, has_weights = has_weights, - engine = engine, - bootstrap_type = bootstrap_type, + engine = engine, + bootstrap_type = bootstrap_type, clustid_char = clustid ) @@ -248,7 +439,7 @@ preprocess2.felm <- param, fe, engine, - bootcluster, + bootcluster, bootstrap_type) { #' preprocess data for objects of type felm #' @@ -263,8 +454,8 @@ preprocess2.felm <- #' "WildBootTests.jl" #' @param bootcluster a character string containing the name(s) of the #' bootcluster variables. Alternatively, "min" or "max" - #' @param bootstrap_type Which bootstrap type should be run? - #' '11', '13', '31' or 33'. For more infos, + #' @param bootstrap_type Which bootstrap type should be run? + #' '11', '13', '31' or 33'. For more infos, #' see the paper by MacKinnon, Nielsen & Webb (2022) #' #' @noRd @@ -297,7 +488,7 @@ preprocess2.felm <- rlang::abort( "IV regression is currently not supported by boottest() for objects of type 'felm'. You can either use 'fixest::feols()' - or 'ivreg::ivreg' for IV-regression.", + or 'ivreg::ivreg' for IV-regression.", use_cli_format = TRUE ) is_iv <- TRUE @@ -333,8 +524,8 @@ preprocess2.felm <- fe = fe, N = N, has_weights = has_weights, - engine = engine, - bootstrap_type = bootstrap_type, + engine = engine, + bootstrap_type = bootstrap_type, clustid_char = clustid ) X <- get_fe$X @@ -396,7 +587,7 @@ preprocess2.felm <- R0[, 1:p] <- R } # } - + res <- list( Y = Y, X = X, @@ -421,13 +612,13 @@ preprocess2.felm <- has_fe = has_fe, all_c = all_c ) - + if (is_iv) { class(res) <- c("preprocess", "iv") } else { class(res) <- c("preprocess", "ols") } - + res } @@ -438,7 +629,7 @@ preprocess2.lm <- R, param, engine, - bootcluster, + bootcluster, bootstrap_type = NULL) { #' preprocess data for objects of type lm #' @@ -455,23 +646,23 @@ preprocess2.lm <- #' @noRd #' #' @method preprocess2 lm - + call <- object$call call_env <- environment(formula(object)) fml <- formula(object) fml <- Formula::as.Formula(fml) - + N <- nobs(object) # lm and felm don't drop NAs due to multicollinearity, while fixest does k <- length(na.omit(coef(object))) p <- object$p - + is_iv <- FALSE - + X <- model_matrix(object, collin.rm = TRUE) Y <- model.response(model.frame(object)) has_fe <- FALSE - + weights <- weights(object) if (is.null(weights)) { has_weights <- FALSE @@ -479,7 +670,7 @@ preprocess2.lm <- } else { has_weights <- TRUE } - + # get cluster variable if (!is.null(clustid)) { clustid_list <- get_cluster( @@ -489,13 +680,13 @@ preprocess2.lm <- bootcluster = bootcluster, call_env = call_env ) - + vcov_sign <- clustid_list$vcov_sign clustid <- clustid_list$clustid clustid_dims <- ncol(clustid) N_G <- clustid_list$N_G cluster_names <- clustid_list$cluster_names - + cluster_bootcluster <- clustid_list$cluster_bootcluster bootcluster <- clustid_list$bootcluster all_c <- clustid_list$all_c @@ -505,8 +696,8 @@ preprocess2.lm <- clustid <- bootcluster <- N_G <- cluster_names <- NULL cluster_bootcluster <- bootcluster <- all_c <- NULL } - - + + instruments <- X_exog <- X_endog <- NULL if (!is.matrix(R)) { R0 <- rep(0, length(colnames(X))) @@ -518,9 +709,9 @@ preprocess2.lm <- R0 <- matrix(0, q, ncol(X)) R0[, 1:p] <- R } - - - + + + res <- list( Y = Y, X = X, @@ -545,9 +736,9 @@ preprocess2.lm <- has_fe = has_fe, all_c = all_c ) - + class(res) <- c("preprocess", "ols") - + res } @@ -559,7 +750,7 @@ preprocess2.ivreg <- R, param, engine, - bootcluster, + bootcluster, bootstrap_type = NULL) { #' preprocess data for objects of type ivreg #' @@ -576,21 +767,21 @@ preprocess2.ivreg <- #' @noRd #' #' @method preprocess2 ivreg - + call <- object$call call_env <- environment(formula(object)) fml <- formula(object) fml <- Formula::as.Formula(fml) is_iv <- TRUE - + N <- nobs(object) # lm and felm don't drop NAs due to multicollinearity, while fixest does k <- length(na.omit(coef(object))) p <- object$p - + is_iv <- FALSE has_fe <- FALSE - + X_endog <- model.matrix( object, @@ -607,11 +798,11 @@ preprocess2.ivreg <- component = "instruments", na.rm = TRUE )[, object$instruments, drop = FALSE] Y <- model.response(model.frame(object)) - + n_exog <- length(object$exogenous) n_endog <- length(object$endogenous) n_instruments <- length(object$instruments) - + weights <- weights(object) if (is.null(weights)) { has_weights <- FALSE @@ -619,7 +810,7 @@ preprocess2.ivreg <- } else { has_weights <- TRUE } - + # get cluster variable if (!is.null(clustid)) { clustid_list <- get_cluster( @@ -629,13 +820,13 @@ preprocess2.ivreg <- bootcluster = bootcluster, call_env = call_env ) - + vcov_sign <- clustid_list$vcov_sign clustid <- clustid_list$clustid clustid_dims <- ncol(clustid) N_G <- clustid_list$N_G cluster_names <- clustid_list$cluster_names - + cluster_bootcluster <- clustid_list$cluster_bootcluster bootcluster <- clustid_list$bootcluster all_c <- clustid_list$all_c @@ -645,23 +836,23 @@ preprocess2.ivreg <- clustid <- bootcluster <- N_G <- cluster_names <- NULL cluster_bootcluster <- bootcluster <- all_c <- NULL } - - + + # iv prep R0 <- rep(0, n_exog + n_endog) R0[ match( - param, + param, c( names(object$exogenous), names(object$endogenous) ) ) ] <- R names(R0) <- c( - names(object$exogenous), + names(object$exogenous), names(object$endogenous) ) - + res <- list( Y = Y, X = NULL, @@ -687,10 +878,10 @@ preprocess2.ivreg <- has_fe = has_fe, all_c = all_c ) - - + + class(res) <- c("preprocess", "iv") - + res } @@ -709,13 +900,13 @@ demean_fe <- function(X, Y, fe, has_weights, N) { #' design matrix X and depvar Y #' #' @noRd - + g <- collapse::GRP(fe, call = FALSE) X <- collapse::fwithin(X, g) Y <- collapse::fwithin(Y, g) - + fixed_effect_W <- as.factor(fe[, 1]) - + if (!has_weights) { levels(fixed_effect_W) <- (1 / table(fe)) # because duplicate levels are forbidden @@ -723,15 +914,15 @@ demean_fe <- function(X, Y, fe, has_weights, N) { rlang::abort( "Currently, boottest() does not jointly support regression weights / WLS and fixed effects. If you want to use - boottest() for inference based on WLS, please set fe = NULL.", + boottest() for inference based on WLS, please set fe = NULL.", use_cli_format = TRUE ) # levels(fixed_effect_W) <- 1 / table(fixed_effect) } - + W <- Matrix::Diagonal(N, as.numeric(as.character(fixed_effect_W))) n_fe <- length(unique(fe[, 1])) - + res <- list( X = X, Y = Y, @@ -739,13 +930,13 @@ demean_fe <- function(X, Y, fe, has_weights, N) { W = W, n_fe = n_fe ) - + res } transform_fe <- - function(object, X, Y, fe, has_weights, N, engine, + function(object, X, Y, fe, has_weights, N, engine, bootstrap_type, clustid_char) { #' preprocess the model fixed effects #' @@ -769,29 +960,29 @@ transform_fe <- #' be projected out #' #' @noRd - + all_fe <- model_matrix(object, type = "fixef", collin.rm = TRUE) # make sure all fixed effects variables are characters - + n_fe <- ncol(all_fe) all_fe_names <- names(all_fe) k2 <- Reduce("+", lapply(all_fe, function(x) { length(unique(x)) })) - + fe_df <- W <- n_fe <- NULL - + # if a fe is to be projected out in the bootstrap if (!is.null(fe)) { # add all fe except for fe to data frame fe_df <- all_fe[, fe, drop = FALSE] add_fe <- all_fe[, all_fe_names != fe, drop = FALSE] add_fe_names <- names(add_fe) - + # nothing to add if only one fixed effect in model if (length(add_fe_names) != 0) { fml_fe <- reformulate(add_fe_names, response = NULL) - + add_fe_dummies <- model.matrix( fml_fe, model.frame( @@ -807,18 +998,19 @@ transform_fe <- X <- as.matrix(collapse::add_vars(as.data.frame(X), add_fe_dummies)) } - + # project out fe if (engine == "R") { + if(bootstrap_type != "fnw11"){ if(fe != clustid_char){ - rlang::abort("Only cluster fixed effects are supported for bootstrap_types - '11', '13', '31', '33'.", + rlang::abort("Only cluster fixed effects are supported for bootstrap_types + '11', '13', '31', '33'.", use_cli_format = TRUE ) } } - + # WildBootTests.jl does demeaning internally prep_fe <- demean_fe(X, Y, fe_df, has_weights, N) X <- prep_fe$X @@ -835,7 +1027,7 @@ transform_fe <- X <- as.matrix(collapse::add_vars(as.data.frame(X), add_fe_dummies)) } - + res <- list( X = X, Y = Y, @@ -844,7 +1036,7 @@ transform_fe <- k2 = k2, fixed_effect = fe_df ) - + res } @@ -871,7 +1063,7 @@ get_cluster <- #' cluster variables, #' a data.frame of the bootcluster variable(s), and a helper #' matrix, all_c, used in `engine_julia()` - + # ----------------------------------------------------------------------- # # Note: a large part of the following code was taken and adapted from the # sandwich R package, which is distributed under GPL-2 | GPL-3 @@ -879,17 +1071,17 @@ get_cluster <- # An object-Oriented Implementation of Clustered Covariances in R." # _Journal of Statistical Software_, *95*(1), 1-36. # doi: 10.18637/jss.v095.i01 (URL: https://doi.org/10.18637/jss.v095.i01). - + # changes by Alexander Fischer: # no essential changes, but slight reorganization of pieces of code - + dreamerr::check_arg(clustid_char, "character scalar|charakter vector") dreamerr::check_arg(bootcluster, "character scalar | character vector") - + clustid_fml <- reformulate(clustid_char) # Step 1: create cluster df - - + + manipulate_object <- function(object){ if(inherits(object, "fixest")){ if(!is.null(object$fixef_vars)){ @@ -901,13 +1093,13 @@ get_cluster <- object } } - + cluster_tmp <- if ("Formula" %in% loadedNamespaces()) { ## FIXME to suppress potential warnings due to | in Formula suppressWarnings( expand.model.frame( - model = + model = manipulate_object(object), extras = clustid_fml, na.expand = FALSE, @@ -916,38 +1108,38 @@ get_cluster <- ) } else { expand.model.frame( - model = + model = manipulate_object(object), extras = clustid_fml, na.expand = FALSE, envir = call_env ) } - + # if(inherits(cluster_tmp, "try-error")){ # if(inherits(object, "fixest") || inherits(object, "felm")){ - # if(grepl("non-numeric argument to binary operator$", + # if(grepl("non-numeric argument to binary operator$", # attr(cluster_tmp, "condition")$message)){ # rlang::abort("In your model, you have specified multiple fixed effects, # none of which are of type factor. While `fixest::feols()` and - # `lfe::felm()` handle this case without any troubles, - # `boottest()` currently cannot handle this case - please + # `lfe::felm()` handle this case without any troubles, + # `boottest()` currently cannot handle this case - please # change the type of (at least one) fixed effect(s) to factor. # If this does not solve the error, please report the issue # at https://github.com/s3alfisc/fwildclusterboot.") # } # if(grepl("operations are possible only for numeric, logical - # or complex types$", + # or complex types$", # attr(cluster_tmp, "condition")$message)){ # rlang::abort("Either a fixed effect or a cluster variable in your fixest() - # or felm() model is currently specified as a character. + # or felm() model is currently specified as a character. # 'boottest()' relies on 'expand.model.frame', # which can not handle these variable types in models. # Please change these character variables to factors. ") # } # } # } - + cluster_df <- model.frame(clustid_fml, cluster_tmp, na.action = na.pass) # without cluster intersection @@ -955,9 +1147,9 @@ get_cluster <- vapply(cluster_df, function(x) { length(unique(x)) }, numeric(1)) - + # Step 1: decode bootcluster variable - + # create a bootcluster vector if (length(bootcluster) == 1) { if (bootcluster == "max") { @@ -972,7 +1164,7 @@ get_cluster <- } else { bootcluster_char <- bootcluster } - + # add bootcluster variable to formula of clusters cluster_bootcluster_fml <- update( @@ -983,14 +1175,14 @@ get_cluster <- ) ) ) - - + + cluster_bootcluster_tmp <- if ("Formula" %in% loadedNamespaces()) { ## FIXME to suppress potential warnings due to | in Formula suppressWarnings( expand.model.frame( - model = + model = manipulate_object(object), extras = cluster_bootcluster_fml, na.expand = FALSE, @@ -999,21 +1191,21 @@ get_cluster <- ) } else { expand.model.frame( - model = + model = manipulate_object(object), extras = cluster_bootcluster_fml, na.expand = FALSE, envir = call_env ) } - + # data.frame as needed for WildBootTests.jl cluster_bootcluster_df <- model.frame( cluster_bootcluster_fml, cluster_bootcluster_tmp, na.action = na.pass ) - + # if(inherits(cluster_tmp, "try-error")){ # if(inherits(object, "fixest") || inherits(object, "felm")){ # if( @@ -1025,12 +1217,12 @@ get_cluster <- # ) # ){ # rlang::abort( - # "In your model, you have specified multiple fixed effects, - # none of which are of type factor. While `fixest::feols()` and + # "In your model, you have specified multiple fixed effects, + # none of which are of type factor. While `fixest::feols()` and # `lfe::felm()` handle this case without any troubles, `boottest()` # currently cannot handle this case - please change the type of # (at least one) fixed effect(s) to factor. If this does not solve - # the error, please report the issue at + # the error, please report the issue at # https://github.com/s3alfisc/fwildclusterboot." # ) # } @@ -1044,20 +1236,20 @@ get_cluster <- # ){ # rlang::abort( # "Either a fixed effect or a cluster variable in your fixest() or - # felm() model is currently specified as a character. 'boottest()' + # felm() model is currently specified as a character. 'boottest()' # relies on 'expand.model.frame', which can not handle these variable - # types in models. + # types in models. # Please change these character variables to factors." # ) # } # } # } - + # data.frames with clusters, bootcluster cluster <- cluster_bootcluster_df[, clustid_char, drop = FALSE] bootcluster <- cluster_bootcluster_df[, bootcluster_char, drop = FALSE] - + if (!any(bootcluster_char %in% clustid_char)) { is_subcluster <- TRUE if (!(any(names(bootcluster) %in% c(clustid_char, names(coef( @@ -1065,79 +1257,79 @@ get_cluster <- )))))) { rlang::abort( "A bootcluster variable is neither contained in the cluster - variables nor in the model coefficients.", + variables nor in the model coefficients.", use_cli_format = TRUE ) } } else { is_subcluster <- FALSE } - + ## handle omitted or excluded observations (works for lfe, lm) if ((N != NROW(cluster)) && !is.null(object$na.action) && (class(object$na.action) %in% c("exclude", "omit"))) { cluster <- cluster[-object$na.action, , drop = FALSE] } - + if ((N != NROW(bootcluster)) && !is.null(object$na.action) && (class(object$na.action) %in% c("exclude", "omit"))) { bootcluster <- bootcluster[-object$na.action, , drop = FALSE] } - + if (N != nrow(cluster) && inherits(object, "fixest")) { cluster <- cluster[unlist(object$obs_selection), , drop = FALSE] bootcluster <- bootcluster[unlist(object$obs_selection), , drop = FALSE] } - + if (NROW(cluster) != N) { rlang::abort( - "The number of observations in 'cluster' and 'nobs()' do not match", + "The number of observations in 'cluster' and 'nobs()' do not match", use_cli_format = TRUE ) } if (NROW(bootcluster) != N) { rlang::abort( - "The number of observations in 'bootcluster' and 'nobs()' do not match", + "The number of observations in 'bootcluster' and 'nobs()' do not match", use_cli_format = TRUE ) } - + if (any(is.na(cluster))) { rlang::abort( "`boottest()` cannot handle NAs in `clustid` variables that are not - part of the estimated model object.", + part of the estimated model object.", use_cli_format = TRUE ) } if (any(is.na(bootcluster))) { rlang::abort( "`boottest()` cannot handle NAs in `bootcluster` variables that are - not part of the estimated model object.", + not part of the estimated model object.", use_cli_format = TRUE ) } - - + + clustid_dims <- length(clustid_char) - + i <- !vapply(cluster, is.numeric, logical(1)) cluster[i] <- lapply(cluster[i], as.character) - + # taken from multiwayvcov::cluster.boot acc <- list() for (i in 1:clustid_dims) { acc <- append(acc, utils::combn(1:clustid_dims, i, simplify = FALSE)) } - + vcov_sign <- vapply(acc, function(i) { (-1)^(length(i) + 1) }, numeric(1)) acc <- acc[-1:-clustid_dims] - + if (clustid_dims > 1) { for (i in acc) { cluster <- cbind(cluster, Reduce(paste, cluster[, i])) @@ -1145,11 +1337,11 @@ get_cluster <- Reduce(paste, names(cluster[, i])) } } - + N_G <- vapply(cluster, function(x) { length(unique(x)) }, numeric(1)) - + # now do all the other bootcluster things c1 <- bootcluster_char[which(!(bootcluster_char %in% clustid_char))] @@ -1159,13 +1351,13 @@ get_cluster <- # only error cluster: variables in clustid_char not in c1, c2 c3 <- clustid_char[which(!(clustid_char %in% c(c1, c2)))] all_c <- c(c1, c2, c3) - + if (length(bootcluster_char) > 1) { bootcluster <- as.data.frame(Reduce(paste, bootcluster)) names(bootcluster) <- Reduce(paste, bootcluster_char) } - - + + res <- list( vcov_sign = vcov_sign, clustid_dims = clustid_dims, @@ -1176,7 +1368,7 @@ get_cluster <- bootcluster = bootcluster, cluster_bootcluster = cluster_bootcluster_df ) - + res } diff --git a/man/boottest.plm.Rd b/man/boottest.plm.Rd new file mode 100644 index 00000000..888994cf --- /dev/null +++ b/man/boottest.plm.Rd @@ -0,0 +1,381 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/boottest_plm.r +\name{boottest.plm} +\alias{boottest.plm} +\title{Fast wild cluster bootstrap inference for object of class plm} +\usage{ +\method{boottest}{plm}( + object, + param, + B, + clustid = NULL, + bootcluster = "max", + fe = NULL, + sign_level = 0.05, + conf_int = TRUE, + R = NULL, + r = 0, + beta0 = NULL, + type = "rademacher", + impose_null = TRUE, + bootstrap_type = "fnw11", + p_val_type = "two-tailed", + tol = 1e-06, + maxiter = 10, + sampling = "dqrng", + nthreads = getBoottest_nthreads(), + ssc = boot_ssc(adj = TRUE, fixef.K = "none", cluster.adj = TRUE, cluster.df = + "conventional"), + engine = getBoottest_engine(), + floattype = "Float64", + maxmatsize = FALSE, + bootstrapc = FALSE, + getauxweights = FALSE, + ... +) +} +\arguments{ +\item{object}{An object of class plm and estimated via \code{plm::feols()}. +Non-linear models are not supported.} + +\item{param}{A character vector or rhs formula. The name of the regression +coefficient(s) for which the hypothesis is to be tested} + +\item{B}{Integer. The number of bootstrap iterations. When the number of +clusters is low, +increasing B adds little additional runtime.} + +\item{clustid}{A character vector or rhs formula containing the names of the +cluster variables. If NULL, +a heteroskedasticity-robust (HC1) wild bootstrap is run.} + +\item{bootcluster}{A character vector or rhs formula of length 1. Specifies +the bootstrap clustering variable or variables. If more +than one variable is specified, then bootstrapping is clustered by the +intersections of +clustering implied by the listed variables. To mimic the behavior of +stata's boottest command, +the default is to cluster by the intersection of all the variables +specified via the \code{clustid} argument, +even though that is not necessarily recommended (see the paper by +Roodman et al cited below, section 4.2). +Other options include "min", where bootstrapping is clustered by +the cluster variable with the fewest clusters. +Further, the subcluster bootstrap (MacKinnon & Webb, 2018) is +supported - see the \code{vignette("fwildclusterboot", package = "fwildclusterboot")} for details.} + +\item{fe}{A character vector or rhs formula of length one which contains +the name of the fixed effect to be projected +out in the bootstrap. Note: if regression weights are used, fe +needs to be NULL.} + +\item{sign_level}{A numeric between 0 and 1 which sets the significance level +of the inference procedure. E.g. sign_level = 0.05 +returns 0.95\% confidence intervals. By default, sign_level = 0.05.} + +\item{conf_int}{A logical vector. If TRUE, boottest computes confidence +intervals by test inversion. If FALSE, only the p-value is returned.} + +\item{R}{Hypothesis Vector giving linear combinations of coefficients. +Must be either NULL or a vector of the same length as \code{param}. If NULL, +a vector of ones of length param.} + +\item{r}{A numeric. Shifts the null hypothesis +H0: param = r vs H1: param != r} + +\item{beta0}{Deprecated function argument. Replaced by function argument 'r'.} + +\item{type}{character or function. The character string specifies the type +of boostrap to use: One of "rademacher", "mammen", "norm" +and "webb". Alternatively, type can be a function(n) for drawing +wild bootstrap factors. "rademacher" by default. +For the Rademacher distribution, if the number of replications B +exceeds +the number of possible draw ombinations, 2^(#number of clusters), +then \code{boottest()} +will use each possible combination once (enumeration).} + +\item{impose_null}{Logical. Controls if the null hypothesis is imposed on +the bootstrap dgp or not. Null imposed \code{(WCR)} by default. +If FALSE, the null is not imposed \code{(WCU)}} + +\item{bootstrap_type}{Determines which wild cluster bootstrap type should be +run. Options are "fnw11","11", "13", "31" and "33" for the wild cluster +bootstrap and "11" and "31" for the heteroskedastic bootstrap. +For more information, see the details section. "fnw11" is the default for +the cluster bootstrap, which runs a "11" type +wild cluster bootstrap via the algorithm outlined in "fast and wild" +(Roodman et al (2019)). "11" is the default for the heteroskedastic +bootstrap.} + +\item{p_val_type}{Character vector of length 1. Type of p-value. +By default "two-tailed". Other options include "equal-tailed", +">" and "<".} + +\item{tol}{Numeric vector of length 1. The desired accuracy +(convergence tolerance) used in the root finding procedure to find +the confidence interval. +1e-6 by default.} + +\item{maxiter}{Integer. Maximum number of iterations used in the root +finding procedure to find the confidence interval. +10 by default.} + +\item{sampling}{'dqrng' or 'standard'. If 'dqrng', the 'dqrng' package is +used for random number generation (when available). If 'standard', +functions from the 'stats' package are used when available. +This argument is mostly a convenience to control random number generation in +a wrapper package around \code{fwildclusterboot}, \code{wildrwolf}. +I recommend to use the fast' option.} + +\item{nthreads}{The number of threads. Can be: a) an integer lower than, +or equal to, the maximum number of threads; b) 0: meaning +all available threads will be used; c) a number strictly +between 0 and 1 which represents the fraction of all threads +to use. The default is to use 1 core.} + +\item{ssc}{An object of class \code{boot_ssc.type} obtained with the function +\code{\link[=boot_ssc]{boot_ssc()}}. Represents how the small sample +adjustments are computed. The defaults are \verb{adj = TRUE, fixef.K = "none", cluster.adj = "TRUE", cluster.df = "conventional"}. +You can find more details in the help file for \code{boot_ssc()}. +The function is purposefully designed to mimic plm's +\code{\link[fixest:ssc]{fixest::ssc()}} function.} + +\item{engine}{Character scalar. Either "R", "R-lean" or "WildBootTests.jl". +Controls if \code{boottest()} should run via its native R implementation +or \code{WildBootTests.jl}. +"R" is the default and implements the cluster bootstrap +as in Roodman (2019). "WildBootTests.jl" executes the +wild cluster bootstrap via the WildBootTests.jl +package. For it to run, Julia and WildBootTests.jl need +to be installed. +The "R-lean" algorithm is a memory friendly, but less +performant rcpp-armadillo based implementation of the wild +cluster bootstrap. +Note that if no cluster is provided, boottest() always +defaults to the "lean" algorithm. You can set the employed +algorithm globally by using the +\code{setBoottest_engine()} function.} + +\item{floattype}{Float64 by default. Other option: Float32. Should floating +point numbers in Julia be represented as 32 or 64 bit? Only relevant when +'engine = "WildBootTests.jl"'} + +\item{maxmatsize}{NULL by default = no limit. Else numeric scalar to set +the maximum size of auxilliary weight matrix (v), in gigabytes. Only +relevant when 'engine = "WildBootTests.jl"'} + +\item{bootstrapc}{Logical scalar, FALSE by default. TRUE to request +bootstrap-c instead of bootstrap-t. Only relevant when +'engine = "WildBootTests.jl"'} + +\item{getauxweights}{Logical. Whether to save auxilliary weight matrix (v)} + +\item{...}{Further arguments passed to or from other methods.} +} +\value{ +An object of class \code{boottest} + +\item{p_val}{The bootstrap p-value.} +\item{conf_int}{The bootstrap confidence interval.} +\item{param}{The tested parameter.} +\item{N}{Sample size. Might differ from the regression sample size if +the cluster variables contain NA values.} +\item{boot_iter}{Number of Bootstrap Iterations.} +\item{clustid}{Names of the cluster Variables.} +\item{N_G}{Dimension of the cluster variables as used in boottest.} +\item{sign_level}{Significance level used in boottest.} +\item{type}{Distribution of the bootstrap weights.} +\item{impose_null}{Whether the null was imposed on the bootstrap dgp or not.} +\item{R}{The vector "R" in the null hypothesis of interest Rbeta = r.} +\item{r}{The scalar "r" in the null hypothesis of interest Rbeta = r.} +\item{point_estimate}{R'beta. A scalar: the constraints vector times the +regression coefficients.} +\item{grid_vals}{All t-statistics calculated while calculating the +confidence interval.} +\item{p_grid_vals}{All p-values calculated while calculating the confidence +interval.} +\item{t_stat}{The 'original' regression test statistics.} +\item{t_boot}{All bootstrap t-statistics.} +\item{regression}{The regression object used in boottest.} +\item{call}{Function call of boottest.} +\item{engine}{The employed bootstrap algorithm.} +\item{nthreads}{The number of threads employed.} +} +\description{ +\code{boottest.plm} is a S3 method that allows for fast wild cluster +bootstrap inference for objects of class plm by implementing +fast wild bootstrap algorithms as developed in Roodman et al., 2019 +and MacKinnon, Nielsen & Webb (2022). +} +\section{Setting Seeds}{ + +To guarantee reproducibility, you need to +set a global random seed via +\itemize{ +\item \code{set.seed()} when using +\enumerate{ +\item the lean algorithm (via \code{engine = "R-lean"}) including the +heteroskedastic wild bootstrap +\item the wild cluster bootstrap via \code{engine = "R"} with Mammen weights or +\item \code{engine = "WildBootTests.jl"} +} +\item \code{dqrng::dqset.seed()} when using \code{engine = "R"} for Rademacher, Webb +or Normal weights +} +} + +\section{Confidence Intervals}{ + +\code{boottest} computes confidence intervals by inverting p-values. +In practice, the following procedure is used: +\itemize{ +\item Based on an initial guess for starting values, calculate p-values +for 26 equal spaced points between the starting values. +\item Out of the 26 calculated p-values, find the two pairs of values x +for which the corresponding p-values px cross the significance +sign_level sign_level. +\item Feed the two pairs of x into an numerical root finding procedure +and solve for the root. boottest currently relies on +\code{stats::uniroot} and sets an absolute tolerance of 1e-06 and +stops the procedure after 10 iterations. +} +} + +\section{Standard Errors}{ + +\code{boottest} does not calculate standard errors. +} + +\section{Stata, Julia and Python Implementations}{ + +The fast wild cluster bootstrap algorithms are further implemented in the +following software packages: +\itemize{ +\item Stata:\href{https://github.com/droodman/boottest}{boottest} +\item Julia:\href{https://github.com/droodman/WildBootTests.jl}{WildBootTests.jl} +\item Python:\href{https://github.com/s3alfisc/wildboottest}{wildboottest} +} +} + +\examples{ +\dontrun{ +library("plm") +data("EmplUK", package="plm") + +E <- pdata.frame(EmplUK, index=c("firm","year"), drop.index=TRUE, row.names=TRUE) +head(E) + +plm_fit <- plm(inv~value+capital, data = Grunfeld, model = "within") + +boot1 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = "firm" +) + +boot2 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = c("firm", "year") +) + +boot3 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = c("firm", "year"), + fe = "firm" +) + +boot4 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = c("firm", "year"), + fe = "firm", + sign_level = 0.2, + r = 2 +) +# test treatment + ideology1 = 2 +boot5 <- boottest(plm_fit, + B = 9999, + clustid = c("firm", "year"), + param = c("value", "capital"), + R = c(1, 1), + r = 2 +) +summary(boot1) +print(boot1) +plot(boot1) +nobs(boot1) +pval(boot1) +confint(boot1) +generics::tidy(boot1) + +# run different bootstrap types following MacKinnon, Nielsen & Webb (2022): + +# default: the fnw algorithm +boot_fnw11 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = "firm", + bootstrap_type = "fnw11" +) + +# WCR 31 +boot_WCR31 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = "firm", + bootstrap_type = "31" +) + +# WCU33 +boot_WCR31 <- boottest(plm_fit, + B = 9999, + param = "value", + clustid = "firm", + bootstrap_type = "33", + impose_null = FALSE +) + + +} + +} +\references{ +Roodman et al., 2019, "Fast and wild: Bootstrap inference in +STATA using boottest", The STATA Journal. +(\url{https://ideas.repec.org/p/qed/wpaper/1406.html}) + +MacKinnon, James G., Morten Ørregaard Nielsen, and +Matthew D. Webb. Fast and reliable jackknife and bootstrap +methods for cluster-robust inference. No. 1485. 2022. + +Cameron, A. Colin, Jonah B. Gelbach, and Douglas L. Miller. +"Bootstrap-based improvements for inference with clustered errors." +The Review of Economics and Statistics 90.3 (2008): 414-427. + +Cameron, A.Colin & Douglas L. Miller. +"A practitioner's guide to cluster-robust inference" +Journal of Human Resources (2015) \doi{doi:10.3368/jhr.50.2.317} + +Davidson & MacKinnon. "Wild Bootstrap Tests for IV regression" +Journal of Economics and Business Statistics (2010) +\doi{https://doi.org/10.1198/jbes.2009.07221} + +MacKinnon, James G., and Matthew D. Webb. +"The wild bootstrap for few (treated) clusters. +" The Econometrics Journal 21.2 (2018): 114-135. + +MacKinnon, James G., and Matthew D. Webb. +"Cluster-robust inference: A guide to empirical practice" +Journal of Econometrics (2022) +\doi{https://doi.org/10.1016/j.jeconom.2022.04.001} + +MacKinnon, James. "Wild cluster bootstrap confidence intervals." +L'Actualite economique 91.1-2 (2015): 11-33. + +Webb, Matthew D. Reworking wild bootstrap based inference for +clustered errors. No. 1315. Queen's Economics Department Working Paper, +2013. +} diff --git a/tests/testthat/test_plm.R b/tests/testthat/test_plm.R new file mode 100644 index 00000000..30cfbe6d --- /dev/null +++ b/tests/testthat/test_plm.R @@ -0,0 +1,208 @@ +test_that("test plm method", { + + + library(fwildclusterboot) + library("plm") + data("Grunfeld", package="plm") + + #Grunfeld <- pdata.frame(Grunfeld, index=c("firm","year"), drop.index=TRUE, row.names=TRUE) + #head(Grunfeld) + + + # -------------------------------------------------------------- # + # within estimation + + # model = within, effect = individual + + plm_fit <- plm( + inv~value+capital, + data = Grunfeld, + model = "within", + effect = "individual" + ) + broom::tidy(plm_fit) + lm_fit <- lm(inv~value+capital + as.factor(firm), data = Grunfeld) + broom::tidy(lm_fit)[c(1, 2, 3),] + + dqrng::dqset.seed(233) + lm_boot <- boottest( + lm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + fe = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + + # model = within, effect = time + + plm_fit <- plm( + inv~value+capital, + data = Grunfeld, + model = "within", + effect = "time" + ) + + broom::tidy(plm_fit) + lm_fit <- lm(inv ~ value + capital + as.factor(year), data = Grunfeld) + broom::tidy(lm_fit)[c(1, 2, 3),] + + dqrng::dqset.seed(233) + lm_boot <- boottest( + lm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + fe = "year", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + # model = within, effect = "twoway" + + plm_fit <- plm( + inv~value+capital, + data = Grunfeld, + model = "within", + effect = "twoway" + ) + + broom::tidy(plm_fit) + lm_fit <- lm(inv ~ value + capital + as.factor(firm) + as.factor(year), data = Grunfeld) + broom::tidy(lm_fit)[c(1, 2, 3),] + + dqrng::dqset.seed(233) + lm_boot <- boottest( + lm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = "firm", + fe = "year", + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + + # -------------------------------------------------------------- # + # pooling estimation + + plm_fit <- plm(inv~value+capital, data = Grunfeld, model = "pooling") + summary(plm_fit) + lm_fit <- lm(inv~value+capital , data = Grunfeld) + summary(lm_fit) + + dqrng::dqset.seed(233) + lm_boot <- boottest( + lm_fit, + param = "value", + clustid = c("year", "firm"), + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + dqrng::dqset.seed(233) + plm_boot <- boottest( + plm_fit, + param = "value", + clustid = c("year", "firm"), + B = 999, + ssc = boot_ssc(adj = FALSE, cluster.adj = FALSE) + ) + + expect_equal(coef(lm_fit)[c("(Intercept)","value", "capital")], coef(plm_fit)) + expect_equal(pval(lm_boot), pval(plm_boot)) + expect_equal(teststat(lm_boot), teststat(plm_boot)) + expect_equal(confint(lm_boot), confint(plm_boot)) + + + + + + + + + +})