Copula endogeneity correction of Liengaard et al. (2025) (JAMS)

Generalised Gaussian copula estimator of Liengaard, Becker, Bennedsen, Heiler, Taylor and Ringle (2025). Endogenous regressors are transformed to normal scores and combined with the (also transformed) continuous exogenous regressors through the inverse of their empirical variance-covariance matrix (Equation 17) to form the copula control terms added to the structural model (Equation 19). The estimator's distinctive feature is that this copula structure may be estimated separately within the joint categories formed by the discrete (factor or character) exogenous regressors (Equations 20-21), controlled by conditional. Discrete exogenous regressors always stay in the structural model as ordinary controls; conditional only decides whether they also split the copula terms.

Usage

CopRegJAMS(
  formula,
  data,
  conditional = TRUE,
  cdf = "ecdf.adj",
  ties = "max",
  nboots = 199,
  subset = NULL,
  contrasts = NULL,
  parallel = FALSE,
  ncores = NULL,
  verbose = interactive()
)

Arguments

formula: A two-part Formula of the form y ~ endogenous_1 + endogenous_2 | exogenous_1 + exogenous_2. Terms before the bar are endogenous and receive a copula correction; terms after it are exogenous and receive none. Anything lm() accepts is allowed, including interactions, polynomials, factors and - 1 to drop the intercept. Copula terms are only generated for first-order endogenous terms (matching \(d_P\) in Equation 17 of the source paper); giving an interaction or a power of an endogenous regressor its own position before the bar produces a warning.
data: A data.frame containing the variables named in formula. Missing values are handled by na.omit.
conditional: Controls which exogenous regressors split the copula structure into categories (Equations 20-21). TRUE (the default) uses the joint categories of all factor and character regressors in the exogenous part of formula: two binary variables give four cells, not two. FALSE imposes one common copula structure across the whole sample (Equation 18). A character vector instead names the specific exogenous variables whose categories the structure may vary over, which is useful when the joint number of cells implied by TRUE would be too large for the sample. Under all three settings the discrete regressors remain in the structural model as ordinary controls and are never copula-transformed themselves: \(W\) in Equation 17 is always the continuous exogenous vector.
cdf: Character string naming the estimator of the marginal CDF used in the copula transformation. Defaults to "ecdf.adj", the adjusted empirical CDF proposed in the source paper.
ties: Character string controlling how ties are broken in the CDF estimate: "max" (the default) for the counting function, or "average" for midranks.
nboots: Number of bootstrap replicates used to compute standard errors.
subset: Optional logical vector selecting a subset of data, as in lm().
contrasts: Optional contrasts specification passed to model.matrix.
parallel: FALSE (the default), TRUE, "multicore" or "snow" to spread the bootstrap replicates across workers.
ncores: Number of workers to use when parallel is not FALSE. NULL (the default) uses one less than the number of detected cores.
verbose: Logical; show a progress bar while bootstrapping. Defaults to interactive().

Value

An object of class "copreg": a list including

coefficients, std.error, vcov: Point estimates, bootstrap standard errors and bootstrap covariance matrix for the structural coefficients and the copula terms.
rho, rho.se: Estimated correlation(s) between the structural residual and each endogenous regressor's copula data, with bootstrap standard errors.
boot, boot.rho, boot.ols: The raw bootstrap replicates underlying the above.
ols.coefficients, ols.std.error: The uncorrected OLS benchmark, estimated on the same bootstrap resamples.
icon: Standard error inflation of the corrected estimator relative to uncorrected OLS.
fitted.values, residuals: Fitted values and residuals of the structural model (copula terms excluded).
copula.terms, copula.data, copula.groups: The copula control terms actually added to the regression, the underlying normal scores of the endogenous regressors, and (for JAMS) which endogenous regressor and cell each copula column belongs to.
diagnostics: Identification and specification diagnostics used by validity().
call, method, cdf, ties, n: Book-keeping: the matched call, an estimator label, the CDF and tie method used, and the sample size.

Methods exist for print(), summary, coef(), vcov(), residuals(), fitted(), predict(), confint(), nobs(), formula(), update() and validity(). For JAMS specifically, summary() additionally reports a bootstrap Wald test of the joint significance of the copula terms and, whenever conditional produces more than one cell, a second bootstrap Wald test of whether the copula terms differ across those cells.

References

Liengaard, B. D., J.-M. Becker, M. Bennedsen, P. Heiler, L. N. Taylor, and C. M. Ringle (2025). Dealing with regression models' endogeneity by means of an adjusted estimator for the Gaussian copula approach. Journal of the Academy of Marketing Science 53, 279-299.

Examples

set.seed(1)
n <- 300
z <- factor(sample(c("a", "b"), n, replace = TRUE))
w <- rnorm(n)
u <- rnorm(n)
p <- 0.6 * u + rnorm(n)
y <- 1 + 2 * p + 0.5 * w + (z == "b") + u
dat <- data.frame(y = y, p = p, w = w, z = z)

fit <- CopRegJAMS(y ~ p | w + z, data = dat, nboots = 25)
summary(fit)
#> 
#> Copula endogeneity correction: JAMS (Liengaard et al. 2025) 
#> 
#> Call:
#> CopRegJAMS(formula = y ~ p | w + z, data = dat, nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>      Min       1Q   Median       3Q      Max 
#> -2.33260 -0.54857 -0.06364  0.71824  2.83336 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>       Min        1Q    Median        3Q       Max 
#> -2.964652 -0.743361 -0.002778  0.693468  3.099066 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.13409    0.08556  13.255  < 2e-16 ***
#> p            2.92685    0.37752   7.753 8.98e-15 ***
#> w            0.42374    0.07180   5.902 3.60e-09 ***
#> zb           0.70984    0.11570   6.135 8.50e-10 ***
#> p_cop.za    -0.60715    0.48989  -1.239    0.215    
#> p_cop.zb    -0.46754    0.45689  -1.023    0.306    
#> ---
#> 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.5180     0.3666  -1.413    0.158
#> 
#> Joint significance of the copula terms (bootstrap Wald test):
#>   chi-squared = 2.335 on 2 df, p = 0.3112
#> Copula structure constant across the 2 categories of the discrete
#>   regressors (bootstrap Wald test): chi-squared = 1.751 on 1 df,
#>   p = 0.1858. A small p favours the grouped estimation.
#> 
#> Fit, on 294 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.9062    1.0647    
#> R-squared               0.9286    0.9014    
#> Adjusted R-squared      0.9274    0.8997    
#>   sigma above is the standard error of the structural model, the one 
#>   entering xi* = xi / sigma.
#> Standard errors from 25 bootstrap replicates; cdf = "ecdf.adj", 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 0.1927 0.8943 0.9966
#> 
#> 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.07532    w   0.3869 0.005728  0.4261
#>  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.za     0.9975 0.01504
#> p_cop.zb     0.9883 0.01640
#> 

# \donttest{
fit_unconditional <- CopRegJAMS(y ~ p | w + z, data = dat,
                                conditional = FALSE, nboots = 25)
summary(fit_unconditional)
#> 
#> Copula endogeneity correction: JAMS (Liengaard et al. 2025), unconditional 
#> 
#> Call:
#> CopRegJAMS(formula = y ~ p | w + z, data = dat, conditional = FALSE, 
#>     nboots = 25)
#> 
#> Residuals of the augmented regression (u = xi - C gamma):
#>      Min       1Q   Median       3Q      Max 
#> -2.33667 -0.55235 -0.05705  0.68338  2.79300 
#> 
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#>      Min       1Q   Median       3Q      Max 
#> -5.07144 -1.26545 -0.09888  1.25522  4.95035 
#> 
#> Coefficients:
#>             Estimate Std. Error z value Pr(>|z|)    
#> (Intercept)  1.11151    0.11754   9.456  < 2e-16 ***
#> p            3.68561    0.85953   4.288 1.80e-05 ***
#> w            0.37062    0.09483   3.908 9.29e-05 ***
#> zb           0.83224    0.08889   9.363  < 2e-16 ***
#> p_cop       -1.47954    1.06890  -1.384    0.166    
#> ---
#> 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.8561     0.6395  -1.339    0.181
#> 
#> Fit, on 295 residual degrees of freedom:
#>                         augmented structural
#> Residual standard error 0.9050    1.7632    
#> R-squared               0.9285    0.7287    
#> Adjusted R-squared      0.9276    0.7250    
#>   sigma above is the standard error of the structural model, the one 
#>   entering xi* = xi / sigma.
#> Standard errors from 25 bootstrap replicates; cdf = "ecdf.adj", 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 0.1927 0.8943 0.9966
#> 
#> 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.09902   zb   0.1738  0.01545 0.09911
#>  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.997 0.001903
#> 
# }