Why interactive fixed effects?

Standard two-way fixed effects (TWFE) controls for unobserved heterogeneity through an additive structure:

\[y_{it} = \alpha_i + \xi_t + X_{it}'\beta + u_{it}\]

The unit effect \(\alpha_i\) absorbs any time-invariant unit characteristic; the time effect \(\xi_t\) absorbs any unit-invariant time shock. This works well when heterogeneous responses to common shocks are negligible.

Many empirical settings violate this assumption. Consider:

• Globalisation shocks: countries are exposed differently depending on their trade openness (a unit characteristic that interacts with a common world-demand shock).
• Monetary policy: firms with different financial leverage respond asymmetrically to a common interest-rate change.
• COVID-19: regional outcomes depend on a national shock weighted by pre-existing health-system capacity.

In each case, the unobserved confounder has the product form \(\lambda_i' F_t\): a unit-specific loading \(\lambda_i\) multiplied by a common factor \(F_t\). If \(X_{it}\) is correlated with \(\lambda_i' F_t\), TWFE is inconsistent even as \(N, T \to \infty\).

The IFE model

Bai (2009) proposes replacing the additive error structure with:

\[y_{it} = \alpha_i + \xi_t + X_{it}'\beta + \lambda_i' F_t + u_{it} \qquad (1)\]

where:

• \(\beta \in \mathbb{R}^p\): structural parameters of interest,
• \(F_t \in \mathbb{R}^r\): common factors (\(T \times r\) matrix \(\hat{F}\), normalised \(\hat{F}'\hat{F}/T = I_r\)),
• \(\lambda_i \in \mathbb{R}^r\): unit-specific factor loadings (\(N \times r\) matrix \(\hat{\Lambda}\)),
• \(u_{it}\): idiosyncratic error with \(\mathbb{E}[u_{it} \mid X, F, \Lambda] = 0\).

The \(r\) interactive terms \(\lambda_{ik} F_{kt}\) generalise both additive unit effects (\(r=1\), \(F_t = 1\)) and additive time effects (\(r=1\), \(\lambda_i = 1\)). With \(r\) unrestricted, (1) subsumes TWFE as the special case \(r = 0\).

What xtife provides

xtife is a pure base-R package (no compiled code, no external dependencies) offering:

Data

The package includes the cigar dataset: 46 US states \(\times\) 30 years (1963–1992) of cigarette sales and prices from Baltagi (1995).

library(xtife)
data(cigar)
dim(cigar)
#> [1] 1380    9
head(cigar[, c("state", "year", "sales", "price")], 3)
#>   state year sales price
#> 1     1   63  93.9  28.6
#> 2     1   64  95.4  29.8
#> 3     1   65  98.5  29.8

The panel is balanced: every state appears in every year.

Quick start

fit <- ife(sales ~ price,
           data  = cigar,
           index = c("state", "year"),
           r     = 2,
           force = "two-way",
           se    = "cluster")
print(fit)
#> 
#> Interactive Fixed Effects (Bai 2009, Econometrica)
#> -------------------------------------------------------
#> Panel    : N = 46 units, T = 30 periods
#> Factors  : r = 2
#> Force    : two-way fixed effects
#> SE type  : cluster (by state)
#> Outcome  : sales
#> -------------------------------------------------------
#>       Estimate Std.Error t.value Pr(>|t|)             95% CI    
#> price  -0.5242    0.0802 -6.5333   0.0000 [-0.6816, -0.3667] ***
#> ---
#> Signif. codes: *** <0.01  ** <0.05  * <0.1
#> -------------------------------------------------------
#> sigma^2  : 18.456076 | df = 1157
#> Converged: YES | Iterations: 10
#> -------------------------------------------------------
#> Factor selection criteria at r = 2 [IC1-3: Bai & Ng 2002; IC_bic/PC: Bai 2009]:
#>   IC1 =   3.2347  |  IC2 =   3.2900  |  IC3 =   3.1421
#>   PC  =  36.8093  |  IC (BIC-style) =   4.2661
#>   -> Run ife_select_r() to compare criteria across r = 0, 1, ..., r_max
#>      and identify the IC-minimising number of factors.

The force argument controls additive fixed effects: "none" | "unit" | "time" | "two-way" (default).

Key output fields:

fit$coef          # named coefficient vector
#>      price 
#> -0.5241574
fit$se            # standard errors
#>     price 
#> 0.0802281
fit$ci            # 95% confidence intervals (matrix: coef x 2)
#>         CI.lower   CI.upper
#> price -0.6815663 -0.3667486
fit$converged     # TRUE if outer loop converged
#> [1] TRUE
fit$n_iter        # number of outer iterations
#> [1] 10
dim(fit$F_hat)    # T x r estimated factors
#> [1] 30  2
dim(fit$Lambda_hat) # N x r estimated loadings
#> [1] 46  2

Estimation algorithm

The Bai (2009) estimator alternates between two steps until convergence:

Step 1 — Factor extraction. Given \(\hat\beta\), form \(W_{it} = \ddot y_{it} - \ddot X_{it}'\hat\beta\) (residuals after additive FE demeaning). The top-\(r\) eigenvectors of \(WW'/NT\) (scaled to satisfy \(\hat F' \hat F / T = I_r\)) give \(\hat F\). The loadings follow by regression: \(\hat\Lambda = W'\hat F / T\).

Step 2 — Coefficient update. Project the demeaned regressors onto the orthogonal complement of \(\hat F\):

\[\tilde X_{it} = \ddot X_{it} - \hat F_t (\hat F' \hat F)^{-1} \hat F' \ddot X_i\]

Then \(\hat\beta = (\tilde X'\tilde X)^{-1}\tilde X'\ddot y\) (Frisch-Waugh-Lovell).

The outer loop converges when \(\max|\hat\beta^{(k+1)} - \hat\beta^{(k)}| < \text{tol}\) (default \(10^{-9}\)).

Standard error types

fit_std <- ife(sales ~ price, data = cigar,
               index = c("state", "year"), r = 2, se = "standard")
fit_rob <- ife(sales ~ price, data = cigar,
               index = c("state", "year"), r = 2, se = "robust")
fit_cl  <- ife(sales ~ price, data = cigar,
               index = c("state", "year"), r = 2, se = "cluster")

data.frame(
  se_type = c("standard", "robust (HC1)", "cluster"),
  coef    = round(c(fit_std$coef, fit_rob$coef, fit_cl$coef), 5),
  se      = round(c(fit_std$se,   fit_rob$se,   fit_cl$se),   5),
  t_stat  = round(c(fit_std$tstat, fit_rob$tstat, fit_cl$tstat), 3)
)
#>        se_type     coef      se  t_stat
#> 1     standard -0.52416 0.03987 -13.146
#> 2 robust (HC1) -0.52416 0.04510 -11.623
#> 3      cluster -0.52416 0.08023  -6.533

All three use the sandwich formula \(\widehat{\mathrm{Var}}(\hat\beta) = A^{-1} B A^{-1}\) where \(A = \tilde X'\tilde X\) and:

Rule of thumb: use "cluster" for most macro panels; "robust" when clustering by unit is not appropriate; "standard" only as a benchmark.

Comparison with TWFE

Setting \(r = 0\) reduces ife() to the standard TWFE estimator:

fit0 <- ife(sales ~ price, data = cigar,
            index = c("state", "year"), r = 0, force = "two-way")

# Manual within-transformation
cigar$y_dm <- cigar$sales - ave(cigar$sales, cigar$state) -
                             ave(cigar$sales, cigar$year) + mean(cigar$sales)
cigar$x_dm <- cigar$price - ave(cigar$price, cigar$state) -
                             ave(cigar$price, cigar$year) + mean(cigar$price)
lm0 <- lm(y_dm ~ x_dm - 1, data = cigar)

cat(sprintf("ife (r=0): %.7f\n", fit0$coef["price"]))
#> ife (r=0): -1.0847117
cat(sprintf("lm TWFE:  %.7f\n", coef(lm0)["x_dm"]))
#> lm TWFE:  -1.0847117
cat(sprintf("diff:     %.2e\n",
            abs(fit0$coef["price"] - coef(lm0)["x_dm"])))
#> diff:     4.44e-16

The difference is at machine precision, confirming algebraic equivalence. TWFE gives \(\hat\beta \approx -0.38\); IFE with \(r=2\) gives \(\approx -0.52\): the gap reflects the OVB from ignoring interactive confounders.

Factor number selection

ife_select_r() fits models for \(r = 0, \ldots, r_{\max}\) and reports five information criteria.

# Not run during build (takes ~20 s); set verbose=FALSE in scripts
sel <- ife_select_r(sales ~ price, data = cigar,
                    index = c("state", "year"),
                    r_max = 6, force = "two-way", verbose = FALSE)
attr(sel, "suggested")   # named vector: r chosen by each IC

The five criteria:

Recommendation: IC(BIC) for small-to-moderate panels (\(\min(N,T) < 60\)); IC1 or IC2 for large panels.

Asymptotic bias correction (balanced panel)

fit_bc <- ife(sales ~ price, data = cigar,
              index = c("state", "year"), r = 2,
              se = "standard", bias_corr = TRUE)
cat(sprintf("Raw IFE:  %.5f\n", fit_bc$coef_raw["price"]))
#> Raw IFE:  -0.52416
cat(sprintf("Corrected:%.5f\n", fit_bc$coef["price"]))
#> Corrected:-0.52847
cat(sprintf("B_hat:    %.5f  (B_hat/N = %.5f)\n",
            fit_bc$B_hat["price"], fit_bc$B_hat["price"] / fit_bc$N))
#> B_hat:    0.19466  (B_hat/N = 0.00423)
cat(sprintf("C_hat:    %.5f  (C_hat/T = %.5f)\n",
            fit_bc$C_hat["price"], fit_bc$C_hat["price"] / fit_bc$T))
#> C_hat:    0.00254  (C_hat/T = 0.00008)

fit_dyn <- ife(sales ~ price, data = cigar,
               index = c("state", "year"), r = 2,
               se = "standard", method = "dynamic",
               bias_corr = TRUE, M1 = 1L)
cat(sprintf("Dynamic BC coef: %.5f\n", fit_dyn$coef["price"]))
#> Dynamic BC coef: -0.53161
cat(sprintf("B1/T (Nickell):  %.5f\n",
            fit_dyn$B1_hat["price"] / fit_dyn$T))
#> B1/T (Nickell):  -0.00009

Motivation

Real-world panels are rarely balanced. Survey non-response, merger events, data-collection gaps, and rotating samples all cause missing \((i,t)\) cells. Simply deleting units with incomplete records (list-wise deletion) induces selection bias; imputation by column means distorts the factor structure.

ife_unbalanced() extends the Bai (2009) interactive fixed effects estimator to genuinely unbalanced panels, following the estimation and inference theory of Su, Wang and Wang (2025), and combines three published ideas:

1. Factor analysis with missing data. The factor structure is recovered from the observed cells by an expectation-maximisation (EM) algorithm that treats the unobserved entries as latent, in the spirit of Bai (2009, Appendix B) and the missing-data factor / matrix-completion framework of Bai and Ng (2021).
2. Nuclear-norm initialisation. A nuclear-norm-regularised (soft-impute) estimator (Mazumder, Hastie and Tibshirani 2010) provides a consistent warm start for the factor structure, stabilising the alternating loop when the panel is severely unbalanced.
3. Analytical inference. Standard errors use a sandwich estimator on factor-projected regressors (the unbalanced analogue of the balanced formula of Bai 2009), with heteroskedasticity-, cluster- and autocorrelation-robust variants; an analytical incidental-parameter bias correction is provided for both strictly exogenous and predetermined (weakly exogenous) regressors. The inference and bias-correction formulae for the unbalanced case follow Su, Wang and Wang (2025), extending the bias-correction principles of Bai (2009) and Moon and Weidner (2017).

In addition, additive unit and/or time fixed effects can be removed jointly with the factors through the force argument (Section 3.2).

The unbalanced IFE model

Let \(d_{it} \in \{0,1\}\) indicate whether \((i,t)\) is observed. The model is:

\[y_{it} = X_{it}'\beta + \lambda_i' F_t + u_{it} \qquad \text{for all } (i,t) \text{ with } d_{it}=1 \qquad (2)\]

By default (force = "none") there are no additive fixed effects in (2): individual and time heterogeneity is absorbed entirely by the interactive factor structure. Additive unit and/or time effects can instead be modelled explicitly through the force argument (see below), which is the recommended route when the data contain a strong common level or trend.

Quick start with unbalanced data

# Create a 10% randomly missing version of the cigar panel
set.seed(42)
cigar_unb <- cigar[sample(nrow(cigar), size = floor(0.9 * nrow(cigar))), ]
cat(sprintf("Observations: %d / %d  (%.0f%% fill)\n",
            nrow(cigar_unb), nrow(cigar), 100 * nrow(cigar_unb) / nrow(cigar)))
#> Observations: 1242 / 1380  (90% fill)

fit_unb <- ife_unbalanced(sales ~ price,
                           data  = cigar_unb,
                           index = c("state", "year"),
                           r     = 2L,
                           se    = "cluster")
print(fit_unb)
#> 
#> Unbalanced Panel Interactive Fixed Effects
#> ---------------------------------------------------------- 
#> Panel    : N = 46 units,  TT = 30 periods (max)
#> Observed : n_obs = 1242 cells  (90.0% of N x TT)
#> Factors  : r = 2
#> SE type  : cluster-robust by unit (state)
#> Exog     : strict (exogenous regressors)
#> Init     : nnr  (nu = 46.0000)
#> Outcome  : sales
#> ---------------------------------------------------------- 
#>       Estimate Std.Error t.value Pr(>|t|)            95% CI  
#> price   0.0754    0.1108  0.6801   0.4966 [-0.1421, 0.2928]  
#> ---
#> Signif. codes: *** < 0.01  ** < 0.05  * < 0.10
#> ---------------------------------------------------------- 
#> sigma^2 = 53.1211  |  df = 1093
#> Converged: YES  |  Outer iterations: 360

Key output fields mirror those of ife():

fit_unb$coef       # estimated beta
#>      price 
#> 0.07537269
fit_unb$se         # standard errors
#>    price 
#> 0.110819
fit_unb$ci         # 95% confidence intervals
#>            2.5 %    97.5 %
#> price -0.1420693 0.2928146
fit_unb$converged  # outer-loop convergence
#> [1] TRUE
fit_unb$n_obs      # number of observed (i,t) pairs
#> [1] 1242
fit_unb$N          # number of units
#> [1] 46
fit_unb$TT         # number of distinct time periods
#> [1] 30
dim(fit_unb$F_hat)       # TT x r
#> [1] 30  2
dim(fit_unb$Lambda_hat)  # N  x r
#> [1] 46  2

Estimation algorithm

# Both converge to the same coefficient for mild unbalancedness
fit_ols <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, init = "ols")
fit_nnr <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, init = "nnr")
cat(sprintf("OLS init: coef = %.5f,  %d outer iterations\n",
            fit_ols$coef["price"], fit_ols$n_iter))
#> OLS init: coef = 0.07537,  362 outer iterations
cat(sprintf("NNR init: coef = %.5f,  %d outer iterations\n",
            fit_nnr$coef["price"], fit_nnr$n_iter))
#> NNR init: coef = 0.07537,  360 outer iterations

Analytical standard errors

The sandwich variance of \(\hat\beta\) follows Su, Wang and Wang (2025) and requires projecting each regressor out of both the factor space and the loading space, giving doubly-projected regressors \(\hat x_{itk}\) (the unbalanced analogue of the Frisch–Waugh projection used in the balanced case):

\[\hat x_{itk} = x_{itk} - \hat\delta_{ki}'\hat F_t - \hat\omega_{kt}'\hat\lambda_i\]

where \(\hat\delta_{ki}\) and \(\hat\omega_{kt}\) solve, over the observed cells, the alternating least-squares system

\[\hat\delta_{ki} = \Bigl(\sum_t d_{it}\hat F_t\hat F_t'\Bigr)^{-1} \!\sum_t d_{it}\hat F_t(x_{itk} - \hat\lambda_i'\hat\omega_{kt})\]

\[\hat\omega_{kt} = \Bigl(\sum_i d_{it}\hat\lambda_i\hat\lambda_i'\Bigr)^{-1} \!\sum_i d_{it}\hat\lambda_i(x_{itk} - \hat F_t'\hat\delta_{ki})\]

The sandwich variance is then

\[\widehat{\mathrm{Var}}(\hat\beta) = \hat W_X^{-1} \hat\Omega_X \hat W_X^{-1}\]

with \(\hat W_X = \frac{1}{NT}\sum_{it} d_{it}\hat x_{it}\hat x_{it}'\) and \(\hat\Omega_X\) depending on the se argument:

The cluster-robust variant follows the within-unit clustering of Arellano (1987) and Cameron, Gelbach and Miller (2011); the "hac" variant uses the Bartlett kernel of Newey and West (1987), \(\Gamma(u) = (1-u)\mathbf{1}\{u \leq 1\}\), with bandwidth defaulting to \(L_T = \lfloor 2T^{1/5}\rfloor\).

fit_std <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, se = "standard")
fit_rob <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, se = "robust")
fit_cl  <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, se = "cluster")
fit_hac <- ife_unbalanced(sales ~ price, data = cigar_unb,
                           index = c("state", "year"), r = 2, se = "hac")

data.frame(
  se_type = c("standard", "robust", "cluster", "hac"),
  coef    = round(c(fit_std$coef, fit_rob$coef, fit_cl$coef, fit_hac$coef), 5),
  se      = round(c(fit_std$se,   fit_rob$se,   fit_cl$se,   fit_hac$se),   5)
)
#>    se_type    coef      se
#> 1 standard 0.07537 0.02503
#> 2   robust 0.07537 0.03657
#> 3  cluster 0.07537 0.11082
#> 4      hac 0.07537 0.05090

When to use "hac": when the errors \(u_{it}\) are serially correlated (e.g., in dynamic models or when residuals show autocorrelation). For strictly exogenous regressors with i.i.d. errors, "cluster" is sufficient.

Factor number selection

For unbalanced panels, ife_select_r_unb() uses the singular value thresholding (SVT) rule of Su, Wang and Wang (2025) applied to the nuclear-norm regularised matrix \(\hat\Theta^{(0)}\) — a missing-data counterpart of the eigenvalue-ratio / information-criterion rules of Bai and Ng (2002):

\[\hat r = \sum_{s=1}^{N \wedge T} \mathbf{1}\!\left\{ \sigma_s\!\left(\frac{\hat\Theta^{(0)}}{\sqrt{NT}}\right) \geq c_f\sqrt{c_{NT}^{-1/4}\,\sigma_1\!\left(\frac{\hat\Theta^{(0)}}{\sqrt{NT}}\right)} \right\}\]

where \(c_{NT} = \min(\sqrt{N}, \sqrt{T})\) and \(c_f = 0.6\) is the default threshold constant, which performs well for moderate degrees of unbalancedness.

# Not run during build (NNR penalty selection takes ~10 s)
sel_unb <- ife_select_r_unb(sales ~ price, data = cigar_unb,
                              index = c("state", "year"), verbose = FALSE)
cat(sprintf("SVT selects r_hat = %d\n", sel_unb$r_hat))

The function returns (invisibly):

• r_hat: the selected number of factors (at least 1);
• sv: the normalised singular values used in the threshold comparison;
• threshold: the computed threshold value;
• nu_used: the NNR penalty selected from the grid of Su, Wang and Wang (2025, footnote 12).

Analytical bias correction

Even with the correct \(r\), the unbalanced IFE estimator carries an incidental-parameter bias of order \(O(N^{-1}) + O(T^{-1})\), analogous to the \(B/N + C/T\) bias of the balanced estimator (Bai 2009) and the dynamic bias of Moon and Weidner (2017). Setting bias_corr = TRUE applies the analytical correction of Su, Wang and Wang (2025), \(\hat\beta^{abc} = \hat\beta - (NT)^{-1/2}\hat W_X^{-1}\hat b\), where, for strictly exogenous regressors, \(\hat b = \hat b_3 + \hat b_4 + \hat b_5 + \hat b_6\) with

\[\hat b_{3k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\delta_{ki}'[\hat L_{ff}^{-1}]_t\hat\lambda_i\]

\[\hat b_{4k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\omega_{kt}'[\hat L_{\lambda\lambda}^{-1}]_i\hat F_t\]

\[\hat b_{5k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat\lambda_i'\hat\Xi_{kt}\hat\lambda_i\]

\[\hat b_{6k} = \frac{1}{\sqrt{NT}}\sum_{i,t} d_{it}\hat u_{it}^2 \hat F_t'\hat\Delta_{ki}\hat F_t\]

Here \([\hat L_{ff}^{-1}]_t = (\sum_i d_{it}\hat\lambda_i\hat\lambda_i')^{-1}\), \([\hat L_{\lambda\lambda}^{-1}]_i = (\sum_t d_{it}\hat F_t\hat F_t')^{-1}\), and \(\hat\Delta_{ki}\), \(\hat\Xi_{kt}\) are \(r \times r\) auxiliary matrices formed from \(\hat\delta\), \(\hat\omega\) and the above inverses; \(\hat u_{it}\) are the model residuals. Full expressions are given in the .ife_bias_unb() source documentation.

fit_bc_strict <- ife_unbalanced(sales ~ price, data = cigar_unb,
                                 index = c("state", "year"), r = 2,
                                 se = "standard", bias_corr = TRUE,
                                 exog = "strict")
cat(sprintf("Raw beta:       %.5f\n", fit_bc_strict$coef_raw["price"]))
#> Raw beta:       0.07537
cat(sprintf("Corrected beta: %.5f\n", fit_bc_strict$coef["price"]))
#> Corrected beta: 0.08183
cat(sprintf("b_hat (b3+...+b6): %.5f\n", fit_bc_strict$b_hat["price"]))
#> b_hat (b3+...+b6): NA

# Not run: exog="weak" adds the b2 dynamic term
# (requires HAC SE for valid inference in dynamic models)
fit_bc_weak <- ife_unbalanced(y ~ lag_y + x, data = df_dynamic,
                               index = c("unit", "time"), r = 2,
                               se = "hac", bias_corr = TRUE,
                               exog = "weak")

Balanced vs unbalanced on the same data

As a sanity check, ife_unbalanced() on a complete (balanced) panel should agree with ife() under the same force specification. We use a simulated mean-zero factor panel (on data with a large level, such as cigar, force = "none" is ill-conditioned for both estimators – the factors must absorb the level – so use force = "two-way" there):

set.seed(99)
N <- 30; TT <- 22
F0 <- matrix(rnorm(TT * 2), TT, 2); L0 <- matrix(rnorm(N * 2), N, 2)
X  <- matrix(rnorm(N * TT), TT, N)
Y  <- F0 %*% t(L0) + 0.7 * X + matrix(rnorm(N * TT), TT, N)
sim <- data.frame(unit = rep(seq_len(N), each = TT),
                  time = rep(seq_len(TT), times = N),
                  Y = as.vector(Y), X = as.vector(X))

fit_bal <- ife(Y ~ X, data = sim, index = c("unit", "time"),
               r = 2, force = "none", se = "standard")

fit_unb2 <- ife_unbalanced(Y ~ X, data = sim, index = c("unit", "time"),
                            r = 2, force = "none", se = "standard")

cat(sprintf("ife (balanced):     %.6f\n", fit_bal$coef["X"]))
#> ife (balanced):     0.635347
cat(sprintf("ife_unbalanced:     %.6f\n", fit_unb2$coef["X"]))
#> ife_unbalanced:     0.632448
cat(sprintf("Difference:         %.2e\n",
            abs(fit_bal$coef["X"] - fit_unb2$coef["X"])))
#> Difference:         2.90e-03

The two estimators agree closely (they follow different iteration paths, so agreement is to roughly \(10^{-2}\) rather than machine precision).

ife() — balanced panel

ife_unbalanced() — unbalanced panel