Weighted normal errors regression with censoring

I have some data which I would model via standard multiple regression except:

1. There is censoring (left-censored, fixed but varying censoring points which are known)
2. The errors are assumed independent normal but of non-constant variance. Weights are available.

If it was constant variance, I would use the Tobit model and survreg() function in R. Does anyone know of the/an approach when the variance is not constant (but weights for variances are available)?

• Does the solution need to be in R? Feb 20, 2014 at 1:38
• See the R survival package survreg function. May 1, 2017 at 12:14
• You can use survival::survreg if you ask it for robust standard errors -- it allows for weights, but it doesn't interpret them as precision weights. You'll get the same point estimates as with precision weights and the robust standard errors will give consistent variance estimates. Sep 12, 2021 at 21:45

There must be some weight arguments to the survreg function? Anyhow, this can be solved by setting up a likelihood function from first principles.
You have a normal model (with independent observations) and known weight, the optimal weights are the inverse variances, so write the weight as $w_i$ taken to be the inverse of known variances. Then we can write the density as $$f(x:\mu) = \frac{\sqrt{w_i}}{\sqrt{2\pi}} e^{-\frac12 w_i (x_i-\mu)^2}$$ Assume the censoring points are at $t_i$, the first $r$ obs are fully observed and observations $r+1 \dotsc n$ censored. Then the likelihood becomes $$L(\mu) = \prod_1^r f(x_i; \mu) \prod_{r+1}^n \Phi(\sqrt{w_i}(t_i-\mu))$$ where $\Phi$ denotes the standard normal cdf. The loglikelihood becomes $$l(\mu) = -\frac12\sum_1^r w_i (x_i-\mu)^2 + \sum_{r+1}^n \log \Phi(\sqrt{w_i}(t_i-\mu))$$ where we have left out some terms not influencing the shape of the loglikelihood function. Now this function can be sent to a numerical optimization routine to find the maximum likelihood estimator of $\mu$.
• Indeed. If you are willing to be Bayesian, you of course have a great package that would do all of that for you: brms. Just needs something like brm( y | cens(censor) + weights(weight) ~ 1 + various terms, sigma ~ 1 + offset(log of known_factor affecting SD), family=gaussian()). Feb 18 at 15:01