v1.8.0 — mbco bootstrap checks, exact DOP at |rho| = 1, R 4.1 S7 fix

Pre-release from main (merge commit 7e7528c). CRAN still serves 1.6.1; install with pak::pak("Data-Wise/rmediation@v1.8.0"). The SEMFit mbco() route is experimental.

• S7 property access now works under the declared R 4.1 minimum.

Testing

• Expanded public OpenMx MBCO argument, bootstrap-routing and optimizer-fallback
  tests, including result-field mutation checks and fitted RAM smoke tests.
  Coverage CI now retains line/expression summaries and uncovered-line artifacts.
• The SEMFit MBCO tests and the mbco() example pass fixed_x and
  meanstructure to medfit::fit_sem() explicitly, at today's values, so
  medfit 0.14.0's default-change warning does not appear (medfit flips these
  defaults in 0.15.0). Results are unchanged; the example now needs medfit

  = 0.9.0 (#85).

Removed

• medci_prototype() is removed. It was an exported S7-core prototype of
  medci() (shipped in 1.6.1) that nothing in the package called. Use
  medci(), or ci() on a ProductNormal, which return the same intervals.

Bug fixes

• mbco(type = "parametric") now checks every bootstrap replicate like type = "semi" does.
  It used to hand the bootstrap to OpenMx::mxCompare(boot = TRUE), so a replicate null fit
  that violated its constraint (for example OpenMx/OpenMx#433 dropping a zero-Jacobian
  equality at a replicate's start) entered the null distribution unchecked and understated its
  statistic. The package now simulates each replicate from the fitted h0 with
  OpenMx::mxGenerateData(), fits h0 and h1, refits a violated h0 constraint, and stops
  with an error of class RMediation_bootstrap_infeasible when it cannot be repaired. Behavior
  changes: a replicate with a non-finite statistic or an h1 fit worse than its h0 fit now stops the
  call (mxCompare dropped such replicates; the last check also applies to type = "semi"; a
  replicate that does not converge is dropped within a 5% bound, see the next entry), and a model without a
  single raw data set (a covariance-data or multigroup model) is refused before sampling. The
  observed statistic, df and the result shape are
  unchanged; seeded p-values differ because the replicates are drawn by a different loop.

• The DOP CDF (pprodnormal(), qprodnormal(), medci()) is exact at |rho| = 1: it now uses
  Meeker & Escobar's (1994) quadratic-root solution instead of integrating a step function, which
  was wrong by up to 4e-5 and could stop with "the integral is probably divergent". For
  |rho| < 1 the integrand is evaluated from the factored quadratic (no cancellation at large
  means) with the roots as breakpoints, and a roundoff report from integrate() relaxes the
  tolerance to 1e-8, then 1e-6, instead of failing. All 495 evaluations in the article's own test
  domain (means of +/-3e5, |rho| up to 0.99999) now succeed; 10 used to stop with an integrator
  error.

• mbco_product_null() stops, with class RMediation_h1_nonconverged, when the full model h1 did
  not converge or was never run. mbco() warns, with class RMediation_fit_nonconverged, when a
  supplied h0 or h1 carries an unaccepted OpenMx status; the result is unchanged.

• Point masses now have the right-continuous CDF. pprodnormal3() returned 0.5 at q == 0 for a
  product identically equal to 0 (and for a single zero-variance, zero-mean first component);
  it now returns 1, matching P(XYZ <= q). pprodnormal() with a zero standard error follows
  the same rule (q >= c lower, q < c upper). Continuous zero-mean distributions still give
  0.5 at 0.

• pprodnormal3() rejects a covariance that is not finite and positive semidefinite, including
  negative variances at very different scales, instead of clamping negative variances to zero.

• pprodnormal(lower.tail = FALSE) and qprodnormal(lower.tail = FALSE) with
  type = "dop" now integrate the upper tail directly instead of computing
  1 - P(XY <= q), which cancelled below about 1e-10 (an upper quantile at
  p = 1e-12 round-tripped to a relative error of 7e-5, and smaller tails
  could not be resolved). The quantile root finder's tolerance is tightened
  from the uniroot() default (1.2e-4) to 1e-12, and the integrand handles
  |rho| = 1 (a point-mass conditional) without a 0/0. Probabilities below
  about 1e-22 remain unresolved (the integration window leaves out 1.5e-23).

• medci(type = "mc") reported MC.Error as SE / n.mc; it is now
  SE / sqrt(n.mc), the standard error of the simulated mean, as documented
  elsewhere. Intervals and estimates are unchanged.

• medci() Monte Carlo and plotting no longer stop with a Cholesky error when
  the covariance is singular (rho = +/-1 or a zero standard error); type = "dop"/"all" with a zero standard error now fails with an actionable
  message instead of "non-finite function value".

• plot/plotCI must be a single TRUE/FALSE in medci() and ci()
  (a vector or NA previously failed with an obscure condition error).

• ci(<ProductNormal>, type = "dop") shared the large-mean defect of #60 through
  a second copy of the integral (it stopped with "please use type='MC'" at
  mu/se of about 20 where medci() already worked). The S7 core now calls
  qprodnormalMeeker(), so both paths use one integral.

• pprodnormal3() (and cdf()/dist_quantile() on ProductNormal3) returned 0
  for every q, or an error, when two or three components had zero variance.
  It now returns the exact answer: a point mass at prod(mean) when all three
  are fixed, and a scaled normal when exactly one is random (#64).

• medci() was defined twice, in R/medci.R and R/medci_legacy_backup.R
  (identical code). The backup sorted later and was the live definition, and
  ?medci merged both sets of documentation into duplicated Usage, Value and
  Details sections. The backup file is removed; behavior is unchanged. A test
  now fails if any top-level function is defined in more than one R/ file.

• pprodnormal() and qprodnormal() with type = "dop" no longer collapse at
  large standardized means (pprodnormal(400, 20, 20, 1, 1) returned 8e-7
  instead of about 0.507) and are accurate in far tails (previously wrong by up
  to 4e-6 absolute). Both now share one integral over the support of the
  conditioning variable at a tight tolerance. The CDF is also clamped to
  [0, 1], so lower.tail = FALSE can no longer be slightly negative (#60).

New features

• dprodnormal3(): exact density of the product of three normal variables with any mean vector and
  positive-definite covariance (the counterpart of pprodnormal3()), computed as a
  two-dimensional integral that conditions on the first factor and reuses the dprodnormal()
  integral for the other two. Matches the derivative of pprodnormal3() to about 7 digits.

• dprodnormal() now resolves very small z at nonzero means (it previously could stop with
  "the integral is probably divergent").

• dprodnormal() is the exact density of the product of two normal random
  variables (any means, any -1 < rho < 1), completing the d/p/q family. It
  uses the one-dimensional integral of Rohatgi (Glen et al., 2004) in general
  and the Bessel-K closed form of Nadarajah and Pogany (2016) when both means
  are zero; the two agree to near machine precision and each is the other's
  test oracle. The density is Inf at exactly 0 (a log singularity).

• pprodnormal() and qprodnormal() accept a vector q / p and return one
  probability / quantile per element (a list of two vectors for type = "all").
  For type = "MC" all elements share one simulated sample, so the estimates
  are mutually consistent. Scalar calls are unchanged.

• ml_omitted_confounder() adds a sensitivity analysis for an omitted
  confounder of the mediator-outcome relation in two-level (1-1-1) mediation,
  after Tofighi and Kelley (2016). It returns an MLOmittedConfounder object
  with print, summary, plot and as.data.frame methods. The default
  screen keeps only positive-definite grid points, so it does not reproduce the
  published tables; use screen = "none" for those.

• mbco(h0, h1, type = "asymp") accepts two fitted medfit SEMFit objects
  with medfit >= 0.8.0. It returns MBCOResult using stats::anova() without
  refitting or requiring OpenMx. The restricted fit must come first; df is the
  actual difference. Active tested bounds refuse an invalid chi-square p-value.
  Bootstrap types, mixed classes and unused fitting controls error. Existing
  OpenMx paths and CI integration with medfit >= 0.2.0 are preserved.

• medfit remains optional in Suggests >= 0.2.0. The SEMFit route needs the
  GitHub/r-universe medfit release; CRAN medfit 0.3.2 cannot exercise it. This
  development change does not lift the CRAN submission or development-to-main hold.

• OpenMx mbco() now checks the constraints of h0 and h1 at the fitted
  solution (issue #34). The OpenMx default feasibility tolerance (0.05) can leave
  a constrained fit off the null set and understate the statistic. A returned
  number now changes in that case: a violated h0 (scaled violation above
  1e-6) is refit once at tolerance 1e-8 for all three types, the refit model
  supplies the observed statistic and the bootstrap, a warning of class
  RMediation_h0_refit reports the change, and the new refit slot of
  MBCOResult records it (NULL when no refit was needed). Clean fits and
  constraint-free models are unchanged. If the refit fails, does not converge or leaves
  the constraint violated, mbco() stops with an error of class
  RMediation_h0_refit_failed naming the constraint, instead of returning a
  statistic from a known-invalid fit. A
  constraint that cannot be evaluated now errors. See the wrong-branch warning below.

• OpenMx mbco() now also warns, with class RMediation_constraint_not_enforced and right
  after RMediation_h0_refit, when the refit ran and h0's fit was identical to
  h1's (statistic 0 before the refit) although a constraint of h0 was violated
  (issue #84). The constraint may not have been enforced. OpenMx (2.19.1 and later)
  silently drops an equality constraint whose Jacobian is zero at the starting values
  (OpenMx/OpenMx#433), for example b1 * b2 == 0 started at zero, and a violation
  within the feasibility tolerance looks the same; a fitted h0 keeps no starting
  values, so the warning cannot say which. The statistic and the refit slot are
  unchanged (the seeded example goes from 0 to 0.0878); for a product null use
  mbco_product_null(). Bootstrap replicates are refit and checked the same way
  for both bootstrap types (see the type = "parametric" entry under Bug fixes).

• mbco(type = "semi") and mbco(type = "parametric") now drop a bootstrap replicate whose
  fits did not converge instead of stopping at the first one (issue #92). Up to 5% of the R
  replicates may be dropped; above that the call stops with an error of class
  RMediation_bootstrap_too_many_dropped (with fewer than 20 replicates none may be dropped).
  The new dropped slot of MBCOResult reports the count (NULL for results without a
  bootstrap), and the p-value is the proportion of the kept replicates strictly above the
  observed statistic, with no (1 + k)/(R + 1) correction: the rule of
  OpenMx::mxCompare(boot = TRUE), whose denominator was also the kept replicates. Only
  non-convergence is dropped; an unrepairable constraint violation
  (RMediation_bootstrap_infeasible), a non-finite statistic and an h1 fit worse than its h0
  fit still stop the call. ?mbco states the rule.

• OpenMx mbco() now warns (class RMediation_constraint_branch) when a
  feasible, converged h0 is not the constrained optimum of an equality
  constraint: fixing one of the parameters the constraint depends on at 0 gives a
  converged, feasible fit with -2LL lower by more than 1e-4 (issue #35). It
  catches the single-constraint product null that lands on the larger-statistic
  branch (memory_exp: 221.05 against 0.0831) without parsing the constraint.
  It runs once per call for every type, costs one extra fit per involved
  parameter (a constraint with more than 5 is skipped with a message), changes no
  result, and cannot fail the call. In a simulation (60 seeds at each of three
  effect sizes) it flagged all 28 wrong-branch fits with no false positives.
  mbco_product_null() avoids the problem.

• mbco_product_null(h1, paths, ...) tests a product null (a*b = 0, or a
  serial product) with OpenMx as an intersection-union test: it fixes each path
  at 0 in a copy of the fitted full model, runs mbco() on each, and returns the
  result of the branch with the largest p-value, with the new optional
  MBCOResult@branch slot naming it. Unlike a single nonlinear constraint such
  as ind == 0, it cannot land on the larger-statistic branch (issue #35; on
  memory_exp it returns 0.0831 where one constraint gives 221.05). The
  p-value is conservative when more than one path is 0. A branch whose fit
  does not converge is an error naming the branch.

• For numeric/lavaan formula CIs, MCError now uses SE / sqrt(n.mc),
  the sampling standard error of the simulated mean. The previous SE / n.mc
  understated simulation uncertainty. This is not an endpoint error estimate.

• ci(type = "all", H0 = TRUE) now forwards H0, mu0 and Sigma0 to
  its Monte Carlo component, preserving the alternative/null result pair
  rather than silently ignoring the requested null distribution.

• OpenMx MBCO now rejects equal-df and reversed model pairs, and nonfinite
  or negative observed statistics, before bootstrap sampling. Previously the
  semiparametric route could report a significant result with df 0 or -1.
  The restricted model must be supplied first; substantive nesting is still
  the analyst's responsibility.

• OpenMx bootstrap workers now restore global options on success and error.
  Semiparametric bootstrap rejects failed refits and nonfinite statistics, with
  explicit errors for singular sample or null covariances. Parametric bootstrap
  rejects invalid reported statistics, degrees of freedom and p-values.

• mbco() no longer rejects fitted OpenMx RAM models with "h0 must be an
  MxModel". The check compared class(h0) to "MxModel", which fails for
  subclasses such as MxRAMModel; it now uses methods::is(). Verified for
  type = "asymp" only; the bootstrap types ("parametric", "semi") have
  not been exercised on fitted RAM models (rmediation#35).

• The ?mbco example had shifted path labels (endVar started with "x", so
  the labels a1, a2, cp landed on the wrong paths) and passed a tibble to
  mxData(). Both are fixed. The example now tests a1 * b1 = 0 with one
  linear null per path (a1 == 0, b1 == 0) and reports the smaller statistic
  (0.0831), because a single nonlinear ind1 == 0 null can converge to the
  wrong branch (221.05) (rmediation#35).

Documentation

• ?mbco now has a Details section on testing a product null (a*b = 0): fit
  one linear null per path and report the smaller statistic, because a single
  nonlinear null can converge to the wrong branch. It also says to build h0
  from the unfitted model, since a constraint refit from the fitted model's
  estimates can stop at an infeasible point (diffLL 0) (rmediation#35).

Tests

• Added a known-answer OpenMx fixture on memory_exp (linear nulls
  a1 == 0: 221.045546 and b1 == 0: 0.083100) and an end-to-end mbco()
  test against it.