Unbiased estimation of covariance matrix for multiply censored data

Chemical analyses of environmental samples are often censored below at reporting limits or various detection/quantitation limits. The latter can vary, usually in proportion to the values of other variables. For example, a sample with a high concentration of one compound might need to be diluted for analysis, resulting in proportional inflation of the censoring limits for all other compounds analyzed at the same time in that sample. As another example, sometimes the presence of a compound can alter the response of the test to other compounds (a "matrix interference"); when this is detected by the laboratory, it will inflate its reporting limits accordingly.

I am seeking a practical way to estimate the entire variance-covariance matrix for such datasets, especially when many of the compounds experience more than 50% censoring, which is often the case. A conventional distributional model is that the logarithms of the (true) concentrations are multinormally distributed, and this appears to fit well in practice, so a solution for this situation would be useful.

(By "practical" I mean a method that can reliably be coded in at least one generally available software environment like R, Python, SAS, etc., in a way that executes quickly enough to support iterative recalculations such as occur in multiple imputation, and which is reasonably stable [which is why I am reluctant to explore a BUGS implementation, although Bayesian solutions in general are welcome].)

Many thanks in advance for your thoughts on this matter.

• Just so I understand correctly the censoring issue: When you dilute a sample the concentration of a compound falls so low that the test instrument can fail to detect its presence. Is that an accurate re-phrasing of the censoring problem? – user28 Aug 17 '10 at 14:34
• Yes, that is correct: dilution by a factor of D increases all detection limits by a factor of D as well. (The matrix interference issue is more difficult to quantify and the general situation is extremely complex. To simplify this, the conventional model is that a suite of tests on one sample yields a vector (x, ..., x[k]) where the x[i] are either real numbers or are intervals of reals, typically with left endpoint at -infinity; an interval identifies a set in which the true value is assumed to lie.) – whuber Aug 17 '10 at 18:20
• Why would the detection limits go up? Are they not a feature of the test instrument rather than that of the sample being tested? – user28 Aug 17 '10 at 21:33
• As an example, suppose an instrument's detection limit is 1 microgram per Liter (ug/L). A sample is diluted 10:1 (with great precision, so we don't worry about error here) and the instrument reads "<1"; that is, nondetectable, for the diluted sample. The laboratory infers that the concentration in the sample is less than 10*1 = 10 ug/L and reports it as such; that is, as "<10". – whuber Aug 18 '10 at 14:06
• @amoeba I see I should have explained those things in the question itself. The answers are: PCA; dimensionality will vary from 3 to a few hundred; sample sizes always greatly exceed dimensionality but censoring rates may be very high (being able to handle up to 50% is needed and up to 95% is desirable). – whuber Dec 1 '17 at 21:40

I have not full internalized the issue of matrix interference but here is one approach. Let:

$Y$ be a vector that represents the concentration of all the target compounds in the undiluted sample.

$Z$ be the corresponding vector in the diluted sample.

$d$ be the dilution factor i.e., the sample is diluted $d$:1.

Our model is:

$Y \sim N(\mu,\Sigma)$

$Z = \frac{Y}{d} + \epsilon$

where $\epsilon \sim N(0,\sigma^2\ I)$ represents the error due to dilution errors.

Therefore, it follows that:

$Z \sim N(\frac{\mu}{d}, \Sigma + \sigma^2\ I)$

Denote the above distribution of $Z$ by $f_Z(.)$.

Let $O$ be the observed concentrations and $\tau$ represent the test instrument's threshold below which it cannot detect a compound. Then, for the $i^{th}$ compound we have:

$O_i = Z_i I(Z_i > \tau) + 0 I(Z_i \le \tau)$

Without loss of generality let the first $k$ compounds be such that they are below the threshold. Then the likelihood function can be written as:

$L(O_1, ... O_k, O_{k+1},...O_n |- ) = [\prod_{i=1}^{i=k}{Pr(Z_i \le \tau)}] [\prod_{i=k+1}^{i=n}{f(O_i |-)}]$

where

$f(O_i |-) = \int_{j\neq i}{f_Z(O_i|-) I(O_i > \tau)}$

Estimation is then a matter of using either maximum likelihood or bayesian ideas. I am not sure how tractable the above is but I hope it gives you some ideas.

• Thank you very much for this thought. Indeed, this is a standard and well-documented approach to multiple censoring. One difficulty lies in its intractability: those integrals are notoriously difficult to compute. There's a modeling problem lurking here, too: the value of d is usually positively correlated with Y, as implied by the first paragraph of my description. – whuber Aug 20 '10 at 14:22

Another more computationally efficient option would be to fit the covariance matrix by moment matching using a model that has been called the "dichomized Gaussian", really just a Gaussian copula model.

A recent paper from Macke et al 2010 describes a closed form procedure for fitting this model which involves only the (censored) empirical covariance matrix and the calculation of some bivariate normal probabilities. The same group (Bethge lab at MPI Tuebingen) has also described hybrid discrete/continuous Gaussian models which are probably what you want here (i.e., since the Gaussian RVs aren't fully "dichotomized" -- only those below threshold).

Critically, this is not an ML estimator, and I'm afraid I don't know what its bias properties are.

• @jp Thank you: I will look into this. (It may take some time...) – whuber Jul 6 '11 at 14:48

How many compounds are in your sample? (Or, how big is the covariance matrix in question?).

Alan Genz has some very nice code in a variety of languages (R, Matlab, Fortran; see here) for computing integrals of multivariate normal densities over hyper-rectangles (i.e., the kinds of integrals you need to evaluate the likelihood, as noted by user28).

I've used these functions ("ADAPT" and "QSIMVN") for integrals up to around 10-12 dimensions, and several functions on that page advertise integrals (and associated derivatives you might need) for problems up to dimension 100. I don't know if that's enough dimensions for your purposes, but if so it could presumably allow you to find maximum likelihood estimates by gradient ascent.

• Oh, sorry—I'm new here and didn't notice how long ago this was posted—probably too late to be of much help! – jpillow Jul 6 '11 at 0:40
• @jp This is an ongoing important problem, so the elapsed time between the question and answer is of little consequence. Thank you for replying! – whuber Jul 6 '11 at 14:46