fix: restore advisor's quarto book and properly merge feature-cm

DESCRIPTION

@@ -2,7 +2,7 @@ Package: fect
 Type: Package
 Title: Fixed Effects Counterfactual Estimators
 Version: 2.2.0
-Date: 2026-03-21
+Date: 2026-03-20
 Authors@R: 
     c(person("Licheng", "Liu", , "lichengl@stanford.edu", role = c("aut")), 
       person("Ziyi", "Liu", , "zyliu2023@berkeley.edu", role = c("aut")),
@@ -11,9 +11,12 @@ Authors@R:
       person("Tianzhu", "Qin", , "tianzhu@stanford.edu", role = c("aut")),
       person("Shiyun", "Hu", , "hushiyun@pku.edu.cn", role = c("aut")),
       person("Rivka", "Lipkovitz", , "rivkal@mit.edu", role = c("aut")))
+Author: Licheng Liu [aut], Ziyi Liu [aut], Ye Wang [aut], Yiqing Xu [aut, cre],
+    Tianzhu Qin [aut], Shiyun Hu [aut], Rivka Lipkovitz [aut]
 Maintainer: Yiqing Xu <yiqingxu@stanford.edu>
 Description: Provides tools for estimating causal effects in panel data using counterfactual methods, as well as other modern DID estimators. It is designed for causal panel analysis with binary treatments under the parallel trends assumption. The package supports scenarios where treatments can switch on and off and allows for limited carryover effects. It includes several imputation estimators, such as Gsynth (Xu 2017), linear factor models, and the matrix completion method. Detailed methodology is described in Liu, Wang, and Xu (2024) <doi:10.48550/arXiv.2107.00856> and Chiu et al. (2025) <doi:10.48550/arXiv.2309.15983>. Optionally integrates with the "HonestDiDFEct" package for sensitivity analyses compatible with imputation estimators. "HonestDiDFEct" is not on CRAN but can be obtained from <https://github.com/lzy318/HonestDiDFEct>.
-URL: https://yiqingxu.org/packages/fect/
+URL: https://yiqingxu.org/packages/fect/, https://github.com/xuyiqing/fect
+BugReports: https://github.com/xuyiqing/fect/issues
 NeedsCompilation: yes
 License: MIT + file LICENSE
 Imports: 
@@ -25,14 +28,14 @@ Imports:
     foreach (>= 1.4.3), 
     abind (>= 1.4-0),
     codetools,
-    MASS,
+    MASS, 
     gridExtra, 
     grid, 
     fixest, 
     doRNG,
-    future,
-    mvtnorm,
+    future, 
     parallelly,
+    mvtnorm,
     dplyr,
     future.apply,
     reshape2,
@@ -49,6 +52,7 @@ Suggests:
 Depends: R (>= 4.1.0)
 LinkingTo: Rcpp, RcppArmadillo
 RoxygenNote: 7.3.2
-Packaged: 2026-03-20 00:00:00 UTC; yiqingxu
+Packaged: 2024-01-26 03:25:56 UTC; ziyil
 Encoding: UTF-8
 Config/testthat/edition: 3
+LazyData: true

R/boot.R

@@ -66,6 +66,8 @@ fect_boot <- function(
   cl = NULL,
   I,
   II,
+  cm=FALSE,
+  II.cm=NULL,
   T.on,
   T.off = NULL,
   T.on.carry = NULL,
@@ -197,6 +199,8 @@ fect_boot <- function(
         W = W,
         I = I,
         II = II,
+        cm = cm,
+        II.cm = II.cm,
         T.on = T.on,
         T.off = T.off,
         T.on.carry = T.on.carry,
@@ -311,6 +315,8 @@ fect_boot <- function(
           kappaQ.id = kappaQ.id,
           I = I,
           II = II,
+          cm = cm,
+          II.cm = II.cm,
           T.on = T.on,
           T.off = T.off,
           T.on.carry = T.on.carry,
@@ -643,6 +649,8 @@ fect_boot <- function(
           W = W,
           I = I,
           II = II,
+          cm = cm,
+          II.cm = II.cm,
           T.on = T.on,
           T.off = T.off,
           T.on.carry = T.on.carry,
@@ -704,7 +712,11 @@ fect_boot <- function(
           att.carryover.W = NA,
           balance.avg.att = NA,
           balance.time = NA,
-          group.output = list()
+          group.output = list(),
+          eff = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          D = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          I = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          boot.id = NULL
         )
         return(boot0)
       } else {
@@ -1024,6 +1036,7 @@ fect_boot <- function(
       D.boot <- out$D[, id.boot, drop = FALSE]
       I.boot <- out$I[, id.boot, drop = FALSE]
       II.boot <- out$II[, id.boot, drop = FALSE]
+      II.cm.boot <- out$II.cm[, id.boot, drop = FALSE]
       W.boot <- NULL
       if (!is.null(W)) {
         W.boot <- NULL
@@ -1103,7 +1116,11 @@ fect_boot <- function(
           count.off.W = NA,
           time.off.W = NA,
           att.carryover.W = NA,
-          group.output = list()
+          group.output = list(),
+          eff = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          D = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          I = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          boot.id = NULL
         )
         return(boot0)
       } else {
@@ -1161,6 +1178,8 @@ fect_boot <- function(
             W = W,
             I = I,
             II = II,
+            cm = cm,
+            II.cm = II.cm,
             T.on = T.on,
             T.off = T.off,
             T.on.carry = T.on.carry,
@@ -1264,7 +1283,11 @@ fect_boot <- function(
           count.off.W = NA,
           time.off.W = NA,
           att.carryover.W = NA,
-          group.output = list()
+          group.output = list(),
+          eff = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          D = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          I = if (keep.sims) matrix(NA_real_, TT, N) else NULL,
+          boot.id = NULL
         )
         return(boot0)
       } else {
@@ -1403,7 +1426,11 @@ fect_boot <- function(
           count.off.W = NA,
           time.off.W = NA,
           att.carryover.W = NA,
-          group.out = list()
+          group.out = list(),
+          eff = matrix(NA_real_, TT, length(boot.id)),
+          D = matrix(NA_real_, TT, length(boot.id)),
+          I = matrix(NA_real_, TT, length(boot.id)),
+          boot.id = boot.id
         )
         return(boot0)
       } else {
@@ -1470,6 +1497,8 @@ fect_boot <- function(
               W = W.boot,
               I = I.boot,
               II = II[, boot.id],
+              cm = cm,
+              II.cm = II.cm[, boot.id],
               T.on = T.on[, boot.id],
               T.off = T.off.boot,
               T.on.carry = T.on.carry[, boot.id],
@@ -1625,7 +1654,11 @@ fect_boot <- function(
             att.off.W = NA,
             count.off.W = NA,
             time.off.W = NA,
-            att.carryover.W = NA
+            att.carryover.W = NA,
+            eff = matrix(NA_real_, TT, length(boot.id)),
+            D = matrix(NA_real_, TT, length(boot.id)),
+            I = matrix(NA_real_, TT, length(boot.id)),
+            boot.id = boot.id
           )
           return(boot0)
         } else {
@@ -2145,6 +2178,12 @@ fect_boot <- function(
         }
       }
     }
+    if (keep.sims) {
+      eff.boot <- eff.boot[, , -boot.rm, drop = FALSE]
+      D.boot <- D.boot[, , -boot.rm, drop = FALSE]
+      I.boot <- I.boot[, , -boot.rm, drop = FALSE]
+      colnames.boot <- colnames.boot[-boot.rm]
+    }
   }
   if (dis) {
     message(dim(att.boot)[2], " runs\n", sep = "")

R/cfe.R

@@ -17,6 +17,8 @@ fect_cfe <- function(
     kappaQ.id = NULL,
     I,
     II,
+    cm=FALSE,
+    II.cm=NULL,
     T.on,
     T.off = NULL,
     T.on.carry = NULL,

R/default.R

@@ -90,7 +90,8 @@ fect <- function(
     permute = FALSE, ## permutation test
     m = 2, ## block length
     normalize = FALSE, # accelerate option
-    keep.sims = FALSE # keep individual bootstrap/jackknife simulations
+    keep.sims = FALSE, # keep individual bootstrap/jackknife simulations
+    cm = FALSE # causal moderation
 ) {
     UseMethod("fect")
 }
@@ -163,7 +164,8 @@ fect.formula <- function(
     permute = FALSE, ## permutation test
     m = 2, ## block length
     normalize = FALSE,
-    keep.sims = FALSE
+    keep.sims = FALSE,
+    cm = FALSE
 ) {
     ## parsing
     varnames <- all.vars(formula)
@@ -266,7 +268,8 @@ fect.formula <- function(
         permute = permute,
         m = m,
         normalize = normalize,
-        keep.sims = keep.sims
+        keep.sims = keep.sims,
+        cm = cm
     )
 
     out$call <- match.call()
@@ -343,7 +346,8 @@ fect.default <- function(
     permute = FALSE, ## permutation test
     m = 2, ## block length
     normalize = FALSE,
-    keep.sims = FALSE
+    keep.sims = FALSE,
+    cm=FALSE
 ) {
     ## -------------------------------##
     ## Checking Parameters
@@ -381,6 +385,11 @@ fect.default <- function(
     id <- index[1]
     time <- index[2]
 
+
+    if (cm == TRUE & ! method %in% c("fe", "ife")) {
+        stop("\"cm\" option is only available for the \"fe\" and \"ife\" methods.")
+    }
+
     if (se == 1) {
         if (!vartype %in% c("bootstrap", "jackknife", "parametric")) {
             stop("\"vartype\" option misspecified.")
@@ -1483,6 +1492,12 @@ fect.default <- function(
     II <- I
     II[which(D == 1)] <- 0 ## regard treated values as missing
 
+    II.cm <- matrix(0, TT, N)
+    II.cm[which(D == 1)] <- 1
+    II.cm[is.nan(Y.ind)] <- 0
+
+    # Unbalance Check
+    ## 1. remove units that have too control status
     T0 <- apply(II, 2, sum)
     T0.min <- min(T0)
 
@@ -1571,6 +1586,7 @@ fect.default <- function(
         I <- as.matrix(I[, -rm.id]) ## after removing
         I.D <- as.matrix(I.D[, -rm.id])
         II <- as.matrix(II[, -rm.id])
+        II.cm <- as.matrix(II.cm[, -rm.id])
         if (!is.null(group)) {
             G <- as.matrix(G[, -rm.id])
         }
@@ -1610,6 +1626,7 @@ fect.default <- function(
         I <- I[-which(I.use == 0), ] ## remove that period
         I.D <- I.D[-which(I.use == 0), ]
         II <- II[-which(I.use == 0), ] ## remove that period
+        II.cm <- II.cm[-which(I.use == 0), ] ## remove that period
         D <- D[-which(I.use == 0), ] ## remove that period
         Y <- Y[-which(I.use == 0), ] ## remove that period
 
@@ -2048,6 +2065,8 @@ fect.default <- function(
                     W = W,
                     I = I,
                     II = II,
+                    cm=cm,
+                    II.cm = II.cm,
                     T.on = T.on,
                     T.off = T.off,
                     r = r,
@@ -2251,6 +2270,8 @@ fect.default <- function(
             W = W,
             I = I,
             II = II,
+            cm = cm,
+            II.cm = II.cm,
             T.on = T.on,
             T.off = T.off,
             T.on.carry = T.on.carry,

R/effect.R

@@ -82,15 +82,33 @@ effect <- function(x, ## a fect object
   if (is.null(x$est.avg)) {
     cat("No uncertainty estimates.")
   } else {
-    # Perform bootstrap analysis
-    nboots <- length(x$att.avg.boot)
+    # Determine nboots from the 3D arrays
+    if (!is.null(x$eff.boot) && length(dim(x$eff.boot)) == 3) {
+      nboots <- dim(x$eff.boot)[3]
+    } else {
+      stop("eff.boot is not a valid 3D array. Ensure keep.sims = TRUE in fect().")
+    }
+    if (nboots == 0) {
+      stop("All bootstrap iterations failed. Cannot compute cumulative effects.")
+    }
+
+    # Validate D.boot and I.boot exist with matching dimensions
+    has.D.boot <- !is.null(x$D.boot) && length(dim(x$D.boot)) == 3 && dim(x$D.boot)[3] >= nboots
+    has.I.boot <- !is.null(x$I.boot) && length(dim(x$I.boot)) == 3 && dim(x$I.boot)[3] >= nboots
+
     catt.boot <- matrix(NA, period[2] - period[1] + 1, nboots)
 
     # Calculate treatment effect for each bootstrap sample
     for (i in 1:nboots) {
       # Extract bootstrap matrices
-      D.boot <- x$D.boot[, , i]
-      I.boot <- x$I.boot[, , i]
+      if (has.D.boot) {
+        D.boot <- x$D.boot[, , i]
+        I.boot <- x$I.boot[, , i]
+      } else {
+        # Fallback: use original D.dat and I.dat (less accurate but prevents crash)
+        D.boot <- x$D.dat
+        I.boot <- x$I.dat
+      }
       eff.boot <- x$eff.boot[, , i]
 
       # Select treated units in bootstrap sample