Skip to contents

Maximises \(\ell(\boldsymbol{\rho})+\tfrac12\log\det I_{\boldsymbol{\rho}}\) on the ILR chart. The penalty is Jeffreys' prior for the expected Fisher information of one bulk column; extra replicate columns add only a constant, so they do not change the penalty gradient. The returned Wald covariance is that information multiplied by the number of columns and evaluated at the penalised estimate.

The penalty is differentiated numerically. The log-likelihood score stays analytic. This is not an exact removal of the \(O(N^{-1})\) bias: the convolution is a curved exponential family.

Usage

fit_firth_decovart(
  y,
  mean_signature_matrix,
  Sigma,
  p_start = NULL,
  epsilon = 1e-08,
  itmax = 200L
)

Arguments

y

Numeric vector of length \(G\), or \(G\times N\) matrix of replicate bulk profiles that share one composition; their log-likelihoods are summed (the pooled, population-level fit).

mean_signature_matrix

Numeric matrix \(\boldsymbol{\mu}\in\mathcal{M}_{G\times J}\) (plug-in means).

Sigma

Array of cell-type covariances in \(\mathcal{M}_{G\times G\times J}\).

p_start

Start on the simplex. NULL uses the unpenalised MLE, pulled into the open simplex when that MLE lies on a face.

epsilon, itmax

Relative convergence tolerance and iteration budget passed to stats::nlminb().

Value

A list with coefficients, rho, loglik (unpenalised), penalty, penalised_loglik, vcov, vcov_rho and convergence.

References

Firth D (1993). “Bias Reduction of Maximum Likelihood Estimates.” Biometrika, 80(1), 27–38. doi:10.1093/biomet/80.1.27 .

Examples

mu <- matrix(c(20, 40, 15, 40, 20, 25), nrow = 3)
colnames(mu) <- paste0("ct", 1:2)
Sigma <- array(c(diag(3), diag(3)), dim = c(3, 3, 2))
y <- drop(mu %*% c(0.6, 0.4))
fit_firth_decovart(y, mu, Sigma)$coefficients
#>       ct1       ct2 
#> 0.5986526 0.4013474