Copula endogeneity correction of Yang, Qian and Xie (2025) (2sCOPE)

Fits the two-stage copula endogeneity correction of Yang, Qian and Xie (2025). In the first stage the endogenous and exogenous regressors are each transformed to normal scores via their empirical CDFs (P* and W*), every endogenous P*_k is regressed on W* (with an intercept, following the paper's implementation rather than its derivation) and the residual is kept. In the second stage those residuals are added to the structural regression as generated regressors, so the correction is orthogonal to the exogenous regressors by construction. Without an exogenous part the first stage has nothing to project on and the call is redirected to CopRegPG() with a warning. Standard errors are obtained by a nonparametric bootstrap.

Usage

CopReg2sCOPE(
  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 y ~ endogenous_1 + endogenous_2 + ... | exogenous_1 + exogenous_2 + .... Position decides: a term written before the | is endogenous and receives a copula correction, the same term written after it is exogenous and receives none. Everything lm() accepts is allowed. Endogenous regressors must be numeric. Without an exogenous part 2sCOPE is identical to Park and Gupta (2012) and the call is redirected there with a warning.
data: A data.frame containing the variables referenced in formula.
cdf: CDF estimator used for the copula transformation of both the endogenous and the exogenous regressors. One of "kde.silverman" (integral of a Gaussian kernel density with Silverman's bandwidth), "kde.cv" (same, cross-validated bandwidth), "kde.plugin" (same, plug-in bandwidth), "ecdf.fixed" (empirical CDF with a replaced boundary), "ecdf.adj" (adjusted empirical CDF), "rank.n" (default; rescaled empirical CDF, Rank/n with a correction, the algorithm of Equation 9 in Qian, Koschmann and Xie (2025), which is what 2sCOPE uses), or "rank.n1" (Rank/(n+1)).
ties: How tied values are ranked when computing the empirical CDF: "max" (default, the counting function) or "average" (midranks). This matters more here than for Park and Gupta, because the exogenous regressors are transformed as well and dummies and other discrete controls carry many ties.
nboots: Number of bootstrap replicates used for the standard errors.
subset: An optional logical vector specifying a subset of observations to use, as in lm(). NULL (default) uses all rows.
contrasts: An optional list of contrasts passed to model.matrix(). NULL (default) uses the defaults in options().
parallel: FALSE (default), TRUE, "multicore" or "snow" to spread the bootstrap over ncores workers. No seed is set internally either way; call set.seed() beforehand for reproducible standard errors.
ncores: Number of workers used when parallel is not FALSE. NULL (default) uses one less than the number of detected cores.
verbose: Show a progress bar during the bootstrap. Defaults to interactive().

Value

An object of class "copreg", a list including the structural coefficients and their bootstrap standard errors, the matched call, the uncorrected OLS benchmark from the same bootstrap resamples (ols.coefficients, ols.std.error), the raw bootstrap draws (boot), identification diagnostics (diagnostics), and the standard error inflation relative to OLS (icon). Methods exist for coef(), vcov(), confint(), residuals(), fitted(), predict(), print() and summary(). Use validity() on the result for the decision tree of Yang, Qian and Xie (2025, Figure 2) and the ICON statistic.

References

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.

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

Examples

set.seed(42)
n <- 300
w <- rnorm(n)
u <- rnorm(n)
p <- rchisq(n, df = 3) + 0.3 * u   # non-normal endogenous regressor
y <- 1 + 2 * p + 0.5 * w + u
dat <- data.frame(y = y, p = p, w = w)

fit <- CopReg2sCOPE(y ~ p | w, data = dat, nboots = 25)
summary(fit)
#> 
#> Copula endogeneity correction: 2sCOPE (Yang, Qian & Xie 2025) 
#> 
#> Call:
#> CopReg2sCOPE(formula = y ~ p | w, data = dat, nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>      Min       1Q   Median       3Q      Max 
#> -2.69286 -0.65018 -0.09691  0.62174  3.07035 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>     Min      1Q  Median      3Q     Max 
#> -2.5558 -0.6215 -0.0672  0.6113  3.3038 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.36456    0.41410   3.295 0.000983 ***
#> p            1.86280    0.15360  12.127  < 2e-16 ***
#> w            0.54516    0.06464   8.434  < 2e-16 ***
#> p_cop        0.50003    0.30698   1.629 0.103337    
#> ---
#> 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(p*, xi*)   0.4604     0.1656    2.78  0.00544 **
#> ---
#> Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
#> 
#> Fit, on 296 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.9646    1.0867    
#> R-squared               0.9647    0.9552    
#> Adjusted R-squared      0.9643    0.9547    
#>   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
#> p 8.608 5.845e-21 0.001517
#> 
#> 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
#> p    0.02635    w   0.6496 0.0006945  0.6494
#>  full matrix in summary(object)$diagnostics$exog.correlation.matrix
#> 
#> Collinearity of the copula terms (omega near 0 = weakly identified):
#>       corr(P, C)  omega
#> p_cop     0.9147 0.1603
#>