ppi: Estimation of Polynomially-Tilted Pairwise Interaction (PPI) Model

Description

Estimates the parameters of the Polynomially-Tilted Pairwise Interaction (PPI) model scealy2023scscorematchingad for compositional data. By default ppi() uses cppad_closed() to find estimate. For many situations a hard-coded implementation of the score matching estimator is also available.

For a given parameter vector evalparam, ppi_smvalues() computes the score matching discrepancy, the gradient and the Hessian of the score matching discrepancy (see smvalues()) and the gradient offset of the score matching discrepancy (see quadratictape_parts() and tape_gradoffset()).

Usage

ppi(
  Y,
  paramvec = NULL,
  trans,
  method = "closed",
  w = rep(1, nrow(Y)),
  constrainbeta = FALSE,
  bdryw = "ones",
  acut = NULL,
  bdrythreshold = 1e-10,
  shiftsize = bdrythreshold,
  approxorder = 10,
  control = list(tol = 1e-15, checkgrad = TRUE),
  paramvec_start = NULL
)
ppi_smvalues(
  Y,
  paramvec = NULL,
  evalparam,
  trans,
  method = "closed",
  w = rep(1, nrow(Y)),
  bdryw = "ones",
  acut = NULL,
  bdrythreshold = 1e-10,
  shiftsize = bdrythreshold,
  approxorder = 10,
  average = TRUE
)

Value

ppi() returns: A list of est, SE and info.

est contains the estimates in vector form, paramvec, and as $A_L$, $b_L$ and $\beta$.
SE contains estimates of the standard errors if computed. See cppad_closed().
info contains a variety of information about the model fitting procedure and results.

ppi_smvalues() returns a list of

obj the score matching discrepancy value
grad the gradient of the score matching discrepancy
hess the Hessian of the score matching discrepancy
offset gradient offset (see quadratictape_parts())

Arguments

Y: A matrix of measurements. Each row is a compositional measurement (i.e. each row sums to 1 and has non-negative elements).
paramvec: Optionally a vector of the PPI models parameters. NA-valued elements of this vector are estimated and non-NA values are fixed. Generate paramvec easily using ppi_paramvec(). If NULL then all elements of $A_L$, $b_L$ and $\beta$ are estimated.
trans: The name of the transformation of the manifold in Hyvärinen divergence (See scorematchingtheory): "clr" (centred log ratio), "alr" (additive log ratio), "sqrt" or "none".
method: "closed" uses CppAD to solve in closed form the a quadratic score matching discrepancy using cppad_closed(). "hardcoded" uses hardcoded implementations. "iterative" uses cppad_search() (which uses CppAD and optimx::Rcgmin()) to iteratively find the minimum of the weighted Hyvärinen divergence.
w: Weights for each observation, if different observations have different importance. Used by Windham() and ppi_robust() for robust estimation.
constrainbeta: If TRUE, elements of $\beta$ that are less than -1 are converted to -1 + 1E-7.
bdryw: The boundary weight function for down weighting measurements as they approach the manifold boundary. Either "ones", "minsq" or "prodsq". See details.
acut: The threshold $a_c$ in bdryw to avoid over-weighting measurements interior to the simplex
bdrythreshold: iterative or closed methods only. For measurements within bdrythreshold of the simplex boundary a Taylor approximation is applied by shifting the measurement shiftsize towards the center of the simplex.
shiftsize: iterative or closed methods only. For measurements within bdrythreshold of the simplex boundary a Taylor approximation is applied by shifting the measurement shiftsize towards the center of the simplex.
approxorder: iterative or closed methods only. Order of the Taylor approximation for measurements on the boundary of the simplex.
control: iterative only. Passed to optimx::Rcgmin() to control the iterative solver.
paramvec_start: iterative method only. The starting guess for Rcgmin. Generate paramvec_start easily using ppi_paramvec().
evalparam: The parameter set to evaluate the score matching values. This is different to paramvec, which specifies which parameters to estimate. All elements of evalparam must be non-NA, and any parameters fixed by paramvec must have the same value in evalparam.
average: If TRUE return the (weighted average) of the measurements, otherwise return the values for each measurement.

PPI Model

The PPI model density is proportional to $$\exp(z_L^TA_Lz_L + b_L^Tz_L)\prod_{i=1}^p z_i^{\beta_i},$$ where $p$ is the dimension of a compositional measurement $z$, and $z_L$ is $z$ without the final ($p$th) component. $A_L$ is a $p-1 \times p-1$ symmetric matrix that controls the covariance between components. $b_L$ is a $p-1$ vector that controls the location within the simplex. The $i$th component $\beta_i$ of $\beta$ controls the concentration of density when $z_i$ is close to zero: when $\beta_i \geq 0$ there is no concentration and $\beta_i$ is hard to identify; when $\beta_i < 0$ then the probability density of the PPI model increases unboundedly as $z_i$ approaches zero, with the increasing occurring more rapidly and sharply the closer $\beta_i$ is to $-1$.

Details

Estimation may be performed via transformation of the measure in Hyvärinen divergence from Euclidean space to the simplex (inverse of the additive log ratio transform), from a hyperplane to the simplex (inverse of the centred log ratio transform), from the positive quadrant of the sphere to the simplex (inverse of the square root transform), or without any transformation. In the latter two situations there is a boundary and weighted Hyvärinen divergence @Equation 7, @scealy2023scscorematchingad is used. Properties of the estimator using the square root transform were studied by scealy2023sc;textualscorematchingad. Properties of the estimator using the additive log ratio transform were studied by scealy2024ro;textualscorematchingad.

There are three boundary weight functions available:

The function "ones" applies no weights and should be used whenever the manifold does not have a boundary.
The function "minsq" is the minima-based boundary weight function for the PPI model @Equation 12, @scealy2023scscorematchingad $$\tilde{h}(z)^2 = \min(z_1^2, z_2^2, ..., z_p^2, a_c^2).$$ where $z$ is a point in the positive orthant of the p-dimensional unit sphere and $z_j$ is the jth component of z.
The function "prodsq" is the product-based @Equation 9, @scealy2023scscorematchingad $$\tilde{h}(z)^2 = \min(\prod_{j=1}^{p} z_j^2, a_c^2).$$ where $z$ is a point in the positive orthant of the p-dimensional unit sphere and $z_j$ is the jth component of z.

Scealy and Wood @Theorem 1, @scealy2023scscorematchingad prove that minimising the weighted Hyvärinen Divergence is equivalent to minimising $\psi(f, f_0)$ (See scorematchingtheory) when the boundary weight function is smooth or for the functions "minsq" and "prodsq" above when the manifold is the simplex or positive orthant of a sphere.

Hard-coded estimators are available for the following situations:

Square root transformation ("sqrt") with the "minsq" boundary weight function:
- full parameter vector (paramvec not provided)
- paramvec fixes only the final element of $\beta$
- paramvec fixes all elements of $\beta$
- paramvec fixes $b_L = 0$ and provides fixed values of $\beta$
- paramvec fixes $A_L=0$ and $b_L=0$, leaving $\beta$ to be fitted.
Square root transformation ("sqrt") with the "prodsq" boundary weight function:
- paramvec fixes all elements of $\beta$
- paramvec fixes $b_L = 0$ and provides fixed values of $\beta$
- paramvec fixes $A_L=0$ and $b_L=0$, leaving $\beta$ to be fitted.
The additive log ratio transformation ("alr") using the final component on the denominator, with $b_L=0$ and fixed final component of $\beta$.

References

Examples

Run this code

model <- rppi_egmodel(100)
estalr <- ppi(model$sample,
              paramvec = ppi_paramvec(betap = -0.5, p = ncol(model$sample)),
              trans = "alr")
estsqrt <- ppi(model$sample,
              trans = "sqrt",
              bdryw = "minsq", acut = 0.1)

Run the code above in your browser using DataLab