Added conf_int() to calculate confidence intervals for individual reg…

…ression coefficients.

DESCRIPTION

@@ -41,4 +41,4 @@ Suggests:
     plm (>= 1.6-4),
     testthat,
     rmarkdown
-RoxygenNote: 5.0.1
+RoxygenNote: 6.0.1

NAMESPACE

@@ -7,6 +7,7 @@ S3method(bread,mlm)
 S3method(print,Wald_test_clubSandwich)
 S3method(print,clubSandwich)
 S3method(print,coef_test_clubSandwich)
+S3method(print,conf_int_clubSandwich)
 S3method(vcovCR,default)
 S3method(vcovCR,glm)
 S3method(vcovCR,gls)
@@ -20,6 +21,7 @@ S3method(vcovCR,rma.uni)
 S3method(vcovCR,robu)
 export(Wald_test)
 export(coef_test)
+export(conf_int)
 export(impute_covariance_matrix)
 export(vcovCR)
 import(stats)

NEWS

@@ -1,3 +1,8 @@
+clubSandwich 0.3.1.999
+=======================
+* Added conf_int() to provide easy cluster-robust confidence intervals.
+* Added examples to documentaiton for conf_int() and coef_test().
+
 clubSandwich 0.3.1
 =======================
 * Added "coefs" option to coef_test() to allow testing of subsets of coefficients.

R/coef_test.R

@@ -78,7 +78,7 @@ get_which_coef <- function(beta, coefs) {
 # coeftest for all model coefficients
 #---------------------------------------------
 
-#' Test all regression coefficients in a fitted model
+#' Test all or selected regression coefficients in a fitted model
 #'
 #' \code{coef_test} reports t-tests for each coefficient estimate in a fitted
 #' linear regression model, using a sandwich estimator for the standard errors
@@ -108,6 +108,15 @@ get_which_coef <- function(beta, coefs) {
 #'
 #' @seealso \code{\link{vcovCR}}
 #'
+#' @examples 
+#' data("Produc", package = "plm")
+#' lm_individual <- lm(log(gsp) ~ 0 + state + log(pcap) + log(pc) + log(emp) + unemp, data = Produc)
+#' individual_index <- !grepl("state", names(coef(lm_individual)))
+#' coef_test(lm_individual, vcov = "CR2", cluster = Produc$state, coefs = individual_index)
+#' 
+#' V_CR2 <- vcovCR(lm_individual, cluster = Produc$state, type = "CR2")
+#' coef_test(lm_individual, vcov = V_CR2, coefs = individual_index)
+#' 
 #' @export
 
 coef_test <- function(obj, vcov, test = "Satterthwaite", coefs = "All", ...) {

R/conf_int.R

@@ -0,0 +1,89 @@
+
+#--------------------------------------------------
+# confidence intervals for all model coefficients
+#---------------------------------------------------
+
+#' Calculate confidence intervals for all or selected regression coefficients in a fitted model
+#'
+#' \code{conf_int} reports confidence intervals for each coefficient estimate in a fitted
+#' linear regression model, using a sandwich estimator for the standard errors
+#' and a small sample correction for the critical values. The small-sample correction is
+#' based on a Satterthwaite approximation.
+#'
+#' @param obj Fitted model for which to calculate confidence intervals.
+#' @param level Desired coverage level for confidence intervals. 
+#' @inheritParams coef_test
+#' 
+#' @return A data frame containing estimated regression coefficients, standard errors, and confidence intervals. 
+#'
+#' @seealso \code{\link{vcovCR}}
+#' 
+#' @examples 
+#' data("Produc", package = "plm")
+#' lm_individual <- lm(log(gsp) ~ 0 + state + log(pcap) + log(pc) + log(emp) + unemp, data = Produc)
+#' individual_index <- !grepl("state", names(coef(lm_individual)))
+#' conf_int(lm_individual, vcov = "CR2", cluster = Produc$state, coefs = individual_index)
+#' 
+#' V_CR2 <- vcovCR(lm_individual, cluster = Produc$state, type = "CR2")
+#' conf_int(lm_individual, vcov = V_CR2, level = .99, coefs = individual_index)
+#'
+#' @export
+
+conf_int <- function(obj, vcov, level = .95, test = "Satterthwaite", coefs = "All", ...) {
+
+  if (level <= 0 | level >= 1) stop("Confidence level must be between 0 and 1.")
+
+  beta_full <- coef_CS(obj)
+  beta_NA <- is.na(beta_full)
+
+  which_beta <- get_which_coef(beta_full, coefs)
+
+  beta <- beta_full[which_beta & !beta_NA]
+
+  if (is.character(vcov)) vcov <- vcovCR(obj, type = vcov, ...)
+  if (!("clubSandwich" %in% class(vcov))) stop("Variance-covariance matrix must be a clubSandwich.")
+
+  all_tests <- c("z","naive-t","Satterthwaite")
+  test <- match.arg(test, all_tests, several.ok = FALSE)
+
+  SE <- sqrt(diag(vcov))[which_beta[!beta_NA]]
+
+  if (test=="Satterthwaite") {
+    P_array <- get_P_array(get_GH(obj, vcov))[,,which_beta[!beta_NA],drop=FALSE]
+  }
+
+  df <- switch(test, 
+               z = Inf,
+               `naive-t` = nlevels(attr(vcov, "cluster")) - 1,
+               `Satterthwaite` = Satterthwaite(beta = beta, SE = SE, P_array = P_array)$df
+  )
+
+  crit <- qt(1 - (1 - level) / 2, df = df)
+
+  result <- data.frame(
+    beta = beta, 
+    SE = SE,
+    CI_L = beta - SE * crit,
+    CI_U = beta + SE * crit
+  )
+
+  class(result) <- c("conf_int_clubSandwich", class(result))
+  attr(result, "type") <- attr(vcov, "type")
+  attr(result, "level") <- level
+  result
+}
+
+#---------------------------------------------
+# print method for conf_int
+#---------------------------------------------
+
+#' @export
+
+print.conf_int_clubSandwich <- function(x, digits = 3, ...) {
+  lev <- paste0(100 * attr(x, "level"), "%")
+  res <- data.frame("Coef" = rownames(x), x)
+  rownames(res) <- NULL
+  names(res) <- c("Coef", "Estimate", "SE", paste(c("Lower", "Upper"), lev, "CI"))
+  print(format(res, digits = 3))
+}
+

tests/testthat/test_conf_int.R

@@ -0,0 +1,62 @@
+context("confidence intervals")
+
+library(nlme, quietly=TRUE, warn.conflicts=FALSE)
+
+data(Ovary, package = "nlme")
+
+Ovary$time_int <- 1:nrow(Ovary)
+
+gls_fit <- gls(follicles ~ sin(2*pi*Time) + cos(2*pi*Time), data = Ovary,
+              correlation = corAR1(form = ~ time_int | Mare), 
+              weights = varPower())
+
+CRs <- paste0("CR", 0:4)
+
+test_that("vcov arguments work", {
+  VCR <- lapply(CRs, function(t) vcovCR(gls_fit, type = t))
+  CI_A <- lapply(VCR, function(v) conf_int(gls_fit, vcov = v, level = .98))
+  CI_B <- lapply(CRs, function(t) conf_int(gls_fit, vcov = t, level = .98))
+  expect_identical(CI_A, CI_B)
+})
+
+test_that("coefs argument works", {
+  which_grid <- expand.grid(rep(list(c(FALSE,TRUE)), length(coef(gls_fit))))
+  tests_all <- conf_int(gls_fit, vcov = "CR0", coefs = "All")
+
+  CI_A <- apply(which_grid[-1,], 1, function(x) tests_all[x,])
+  CI_B <- apply(which_grid[-1,], 1, function(x) conf_int(gls_fit, vcov = "CR0", coefs = x))
+  expect_identical(CI_A, CI_B)
+})
+
+test_that("printing works", {
+  CIs <- conf_int(gls_fit, vcov = "CR0")
+  expect_output(print(CIs))
+})
+
+test_that("level checks work", {
+  expect_error(conf_int(gls_fit, vcov = "CR0", level = -0.01))
+  expect_error(conf_int(gls_fit, vcov = "CR0", level = 95))
+  expect_output(print(conf_int(gls_fit, vcov = "CR0", level = runif(1))))
+})
+
+test_that("CI boundaries are ordered", {
+  lev <- runif(1)
+  CI_z <- conf_int(gls_fit, vcov = "CR0", test = "z", level = lev)
+  CI_t <- conf_int(gls_fit, vcov = "CR0", test = "naive-t", level = lev)
+  CI_Satt <- conf_int(gls_fit, vcov = "CR0", test = "Satterthwaite", level = lev)
+  expect_true(all(CI_t$CI_L < CI_z$CI_L))
+  expect_true(all(CI_t$CI_U > CI_z$CI_U))
+  expect_true(all(CI_Satt$CI_L < CI_z$CI_L))
+  expect_true(all(CI_Satt$CI_U > CI_z$CI_U))
+})
+
+test_that("conf_int() is consistent with coef_test()", {
+  lev <- runif(1)
+  CIs <- lapply(CRs, function(v) conf_int(gls_fit, vcov = v, test = "Satterthwaite", level = lev))
+  ttests <- lapply(CRs, function(v) coef_test(gls_fit, vcov = v, test = "Satterthwaite"))
+  CI_L <- lapply(ttests, function(x) x$beta - x$SE * qt(1 - (1 - lev) / 2, df = x$df))
+  CI_U <- lapply(ttests, function(x) x$beta + x$SE * qt(1 - (1 - lev) / 2, df = x$df))
+  expect_identical(lapply(CIs, function(x) x$CI_L), CI_L)
+  expect_identical(lapply(CIs, function(x) x$CI_U), CI_U)
+})
+