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.