Assume that a $n$-dimensional vector of real valued variables $\mathbf{Y}=(Y_1,...,Y_n)^T$ is given. We know that they are jointly normal which expectation vector $\boldsymbol{\mu}=(\mu_1,...,\mu_n)^T$ and Variance-Covariance Matrix $\boldsymbol{\Sigma}\neq I_{n}\sigma^2$ (i.e. the observations are dependent). So we have $$ \mathbf{Y}\sim N_n(\boldsymbol{\mu},\boldsymbol{\Sigma})$$.

Now assume furthermore that there is a share of censored observations (i.e. we know they are smaller than zero but only observe a zero) $\mathbf{Y}_m=(0,..,0)^T\subset \mathbf{Y}$ and a share of observed observations $\mathbf{Y}_o$. Hence we have $\mathbf{Y}=(\mathbf{Y}_o,\mathbf{Y}_m)^T$

Now if i know that the expectation is a function of some parameters, like $\boldsymbol{\mu}=(X_1^T\beta,...,X_n^T\beta)^T$ and also know the structure of the covariance $\boldsymbol{\Sigma}=\boldsymbol{\Sigma}(\lambda)$ how would I estimate those parameters $\theta=(\beta,\lambda)^T$?

I thought about the following in order to form the likelihood:

$$ f(\mathbf{Y},X,\theta)= f(\mathbf{Y}_o,X_o,\theta)f(\mathbf{Y}_m,X_m,\theta|\mathbf{Y}_o,\mathbf{Y}_m<0)$$

i.e. I would take a multivariate normal density for the uncensored part ($f(\mathbf{Y}_o,X_o,\theta)$) and then estimate the probability of a multivariate normal cdf with conditional mean ($E[\mathbf{Y}_m|\mathbf{Y}_o,X]$) and variance ($V[\mathbf{Y}_m|\mathbf{Y}_o,X]$) that all $\mathbf{Y}_m$ are smaller than zero. (Conditional Mean and Variance are easily found by applying a partition as in http://www.maths.manchester.ac.uk/~mkt/MT3732%20(MVA)/Notes/MVA_Section3.pdf)

Although I thought that this is a straight generalization of the standard censored regression model with independent observations (see for example https://cran.r-project.org/web/packages/censReg/vignettes/censReg.pdf), I have not found any paper or implementation.

Hence my questions: Do you know a solid way to solve this problem? Is there even some implementation, for example in R?

  • $\begingroup$ btw, in the univariate case, the model you are referring to is the Tobit model $\endgroup$
    – Cliff AB
    Oct 2, 2018 at 23:57

1 Answer 1


It is not so straightforward when responses are correlated, it turns out.

The method described in this paper: https://arxiv.org/pdf/1808.10541.pdf is effectively what you are looking for (except they deal with right-censoring, but the modification is straightforward, they suggest). The citation for this paper (per gung's request) is: Molstad, A. J., Hsu, L., and Sun, W. (2018). Gaussian process regression for survival time prediction with genome-wide gene expression. arXiv preprint, arXiv:1808.10541.

The authors show that you can fit such a model using a Monte-Carlo expectation conditional maximization algorithm. The complication is that in performing the $E$-step, the distribution of the censored outcomes is intractable due to the dependence across censored responses, so they use a Gibbs sampler to approximate the $Q$-function (that is, the expected value of the complete data log-likelihood). To modify their method for your setup, you would need to (i) change the hypercube over which they draw samples for the approximate $E$-step, and (ii) change their conditional $M$-step for the variance components according to whatever parameterization you use for your covariance matrix.

  • $\begingroup$ Welcome to the site. We are trying to build a permanent repository of high-quality statistical information in the form of questions & answers. Thus, we're wary of link-only answers, due to linkrot. Can you post a full citation & a summary of the information at the link, in case it goes dead? $\endgroup$ Oct 3, 2018 at 1:09

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.