R/nmf.ffb.test.R

In nmfkc: Non-Negative Matrix Factorization with Kernel Covariates

#' Test the feed-forward null of an NMF-FFB fit
#'
#' Calibrated test of the hypothesis that the outcomes need no indicator-level
#' feedback: the feed-forward model \eqn{Y_1 = X(\Theta_2 Y_2 + U) + \mathcal{E}}
#' with a full random-effect covariance against the feedback model with
#' \eqn{\Theta_1 \ne 0}.  This is the only inferential step of the procedure
#' (see the workflow below); \code{\link{nmf.ffb.inference}} is for the
#' uncertainty of the coefficients \emph{after} the null has been rejected, and
#' \code{\link{nmf.ffb.diagnostics}} for what the selected feedback looks like.
#'
#' The likelihood ratio does not have a \eqn{\chi^2} null distribution: the
#' constraint \eqn{\Theta_1 \ge 0} puts the null on the boundary of the
#' parameter space, and the support of \eqn{\Theta_1} is chosen from the same
#' data.  It is therefore calibrated by a parametric bootstrap from the fitted
#' feed-forward null,
#' \deqn{Y_1^{*} = (\hat X(\hat\Theta_2 Y_2 + U^{*}) + \mathcal{E}^{*})_{+},
#'       \quad U^{*} \sim N(0, \hat\Phi), \; \mathcal{E}^{*} \sim N(0, \hat\Psi),}
#' with negative draws clipped at zero, and the reported \eqn{p}-value is
#' \eqn{(1 + \#\{LR^{*} \ge LR\}) / (1 + B_{ok})}.
#'
#' \subsection{What is re-applied to each replicate}{
#' \code{calibration = "procedure"} (the default) re-applies the \emph{whole
#' procedure}: stage 1 is re-run on each \eqn{Y_1^{*}}, the basis is
#' re-estimated, the exclusion restriction is re-derived from that basis, and
#' stage 2 follows.  The \eqn{p}-value is then the operating characteristic of
#' the complete exploratory procedure, which is what is applied to a fresh data
#' set.  This is the calibration to use.
#'
#' \code{"conditional"} holds \eqn{\hat X} and the exclusion restriction at
#' their fitted values and re-runs stage 2 only.  It is valid only if basis and
#' restriction are independent of the \eqn{Y_1} being tested, which they are not when
#' all three come from the same sample; it is provided for comparison with that
#' practice, and it is anti-conservative.
#'
#' Sample splitting -- basis and restriction from one half of the units, test on
#' other -- also removes the dependence, and was offered here in 0.9.8.  It was
#' withdrawn in 0.9.8: measured against \code{"procedure"} it has the same size and
#' lower power, because the test uses \eqn{N/2} units, and it cannot be run at
#' all when the halves are too small for stage 1.
#' }
#'
#' \subsection{Reading the result}{
#' Report \code{LR.p.boot} together with \code{prob.select.null} and, for
#' \code{calibration = "procedure"}, \code{C1.restriction.change.rate}.  The last two say
#' whether the exclusion restriction is determined well enough for the test to
#' mean anything: if the dominant factor of an outcome moves from one null
#' replicate to the next, a previously blocked entry becomes free and a feedback
#' coefficient can absorb loading structure that the feed-forward model
#' attributes to \eqn{X}.  A large \code{prob.select.null} with a large null
#' quantile is the signature of a data set in which the procedure finds feedback
#' whether or not there is any.
#' }
#'
#' @param object A fit from \code{\link{nmf.ffb}} with \code{method = "fiml"}.
#' @param Y1,Y2 The endogenous and exogenous matrices the fit was computed from
#'   (variables in rows, units in columns).
#' @param B Number of bootstrap replicates (default 1000).  The smallest
#'   reportable \eqn{p}-value is \eqn{1/(1+B)}.
#' @param calibration What is re-applied to each null replicate; see above.
#'   \code{"conditional"} is selected automatically, with a warning, when the
#'   fit used a basis or an exclusion restriction supplied by the caller, since there is then
#'   nothing to re-estimate.
#' @param seed Base seed; replicate \code{b} uses \code{seed + b}.  The caller's
#'   random stream is restored on exit.
#' @param ... Passed to the bootstrap: \code{cores} (or \code{ncores}) to run
#'   the replicates in parallel, defaulting to \code{getOption("mc.cores", 1L)}
#'   as elsewhere in the package; \code{fiml.maxit} and \code{factr} to control
#'   the optimizer in the replicates; \code{print.trace = TRUE} for progress.
#'
#' @return An object of class \code{"nmf.ffb.test"}:
#' \describe{
#' \item{LR, LR.df}{Observed \code{c(full, selected)} likelihood ratios and the
#'   number of free feedback entries left by the exclusion restriction.}
#' \item{LR.p.boot}{Calibrated \eqn{p}-values of the two statistics.}
#' \item{LR.null.quantile}{The 0.95 quantiles of the null distributions.}
#' \item{LR.boot}{\eqn{B \times 2} matrix of the null replicates.}
#' \item{prob.select.null}{Null false-selection rate: the probability that BIC
#'   retains at least one feedback entry when the feed-forward model is true.}
#' \item{C1.restriction.change.rate}{(\code{"procedure"} only) share of null replicates in which
#'   the exclusion restriction moved.}
#' \item{LR.boot.df}{Free entries per null replicate (varies under \code{"procedure"}).}
#' \item{LR.boot.n.ok, LR.boot.n.nonconv}{Usable replicates, and how many did not
#'   meet the optimizer tolerance.}
#' }
#'
#' @section Workflow:
#' \preformatted{
#' ecv <- nmf.ffb.ecv(Y1, Y2, rank = 1:5)   # 1. choose Q
#' fit <- nmf.ffb(Y1, Y2, rank = ecv$rank)  # 2. estimate; BIC selects the support
#' tst <- nmf.ffb.test(fit, Y1, Y2)         # 3. test the feed-forward null  <- only inference
#' dgn <- nmf.ffb.diagnostics(fit)          # 4. cycles, spectral radius, identifiability
#' inf <- nmf.ffb.inference(fit, Y1, Y2)    # 5. intervals for the retained entries
#' }
#'
#' @seealso \code{\link{nmf.ffb}}, \code{\link{nmf.ffb.ecv}},
#'   \code{\link{nmf.ffb.inference}}, \code{\link{nmf.ffb.diagnostics}}
#' @export
nmf.ffb.test <- function(object, Y1, Y2,
                         B = 1000L,
                         calibration = c("procedure", "conditional"),
                         seed = 123L,
                         ...) {
  if (base::is.null(object$X) || base::is.null(object$C1) || base::is.null(object$C2))
    base::stop("object must be a fit from nmf.ffb().")
  if (!base::identical(object$method, "fiml"))
    base::stop("nmf.ffb.test() applies to likelihood-based fits: use nmf.ffb(..., method = \"fiml\") ",
               "(the default). The multiplicative-update estimator has no likelihood to compare.")
  if ("ci.level" %in% base::names(base::match.call()))
    base::stop("`ci.level` was renamed to `boot.level` in 0.9.8.", call. = FALSE)
  calibration <- base::match.arg(calibration)
  .nmf.ffb.inference.fiml(object, Y1, Y2, B = B, seed = seed,
                          calibration = calibration, what = "test", ...)
}