Copula endogeneity correction of Haschka (2025) (IMA)

Copula correction exploiting between-regressor correlation, following Haschka (2025). Every explanatory variable – endogenous and exogenous – is first transformed to a normal score; each endogenous normal score is then regressed, without an intercept, on the exogenous normal scores, and the residual of that regression is added to the structural regression as the copula control function. This is the same two-stage procedure as 2sCOPE (CopReg2sCOPE(), Yang, Qian & Xie 2025) with one difference: the first-stage (auxiliary) regression here carries no intercept, which is what makes the slope of P* on X* equal to the Pearson correlation used in Haschka's derivation of rho. Without an exogenous part in formula the first stage has nothing to project on, so the call is redirected to CopRegPG() with a warning.

Usage

CopRegIMA(
  formula,
  data,
  cdf = "rank.n",
  ties = "max",
  nboots = 199,
  subset = NULL,
  contrasts = NULL,
  parallel = FALSE,
  ncores = NULL,
  verbose = interactive()
)

Arguments

formula: A two-part formula of the form y ~ endog_1 + endog_2 + ... | exog_1 + exog_2 + .... Position decides: terms written before the | are endogenous and receive a copula correction, terms written after it are exogenous and receive none. Everything lm() accepts is allowed: transformations, e.g. log(price); interactions, e.g. price:feat; polynomials, e.g. I(x^2); factors, e.g. as.factor(store); and -1 to drop the intercept. Endogenous regressors must be numeric. Missing values are handled by na.omit. Writing an interaction or a power of an endogenous regressor before the | gives it its own copula term and produces a warning, because such a term is strongly collinear with the regressor it is built from. Without an exogenous part (no |) the call is redirected to CopRegPG().
data: A data.frame containing the variables referenced in formula.
cdf: Character string naming the estimator of the marginal CDF used in the normal-score transformation. Haschka (2025) leaves this choice open, so the default follows 2sCOPE: "rank.n", the rescaled ECDF of Qian, Koschmann & Xie (2025, Eq. 9). Other options are "kde.silverman" (Park & Gupta 2012), "kde.cv" (Li, Li & Racine 2017), "kde.plugin" (Polansky & Baker 2000), "ecdf.fixed" (Becker, Proksch & Ringle 2022), "ecdf.adj" (Liengaard et al. 2025), and "rank.n1" (Breitung, Mayer & Wied 2024).
ties: Character string, "max" (default) for the counting function or "average" for midranks. Relevant here because the exogenous regressors are transformed to normal scores as well.
nboots: Number of bootstrap replicates used to compute the standard errors. Defaults to 199.
subset: An optional logical vector selecting the rows to use, as in lm(). Defaults to NULL (use every row).
contrasts: An optional list of contrasts for factor variables, passed on to model.matrix(). Defaults to NULL.
parallel: FALSE (default), TRUE, "multicore", or "snow" to spread the bootstrap replicates across worker processes.
ncores: Number of worker processes to use when parallel is not FALSE. NULL (default) uses one less than the number of detected cores.
verbose: Logical; show a progress message and progress bar while the bootstrap runs. Defaults to interactive().

Value

An object of class "copreg": a list including coefficients and std.error (point estimates and bootstrap standard errors for the structural and copula coefficients), vcov (the bootstrap covariance matrix of coefficients), rho and rho.se (the endogeneity measure corr(structural residual, copula data) for each endogenous regressor and its bootstrap standard error), residuals and fitted.values (structural model, copula terms excluded), residuals.augmented and fitted.augmented (including the copula terms), boot (raw bootstrap coefficient draws), ols.coefficients and ols.std.error (the uncorrected OLS benchmark from the same bootstrap resamples), icon (the resulting standard error inflation), diagnostics (identification diagnostics), and call, method, cdf, ties. Methods exist for print(), summary(), coef(), vcov(), nobs(), residuals(), fitted(), predict(), confint(), formula(), and update().

References

Haschka, R. E. (2025). Robustness of copula-correction models in causal analysis: Exploiting between-regressor correlation. IMA Journal of Management Mathematics 36(1), 161-180.

Qian, Y., A. Koschmann, and H. Xie (2025). A practical guide to endogeneity correction using copulas. Journal of Marketing.

Yang, F., Y. Qian, and H. Xie (2025). Addressing endogeneity using a two-stage copula generated regressor approach. Journal of Marketing Research 62(4), 601-623.

Examples

set.seed(1)
n <- 200
w <- rnorm(n)
x <- rnorm(n)
e <- rexp(n, rate = 1) - 1               # nonnormal endogeneity driver
z <- 0.5 * w + 0.5 * x + e               # endogenous regressor
y <- 1 + 2 * z + 0.5 * x + 0.5 * e + rnorm(n, sd = 0.5)
dat <- data.frame(y = y, z = z, x = x, w = w)

fit <- CopRegIMA(y ~ z | x + w, data = dat, nboots = 25)
summary(fit)
#> 
#> Copula endogeneity correction: IMA (Haschka 2025) 
#> 
#> Call:
#> CopRegIMA(formula = y ~ z | x + w, data = dat, nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>      Min       1Q   Median       3Q      Max 
#> -1.41619 -0.43098 -0.05137  0.39325  1.37631 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>      Min       1Q   Median       3Q      Max 
#> -1.40272 -0.42622 -0.04259  0.40404  1.31991 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  0.95156    0.04273  22.268  < 2e-16 ***
#> z            2.51239    0.13621  18.445  < 2e-16 ***
#> x            0.23027    0.07868   2.927  0.00343 ** 
#> w           -0.25856    0.08348  -3.097  0.00195 ** 
#> z_cop       -0.02389    0.19240  -0.124  0.90118    
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Endogeneity: rho(P*, xi*) is the correlation between the normal score 
#>   of an endogenous regressor and that of the structural error, xi* = xi / sigma.
#>              Estimate Std. Error z value Pr(>|z|)
#> rho(z*, xi*) -0.02696    0.20267  -0.133    0.894
#> 
#> Fit, on 195 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.5520    0.5524    
#> R-squared               0.9734    0.9734    
#> Adjusted R-squared      0.9729    0.9728    
#>   sigma above is the standard error of the structural model, the one 
#>   entering xi* = xi / sigma.
#> Standard errors from 25 bootstrap replicates; cdf = "rank.n", ties = "max".
#> Pr(>|z|) in both tables: Wald test using the normal approximation,
#>   z = Estimate / Std. Error, with the bootstrap standard error.
#>   See confint(object, type = "percentile") for bootstrap percentile intervals.
#> 
#> --- Identification diagnostics ------------------------------------
#> 
#> Non-normality of the endogenous regressors (small p = non-normal, good):
#>      AD     AD p    KS p
#> z 5.079 1.35e-12 0.01434
#> 
#> Correlation of the copula terms with the exogenous regressors
#> (this estimator does not assume it is zero; it projects the correlation
#>  out in the first stage):
#>   max |corr| with  p (Holm) joint R2   joint p
#> z     0.4278    w 2.774e-10   0.3369 2.674e-18
#>  full matrix in summary(object)$diagnostics$exog.correlation.matrix
#> 
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#>       corr(P, C)  omega
#> z_cop     0.8238 0.1111
#>