copreg() only rewrites the call and hands it to the estimator's own
function, so both routes behave identically down to the warnings. This
matters because each wrapper carries rules of its own: BMW warns when the
CDF departs from the one its asymptotics are derived for, the two-stage
estimators redirect to Park and Gupta when the formula has no exogenous
part, and each tells summary() which assumptions it makes. CopRegPANEL()
and CopRegBAYES() are deliberately not reachable through copreg(): the
former needs a panel index and a different likelihood, the latter replaces
bootstrap inference by a posterior.
Usage
copreg(
formula,
data,
method = c("pg", "2scope", "ima", "bmw", "jams", "np"),
cdf = NULL,
...
)Arguments
- formula
A two-part formula
y ~ endogenous_1 + endogenous_2 + ... | exogenous_1 + exogenous_2 + ..., seeCopRegPG()for the full syntax.- data
A
data.framecontaining the variables referenced informula.- method
Estimator to use: one of
"pg"(CopRegPG()),"2scope"(CopReg2sCOPE()),"ima","bmw","jams"or"np".- cdf
CDF estimator for the copula transformation; see
CopRegPG()for the seven available options.NULL(default) is dropped from the call so that the selected estimator's own default applies.- ...
Further arguments passed on to the selected estimator's function, e.g.
ties,nboots,subset,contrasts,parallel,ncores,verbose.
Value
An object of class "copreg", identical to what the estimator
named by method would return had it been called directly; see
CopRegPG() and CopReg2sCOPE().
References
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 <- copreg(y ~ p | w, data = dat, method = "pg", nboots = 25)
summary(fit)
#>
#> Copula endogeneity correction: PG (Park & Gupta 2012)
#>
#> Call:
#> copreg(formula = y ~ p | w, data = dat, method = "pg", nboots = 25)
#>
#> Residuals of the augmented regression (u = xi - C gamma):
#> Min 1Q Median 3Q Max
#> -2.6892 -0.6745 -0.0929 0.6247 3.0787
#>
#> Residuals of the structural model (xi = y - mu - P alpha - W beta):
#> Min 1Q Median 3Q Max
#> -2.53982 -0.62128 -0.06731 0.61248 3.31175
#>
#> Coefficients:
#> Estimate Std. Error z value Pr(>|z|)
#> (Intercept) 1.35084 0.43959 3.073 0.00212 **
#> p 1.86613 0.16466 11.333 < 2e-16 ***
#> w 0.53260 0.05763 9.242 < 2e-16 ***
#> p_cop 0.49819 0.33814 1.473 0.14066
#> ---
#> 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.4512 0.1810 2.493 0.0127 *
#> ---
#> 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.9668 1.0834
#> R-squared 0.9645 0.9555
#> Adjusted R-squared 0.9642 0.9550
#> sigma above is the standard error of the structural model, the one
#> entering xi* = xi / sigma.
#> Standard errors from 25 bootstrap replicates; cdf = "kde.silverman", 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
#> (Park & Gupta assume this is zero; 'joint' tests all of them at once,
#> the Holm p value refers to the single largest correlation):
#> max |corr| with p (Holm) joint R2 joint p
#> p 0.02702 w 0.6414 0.0007299 0.6412
#> 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.9214 0.15
#>