CRAN v.1.0.6

Fix bugs with clustering in main function and AR test.

DESCRIPTION

@@ -1,7 +1,7 @@
 Package: ivDiag
 Title: Estimation and Diagnostic Tools for Instrumental Variables Designs
-Version: 1.0.5
-Date: 2023-05-01
+Version: 1.0.6
+Date: 2023-09-16
 Authors@R: c(
   person("Apoorva", "Lal", , "[email protected]", c("aut"), comment = c(ORCID = "0000-0002-3697-614X")),
   person("Yiqing", "Xu", , "[email protected]", c("aut", "cre"), comment = c(ORCID = "0000-0003-2041-6671"))  

NAMESPACE

@@ -19,6 +19,6 @@ importFrom(patchwork, wrap_plots)
 importFrom("graphics", "abline", "axTicks", "axis", "box", "legend",
            "mtext", "points", "rect", "segments", "text")
 importFrom("stats", "as.formula", "complete.cases", "cor", "dnorm",
-           "pnorm", "qnorm", "quantile", "sd")
+           "pnorm", "qnorm", "quantile", "sd", "pf")
 importFrom("utils", "tail")
 importFrom("testthat", "test_that", "expect_equal")

NEWS.md

@@ -1,3 +1,7 @@
+# ivDiag 1.0.6
+Fix clustering issues with both the main function and AR test (thanks to Michael
+Allen @mmooddaa for reporting this bug).
+
 # ivDiag 1.0.5
 Fix bugs with `plot_coef()`.
 

R/AR_test.R

@@ -57,31 +57,37 @@ AR_test = function(
   colnames(d) <- c(Y, D, Z, cl, weights)
 
   # reduced form
-  fmla = formula_lfe(Y = Y, W = Z, X = NULL, FE = NULL, Cl = cl)
-  if (is.null(weights) == TRUE) {
-    m1 = lfe::felm(fmla, data = d)
-  } else {
-    m1 = lfe::felm(fmla, data = d, weights = d[, weights])
-  }
-  m1 <- robustify(m1)
-  s <- summary(m1)
-  Fstat <- s$F.fstat
+  fmla = formula_lfe(Y = Y, D = Z, X = NULL, FE = NULL, cl = cl)
+
+
+  # function to test each value on the real line
+  one.AR <- function(beta, d, Y, D, Z, cl, weights) {
+    d[, Y] <- d[, Y] - beta * d[, D]
+    p_iv <- length(Z)
+    reg <- OLS(data = d, Y = Y, D = Z, X = NULL, FE = NULL, cl = cl, weights = weights)
+    coef <- reg$coef
+    vcov <- reg$vcv
+    df2 <- reg$df
+    df1 <- p_iv
+    Fstat <- c((t(coef) %*% solve(vcov) %*% coef) / p_iv)
+    pval <- pf(Fstat, df1, df2, lower.tail = FALSE)
+    output <- c(Fstat, df1, df2, pval)
+    names(output) <- c("F", "df1", "df2", "p")
+    return(output)
+  } 
 
+  # AR test
+  Fstat <- one.AR(0, d = d, Y = Y, D = D, Z = Z, cl = cl, weights = weights)
+
+
   # Confidence intervals
   if (CI == TRUE) {
     message("AR Test Inversion...\n")
     if (parallel == FALSE) {
       accept <- rep(NA, ngrid)
-      for (i in 1:ngrid) {
-        d[, Y] <- Ytil - beta_seq[i] * Dtil
-        if (is.null(weights) == TRUE) {
-          m2 = lfe::felm(fmla, data = d)
-        } else {
-          m2 = lfe::felm(fmla, data = d, weights = d[, weights])
-        }
-        m2 <- robustify(m2)
-        s2 <- summary(m2)
-        accept[i] <- ifelse(s2$pval >= alpha, 1, 0)
+      for (i in 1:ngrid) {    
+        test.out <-  one.AR(beta_seq[i], d = d, Y = Y, D = D, Z = Z, cl = cl, weights = weights) 
+        accept[i] <- ifelse(test.out[4] >= alpha, 1, 0)
       }
     } else {
       # parallel computing
@@ -92,33 +98,26 @@ AR_test = function(
       # register
       cl.parallel <- future::makeClusterPSOCK(cores, verbose = FALSE)
       doParallel::registerDoParallel(cl.parallel)
-      expfun <- c("OLS", "IV", "formula_lfe", "robustify")
+      expfun <- c("OLS", "IV", "formula_lfe", "OLS")
       accept <- foreach(
         i = 1:ngrid, .combine = c, .inorder = FALSE,
         .export = expfun,
         .packages = c("lfe")
       ) %dopar% {
-        d[, Y] <- Ytil - beta_seq[i] * Dtil
-        if (is.null(weights) == TRUE) {
-          m2 = lfe::felm(fmla, data = d)
-        } else {
-          m2 = lfe::felm(fmla, data = d, weights = d[, weights])
-        }
-        m2 <- robustify(m2)
-        s2 <- summary(m2)
-        return(ifelse(s2$pval >= alpha, 1, 0))
+        test.out <-  one.AR(beta_seq[i], d = d, Y = Y, D = D, Z = Z, cl = cl, weights = weights) 
+        return(ifelse(test.out[4] >= alpha, 1, 0))
       }
       doParallel::stopImplicitCluster()
     }
-    # summarize
+    # summarize ("accept" means reject the null)
     bounded <- FALSE
     if (sum(accept) == ngrid) {
       ci <- c(-Inf, Inf)
       ci.print <- "(-Inf, Inf)" # all accepted
     } else if (sum(accept) == 0) {
       ci <- NA
       ci.print <- "empty"
-    } else if (accept[1] == 0 && accept[ngrid] == 0) {
+    } else if (accept[1] == 0 && accept[ngrid] == 0) { # e.g. 0 0 0 1 1 1 0 0 0
       betas <- range(beta_seq[accept == 1])
       ci <- round(betas, prec)
       ci.print <- paste0("[", sprintf(paste0("%.", prec, "f"), betas[1]), ", ", sprintf(paste0("%.", prec, "f"), betas[2]), "]") # bounded interval
@@ -137,6 +136,11 @@ AR_test = function(
       ci.print <- paste0("(-Inf, ", sprintf(paste0("%.", prec, "f"), betas[2]), "]")
     }
   }
+
+  ## Print acceptance ##
+  # accept <- cbind(beta_seq, accept)
+  # rownames(accept) <- NULL
+  # colnames(accept) <- c("beta", "accept")
 
   # output
   if (CI == TRUE) {

R/effF.R

@@ -12,7 +12,7 @@
 eff_F = function(data, Y, D, Z, controls = NULL, FE = NULL, 
   cl = NULL, weights = NULL, prec = 4
   ) {
-  fmla = formula_lfe(Y = Y, W = D, Z = Z, X = controls, FE = FE, Cl = cl)
+  fmla = formula_lfe(Y = Y, D = D, Z = Z, X = controls, FE = FE, cl = cl)
   if (is.null(weights)) {
     ivmod = lfe::felm(fmla, data = data)
   } else {

R/ivDiag.R

@@ -161,7 +161,7 @@ ivDiag <- function(
       # register
       cl.parallel <- future::makeClusterPSOCK(cores, verbose = FALSE)
       doParallel::registerDoParallel(cl.parallel)
-      expfun <- c("OLS", "IV", "formula_lfe", "robustify")
+      expfun <- c("OLS", "IV", "formula_lfe")
       boot.out <- foreach::foreach(
         i = 1:nboots, .combine = rbind, .inorder = FALSE,
         .export = expfun,

R/ltz.R

@@ -23,15 +23,19 @@ ltz = function(data, Y, D, Z, controls, FE = NULL, cl = NULL, weights = NULL, pr
   Sig <- prior[2]^2
 
   # IV fit
-  fmla = formula_lfe(Y = Y, W = D, Z = Z, X = controls, FE = FE, Cl = cl)
+  fmla = formula_lfe(Y = Y, D = D, Z = Z, X = controls, FE = FE, cl = cl)
   if (is.null(weights)) {
-    m2 = robustify(lfe::felm(fmla, data = data))
+    m2 = lfe::felm(fmla, data = data)
   } else {
-    m2 = robustify(lfe::felm(fmla, data = data, weights = data[, weights]))
+    m2 = lfe::felm(fmla, data = data, weights = data[, weights])
   }
   iv_beta <- c(tail(m2$coefficients, n = 1)) # felm IV fits have endog coef at the tail
-  iv_Var <- c(tail(diag(m2$robustvcv), n = 1))
-  iv_se <- sqrt(iv_Var)
+  if (is.null(cl) == TRUE) {
+    iv_se <- tail(m2$rse, n = 1)
+  } else {
+    iv_se <- tail(m2$cse, n = 1)
+  }
+  iv_Var <- iv_se^2
   iv_t <- iv_beta / iv_se
   iv_ci <- qnorm(c(0.025, 0.975), iv_beta, iv_se)
   iv_p <- (1 - pnorm(abs(iv_t))) * 2

R/plot_coef.R

@@ -138,7 +138,7 @@ plot_coef <- function(out,
   abline(h = 0, col = "red", lwd = 2, lty = "solid")
   segments(y0 = data$lower_ci, x0 = c(1: ncoefs), y1 = data$upper_ci, x1 = c(1: ncoefs), lwd = 2) #CI
   points(1: ncoefs, data$coef, pch = c(rep(16, ncoefs.ols), rep(17, ncoefs.iv)), col = "blue") #point coefs
-  if ("ar" %in% iv.methods) {
+  if ("iv-ar" %in% iv.methods) {
     if (out$AR$bounded == FALSE) {
       if (length(out$AR$ci) == 1) {
         ar.lab <- "Empty CI"