Package {BiFactorCopulaGAMM}

(internal) E[Y_j] numerator and B(x,v), for ONE observation, ONE Y-margin, fixed v

Description

Returns 'list(Bxv, Enum)' such that 'E[Y_j|X=x,V=v] = Enum/Bxv' (eq. 2.5 for continuous Y_j, eq. 2.6 for discrete Y_j). For the integrated-over-V functions, callers combine several v-nodes' 'Bxv'/'Enum' via a further v-quadrature (eqs. 2.10/2.11) before dividing, rather than dividing here and re-weighting afterward.

Usage

.EYj_terms_given_Xv(
  v,
  x_F1,
  x_F2,
  cont_X,
  fam1_X,
  rot1_X,
  cpar1_X,
  dfC1_X,
  fam2_X,
  rot2_X,
  cpar2_X,
  dfC2_X,
  margin_j,
  param_j,
  dfM_j,
  fam1_j,
  rot1_j,
  cpar1_j,
  dfC1_j,
  fam2_j,
  rot2_j,
  cpar2_j,
  dfC2_j,
  is_disc_j,
  nl,
  wl
)

(internal) X-side w-integrand, one observation, at a fixed v

Description

Computes, for ONE observation's X covariates (already transformed to F_l(x_l)/F_l(x_l-)) and a fixed v, the vector (over 'wgrid') of 'product over l of [c_l,1[F_l(x_l),v] * c_l,2[C_l,1[F_l(x_l),v],w]]' (continuous X_l) or '[C_l,2[...,w] - C_l,2[...,w]]' (discrete X_l) – i.e. the integrand of B(x,v) (eqs. 2.3/2.8) before the w-quadrature sum.

Usage

.Xw_terms(
  x_F1,
  x_F2,
  cont_X,
  fam1_X,
  rot1_X,
  cpar1_X,
  dfC1_X,
  fam2_X,
  rot2_X,
  cpar2_X,
  dfC2_X,
  v,
  wgrid
)

Details

condMeanY_bifactor.R

Computes the conditional MEAN of each Y-component Y_j (j=1..d), given the observed X covariates of each observation, following "Bi-factor copulas for clustered data":

- condMeanY_given_XV(): E[Y_j | X=x, V=v], plugging in a PREDICTED value of the cluster-level latent factor V (e.g. from condMeansBifactor()/condMediansBifactor()) – continuous Y_j uses eq. (2.5), discrete Y_j uses eq. (2.6). - condMeanY_given_X(): E[Y_j | X=x], integrating V out entirely (no latent-variable estimate needed at all) – continuous Y_j uses eq. (2.10), discrete Y_j uses eq. (2.11).

Both follow the 'newdata = NULL' convention already used by 'predictContinuous.R': 'newdata = NULL' re-predicts the original training observations; otherwise 'newdata' is a list with an 'x' element (and, for 'condMeanY_given_XV()', a 'clu' element saying which fitted cluster's V to use for each new row).

Source 'condMedianV_bifactor.R' first – this file reuses ‘.prepare_bifactor_data()', '.pmargin_natural()', and the package’s own ‘pcond()'/'dcop()'/'qmargin()' (make sure you’re using the FIXED 'qmargin.R', which this file also depends on for continuous Y margins other than normal/t/laplace – the original had a bug in five of its branches; see that file's own NOTE).

MATH RECAP (see the paper for the full derivation): both eq. (2.5)/ (2.6) and (2.10)/(2.11) share the same "X-side" building block – B(x,v) = integral_0^1 [product over X-margins l of the copula-1 x copula-2 contribution] dw – which is EXACTLY the same quantity our 'cond_median_logh()'/'.lden_grid_bifactor()' machinery already computes as "h_obs" for one observation, just restricted to the X-margins only (no Y-term). '.Xw_terms()'/'.EYj_terms_given_Xv()' below recompute this directly in R (rather than reusing that C-kernel path), since here we additionally need the RAW per-w-node values (not already collapsed by quadrature) to combine with each Y-quantile/Y-mass term – this is a one-off per-row prediction, not a hot loop inside estimation, so plain (not C-accelerated) R is used throughout.

(internal) Bound-aware central-difference gradient

Description

Computes a numerical gradient of 'f' at 'x' using central differences, with a per-parameter step size that is automatically shrunk so the probe points 'x[i] +/- h[i]' never cross 'lower[i]'/ 'upper[i]'. This is what actually fixes (rather than papers over) the "exploding gradient" failure mode: a plain finite-difference gradient has no notion of a feasible region and will happily probe past a tight bound, hitting whatever infeasibility penalty the objective returns there and producing a wildly inflated derivative estimate for that coordinate.

Usage

.central_grad(f, x, lower, upper, h0 = 1e-04)

(internal) Numeric derivative of linkCop()'s cpar map

Description

Returns d(cpar)/d(param), via central difference of 'linkCop(param, family)$cpar' – used by the delta method to convert a raw-scale (pre-link) copula-parameter standard error into the natural-units scale 'coef$copulas' is reported in.

Local copy of the helper from EstBifactorCopulaGAMMFast.R – see '.numerical_hessian_from_grad()' above for why this file now carries its own copy rather than relying on that file's. Returns d(cpar)/d(param), via central difference of 'linkCop(param, family)$cpar'.

Usage

.cop_link_deriv(param_raw, family, h = 1e-06)

.cop_link_deriv(param_raw, family, h = 1e-06)

(internal) Copula family name -> integer code, for the C kernel

Description

Local copy of the same mapping 'family_code()' uses in 'EstMargin1Factor.R'/'EstBifactorCopulaGAMMFast.R' – must match ‘bifactor_kernel.c'’s family codes exactly.

Usage

.family_code(fam)

Conditional (posterior) median of the cluster-level factor V, for ONE cluster, all margins jointly

Description

Implements (2.9)-(2.11) for a single cluster: for each candidate v, computes every member's own per-observation density h(y_i,x_i,v) by integrating out that member's own nuisance w (the inner loop, over ALL dp margins jointly), multiplies across the n_k members (the outer product in 2.10), then finds the v at which the resulting posterior CDF equals 0.5 – using the same quadrature-rescaling + 'invfunc()'-bisection trick as 'MAP.continuous()'/'MAP.discrete()'.

Shared by 'condMedianV_bifactor()' and 'condMeanV_bifactor()': tries the compiled 'cond_median_logh()' C kernel first, falling back to the vectorized pure-R computation (tiling every (member, v-node) pair, summing logs across margins, log-sum-exp across the w-quadrature nodes) if the kernel isn't available or a family isn't covered by it.

Usage

.lden_grid_bifactor(
  vgrid,
  F1,
  F2,
  n_k,
  dp,
  cont,
  fam1,
  rot1,
  cpar1,
  dfC1,
  fam2,
  rot2,
  cpar2,
  dfC2,
  nl,
  wl,
  use_c,
  kernel_path
)

(internal) Load or compile the bifactor C kernel

Description

Loads 'bifactor_kernel.so'/'.dll' if already built and up to date; otherwise compiles 'bifactor_kernel.c' with 'R CMD SHLIB' and loads the result. Safe to call repeatedly (no-op if already loaded).

The clusters processed by 'bifactor_negloglik()' are statistically independent, so the kernel parallelizes the cluster loop with OpenMP when available. This function first tries to compile with 'PKG_CFLAGS=-fopenmp PKG_LIBS=-fopenmp'; if that fails (e.g. plain Xcode clang on macOS with no 'libomp' installed, or a compiler that otherwise doesn't support '-fopenmp'), it automatically falls back to a plain (single-threaded, but still fully correct) build rather than erroring out. Set 'options(bifactor.use_openmp = FALSE)' before calling this to skip the OpenMP attempt entirely and go straight to the plain build.

Usage

.load_bifactor_kernel(kernel_path)

(internal) Load or compile the cond_median_logh() C kernel

Description

Same pattern as EstBifactorCopulaGAMMFast()'s '.load_bifactor_kernel()': loads 'bifactor_kernel.so'/'.dll' if already built and up to date, otherwise compiles 'bifactor_kernel.c' with 'R CMD SHLIB' (trying '-fopenmp' first, falling back to a single-threaded build). This is the SAME source file/compiled object 'EstBifactorCopulaGAMMFast()' already uses for 'bifactor_negloglik()' – 'cond_median_logh()' was added alongside it, reusing the same 'pcond_c()'/'dcop_c()' static helpers, so loading one gives you both. If you've already called '.load_bifactor_kernel()' (e.g. via 'EstBifactorCopulaGAMMFast()'/'EstBifactorTwoStep()') in this session, this is a near-instant no-op ('is.loaded()' short-circuits it).

Usage

.load_cond_median_kernel(kernel_path = "bifactor_kernel.c")

Details

condMedianV_bifactor.R

Computes the conditional (posterior) MEDIAN of the cluster-level latent factor V_k for each cluster, combining all d Y-margins and p X-margins jointly, following eqs. (2.9)-(2.11) of "Bi-factor copulas for clustered data": for cluster k with members i=1..n_k,

h_k(y_k,x_k,v) = prod_[i=1..n_k] h(y_ki, x_ki, v) (built from 2.9 via 2.5) h_k(y_k,x_k) = integral_0^1 h_k(y_k,x_k,v) dv (2.10) posterior of V_k at v = h_k(y_k,x_k,v) / h_k(y_k,x_k) (2.11)

IMPORTANT (this is the correction from the previous version of this file): 'w' in (2.9) is a PER-OBSERVATION nuisance variable, not a second cluster-level latent factor – see section 2.2 of the paper, where the auxiliary variables are written W_[k,i] (both a cluster index k AND a member index i). It is integrated out separately for every single cluster member i (to combine that member's own Y/X into one joint density given V=v), then the n_k per-member results are multiplied together, and ONLY THEN is the integral over v=V_k taken. There is therefore only ONE latent variable with a meaningful cluster-level posterior to report: V_k. This is exactly what ‘bifactor_kernel.c'’s bifactor_negloglik() computes too (its q1 loop is the outer integral over v; for each v, and EACH cluster member row separately, it re-runs its own inner q2 loop over w and multiplies the resulting per-member h_val's together via 'sum_log_h') – this file mirrors that same nested structure, but inverts the resulting posterior CDF of V (via the same quadrature + invfunc-bisection pattern MAP.continuous()/MAP.discrete() use) instead of returning a log-likelihood.

ASSUMPTIONS (please double-check against your fitted object – I could not execute R to verify these directly): - 'coef$copulas' is a dp x 2 matrix of NATURAL copula parameters (cpar1, cpar2), i.e. already passed through linkCop(), and cpar1/ cpar2 are constant per margin (no copula-level covariates), matching how EstBifactorTwoStep() itself fits them. - 'coef$margins' is a dp x 2 matrix in NATURAL units (via link2parM()), NA-padded in column 2 for 1-parameter margins. - 'model' has the same shape used throughout (Sim2Cluster/ EstBifactorTwoStep): model$Y[[j]]$margin$family/dfM, model$Y[[j]]$copula[[1|2]]$family/rot/dfC, and likewise model$X[[j]].

(internal) Load or compile the Step-1 margin/copula-1 C kernel

Description

Same pattern as ‘EstBifactorCopulaGAMMFast()'’s '.load_bifactor_kernel()': loads 'margin1factor_kernel.so'/'.dll' if already built and up to date, otherwise compiles 'margin1factor_kernel.c' with 'R CMD SHLIB', trying '-fopenmp' first (clusters are independent, same as the joint model's kernel) and falling back to a single-threaded build if that fails.

Usage

.load_margin1factor_kernel(kernel_path)

(internal) Numeric derivative of link2parM()'s natural-units map

Description

Returns d(natural units)/d(raw theta) elementwise, via central difference of 'link2parM()' itself – used by the delta method to convert a raw-scale standard error into a natural-units one. Only the first 'nth' entries are meaningful.

Returns d(natural units)/d(raw theta) elementwise, via central difference of 'link2parM()' itself – used by the delta method to convert a raw-scale standard error into a natural-units one. Only the first 'nth' entries are meaningful (matching how many real parameters this margin actually has).

Usage

.margin_link_deriv(theta_raw, marg, nth, h = 1e-06)

.margin_link_deriv(theta_raw, marg, nth, h = 1e-06)

(internal) Numerical Hessian from a gradient function

Description

Forms the Hessian of a scalar function at 'x' via central differences of its gradient (rather than second differences of the function value directly, which is noisier) – 'H[,i]' is '(grad_fn(x + h*e_i) - grad_fn(x - h*e_i)) / (2h)', symmetrized afterward to cancel small numerical asymmetries. Used to get standard errors (via the inverse Hessian) for 'se'/'tval'/'pval'.

Local copy of the helper from EstBifactorCopulaGAMMFast.R (this file previously relied on that one being sourced first, which broke – "could not find function" – whenever an older copy of that file, from before these helpers were added, was still loaded). Same implementation as there: forms the Hessian of a scalar function at 'x' via central differences of its gradient, symmetrized afterward.

Forms the Hessian of a scalar function at 'x' via central differences of its gradient (rather than second differences of the function value directly, which is noisier) – 'H[,i]' is '(grad_fn(x + h*e_i) - grad_fn(x - h*e_i)) / (2h)', symmetrized afterward to cancel small numerical asymmetries.

Usage

.numerical_hessian_from_grad(grad_fn, x, lower, upper, h0 = 1e-04)

.numerical_hessian_from_grad(grad_fn, x, lower, upper, h0 = 1e-04)

.numerical_hessian_from_grad(grad_fn, x, lower, upper, h0 = 1e-04)

(internal) CDF (and, for discrete families, F(y-)) from natural-unit margin parameters

Description

Mirrors qmargin.R's family dispatch/parameterization, but returns the CDF instead of the quantile – needed because the fitted margin parameters are already in natural units (post-link2parM()), whereas margins0() expects the raw pre-link theta.

Usage

.pmargin_natural(y, family, param, dfM = 0)

(internal) Batched, all-clusters-at-once setup for condMediansBifactor()/condMeansBifactor()

Description

Tries to load the 'cond_median_logh_multi()' C kernel (batched across ALL K clusters, with per-cluster v-grids, OpenMP- parallelized across clusters – the multi-cluster companion to 'cond_median_logh()') and, if available, runs the shared "evaluate at the base grid 'nl' for every cluster in ONE call" step both the mean and the median need. This is what lets 'condMediansBifactor()' call 'invfunc()' a single time across all K clusters at once (using ‘invfunc()'’s own built-in vectorization – it subsets 'q'/'v0'/'v1'/ 'th' together by the same convergence mask, so passing cluster indices as 'th' lets our 'func' recover which clusters are still active at each bisection step) instead of K separate per-cluster 'invfunc()' calls. Falls back to 'list(use_multi = FALSE, ...)' if the multi-cluster kernel isn't available (older compiled kernel, no C compiler, an uncovered family, etc.) – callers should fall back to the original per-cluster loop over 'condMedianV_bifactor()'/'condMeanV_bifactor()' in that case.

Usage

.prepare_bifactor_batch(prep, nl, wl, nq, use_c, kernel_path)

(internal) Shared bookkeeping for condMediansBifactor()/condMeansBifactor()

Description

Reconstructs, from 'Y'/'X'/'model'/'coef', everything both wrappers need: per-margin family/rotation/df arrays (Y then X, matching ‘EstBifactorTwoStep()'’s own convention), the natural-unit copula parameters, and each margin's CDF ('F1'/'F2') evaluated at every observation. Factored out so the median and mean wrappers don't duplicate this setup.

Usage

.prepare_bifactor_data(Y, X, clu, model, coef)

Bi-factor copula estimation for clustered mixed data (fast version)

Description

This is a faster replacement for 'EstBifactorCopulaGAMM()'. It is NOT numerically identical to the original: the continuous-margin branch of the likelihood has been corrected to match eq. (2.8) of the bi-factor copula paper (see the note below), so expect the estimates to differ from 'EstBifactorCopulaGAMM()' whenever the model includes a continuous margin. Everything else about the objective function is unchanged; the only remaining differences are *how* it is evaluated.

Why the original is slow: 'nlm()' evaluates the objective function many times (it estimates the gradient numerically, so every iteration costs roughly 'npar+1' evaluations of 'likf()', and 'iterlim=200' allows many iterations). Each call to the original 'likf()' loops over every cluster in R and, for every cluster, calls 'dcop()'/'pcond()' – which each re-dispatch through a 'switch()' on the family name – on freshly allocated vectors. That per-call R-interpreter overhead, multiplied by (number of clusters) x (number of 'nlm' evaluations), is what dominates the runtime; the actual arithmetic is tiny.

This version:

1. Precomputes, once (outside the objective function), everything that does not change across 'nlm' evaluations: the cluster grouping (row order + boundaries), and the integer family/ rotation/df codes for every copula slot.

2. Moves the per-cluster double-Gauss-quadrature loop (previously re-run, in R, on every objective evaluation) into a single call to a compiled C routine, 'bifactor_negloglik()', defined in 'bifactor_kernel.c'. That C routine reproduces the arithmetic of 'dcop()'/'pcond()' exactly (including the 'pcond.R' "fgm" rotation case), with one intentional correctness fix: the continuous-margin branch of the original 'likf()' computed 'temp1 = matff[,j]*dcop(u1,vv,...)' – the density prefactor 'f_j(x_j)*c_[j,1][F_j(x_j),v]' required by eq. (2.8) of the bi-factor copula paper – but then discarded it and multiplied by the conditional cdf 'temp2' in its place. Both this C kernel and the pure-R fallback below now correctly use 'temp1'/'f_j*c_[j,1]' as the multiplier instead.

3. If the C library cannot be found/compiled/loaded, this function automatically falls back to a pure-R evaluator (still with the precomputed cluster grouping, which is a safe, free win on its own) and prints a one-time message explaining how to compile the C kernel for the full speedup.

4. Parallelizes the C kernel's cluster loop with OpenMP (clusters are statistically independent, so this is embarrassingly parallel and exact – not an approximation). This is the single biggest remaining lever if estimation is still slow: with 'method="nlminb"', the numerical gradient costs '2*length(start)' evaluations of the C kernel per outer iteration, so on a machine with several cores this can give a near-linear additional speedup on top of the C rewrite itself. Compiled automatically by '.load_bifactor_kernel()'; falls back to a single-threaded build if the compiler doesn't support '-fopenmp'.

Usage

EstBifactorCopulaGAMMFast(
  Y,
  X = NULL,
  clu,
  model,
  start = NULL,
  LB = NULL,
  UB = NULL,
  nq = 25,
  prediction = TRUE,
  kernel_path = "bifactor_kernel.c",
  method = c("nlminb", "nlm"),
  compute_se = TRUE
)

Arguments

Value

list with 'coef' (list with 'copulas'/'margins' matrices, same shape as ‘EstBifactorTwoStep()'’s 'coef' – 'margins' is in natural units via ‘link2parM()', matching that function’s own convention, not the raw pre-link scale Step 1 optimizes over) and 'gradient' (shaped like 'coef': a list with 'copula' and 'margin' matrices, each 'dp x 2' with the same row order as 'coef$copulas'/ 'coef$margins', rather than a flat vector. Same parameter scale as ‘EstBifactorTwoStep()'’s ‘gradient' – Step 1’s copula-1 gradient is chain-rule-converted from the raw copula-parameter scale it's actually computed on), plus 'se'/'tval'/'pval' (same 'list(copula=,margin=)' shape as 'coef', standard errors and asymptotic Wald tests – see 'se_note' for whether these are the naive two-step version or the full-MLE version), 'se_note' (a short string explaining which), 'loglik', 'aic', 'bic'. These all come from 'final' (the refined joint fit) when 'refine=TRUE'; otherwise they're assembled directly from the Step 1 + Step 2 results – in that case 'se'/'tval'/'pval' are the NAIVE two-step version (Step 1's margin/cpar1 standard errors reused from each ‘EstMargin1Factor()' fit in isolation; Step 2’s cpar2 standard error treats margin/cpar1 as fixed/known), which understate the true uncertainty – use 'refine=TRUE' for valid inference. Also returns 'gradient_copulas', a 'dp x 2' matrix (same row order as 'coef$copulas') with columns '"step1"' (the chain-rule-converted cpar1 gradient from 'EstMargin1Factor()') and '"step2"' (the cpar2 gradient from the Step 2 fit) – always the two-step breakdown, regardless of 'refine', since that distinction doesn't exist once/if the joint refinement runs. Also returns 'step1' (per-margin 'EstMargin1Factor()' results), 'step2' (the copula-2-only fit), 'warm_start' (the assembled full parameter vector – also used to warm-start the refinement step when 'refine=TRUE'), 'final' (the refined joint fit if 'refine=TRUE', else 'NULL'), and 'calc_time' (a named vector with the wall-clock time for each phase).

Two-step ("IFM"-style) estimation for the bi-factor copula GAMM

Description

Implements the two-step estimator suggested by section 2.1 of the bi-factor copula paper: "the parameters of Gj and Dj,1 can be estimated using the CRAN package CopulaGAMM ... with no covariates," applied separately to every margin (Step 1, via 'EstMargin1Factor()', with an analytic gradient built from 'margins()'/'coplik()'), followed by plugging those fixed values into the full bi-factor likelihood to estimate only the copula-2 parameters (Step 2, reusing the existing C kernel, 'bifactor_negloglik()', but now optimizing over just 'dp' parameters instead of the whole model).

This changes the *shape* of the optimization, not just its constant: Step 1 replaces one big, badly-scaled joint optimization with 'dp' small, cheap, single-quadrature ones (no integral over the second latent factor W at all); Step 2 keeps the double quadrature but only has to differentiate 'dp' parameters instead of the full parameter count. Both are typically much cheaper than one joint fit over everything at once.

IMPORTANT CAVEAT: two-step (IFM-style) estimators are consistent but not fully efficient – standard errors read off Step 2 alone understate the true uncertainty, because they treat the Step 1 estimates as known constants rather than propagating their sampling variability. If you need valid inference (not just point estimates), either use a sandwich/Godambe covariance, bootstrap over clusters, or (recommended, and the default here) use these estimates purely as a warm start for one final joint optimization over all parameters ('refine=TRUE'): since it starts essentially at the optimum, that final pass typically converges in a handful of iterations rather than a full cold-start search, so you get IFM-level speed for the hard part and exact joint- MLE point estimates/standard errors at the end.

Usage

EstBifactorTwoStep(
  Y,
  X = NULL,
  clu,
  model,
  start = NULL,
  LB = NULL,
  UB = NULL,
  nq = 25,
  prediction = TRUE,
  kernel_path = "bifactor_kernel.c",
  refine = TRUE,
  compute_se = TRUE
)

Arguments

Value

list with 'coef' (list with 'copulas'/'margins' matrices, same shape as ‘EstBifactorCopulaGAMMFast()'’s 'coef' – 'margins' is in natural units via ‘link2parM()', matching that function’s own convention, not the raw pre-link scale Step 1 optimizes over) and 'gradient' (shaped like 'coef': a list with 'copula' and 'margin' matrices, each 'dp x 2' with the same row order as 'coef$copulas'/ 'coef$margins', rather than a flat vector. Same parameter scale as ‘EstBifactorCopulaGAMMFast()'’s ‘gradient' – Step 1’s copula-1 gradient is chain-rule-converted from the raw copula-parameter scale it's actually computed on), plus 'se'/'tval'/'pval' (same 'list(copula=,margin=)' shape as 'coef', standard errors and asymptotic Wald tests – see 'se_note' for whether these are the naive two-step version or the full-MLE version), 'se_note' (a short string explaining which), 'loglik', 'aic', 'bic'. These all come from 'final' (the refined joint fit) when 'refine=TRUE'; otherwise they're assembled directly from the Step 1 + Step 2 results – in that case 'se'/'tval'/'pval' are the NAIVE two-step version (Step 1's margin/cpar1 standard errors reused from each ‘EstMargin1Factor()' fit in isolation; Step 2’s cpar2 standard error treats margin/cpar1 as fixed/known), which understate the true uncertainty – use 'refine=TRUE' for valid inference. Also returns 'gradient_copulas', a 'dp x 2' matrix (same row order as 'coef$copulas') with columns '"step1"' (the chain-rule-converted cpar1 gradient from 'EstMargin1Factor()') and '"step2"' (the cpar2 gradient from the Step 2 fit) – always the two-step breakdown, regardless of 'refine', since that distinction doesn't exist once/if the joint refinement runs. Also returns 'step1' (per-margin 'EstMargin1Factor()' results), 'step2' (the copula-2-only fit), 'warm_start' (the assembled full parameter vector – also used to warm-start the refinement step when 'refine=TRUE'), 'final' (the refined joint fit if 'refine=TRUE', else 'NULL'), and 'calc_time' (a named vector with the wall-clock time for each phase).

Step 1 of the two-step bi-factor estimator: single-margin, 1-factor fit

Description

Estimates, for one response or covariate margin j in isolation, its marginal-distribution parameters and the parameter of the copula linking it to the cluster-level latent factor V, using only the univariate 1-factor clustered-data likelihood of eq. (2.6)/(2.7) of the bi-factor copula paper – i.e. every other margin and the second-level latent factor W are ignored entirely. This is exactly Step 1 of the two-step estimator described in that section: "the parameters of Gj and Dj,1 can be estimated using the CRAN package CopulaGAMM ... with no covariates," applied separately to each margin.

Compared to plugging this margin into the full bi-factor joint likelihood, this is cheap in two independent ways: (1) it needs only a *single* Gauss-Legendre quadrature over V (no integral over W at all), and (2) it only has ‘nth+1' parameters (the margin’s own 1-2 parameters plus one copula parameter) instead of the full model's parameter count, so even a finite-difference gradient over just this piece is cheap.

Where possible, the gradient used by 'nlm()' is *analytic*, built from the closed-form derivatives already implemented in 'margins()' (margin cdf/pdf and their parameter derivatives) and 'coplik()' (copula conditional-cdf/density and their derivatives), combined via the chain rule – the same style 'mlecop()' in 'coplik.R' already uses for the simple i.i.d. case, extended here to the clustered/latent-factor likelihood (a quadrature-weighted mixture over V).

As with 'EstBifactorCopulaGAMMFast()', the per-cluster/quadrature reduction (previously an R-level loop) is offloaded to a compiled C kernel, 'margin1factor_negloglik_grad()' in 'margin1factor_kernel.c', whenever 'cop1' is one of the families 'coplik()' covers. That kernel also parallelizes across clusters with OpenMP (clusters are independent here too), the same way the joint model's kernel does. 'margins()' itself still runs in R (a single vectorized call per evaluation, not a per-cluster loop, so porting every margin family to C would cost more than it saves); only the reduction and the copula-side derivatives are in C. Falls back automatically to the pure-R evaluator below if the kernel can't be found/compiled/loaded.

Coverage note: 'coplik()' provides analytic derivatives for copula families "gaussian"/"normal", "t", "clayton", "frank", "gumbel", "joe", "plackett", "fgm". It does not cover "galambos", "huesler-reiss", or "ind". For a continuous margin with one of the covered families, the gradient below is fully analytic (both the margin-parameter and the copula-parameter pieces). For a discrete margin, 'coplik()' does not expose an analytic d/du of the conditional cdf C(u|v) (needed for the margin-parameter piece via the chain rule through u=G(y;theta)), so that piece falls back to a cheap central difference – still only 1-2 extra copula evaluations, since this is a 2-3 parameter problem, not a finite-difference gradient over the whole model. For a copula family not covered by 'coplik()' at all, the whole gradient falls back to ‘nlm()'’s own finite-difference gradient (via 'dcop()'/'pcond()'), which is still cheap here for the same reason.

Usage

EstMargin1Factor(
  z,
  clu,
  margin,
  cop1,
  rot1 = 0,
  dfM = 0,
  dfC1 = 0,
  nq = 25,
  start,
  LB,
  UB,
  kernel_path = "margin1factor_kernel.c",
  compute_se = TRUE
)

Arguments

Value

list with 'param' ('c(margin params, cpar1)' at the optimum, raw/pre-link2parM scale), 'coef' ('margin' converted to natural units via 'link2parM()', and 'cpar1' – already natural units, no conversion needed since this function optimizes cpar1 directly, not the joint model's linked "param" scale), 'se'/'tval'/'pval' (same shape as 'coef', standard errors from a numerical Hessian of the negative log-likelihood with the margin part delta-method-converted to match coef's natural units; 'tval' is an asymptotic Wald statistic, not a finite-sample Student-t, so 'pval' uses the normal approximation), 'loglik', 'gradient' (shaped like 'coef': a list with 'margin' and 'cpar1', on the RAW optimization scale – a convergence diagnostic, not meant to pair with ‘coef'’s units the way 'se'/'tval'/'pval' do), 'mle' (the raw 'nlm()' object), and 'used_c_kernel' (whether the compiled kernel was actually used). NOTE: the 'se'/'tval'/'pval' fallback path (for objectives without an analytic gradient attribute) calls '.central_grad()', defined in 'EstBifactorCopulaGAMMFast.R' – source that file too if calling 'EstMargin1Factor()' on its own.

Simulation of clustered data using bi-factor copulas

Description

Generate a random sample of observations from a bi-factor copula mixed regression model.

Usage

Sim2Cluster(model, nclu)

Arguments

Value

Author(s)

Bruno N. Remillard

Examples

nclu=rep(5,25) # 25 clusters of size 5
ccY=list(list(family="gumbel",param=2),list(family="gumbel",param=1.4,rot=90) )
modY1 = list(margin=list(family="bernoulli", param=c(0.6)), copula=ccY)
modY=list(modY1)
ccX1 = list(list(family="clayton",param=0.86,rot=180),list(family="clayton",param=0.2,rot=180))
ccX2 = list(list(family="gaussian",param=0.6),list(family="clayton",param=2))
modX1 = list(margin=list(family="bernoulli", param=c(0.6)),  copula=ccX1)
modX2 = list(margin=list(family="gaussian" , param=c(10,0.2)), copula=ccX2)
modX=list(modX1,modX2)
model=list(Y=modY,X=modX)
set.seed(123)
out=Sim2Cluster(model,nclu)

Starting values, upper and lower bounds for parameter estimation of bi-factor copulas

Description

The starting point is a value giving approximately a Kendall's tau of 0.5, and the upper and lower bounds are values giving a tau of 0.99 (resp. 0.1 or -0.99) for the copula families

Usage

bounds.cop(cop)

Arguments

Value

out List of starting points (start), upper bounds (UB) and lower bounds (LB)

Examples

cop=c("gumbel","gumbel","normal","clayton","frank","clayton")
out=bounds.cop

Starting values, upper and lower bounds for parameter estimation of margins

Description

The starting point is an estimation of the parameters and the upper and lower bounds are approximate 99.7

Usage

bounds.mar(X, margin)

Arguments

Value

out List of starting points (start), upper bounds (UB) and lower bounds (LB)

Examples

 margin=c("bernoulli","poisson","normal")
 U=matrix(runif(300),ncol=3)
 X = cbind(as.double(U[,1]>.3),qpois(U[,2],5),10+2*qnorm(U[,3]))
 out=bounds.mar(X,margin)

Conditional (posterior) MEAN of the cluster-level factor V, for ONE cluster, all margins jointly

Description

Companion to 'condMedianV_bifactor()', computing the conditional posterior MEAN of V_k instead of its median: 'E[V_k | data] = integral v * h_k(v) dv / integral h_k(v) dv'. This is substantially cheaper than the median, since it only needs 'h_k(v)' evaluated ONCE, at the base quadrature grid – unlike the median, which needs ‘invfunc()'’s ~40 bisection evaluations, each at a different rescaled grid. Reuses the exact same '.lden_grid_bifactor()' machinery (and, when available, the compiled 'cond_median_logh()' C kernel) as 'condMedianV_bifactor()' – source that part of the file first.

Usage

condMeanV_bifactor(
  F1,
  F2 = NULL,
  cont,
  fam1,
  rot1,
  cpar1,
  dfC1 = NULL,
  fam2,
  rot2,
  cpar2,
  dfC2 = NULL,
  nq = 31,
  gl = NULL,
  use_c = TRUE,
  kernel_path = "bifactor_kernel.c"
)

Arguments

Value

scalar: the conditional posterior mean of V_k for this cluster

Conditional mean of every Y-component given X, with or without a predicted V

Description

Single entry point combining 'condMeanY_given_XV()' and 'condMeanY_given_X()': if 'V' is supplied, computes E[Y_j|X=x,V=v] (eqs. 2.5/2.6), plugging in each row's cluster's predicted latent factor; if 'V=NULL' (default), computes E[Y_j|X=x] instead (eqs. 2.10/2.11), integrating the latent factor out entirely – in that case ‘clu' isn’t needed either (silently ignored if supplied).

Usage

condMeanY(
  Y = NULL,
  X = NULL,
  clu = NULL,
  model,
  coef,
  V = NULL,
  newdata = NULL,
  nq = NULL
)

Arguments

Value

data.frame, one row per observation, columns 'EY1',...,'EYd'

Conditional mean of every Y-component given X only (V integrated out), for each observation

Description

E[Y_j | X=x] for j=1..d, for every observation, using each observation's OWN X covariates and integrating the cluster-level latent factor V out entirely – no predicted V needed at all (and hence no 'clu'/'Vhat' either). Eq. (2.10) for continuous Y_j, eq. (2.11) for discrete Y_j.

Usage

condMeanY_given_X(X = NULL, model, coef, newdata = NULL, nq = 25)

Arguments

Value

data.frame, one row per observation, columns 'EY1',...,'EYd'

Conditional mean of every Y-component given X and a predicted V, for each observation

Description

E[Y_j | X=x, V=v] for j=1..d, for every observation, using each observation's OWN X covariates and its CLUSTER's predicted latent factor V (from 'condMeansBifactor()'/'condMediansBifactor()') – eq. (2.5) for continuous Y_j, eq. (2.6) for discrete Y_j.

Usage

condMeanY_given_XV(
  Y = NULL,
  X = NULL,
  clu = NULL,
  model,
  coef,
  Vhat,
  newdata = NULL,
  nq = 31
)

Arguments

Value

data.frame, one row per observation, columns 'EY1',...,'EYd'

Conditional mean of V for every cluster, all Y/X components jointly

Description

Same interface/bookkeeping as 'condMediansBifactor()'. When the batched 'cond_median_logh_multi()' kernel is available, computes every cluster's mean from the SAME single, all-K-clusters-at-once evaluation '.prepare_bifactor_batch()' performs (the mean needs no bisection, so this is just one weighted average per cluster after that one shared call) – otherwise falls back to the original per-cluster loop over 'condMeanV_bifactor()'.

Usage

condMeansBifactor(
  Y,
  X = NULL,
  clu,
  model,
  coef,
  nq = 31,
  use_c = TRUE,
  kernel_path = "bifactor_kernel.c"
)

Arguments

Value

data.frame with one row per cluster: 'cluster', 'V' (conditional posterior mean)

Conditional (posterior) median of the cluster-level factor V, for ONE cluster, all margins jointly

Description

Implements (2.9)-(2.11) for a single cluster: for each candidate v, computes every member's own per-observation density h(y_i,x_i,v) by integrating out that member's own nuisance w (the inner loop, over ALL dp margins jointly), multiplies across the n_k members (the outer product in 2.10), then finds the v at which the resulting posterior CDF equals 0.5 – using the same quadrature-rescaling + 'invfunc()'-bisection trick as 'MAP.continuous()'/'MAP.discrete()'.

Usage

condMedianV_bifactor(
  F1,
  F2 = NULL,
  cont,
  fam1,
  rot1,
  cpar1,
  dfC1 = NULL,
  fam2,
  rot2,
  cpar2,
  dfC2 = NULL,
  nq = 31,
  gl = NULL,
  use_c = TRUE,
  kernel_path = "bifactor_kernel.c"
)

Arguments

Value

scalar: the conditional posterior median of V_k for this cluster

Conditional median of V for every cluster, all Y/X components jointly

Description

Top-level wrapper: given the original data and a fitted bi-factor model (as returned by 'EstBifactorTwoStep()'/ ‘EstBifactorCopulaGAMMFast()'), reconstructs each margin’s CDF from its fitted (natural-unit) parameters, then estimates every cluster's V_k.

NOTE (perf): when the batched 'cond_median_logh_multi()' C kernel is available, this calls 'invfunc()' ONCE across ALL K clusters simultaneously (see '.prepare_bifactor_batch()'), so each of its ~40 bisection steps costs a single '.C()' call spanning whichever clusters haven't converged yet – instead of the ~40*K small '.C()' calls the original per-cluster loop made. Falls back automatically to that original per-cluster loop (via 'condMedianV_bifactor()', still using the single-cluster kernel if available) when the multi-cluster kernel can't be used.

Usage

condMediansBifactor(
  Y,
  X = NULL,
  clu,
  model,
  coef,
  nq = 31,
  use_c = TRUE,
  kernel_path = "bifactor_kernel.c"
)

Arguments

Value

data.frame with one row per cluster: 'cluster', 'V' (conditional posterior median)

Link to copula parameter

Description

Computes the copula parameters given a linear combination of covariates.

Usage

cpar2link(cpar, family = "clayton")

Arguments

Details

NOTE (fix): the previous version grouped '"clayton"' together with ‘"huesler-reiss"'/'"galambos"' under 'th = log(cpar)'. That’s the correct inverse for huesler-reiss/galambos (whose forward link in 'linkCop()' is 'cpar = 1/expf = exp(th)'), but NOT for clayton, whose forward link is instead 'cpar = 2*tau/(1-tau)' with 'tau = 1/(1+exp(-th))' (the SAME tau clayton shares with gumbel/joe/ plackett, just combined into a different final formula) – inverting that gives 'th = log(cpar/2)', not 'log(cpar)'. The old formula was systematically wrong by a factor of exactly 2 in 'cpar' (equivalently 'log(2)' in 'th') for every clayton value: 'linkCop(log(cpar), "clayton")$cpar' evaluates to '2*cpar', not 'cpar'. Also added the '"fgm"' case, previously missing entirely even though 'linkCop()' supports it ('cpar = tanh(th)', the same link as gaussian/t).

Value

Examples

th = cpar2link(0.5,"gaussian")

Laplace pdf

Description

This function computes the pdf of Laplace distribution

Usage

dlap(x)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Exponential cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

expcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Geometric with p = 1/(1+exp(-th)) cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

geomcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Inverse function

Description

This function is used to get the inverse of a monotonic function on (0,1), depending on parameters, and using the bisection method

Usage

invfuncNew(q, func, th, lb = 1e-12, ub = 1 - 1e-12, tol = 1e-08, nbreak = 40)

Arguments

Details

NOTE: renamed from 'invfunc()' to 'invfuncNew()'. This package Imports 'CopulaGAMM', which also exports a function named 'invfunc()' – and unlike 'pcond()'/'linkCop()' (which turned out to be identical to CopulaGAMM's own and were dropped in favor of importing them directly), THIS function contains a real, necessary bug fix CopulaGAMM's release doesn't have (see the NOTE (fix) comment below: the three 'func0()' calls used to pass the internal 'iconv' convergence-tracking vector in place of ‘th', which is what caused the "object ’hder' not found" / node-stack-overflow bugs earlier in this project's history). Every caller in this package (‘qcondNew()'’s family-specific helpers, 'condMedianV_bifactor()', 'condMediansBifactor()', 'EstBifactorTwoStep()') has been updated to call 'invfuncNew()' explicitly, so none of them can silently fall back to CopulaGAMM's still-buggy original.

Value

Author(s)

Pavel Krupskii

Laplace cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

lapcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Parameter values from link parameter for margins

Description

Computation of parameter values for margins from the associated link parameters

Usage

link2parM(th, marg)

Arguments

Value

parM Associated margin parameters

Examples

p=link2parM(.43,"bernoulli")

Margins cdf/pdf

Description

This function computes the cdf and pdf

Usage

margins0(z, th, model, dfM)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Examples

out = margins0(c(1:10),c(2.5,0),"poisson")

Negative binomial cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

nbinom1cpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Negative binomial cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

nbinomcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

normal cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

normcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Number of parameters

Description

Gives the number of paramters for a given idstribution

Usage

param.dim(family)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, March 9, 2026

Examples

dd = param.dim("normal")

Laplace cdf

Description

This function computes the cdf of Laplace distribution

Usage

plap(x)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Poisson cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

poiscpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Inverse conditional cdf

Description

This function computes the quantile of conditional cdf C(U|v) for a copula C

Usage

qcondNew(w, v, family, cpar, rot = 0)

Arguments

Details

NOTE: renamed from 'qcond()' to 'qcondNew()'. This package Imports 'CopulaGAMM', which also exports a function named 'qcond()' – keeping the same name here would create a real namespace collision (not just a load-order ambiguity, the way it was when these were loose sourced scripts), and this version fixes real bugs the 'CopulaGAMM' release doesn't have (see 'qcondgum()'/'qcondjoe()'/'qcondgal()'/ ‘qcondhr()'’s own NOTE (fix) comments below), so silently falling back to 'CopulaGAMM::qcond()' anywhere this package forgets to qualify the call would reintroduce those bugs. Every caller in this package (e.g. 'Sim2Cluster()') has been updated to call 'qcondNew()' explicitly.

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Examples

U = qcondNew(0.1,0.2,"gaussian",0.87)

Inverse clayton

Description

Inverse clayton

Usage

qcondcla(w, v, th)

Arguments

Value

Inverse FGM (B10)

Description

Inverse FGM (B10)

Usage

qcondfgm(w, v, th)

Arguments

Value

Inverse Frank

Description

Inverse Frank

Usage

qcondfra(w, v, th)

Arguments

Value

Inverse Galambos

Description

Inverse Galambos

Usage

qcondgal(w, v, th)

Arguments

Value

Inverse Gumbel

Description

Inverse Gumbel

Usage

qcondgum(w, v, th)

Arguments

Value

Inverse Huesler-Reiss

Description

Inverse Huesler-Reiss

Usage

qcondhr(w, v, th)

Arguments

Value

Inverse Joe

Description

Inverse Joe

Usage

qcondjoe(w, v, th)

Arguments

Value

Inverse Gaussian

Description

Inverse Gaussian

Usage

qcondnor(w, v, th)

Arguments

Value

Inverse Plackett

Description

Inverse Plackett

Usage

qcondpla(w, v, th)

Arguments

Value

Inverse Student

Description

Inverse Student

Usage

qcondt(w, v, th)

Arguments

Value

Laplace quantile

Description

This function computes the qunatiles of Laplace distribution

Usage

qlap(u)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Quantiles of margins

Description

This function computes the quantiles of a given margin family

Usage

qmargin(U, family, param, size = NULL, shape = NULL, dfM = NULL)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Examples

out = qmargin(0.5,"normal", param=c(10,2))

Student cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

tcpdf(z, th, df)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022

Weibul cdf/pdf and ders

Description

This function computes the cdf, pdf, and associated derivatives

Usage

weibcpdf(z, th)

Arguments

Value

Author(s)

Pavel Krupskii and Bruno N. Remillard, January 20, 2022