Skip to main content

Modeling observations with a detection limit using a truncated normal distribution with censoring

Abstract

Background

When data are collected subject to a detection limit, observations below the detection limit may be considered censored. In addition, the domain of such observations may be restricted; for example, values may be required to be non-negative.

Methods

We propose a method for estimating population mean and variance from censored observations that accounts for known domain restriction. The method finds maximum likelihood estimates assuming an underlying truncated normal distribution.

Results

We show that our method, tcensReg, has lower bias, Type I error rates, and mean squared error than other methods commonly used for data with detection limits such as Tobit regression and single imputation under a range of simulation settings from mild to heavy censoring and truncation. We further demonstrate the consistency of the maximum likelihood estimators. We apply our method to analyze vision quality data collected from ophthalmology clinical trials comparing different types of intraocular lenses implanted during cataract surgery. All of the methods yield similar conclusions regarding non-inferiority, but estimates from the tcensReg method suggest that there may be greater mean differences and overall variability.

Conclusions

In the presence of detection limits, our new method tcensReg provides a way to incorporate known domain restrictions in dependent variables that substantially improves inferences.

Peer Review reports

Background

Censoring, in which the value of an observation is not known exactly but rather is only known to be above or below a specific value, is prevalent in many data settings. Censoring occurs with time-to-event data, but can also occur when measurements are subject to a detection limit (DL). A detection limit is defined as the lowest quantity or concentration of a compound that can be reliably detected with a given analytical method [1]. Quantities below the DL can be considered censored. Detection limits and the censored observations associated with them are encountered in epidemiology [2, 3], hydrology [4], chemistry [5], toxicology [6], and economics [7, 8].

Estimation of the parameters of a normal distribution based on a sample with censored observations has a long history of investigation. [9] was one of the first authors to develop maximum likelihood estimation methods for this data setting. Much of the early work focused on single mean models [912] or estimation using order statistics [13, 14]. Regression with a dependent variable subject to known censoring gained prominence with the work of [15], who developed the Tobit model, also called a censored regression or Tobit regression model. Linear regression methods with an unspecified censored distribution were later developed by [16].

In settings in which censored data arise due to a DL, estimation is sometimes performed by singly imputing the DL or 1/2 DL for censored observations. While these methods are known to yield biased estimates of the mean and standard deviation, they are still routinely applied due to convenience. Single imputation of the DL yields upwardly biased estimates of the mean [4]. There can also be substantial bias using 1/2 DL imputed values, with the direction of the bias depending on the underlying data mechanism [1, 2, 4]. [2] noted that the bias of parameter estimates when using 1/2 DL is substantial unless the proportion of censored observations is small, defined as less than 10%.

Related to but distinct from censoring is the concept of truncation. A truncated distribution is a conditional distribution that results from restricting the domain of some other probability distribution. For example, the zero-truncated Poisson distribution is the distribution of a Poisson random variable conditional on the value of the random variable being non-zero. Truncation is a strict restriction of the domain of the random variable; observations outside the domain cannot occur. In contrast, when an observation is censored, its true value is known to lie beyond the censoring threshold, and such true values are permitted to occur. Truncation is a property of the population, while censoring is a result of the sampling mechanism.

Methods have been developed for estimating the parameters of a truncated normal distribution. [9] developed methods for the single mean model, while later work by [17] provided solutions for settings involving normal distributions with known or unknown truncation points and single or double truncation. [18] extended the theoretical framework to settings in which the distribution satisfies regularity conditions but is not necessarily normal. Work by [19] expanded the truncated regression framework by demonstrating the consistency and asymptotic normality properties of the maximum likelihood estimator and identifying consistent initial estimators.

Both censored and truncated data are often referred to as limited dependent variables in the economics literature [8, 15]. There are currently methods of parameter estimation applicable to censored observations from normally distributed random variables and uncensored observations from truncated normal random variables. However, estimation for the case in which observations are censored from a truncated normal random variable, are lacking. In this manuscript we propose maximum likelihood estimation techniques for such data. The methods can be seen as analogous to those of [15] with the latent normal distribution replaced with a latent truncated normal distribution.

The paper is organized as follows. In “Methods” we introduce a motivating example from an analysis of visual quality data from clinical trials of intraocular lenses implanted during cataract surgery and the maximum likelihood estimation procedure. The “Results” compares the performance of the proposed new methodology to other common methods using both simulations and our motivating example. In the “Discussion” we reflect on the significance of the results, discusses limitations and note possible areas for extending the work. Finally the “Conclusion” sections highlights the importance of the method when modeling observations subject to a detection limit.

Methods

Problem motivation

Our application concerns left censored non-negative observations arising from contrast sensitivity testing in clinical trials for intraocular lenses implanted during cataract surgeries.

Contrast sensitivity measures the visual quality experienced by a subject by testing his or her ability to distinguish increasingly finer increments of light versus dark. Being unable to distinguish objects when contrast is low, i.e., when there is little difference between light and dark, can make everyday tasks such as night driving, navigating new settings, or perceiving distances difficult [20]. Cataract patients tend to have especially poor contrast sensitivity due to the clouding of the natural lens. The primary treatment for cataracts is surgery, during which the natural clouded lens is removed and replaced with a new synthetical intraocular lens (IOL). The patient should see improvement in visual acuity and visual quality following IOL implantation. As such, contrast sensitivity is an important clinical outcome for patients who receive IOLs during cataract surgery.

Contrast sensitivity testing is performed using standardized charts with alternating light and dark bars, referred to as gratings. To determine the contrast sensitivity of a patient, the intensity of the contrast between the bars as well as the spacing is reduced until the patient is no longer able to perceive separate bars. The contrast is defined as the relative difference in luminance of the bars from the background and may be calculated using the Weber contrast, \(\frac {L_{{max}}-L_{{min}}}{L_{{background}}}\), Michelson contrast, \(\frac {L_{{max}}-L_{{min}}}{L_{{max}}+L_{{min}}}\), or RMS contrast, \(\frac {L_{\sigma }}{L_{\mu }}\), where Lmax,Lmin,Lbackground,Lμ, and Lσ are luminance maximum, minimum, background, mean and standard deviation respectively [21]. Typically when using gratings to test contrast, the Michelson contrast is preferred. Contrast sensitivity is defined as the reciprocal of the threshold contrast, which is the lowest contrast that the patient can identify the grating. The spacing between the bars is measured in cycles per degree (CPD) with higher values of CPD indicating less space between bars. The testing is performed across a variety of different CPD levels, under either bright or dim lighting conditions, and with or without glare. In general, contrast sensitivity scores are lower when testing is performed under dim lighting with glare.

At each CPD level, the subject is presented with a sample grating followed by 8 gratings that progressively decrease the intensity of the image contrast. Figure 1 shows the testing setup used for 12.0 CPD, which is one of the most common visual quality outcomes analyzed in ophthalmic clinical trials. The subject is first asked if they can identify the striped pattern in the sample image, and is subsequently shown each column starting from column 1 and asked to identify whether the striped pattern is in the top, bottom, or neither grating. A contrast sensitivity score is recorded as the lowest level of identifiable contrast, i.e., the lowest intensity contrast for which the patient is able to correctly identify the striped pattern. The contrast sensitivity score ranges from 0 to 8, with 0 representing the ability to identify only the sample grating and 8 corresponding to the last column with the lowest image contrast. If a subject is unable to identify the sample grating, the value is recorded as -1. The scores of 0-8 are converted to continuous log contrast sensitivity values based on the manufacturers’ recommendations as shown in Table 1. Note that by definition, the contrast sensitivity values using the Michelson formula range from [1,), meaning that the log contrast sensitivity scores must be non-negative, i.e., [0,).

Fig. 1
figure1

CSV-1000E1 Contrast Sensitivity Chart for 12.0 CPD. 1This testing chart is distributed by Vector Vision and was accessed from http://www.vectorvision.com/csv1000-contrast-sensitivity/ on 29NOV2018

Table 1 Manufacturer Log Contrast Sensitivity Value 1

Problems arise using this scoring approach when patients are unable to identify the sample grating, i.e., have contrast sensitivity threshold score of -1. The true log contrast sensitivity values for these patients are only known to be below the sample threshold. This is equivalent to having left censored observations. However, we also know that log contrast sensitivity values must be non-negative, implying that the distribution of log contrast sensitivity scores is left truncated at zero. Therefore, to accurately estimate the mean and standard deviation of log contrast sensitivity for a particular IOL, both the censoring and truncation need to be accounted for.

Our data come from two prospective clinical trials for IOLs implanted during cataract surgery. The first study, ClinicalTrials.gov Identifier NCT01510717, compared a monofocal lens and a multifocal lens in a double blind randomized parallel group study with bilateral IOL implantation, and the second study, ClinicalTrials.gov Identifier NCT01424189, compared two different multifocal lenses in a nonrandomized parallel assignment multi-center study again with bilateral IOL implantation. Our analysis will be restricted to data reported for binocular (both eyes open) testing, under dim lighting with and without glare. All observations were taken 6 months after surgery.

Figure 2 shows the marginal histograms with kernel density smoothers of each IOL by glare condition. Visual inspection of the data suggest that it may be reasonable to assume an underlying normal distribution for the log contrast sensitivity scores, with scores being censored at the lower-bound detection limit.

Fig. 2
figure2

Marginal Histograms for Monofocal and Multifocal Lens at 12 CPD under Dim Lighting. Log contrast sensitivity scores are converted from contrast sensitivity threshold scores via Table 1. Detection limit for 12 CPD occurs at ν=0.61. Histogram is shown in the background with Gaussian kernel density estimate in the foreground with bandwidth set to 0.2

The goal of the analysis was to estimate the difference in mean log contrast sensitivity between the monofocal and each multifocal lenses. Historically, monofocal lenses have provided patients with better contrast sensitivity, but multifocal lenses are often preferred by patients for visual acuity as they provide both distance and near vision with increased spectacle independence. A common clinical trial hypothesis is to test whether the multifocal lens is non-inferior to the monofocal lens with respect to contrast sensitivity. Based on regulatory guidelines, losses of 0.3 log units are considered to be clinically significant when they occur at two or more spatial frequencies [22]. The non-inferiority margin for log contrast sensitivity is set at half of this clinically significant magnitude, i.e., loss of 0.15 log units. Our goal was thus to estimate pairwise differences in mean log contrast sensitivity between the monofocal lens and each multifocal lens to test for non-inferiority of visual quality.

Statistical methods

This section first develops maximum likelihood estimation for a sample from a truncated normal distribution with censored observations, i.e., a single mean model. The methods are then extended to include a linear predictor for the mean. We then propose an extension to handle heteroskedastic variances. We focus on the setting in which the truncation and censoring occur only on the lower, or left, end of the distribution, i.e., left truncation and left censoring. The results can be generalized to right truncated and right censoring. Finally, the process of finding the maximum likelihood estimates is described using different optimization techniques.

Single mean model

As a first step, we define the likelihood function for censored observations from a normal distribution. Assume a latent normal distribution with mean μ and variance σ2 for the random variable X. Following [15], assume that the values of X are left-censored at \(\nu, \nu \in \mathbb {R}\), to produce the random variable X defined as

$$ X_{i}=\left\{\begin{array}{lr} \nu & {if} \ X_{i}^{*}\le\nu,\\ X_{i}^{*} &{if} \ X_{i}^{*}>\nu. \end{array}\right. $$
(1)

Here, censored observations are reported as the DL ν. The values of X represent the observed values in the sample, which are a partial representation of the values of X. Assume that a total of n observations are independently drawn with n0 observations censored, i.e., \(n_{0}=\sum _{i=1}^{n}1_{\{x_{i}=\nu \}}\), and n1 observations uncensored, i.e., \(n_{1}=\sum _{i=1}^{n}1_{\{x_{i}>\nu \}}\). The likelihood function for such data is

$$L(\mu,\sigma)=\left[\Phi\left(\frac{\nu-\mu}{\sigma}\right)\right]^{n_{0}}\left[\frac{1}{\sigma}\right]^{n_{1}}\prod_{i\in S_{1}}\phi\left(\frac{x_{i}-\mu}{\sigma}\right), $$

where ϕ(·) and Φ(·) denote the standard normal pdf and cdf and S1 is the set of all uncensored observations. Therefore the log-likelihood is

$$\begin{aligned} l(\mu,\sigma)= n_{0}\ln\left[\Phi\left(\frac{\nu-\mu}{\sigma}\right)\right]-n_{1}\ln(\sigma)+\sum\limits_{i\in S_{1}}\ln\left[\phi\left(\frac{x_{i}-\mu}{\sigma}\right)\right]. \end{aligned} $$

Maximum likelihood estimates \(\hat {\mu }\) and \(\hat {\sigma }\) can be found using iterative optimization techniques, such as the Newton-Raphson algorithm, as discussed in [15].

Now we replace the latent normal distribution above with a latent truncated normal distribution, truncated from the left at the value a. Call this random variable Y, which has the following pdf and cdf from [17]:

$$\begin{aligned} f_{Y_{i}^{*}}\left(y_{i}^{*}\right)=&\frac{1}{1-\Phi\left(\frac{a-\mu}{\sigma}\right)}\left[\frac{1}{\sigma}\phi\left(\frac{y_{i}^{*}-\mu}{\sigma}\right)\right],\\ F_{Y_{i}^{*}}\left(y_{i}^{*}\right)=&\frac{\Phi\left(\frac{y_{i}^{*}-\mu}{\sigma}\right)-\Phi\left(\frac{a-\mu}{\sigma}\right)}{1-\Phi\left(\frac{a-\mu}{\sigma}\right)}. \end{aligned} $$

The distribution of Y is a scaled version of a normally distributed random variable, obtained by dividing the pdf by the constant \(1-\Phi \left (\frac {a-\mu }{\sigma }\right)\) to obtain a proper probability density function that integrates to one. We will denote a latent truncated normal random variable with left truncation at a constant \(a\in \mathbb {R}\) as

$$Y^{*}\sim{TN}\left(\mu, \sigma^{2}, a\right), $$

where we assume that the truncation value a is known and therefore fixed. For non-negative variables, a=0.

Note that the parameter μ denotes the mean of the underlying normal distribution prior to truncation, rather than the mean of the truncated normal distribution, which is \(\mu _{{TN}}=\mu +\frac {\phi \left (\frac {a-\mu }{\sigma }\right)}{1-\Phi \left (\frac {a-\mu }{\sigma }\right)}\sigma \) [8]. Throughout the paper we focus on the estimation of this underlying central tendency parameter μ rather than μTN.

The log-likelihood for the truncated normal distribution with n independent observations drawn from this distribution is

$$\begin{aligned} l(\mu,\sigma)=-n\ln\left[1-\Phi\left(\frac{a-\mu}{\sigma}\right)\right]-n\ln(\sigma)+\sum_{i=1}^{n}\ln\left[\phi\left(\frac{y_{i}^{*}-\mu}{\sigma}\right)\right]. \end{aligned} $$

There is no closed form for the maximum likelihood estimates for μ and σ for the truncated normal distribution. However, [17] described methods for iteratively solving for the maximum likelihood estimates and [19] showed that the estimates were consistent and asymptotically normal.

Further suppose that the values of Y are censored at \(\nu >a, \nu \in \mathbb {R}\), to produce the random variable Y, defined as

$$ Y_{i}=\left\{\begin{array}{lr} \nu &{if} \ Y_{i}^{*}\le\nu,\\ Y_{i}^{*} &{if} \ Y_{i}^{*}>\nu, \end{array}\right. $$
(2)

where ν is a known constant. For example in the contrast sensitivity problem, the testing procedure has an inherent detection limit and cannot detect values below 0.61 log units. The log contrast sensitivity scores must also be non-negative as discussed earlier in “Problem motivation” section, signifying an implicit truncation value of a=0.

The pdf for the truncated random variable with censoring can be expressed as

$$ \begin{aligned} f_{Y_{i}}(y_{i})&=1_{\{y_{i}=\nu\}}\left[\frac{\Phi\left(\frac{\nu-\mu}{\sigma}\right)-\Phi\left(\frac{a-\mu}{\sigma}\right)}{1-\Phi\left(\frac{a-\mu}{\sigma}\right)}\right]\\&+1_{\{y_{i}>\nu\}}\left[\frac{1}{\sigma\left(1-\Phi\left(\frac{a-\mu}{\sigma}\right)\right)}\phi\left(\frac{y_{i}-\mu}{\sigma}\right)\right]. \end{aligned} $$
(3)

The first term of this equation captures censored observations, in which case the data are reported as the detection limit value ν and the cdf of the truncated normal is used to provide information for the likelihood. In the second term, the observation is not censored and we simply use the pdf for the truncated normal. Figure 3 shows an example of the probability density function for such a random variable.

Fig. 3
figure3

Truncated Normal Distribution with Censoring. Potential density for a left truncated normal distribution with left censoring. The density above was created with μ=0.8,σ=0.5,a=0, and ν=0.61. In our application, a=0 and ν=0.61

Out of a total of n observations, again let n0 be censored such that \(n_{0}=\sum _{i=1}^{n}1_{\{y_{i}=\nu \}}\), and n1 be uncensored, \(n_{1}=\sum _{i}^{n}1_{\{y_{i}>\nu \}}\). Let S define the set of all observations, S1 be the set of uncensored observations, and S0 be the set of censored observations, i.e., S0S1=S. Assuming that the observations are drawn independently, the likelihood is

$$ \begin{aligned} L(\mu,\sigma)&=\left(\frac{1}{1-\Phi\left(\frac{a-\mu}{\sigma}\right)}\right)^{n}\left[\Phi\left(\frac{\nu-\mu}{\sigma}\right)-\Phi\left(\frac{a-\mu}{\sigma}\right)\right]^{n_{0}}\\&\quad\left(\frac{1}{\sigma}\right)^{n_{1}}\prod_{i\in S_{1}}\phi\left(\frac{y_{i}-\mu}{\sigma}\right). \end{aligned} $$
(4)

Taking the log of the likelihood, we have

$$ \begin{aligned} l(\mu,\sigma)=&-n\ln\left[1-\Phi\left(\frac{a-\mu}{\sigma}\right)\right] \\&\quad+n_{0}\ln\left[\Phi\left(\frac{\nu-\mu}{\sigma}\right)-\Phi\left(\frac{a-\mu}{\sigma}\right)\right]-n_{1}\ln(\sigma)\\ &+\sum\limits_{i\in S_{1}}\ln\left[\phi\left(\frac{y_{i}-\mu}{\sigma}\right)\right]. \end{aligned} $$
(5)

Similar to the censored only and truncated only log-likelihoods, the maximum likelihood estimators for μ and σ for the truncated normal distribution with censoring do not have a closed form but can be estimated using the iterative process discussed in “Obtaining maximum likelihood estimates” section.

Linear predictor for the mean

The goal of many applications is to understand how certain predictors influence the mean or to compare the means of different populations. This can be accomplished by using the linear predictor \(\mathbf {X}_{i}^{T}\boldsymbol {\beta }\), i.e.,

$$Y_{i}^{*}\sim{TN}\left(\mathbf{X}_{i}^{T}\boldsymbol{\beta}, \sigma^{2}, a\right), $$

where \(\boldsymbol {\beta }=(\beta _{1},\dots,\beta _{p-1})^{T}\) is a (p−1)×1 vector of parameters with p≥2. Again, note that \(\mathbf {X}_{i}^{T}\boldsymbol {\beta }\) is the mean of the underlying normal distribution rather than the mean of the truncated normal distribution, i.e., \(\mathrm {E}[X_{i}^{*}]\) rather than \(\mathrm {E}[Y_{i}^{*}]\). All of the unknown parameters can be collected into the vector θ=(β,σ2)T which has length p. The corresponding pdf and cdf are

$$\begin{aligned} f_{Y_{i}^{*}}\left(y_{i}^{*}\right)=&\frac{1}{1-\Phi\left(\frac{a-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)}\left[\frac{1}{\sigma}\phi\left(\frac{y_{i}^{*}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)\right]\\ F_{Y_{i}^{*}}\left(y_{i}^{*}\right)=&\frac{\Phi\left(\frac{y_{i}^{*}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)-\Phi\left(\frac{a-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)}{1-\Phi\left(\frac{a-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)}. \end{aligned} $$

Suppose that this truncated normal distribution is then censored at the value ν, where ν>a. Let the notation \(a_{i}^{*}=\frac {a-\mathbf {x}_{i}^{T}\boldsymbol {\beta }}{\sigma }\) denote a standardized version of the constant a. The likelihood of the truncated normal distribution with censoring can be expressed as

$$ \begin{aligned} l\left(\boldsymbol{\beta},\sigma^{2}\right)=&-\sum\limits_{i=1}^{n}\ln\left[1-\Phi\left(a_{i}^{*}\right)\right]\\&\quad+\sum\limits_{i\in S_{0}}\ln\left[\Phi\left(\nu_{i}^{*})-\Phi(a_{i}^{*}\right)\right] - n_{1}\ln(\sigma)\\ &+\sum\limits_{i\in S_{1}}\ln\left[\phi\left(\frac{y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}}{\sigma}\right)\right]. \end{aligned} $$
(6)

Heteroskedastic variances

We now relax the assumption of homogeneous variance. We consider the case of independent groups with different variances. Assume we have samples drawn independently from J truncated normals, with each population having a common truncation value but possibly different variance, according to the model

$${}Y^{*}_{{ij}}\sim{TN}\left(\mathbf{X}_{{ij}}^{T}\boldsymbol{\beta}, \sigma^{2}_{j}, a\right), \quad i = 1,\dots,n_{j} \ {and} \ j = 1,\dots,J, $$

where nj is the number of observations in group j and \(Y^{*}_{{ij}} \perp\perp Y^{*}_{i'j'}\) for all ii and jj. Assume that observations are censored at the value ν to create a sample of independent random variables with pdf

$$\begin{aligned} f_{Y_{{ij}}}(y_{{ij}})=&1_{\{y_{{ij}}=\nu\}}\left[\frac{\Phi\left(\frac{\nu - \mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)-\Phi\left(\frac{a - \mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)}{1-\Phi\left(\frac{a - \mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)}\right]\\ &+1_{\{y_{{ij}}>\nu\}}\left[\frac{1}{\sigma_{j}\left(1-\Phi\left(\frac{a - \mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)\right)}\phi\left(\frac{y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)\right]. \end{aligned} $$

Because groups are independent, the log likelihood becomes

$$ \begin{aligned} l\left(\boldsymbol{\beta},\sigma^{2}_{1},\dots,\sigma^{2}_{J}\right)=&\sum\limits_{j=1}^{J}\sum\limits_{i=1}^{n_{j}}-\ln\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]\\&+\sum\limits_{i\in S_{0j}}\ln\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right] - n_{1j}\ln(\sigma_{j})\\ &+\sum\limits_{i\in S_{1j}}\ln\left[\phi\left(\frac{y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}}{\sigma_{j}}\right)\right], \end{aligned} $$
(7)

where \(a_{{ij}}^{*}=\frac {a-\mathbf {x}_{{ij}}^{T}\boldsymbol {\beta }}{\sigma _{j}}, S_{0j}\) and S1j are the sets of censored observations and uncensored observations respectively in the jth group, and n1j is the number of uncensored observations in the jth group.

Obtaining maximum likelihood estimates

Our goal is to find the values \(\hat {\boldsymbol {\theta }}\) that maximize the log-likelihoods of Eqs. 5, 6, and 7, where θ corresponds to the appropriate mean and standard deviation parameters, i.e., θ=(β,σ). However, closed form solutions do not exist. To obtain maximum likelihood estimates, an iterative procedure is required. One approach is to use the Newton-Raphson algorithm using Taylor series expansion, as discussed in Chapter 14 in [23]. While the Newton-Raphson method has attractive local convergence guarantees and reliable performance [24], each step requires evaluation of the Hessian matrix, which can become computationally expensive for a large set of predictors. Alternative optimization routines, such as the quasi-Newton BFGS [2528] or the conjugate gradient [29], which require only evaluation of the likelihood function and corresponding gradient, often require additional evaluations but have reduced memory and faster computing time.

Within R, other optimization packages such as the maxLik package from [30] can also be used to find maximum likelihood estimates. This package is called in the censored only and truncated only maximum likelihood estimation packages in R such as censReg by [31] and truncreg by [32].

We developed the standalone R package tcensReg available in CRAN to solve the novel likelihood equation of the truncated normal distribution with censoring. This software package uses analytic results of the gradient and Hessian via Newton-Raphson optimization for the corresponding model of interest in either Eq. 6 or Eq. 7 derived in Appendix A. Several other optimization routines are available within the software including conjugate gradient, maxLik, and BFGS. This package uses familiar model syntax and has additional functionality to estimate parameters for the censored only or truncated only settings similar to the censReg and truncreg packages, respectively.

Results

Simulation study

We conducted a simulation study to compare the performance of our method to that of five methods of estimating the mean and standard deviation from a truncated normal distribution with censored observations. Method 1 is the gold standard method which uses the true uncensored observations, and accounts for the truncation in the estimation procedure by using the appropriate truncated log-likelihood. Method 2 uses the same uncensored truncated data but does not adjust for truncation in the normal distribution likelihood function. Methods 3-6 use the censored truncated observations but differ in how they treat censoring and truncation. Method 3 imputes all censored values with the detection limit and uses maximum likelihood estimation with a normal distribution likelihood, while Method 4 imputes all values as half of the detection limit and also uses normal maximum likelihood estimation. Both the DL and 1/2 DL imputation methods have been shown to perform poorly in cases with censoring [1, 2, 4], and are expected to show even worse performance in this setting since they do not account for censoring or truncation. However, they are still sometimes used in practice. Method 5 uses Tobit regression, which incorporates censoring into the likelihood. This method is often recommended when the assumption of normality seems reasonable. Finally, Method 6 is our proposed maximum likelihood estimation procedure described in the “Methods” section which takes into account not only the censoring but also the truncation of the underlying distribution. These six methods will be referred to as: 1) Gold Standard (GS), 2) Uncensored with no truncation adjustment (Uncens NT), 3) Detection Limit (DL), 4) 1/2 Detection Limit (1/2 DL), 5) Tobit regression (Tobit), and 6) censored regression with truncation adjustment (tcensReg), our proposed method.

Each of these methods makes different assumptions about the underlying data generating process. Both the GS and Uncens NT methods observe \(Y_{i}^{*}\). The GS method uses the true generating distribution, i.e., \(Y_{i}^{*}\sim {TN}\left (\mu, \sigma ^{2}, a\right)\). With Uncens NT, the distribution is incorrectly assumed to be not truncated, i.e., \(Y_{i}^{*}\sim \mathrm {N}\left (\mu, \sigma ^{2}\right)\). The remaining methods observe only the censored data, Yi, along with the detection limit ν. For DL, 1/2 DL, and Tobit methods, the observed data are assumed to be generated from a normal distribution YiN(μ,σ2), but handle the censored observations differently. In DL and 1/2 DL, point masses are placed at either ν or ν/2 before estimating the parameters μ and σ2. The Tobit method instead assumes that the censored observations fall within the tail region below the detection limit, and an appropriate term is incorporated into the likelihood function. The tcensReg method correctly assumes that the censored observations are generated from an underlying truncated normal distribution, YiTN(μ,σ2,a), with information in the censored observations incorporated using the truncated normal tail probabilities. In simulation studies, we can test each of these methods against the gold standard, since the complete data generating process is known. However, in practice the true underlying distribution is unknown, meaning only the DL, 1/2 DL, Tobit, and tcensReg methods are appropriate to implement.

Set-up

Four different sets of simulation studies were conducted to compare the six methods. In the first simulation study, values from a single mean model were drawn to compare performance of the methods in terms of bias and mean squared error (MSE) for estimating the mean and standard deviation of the underlying latent distribution. The second simulation study focused on estimating the difference of the means of two independent populations and their common variance. In the third simulation study, the performance of the methods in a non-inferiority test setting similar to that of our motivating example was assessed. Finally, the last simulation study examined the consistency of the maximum likelihood estimation procedures as a function of sample size in a single mean model. The simulations were conducted in R version 3.6.2[33].

For the first simulation study, values were simulated from a truncated normal distribution,

$$\begin{array}{*{20}l} Y_{i}^{*}\sim{TN}\left(\mu, \sigma^{2}, a\right), \end{array} $$

and then censored. A constant value of a=0 was used to represent zero-truncation. The values of μ and σ2 were chosen to approximate the marginal distributions of the log contrast sensitivity scores from the application introduced in “Problem motivation” section. Typical values of mean log contrast sensitivity ranged from 0.7 to 1.1, and thus we used μ{0.7,0.8,0.9,1.0,1.1}. Based on standard deviations observed in the contrast sensitivity data, for the simulation we used σ{0.40,0.45,0.50}. This created a total of 15 (5×3) parameter combinations for μ and σ2.

Observations from a truncated normal distribution, \(Y_{i}^{*}\sim {TN}\left (\mu,\sigma ^{2}, a\right)\), can be simulated by transforming samples from a uniform distribution on the interval [0,1] using the inverse probability transformation:

$${}Y_{i}^{*}\,=\,\Phi^{-1}\left\{p\times\left[1-\Phi\left(\frac{a-\mu}{\sigma}\right)\right] \,+\, \Phi\left(\frac{a-\mu}{\sigma}\right)\right\}\times\sigma+\mu, $$

where p represents the sample from the uniform distribution [34]. This method of inverse transformation sampling is implemented with the tcensReg software package. After transforming to the appropriate truncated normal distribution, another dataset was generated by censoring the observations at ν=0.61. Values that fell below ν were either replaced with the DL (0.61), 1/2 DL (0.305), or marked as censored for Tobit and tcensReg estimation.

Each of the six estimation methods was used to estimate μ and σ. The data simulation and estimation procedures were repeated for B=10,000 replications with a sample size of n=100. For each of the six methods, two performance metrics were calculated for θ{μ,σ}: average bias (\(\bar {\hat {\theta }}-\theta \)) and mean squared error (MSE; \(1/B\sum _{k=1}^{B}[\hat {\theta }_{k}-\theta ]^{2}\)) where \(\bar {\hat {\theta }}=1/B\sum _{k=1}^{B}\hat {\theta }_{k}\). When reporting results for simulations throughout this manuscript, log MSE rather than MSE was used as the MSE for the various procedures sometimes differed by orders of magnitude. Applying this monotonic transformation allowed for easier comparison among methods.

In the second simulation study, we simulated values from two truncated normal populations with different means but common variance and truncation value, i.e.,

$$\begin{array}{*{20}l} Y_{1}^{*}&\sim{TN}\left(\mu_{1}, \sigma^{2}, a\right) & Y_{2}^{*}&\sim{TN}\left(\mu_{2},\sigma^{2}, a\right). \end{array} $$

A constant value of a=0 was used to represent zero-truncation and a range of values of μ1,μ2, and σ2 were selected to produce data similar to the application. Population 1 approximated a monofocal intraocular lens, while Population 2 approximated a multifocal intraocular lens. Means for monofocal lenses took values μ1{1.0,1.1}, and the difference between multifocal and monofocal lenses was set to range from no difference to a clinically significant difference of 0.3 log units [22], i.e., δ{−0.3,−0.2,−0.1,0}. Therefore the mean of Population 2 was set as μ2=μ1+δ. The common standard deviation was assumed to take values σ{0.40,0.45,0.50}. The two-population simulation had a total of 24 (2×4×3) parameter combinations for μ1,μ2, and σ2.

A total of 100 observations from each population were sampled. Again, for each population a separate censored dataset was generated using the censoring threshold ν=0.61, and values that fell below ν were replaced with the DL (0.61), 1/2 DL (0.305) or marked as censored for Tobit and tcensReg estimation. Each of the six estimation methods was used to estimate the mean difference between Population 1 and 2, δ, and the common standard deviation, σ. Data simulation and estimation were repeated for B=10,000 replications. For each of the six methods, average bias and MSE were calculated.

The third simulation study assessed the performance of the methods in the context of a non-inferiority test. Here we focused on the Type I error rates of the various methods. Type I error rates are particularly important for non-inferiority tests because, if non-inferiority is falsely accepted, patients may decide between products based on non-efficacy factors such as price and side effects assuming the products to be similar when in fact one is truly superior. For the non-inferiority test, data were simulated with a true difference of δ=−0.15 in the two population model, while varying μ1{1.0,1.1} and σ{0.40,0.45,0.50} for a total of 6 different non-inferiority scenarios. Then each of the six methods was used to construct 1−α confidence intervals for δ. The test is specified as one-sided because it is known that multifocal contrast sensitivity is less than monofocal contrast sensitivity. Constructing 1−α confidence intervals and comparing the lower bound to the non-inferiority margin will result in a one-sided hypothesis test at the α level. If the lower bound of the confidence interval does not cover the true value of δ, then the non-inferiority hypothesis would be falsely accepted. A total of 100 observations were drawn from each population to construct the confidence intervals. The hypothesis test was repeated for B=10,000 replications and the Type I error rate was calculated as the percent of replications where the lower bound of the 1−α confidence interval was greater than −0.15.

The final simulation study compared the performance of the methods as a function of the number of observations (n). The procedure from the first simulation was repeated fixing μ=0.9 and σ=0.45, but sample size varied from n{100,200,400,800,1600}. Average bias and MSE were calculated for estimation of μ and σ. To compare the rate of Type I error as a function of sample size, the procedure for the third simulation was repeated with μ1=1.0 and σ=0.45 while allowing the sample size per population to vary as above, i.e., n=n1=n2{100,200,400,800,1600}.

The choice of initial starting values θ(0) for the tcensReg method is important since optimization algorithms can provide local rather than global convergence. To ensure that the starting values are reasonable, we recommend using initial estimates from a censored regression, i.e., Tobit, model. These estimates showed excellent rates of convergence for our simulation settings.

Simulation results

Table 2 shows expected censoring and truncation percentages and the ratio of censoring to truncation for each of the 15 parameter value scenarios in the single mean model. The expected percent of censoring ranges from 10.8%-38.6% and truncation was typically ≤5% but was slightly higher when the mean was closer to the truncation value of 0, i.e., μ{0.7,0.8}. The ratio of censoring to truncation varies from 4.68 to 35.87 where values greater than 1 indicate more censoring than truncation. In general, for a fixed μ, as σ increases the expected censoring and truncation increase while the ratio of censoring to truncation decreases.

Table 2 Expected Percentages of Truncation and Censoring in Single Mean Simulation Study

Performance metrics of the six methods for μ are shown in Fig. 4. The left panel shows results with respect to average bias. The figure shows that the gold standard was essentially unbiased with only slight negative bias for low values of μ. The average bias of tcensReg was slightly negative for all scenarios, meaning that the estimated mean was smaller than the true value. The slight negative bias increased as the amount of censoring increased, i.e., as μ decreased for a fixed σ and as σ increased. However, the absolute bias remained small and was always below 2% of the true μ value. The Tobit, Uncens NT and DL methods consistently had positive bias, with the estimated mean consistently larger than the true value. The amount of bias also increased more rapidly for these other methods as censoring increased, compared to tcensReg. The 1/2 DL method had variable trends in bias with positive bias in some settings and negative bias in others. In particular, when σ=0.40, the 1/2 DL method overestimated μ when μ≤0.8 and underestimated it when μ>0.8.

Fig. 4
figure4

Performance Metrics for μ from Six Different Estimation Methods in Single Mean Model. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

The tcensReg method was consistently closer to the gold standard compared to the Tobit method, especially in situations with high censoring, i.e., σ=0.50. There were high censoring settings in which the tcensReg method and Tobit method performed similarly; for example, absolute bias was similar when μ=0.9 and σ=0.40, which had 22.47% expected censoring. For this scenario, the ratio of censoring to truncation is 18.42 meaning that there is more than 18 times more expected censoring than truncation. For cases where the ratio of censoring to truncation is low, the tcensReg method tends to outperform the Tobit method. When the ratio of censoring to truncation is greater, the Tobit method is comparable in terms of absolute bias to the tcensReg method, but with different directions of bias.

The precision of the estimates of the mean as reflected by their log MSE is shown in the right panel of Fig. 4. MSE was transformed to the log scale as there was a great disparity between the values for the DL point imputation method and the other five methods. Here, the average log MSE of each method should be compared to that of the gold standard. The DL method consistently had the highest log MSE of all of the methods. The tcensReg method had higher log MSE than the Tobit, 1/2 DL, and Uncens NT methods, suggesting that these methods were outperforming tcensReg. However, the average log MSE for these three methods often fell below that of the gold standard, which reflects an underestimation of the true sampling variance. In contrast, the tcensReg method always maintained a log MSE greater than or about equal to the gold standard. Overall, the tcensReg method avoided false precision.

The left panel of Fig. 5 provides evidence that other methods such as Tobit, 1/2 DL, and Uncens NT also systematically underestimated the variance of the underlying distribution. The Tobit and Uncens NT always have negative bias, meaning that these methods underestimated σ, with this bias being especially pronounced when σ is large and μ is small. Similar to the bias for μ, the bias for the 1/2 DL method when estimating σ is variable with positive bias in low censoring scenarios and negative bias in higher censoring scenarios. Overall the average bias for the tcensReg method was generally closest to the gold standard with only slight overestimates of σ. However, we note that again the ratio of censoring to truncation appears to play an important role when comparing the Tobit and tcensReg methods. The Tobit method slightly outperforms the tcensReg method for scenarios with σ=0.4. In these settings the censoring to truncation ratio is higher as shown in Table 2, meaning that censoring is expected to occur much more frequently than truncation.

Fig. 5
figure5

Performance Metrics for σ from Six Different Estimation Methods in Single Mean Model. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

On the right panel of Fig. 5 we see results for the log MSE for each method. Again, the tcensReg method protected against false precision, with log MSE slightly greater than that of the gold standard. The 1/2 DL method had artificially low MSE for almost all parameter scenarios, while the Uncens NT alternated between over- and underestimation. The Tobit method was often the closest method to the gold standard.

In the two population simulation study, the parameter of primary interest was the difference of the means, δ. Figure 6 shows the performance of each method with respect to average bias for δ. For all values of δ, the tcensReg and gold standard methods were essentially unbiased. The other four methods all had positive bias, corresponding to underestimation of the true difference in means. Similar to the single mean model, the greatest difference between the tcensReg method and the other methods occurred when the amount of censoring was greater, i.e., when μ1=1.0 and |δ| is large for a fixed σ.

Fig. 6
figure6

Average Bias for δ from Six Different Estimation Methods in Two Population Model. The vertical dashed black line corresponds to the case when δ=0, i.e., μ1=μ2. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

For all six methods, the average bias of \(\hat {\delta }\) when δ=0 was approximately 0. As our simulation study for the single mean model showed, bias generally increases as the mean decreases for fixed σ; see Fig. 4. When the two means were equal, i.e., δ=0, the biases for estimating each mean were also equal and thus there was no bias for estimating δ.

Additional performance results for log MSE of δ and common standard deviation σ in the two-population simulation are available in Appendix B, Figs. 11 and 12 respectively. These results were similar to those from the single mean model with tcensReg avoiding errors of false precision and having optimal average bias compared to the gold standard.

Results for testing the non-inferiority hypothesis with α=0.05 are shown in Table 3. For all scenarios, the tcensReg method had slightly higher Type I error rates than the gold standard and the Tobit method had slightly higher Type I error rates than the tcensReg method. The DL method had Type I error rates 2.5-4 times higher than the nominal value. The other single imputation method, 1/2 DL, had inflated Type I errors in the medium and high variance settings. Especially when censoring was higher, i.e., μ1=1.0 and μ2=0.85, the other methods had significantly higher Type I error rates compared to tcensReg.

Table 3 Type I Error Rates for Non-inferiority Test in Simulation Study

Performance results for the numerical consistency of the maximum likelihood estimates with respect to μ in a single mean model are displayed in Fig. 7. In the left panel average bias for the tcensReg and gold standard methods decreased as sample size increased, and is near zero for sample sizes above 200. All other methods show no noticeable differences in average bias as the sample size increased. As shown in the right panel of Fig. 7, the average log MSE remains approximately constant for the tcensReg method as sample size increases, while Tobit, Uncens NT, and DL begin to perform noticeably worse. The 1/2 DL method still shows evidence of false precision as the average log MSE is below the gold standard for all values of n. Similar results for the estimates of σ, shown in Appendix B Fig. 13, demonstrate that the bias of the tcensReg method decreases as a function of sample size with consistent MSE performance. Finally, Table 4 shows that Type I error rates for GS and tcensReg generally decrease as a function of sample size. Other methods have Type I error rates that escalate with sample size, with Uncens NT, DL, 1/2 DL, and Tobit all having more than 2.9 times the nominal value with a sample size of n=1600.

Fig. 7
figure7

Performance of Maximum Likelihood Estimate for μ as Function of Sample Size. The vertical dashed black line on the left figure corresponds to zero bias. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

Table 4 Type 1 Error Rates as Sample Size Increases for Simulation Study

Application results

We now apply the methods to our contrast sensitivity application introduced in “Problem motivation” section. The goal was to compare visual quality, measured as contrast sensitivity, for monofocal vs multifocal lenses implanted following cataract surgery, using data collected from two clinical trials. Table 5 shows the number of participants who received each type of IOL and the percent of patients whose observation was censored under each glare condition. Each group has over 150 patients; Multifocal 2 is larger than the others due to a 2:1 allocation ratio in that trial. As mentioned previously, the amount of censoring is greater with glare than without glare, with an average difference of 9% censoring. However, all groups have at least moderate levels of censoring, ranging from 11%-29%.

Table 5 Observed Censoring Percentage for Intraocular Lenses (IOLs) under Dim Lighting at 12.0 CPD

Separate models were fit for each pairwise comparison between the monofocal lens and each multifocal lens with and without glare. Initially, each pairwise comparison was fit assuming the standard deviation for lens type was heteroskedastic. Figure 8 displays result for the 95% confidence interval of each groups standard deviation, \(\hat {\sigma }_{j}\). The standard deviations for the multifocal lenses are sufficiently close to the monofocal lens for both glare conditions with a high level of overlap in the confidence intervals. This suggested that a common standard deviation for lens type is appropriate. In the subsequent analysis, each model was re-run assuming a common standard deviation in the two groups. The parameters of interest included the mean of each population (monofocal or multifocal), the mean difference between populations, and the common standard deviation. Unlike the simulation design, we observe only the censored observations meaning we are restricted to four methods to estimate the parameters. DL imputation for censored observations, 1/2 DL imputation, Tobit regression and our censored truncated method, tcensReg.

Fig. 8
figure8

95% Confidence Intervals for Separate Standard Deviation for Monofocal vs Multifocal Lens at 12 CPD

We also tested whether the multifocal lens is non-inferior to the monofocal lenses with respect to contrast sensitivity. Based on regulatory guidelines, the non-inferiority margin for contrast sensitivity is set at -0.15 log units [22]. The non-inferiority test was one-sided test with α=0.05, conducted by constructing a two-sided 90% confidence interval (CI) for δ and comparing the lower bound of the CI to the non-inferiority margin. To establish non-inferiority, the lower bound for the CI must be above the non-inferiority margin.

Figure 9 shows 90% CIs for δ using the four methods. For all of the methods, the lower bound of the CIs are below the non-inferiority margin, for all three lens comparisons and both glare conditions. Thus using any of these methods, we would be unable to conclude that a particular multifocal lens is non-inferior to the monofocal lens in terms of contrast sensitivity. For each comparison, the tcensReg method had the largest estimate of |δ| and the DL method had the smallest estimate. The point estimates of δ from the Tobit and 1/2 DL methods were intermediate and similar.

Fig. 9
figure9

90% Confidence Intervals for Difference in Monofocal vs Multifocal Lens at 12 CPD. The horizontal dashed line at δ=−0.15 indicates the non-inferiority margin. DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

The tcensReg method yielded the longest CIs (mean length = 0.164 log units). Mean CI lengths for the Tobit, 1/2 DL and DL methods were 0.158, 0.156 and 0.133 log units, respectively. The differences in CI length are a reflection of differences in estimates of σ, shown in Fig. 10. The DL method had consistently low estimates of σ, while the tcensReg method had the highest estimated values of σ. The 1/2 DL and Tobit methods gave similar estimates intermediate between the other two methods.

Fig. 10
figure10

Estimate of Common Standard Deviation in Monofocal vs Multifocal Lens at 12 CPD. DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

Interpreting these results in light of our simulation studies, if the latent normality assumption holds for these data, we would expect the estimates from the tcensReg method to have the least bias. The other methods typically underestimated both the difference in means and the population standard deviation, which would tend to lead to confidence intervals that are more likely to falsely exclude the non-inferiority margin and thus to higher Type I error rates, especially when there is a high rate of censoring. The wider confidence intervals observed in Fig. 9 from the tcensReg may better protect against Type I errors.

As an additional way to compare the proposed tcensReg model against potential competitors, we compared goodness of fit. As suggested by [35], we focused on measures of fit based on the Pseudo R2. For limited dependent variables, such as those arising from a censored truncated normal distribution, the authors found that the Pseudo R2 metric from [36] performed best for model comparison. Table 6 shows the results for this Pseudo R2 metric for each pairwise comparison with glare. The tcensReg method consistently had the highest value of the Pseudo R2, providing evidence that this model was the most appropriate for the data.

Table 6 Goodness of Fit Metrics for Intraocular Lenses (IOLs) under Dim Lighting at 12.0 CPD with Glare

Discussion

In this manuscript, we developed a maximum likelihood method for estimating parameters from a truncated normal distribution when observations are subject to censoring. We have also developed the R package tcensReg to implement the method. We showed in simulations that our method has substantially less bias for estimating the mean and standard deviation than other commonly used methods for a range of simulation settings. Our method also had close to the nominal Type I error rate for non-inferiority testing, while single imputation with either the detection limit or half the detection limit and Tobit regression often had inflated Type I error rates when the censoring rate was high.

The single mean simulation study showed that as the levels of censoring and truncation increased, the bias for all other methods generally increased, often dramatically. As expected, point imputation of the detection limit was consistently the worst performing method, but even the Tobit method had large bias under certain conditions. An important factor when comparing the Tobit and tcensReg method in the single mean model was not only the raw levels of censoring and truncation but the ratio of censoring to truncation. For low censoring to truncation ratios, generally below 20, the tcensReg method outperformed Tobit, while high ratios tended to have similar performance in mean estimation and worse performance in standard deviation estimation. This trend may point to the scenarios where the censoring dominates truncation in terms of estimation and thus using the Tobit method can lead to more precise estimates, particularly of the standard deviation. It is important to note that the tcensReg method was the only method to have positive bias in estimation of the standard deviation while other methods were typically negatively biased. Positive bias in this sense leads to conservative hypothesis testing.

In the two population simulation study, we observed trends similar to those for the single mean model. The tcensReg method uniformly outperformed all other methods in estimating the average difference and standard deviation between the populations with respect to average bias. The greatest difference between the methods occurred when the censoring rate was high and the difference between the population means was large. Unlike the single mean model, the ratio of censoring to truncation did not appear to play a significant role in the accuracy of parameter estimation for the two-population method. Similar to the single mean model, the precision of alternative methods was also inaccurate, leading to values of the mean squared error that appeared to outperform the gold standard. Across all parameter scenarios, the tcensReg method had significantly lower average bias and did not fall victim to the false precision fallacy.

Results from the non-inferiority simulation study confirmed the trends observed in the previous simulation studies. Methods other than tcensReg tended to underestimate the true variability of the underlying data generating mechanism, which led to precise but biased results. This bias tended to increase as the censoring increased, and this combined with underestimated values of the true standard deviation led to substantially inflated Type I error rates. The tcensReg method consistently had estimates closest to the gold standard.

In the simulation study, the maximum likelihood estimates were shown to be numerically consistent for the tcensReg method. The average bias for population mean and standard deviation tended towards 0 as the sample size grew, and Type I error rates approached the nominal level.

In the application, the main difference among the methods was a greater estimate of the difference between means when using the tcensReg method. The tcensReg method also had the widest confidence interval length due to higher estimates of the variance. The DL method results suggested that it would be the most likely method to lead to a finding of non-inferiority, but as shown in the simulation study, this method can have highly inflated Type I error rates. The 1/2 DL, Tobit, and tcensReg methods are more likely to capture the true difference and variability in the estimate, with the tcensReg method showing evidence that the true difference may be even greater than previously thought using either of the other two methods. In our particular application, none of these methods differed with respect to the ultimate conclusion of the non-inferiority hypothesis test. However, in cases in which non-inferiority is more marginal, the choice of method could make a difference. Goodness of fit metrics based on Pseudo R2 provided additional evidence that the tcensReg method was appropriate to model the observed contrast sensitivity scores. These goodness of fit metrics can be used in conjunction with subject matter knowledge about the data generating process when selecting a model to be used.

The goal of the estimation methods presented here is to estimate the mean and standard deviation of the latent normal distribution rather than the mean and standard deviation of the truncated normal distribution. This affects the interpretation of the results. When censoring is less than 50 percent, the mean of the latent normal corresponds to the mode of the truncated normal. In many settings the mode of the truncated normal may be a more clinically meaningful measure of central tendency than the mean of the truncated normal. For example, in contrast sensitivity testing, it may be clinically relevant to understand what factors affect the mode rather than the mean in order to target where the greatest proportion of patients lie.

The methodology as presented makes a strong parametric assumption of normality. For studies with small to moderate sample sizes, checking the reasonableness of the normality assumption may be difficult and the assumption may be increasingly tenuous as the amount of censoring and truncation increase. While the unobserved data may not be exactly normal, a normality assumption at least reflects an assumption that the unobserved part of the distribution has a monotone decreasing shape, which is frequently a reasonable assumption. In general, whenever an analysis involves unverifiable assumptions, conducting a sensitivity analysis is prudent. The tcensReg method can be viewed as an option in the toolkit of the statistician that can be used as part of a sensitivity analysis, which might also include Tobit regression.

Another potential limitation is that the detection threshold and truncation value are assumed to be known. While assuming that the detection threshold is known is often reasonable, the truncation value is not observable in the data due to the censoring and so must be incorporated in the analysis based on subject matter expertise.

In this manuscript, we have focused on the setting of left censoring and left truncation. The methods can be easily extended to handle right censoring with right truncation. Future work could extend the methods to handle censoring in both tails of the distribution. In our application, we noted some suggestion of censoring due to an upper detection threshold; see Fig. 2. Conducting the analysis adjusting for left and right censoring could potentially improve the accuracy of the estimates.

Other possible extensions include to other non-normal parametric distributions, or to linear mixed effect models with repeated measurements using the extended Newton-Raphson technique proposed by [37]. An additional area for future work is criteria for model and variable selection. Work by [38] introduced the concept of using Focused Information Criteria when performing variable selection for Tobit models. Based on our work, we suggest using Pseudo R2 from [36] when comparing models and AIC or BIC for variable selection in tcensReg models.

Conclusion

In the presence of detection limits, the values of observations below the detection limit are censored. In addition, in many settings, the response domain is known to be restricted, often with zero truncation. We demonstrated a new estimation method that accounts for this restriction and thereby substantially improves inferences.

Appendix A: Analytic derivations

Formula for the gradient

In the model discussed in “Linear predictor for the mean” section, the parameters to be estimated are θ=(β,σ)T where \(\boldsymbol {\beta }=(\beta _{1},\dots,\beta _{p-1})^{T}\) from Eq. 6. The single mean model is a special case of the linear equation where \(\mathbf {X}_{i}^{T}\boldsymbol {\beta }=\mu \).

Investigating the log likelihood from Eq. 6, it is evident that calculating the partial derivatives with respect to ln(σ) rather than σ is optimal. Then to obtain estimates and standard errors for σ, one can apply the inverse transformation and delta method.

Note that another possible parametrization using one-to-one functions is \(\delta =\frac {\beta }{\sigma }\) and \(\eta =\frac {1}{\sigma }\) used by [15]. However, we present results for the derivation using β and ln(σ) analogous to methods used by [31] in the R package censReg.

Throughout this derivation, let \(c_{i}^{*}=\frac {c-\mathbf {x}_{i}^{T}\boldsymbol {\beta }}{\sigma }\) for some constant c. The following properties will be used to calculate the appropriate derivatives:

  • \(\frac {\partial }{\partial x}\Phi (x)=\phi (x)\)

  • \(\frac {\partial }{\partial x}\phi (x)= -x\phi (x)\)

  • \(\frac {\partial }{\partial \ln (\sigma)}\sigma =\frac {\partial }{\partial \ln (\sigma)}\exp \left [\ln (\sigma)\right ]=\sigma \)

  • \(\frac {\partial }{\partial \ln (\sigma)}c_{i}^{*}=-c_{i}^{*}\)

We define the gradient as

$$\nabla l(\boldsymbol{\theta})=\left(\begin{array}{cccc} \frac{\partial l}{\partial \beta_{1}}\\ \vdots \\ \frac{\partial l}{\partial \beta_{p-1}} \\ \frac{\partial l}{\partial \ln(\sigma)} \end{array}\right), \ {for} \ p\ge2. $$

Taking the derivative of the log-likelihood with respect to βk,

$$ {}\begin{aligned} \frac{\partial l}{\partial \beta_{k}}=&-\sum\limits_{i=1}^{n}\frac{x_{{ik}}\phi\left(a_{i}^{*}\right)} {\sigma\left[1-\Phi\left(a_{i}^{*}\right)\right]} -\sum\limits_{i\in S_{0}}\frac{x_{{ik}}\left[\phi\left(\nu_{i}^{*}\right)-\phi\left(a_{i}^{*}\right)\right]} {\sigma\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]}\\ &+\frac{1}{\sigma^{2}}\sum\limits_{i\in S_{1}}x_{{ik}}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right) \ \ {for} \ k = 1,\dots, p-1. \end{aligned} $$
(8)

Then taking the derivative with respect to ln(σ),

$$ {}\begin{aligned} \frac{\partial l}{\partial \ln(\sigma)}&=-\sum\limits_{i=1}^{n}\frac{a_{i}^{*}\phi\left(a_{i}^{*}\right)}{1-\Phi\left(a_{i}^{*}\right)} -\sum\limits_{i\in S_{0}}\frac{\left[\nu_{i}^{*}\phi\left(\nu_{i}^{*}\right)-a_{i}^{*}\phi\left(a_{i}^{*}\right)\right]}{\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)} \\&\quad-n_{1} +\frac{1}{\sigma^{2}}\sum\limits_{i\in S_{1}}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right)^{2} \end{aligned} $$
(9)

Therefore, the gradient vector is

$$ \begin{aligned} \nabla l(\boldsymbol{\theta})=\left(\begin{array}{cccc} -\sum_{i=1}^{n}\frac{x_{i1}\phi(a_{i}^{*})}{\sigma\left[1-\Phi\left(a_{i}^{*}\right)\right]} -\sum_{i\in S_{0}}\frac{x_{i1}\left[\phi\left(\nu_{i}^{*}\right)-\phi\left(a_{i}^{*}\right)\right]}{\sigma\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]} \\ +\frac{1}{\sigma^{2}}\sum_{i\in S_{1}}x_{i1}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right) \\ \vdots \\ -\sum_{i=1}^{n}\frac{x_{i(p-1)}\phi\left(a_{i}^{*}\right)}{\sigma\left[1-\Phi\left(a_{i}^{*}\right)\right]} -\sum_{i\in S_{0}}\frac{x_{i(p-1)}\left[\phi\left(\nu_{i}^{*}\right)-\phi\left(a_{i}^{*}\right)\right]}{\sigma\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]} \\ +\frac{1}{\sigma^{2}}\sum_{i\in S_{1}}x_{i(p-1)}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right)\\ -\sum_{i=1}^{n}\frac{a_{i}^{*}\phi\left(a_{i}^{*}\right)}{1-\Phi\left(a_{i}^{*}\right)} -\sum_{i\in S_{0}}\frac{\left[\nu_{i}^{*}\phi\left(\nu_{i}^{*}\right)-a_{i}^{*}\phi\left(a_{i}^{*}\right)\right]}{\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)} -n_{1} \\ +\frac{1}{\sigma^{2}}\sum_{i\in S_{1}}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right)^{2} \end{array}\right). \end{aligned} $$
(10)

Formula for the hessian

The Hessian is the matrix of second derivatives,

$$\nabla^{2}l(\boldsymbol{\theta})=\left[\begin{array}{cccc} \frac{\partial^{2} l}{\partial \beta_{1}^{2}} & \cdots & \frac{\partial^{2} l}{\partial \beta_{1} \partial\beta_{p-1}} & \frac{\partial^{2} l}{\partial \beta_{1} \partial\ln(\sigma)}\\ \vdots & \ddots & \vdots & \vdots \\ \frac{\partial^{2} l}{\partial\beta_{p-1} \partial \beta_{1}}& \cdots & \frac{\partial^{2} l}{\partial \beta_{p-1}^{2}} & \frac{\partial^{2} l}{\partial \beta_{p-1} \partial \ln(\sigma)}\\ \frac{\partial^{2} l}{\partial\ln(\sigma) \partial \beta_{1}} & \cdots & \frac{\partial^{2} l}{\partial \ln(\sigma) \partial \beta_{p-1}} & \frac{\partial^{2} l}{\partial \ln^{2}(\sigma)} \end{array}\right]. $$

Note that this Hessian matrix is symmetric so that 2l(θ)ij=2l(θ)ji for ij.

The individual components of this matrix are calculated as

$$ \begin{aligned} \frac{\partial^{2} l}{\partial \beta_{k}\partial\beta_{l}}=&-\sum\limits_{i=1}^{n}\frac{x_{{ik}}x_{{il}}\left\{a_{i}^{*}\left[1-\Phi\left(a_{i}^{*}\right)\right]\phi\left(a_{i}^{*}\right)-\phi^{2}\left(a_{i}^{*}\right)\right\}}{\sigma^{2}\left[1-\Phi\left(a_{i}^{*}\right)\right]^{2}}\\ &-\sum_{i\in S_{0}}\frac{x_{{ik}}x_{{il}}\left\{ \begin{array}{c} \left[\nu_{i}^{*}\phi\left(\nu_{i}^{*}\right)-a_{i}^{*}\phi\left(a_{i}^{*}\right)\right]\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]\\ +\left[\phi\left(\nu_{i}^{*}\right)-\phi(a_{i}^{*})\right]^{2} \end{array} \right\}}{\sigma^{2}[\Phi(\nu_{i}^{*})-\Phi(a_{i}^{*})]^{2}}\\ &-\frac{1}{\sigma^{2}}\sum_{i\in S_{1}}x_{{ik}}x_{{il}}\ {for} \ k= 1,\dots,p-1 \ {and} \ l = k,\dots,p-1, \end{aligned} $$
(11)
$$ {\begin{aligned} \frac{\partial^{2} l}{\partial \beta_{k} \partial \ln(\sigma)}&=\sum_{i=1}^{n}\frac{x_{{ik}}\left\{\left[1-\Phi\left(a_{i}^{*}\right)\right]\phi\left(a_{i}^{*}\right)\left[1-\left(a_{i}^{*}\right)^{2}\right]+a_{i}^{*}\phi^{2}\left(a_{i}^{*}\right)\right\}}{\sigma\left[1-\Phi\left(a_{i}^{*}\right)\right]^{2}}-\\ &\quad\sum\limits_{i\in S_{0}} \frac{x_{{ik}}\left\{ \begin{array}{c} \left(\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right)\left(\phi\left(\nu_{i}^{*}\right)\left[1-\left(\nu_{i}^{*}\right)^{2}\right]-\phi\left(a_{i}^{*}\right)\left[1-\left(a_{i}^{*}\right)^{2}\right]\right)\\ -\left[\phi\left(\nu_{i}^{*}\right)-\phi\left(a_{i}^{*}\right)\right]\left[\phi\left(\nu_{i}^{*}\right)\nu_{i}^{*}-\phi\left(a_{i}^{*}\right)a_{i}^{*}\right] \end{array}\right\}}{\sigma\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]^{2}}\\ &\quad-\frac{2}{\sigma^{2}}\sum\limits_{i\in S_{1}}x_{{ik}}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right) \ {for} \ k = 1,\dots,p-1, \end{aligned}} $$
(12)
$$ {\begin{aligned} \frac{\partial^{2} l}{\partial \ln^{2}(\sigma)}&=\sum\limits_{i=1}^{n}\frac{a_{i}^{*}\left\{\left[1-\Phi\left(a_{i}^{*}\right)\right]\phi\left(a_{i}^{*}\right)\left[1-\left(a_{i}^{*}\right)^{2}\right]+a_{i}^{*}\phi^{2}\left(a_{i}^{*}\right)\right\}}{\left[1-\Phi\left(a_{i}^{*}\right)\right]^{2}}\\ &-\sum\limits_{i\in S_{0}}\frac{\left\{ \begin{array}{c} \left(\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right)\left(\phi\left(\nu_{i}^{*}\right)\left[\left(\nu_{i}^{*}\right)^{3}-\nu_{i}^{*}\right]-\phi\left(a_{i}^{*}\right)\left[\left(a_{i}^{*}\right)^{3}-a_{i}^{*}\right]\right)\\ +\left[\phi\left(\nu_{i}^{*}\right)\nu_{i}^{*}-\phi\left(a_{i}^{*}\right)a_{i}^{*}\right]^{2} \end{array} \right\}} {\left[\Phi\left(\nu_{i}^{*}\right)-\Phi\left(a_{i}^{*}\right)\right]^{2}}\\ &-\frac{2}{\sigma^{2}}\sum\limits_{i\in S_{1}}\left(y_{i}-\mathbf{x}_{i}^{T}\boldsymbol{\beta}\right)^{2}. \end{aligned}} $$
(13)

Gradient for heteroskedastic model

Assuming that there are samples from J independent populations as discussed in “Heteroskedastic variances” section, the parameters to be estimated are \(\boldsymbol {\theta } = \left (\boldsymbol {\beta }, \sigma _{1},\dots,\sigma _{J}\right)^{T}\) where \(\boldsymbol {\beta }=\left (\beta _{1},\dots,\beta _{p-1}\right)^{T}\). Again, the form of log likelihood from Eq. 7 suggests calculating the partial derivatives with respect to ln(σj) rather than σj is optimal. Let \(c_{{ij}}^{*}=\frac {c-\mathbf {x}_{{ij}}^{T}\boldsymbol {\beta }}{\sigma _{j}}\) for some constant c.

We define the gradient as

$$\nabla l(\boldsymbol{\theta})=\left(\begin{array}{cc} \frac{\partial l}{\partial \beta_{1}}\\ \vdots \\ \frac{\partial l}{\partial \beta_{p-1}} \\ \frac{\partial l}{\partial \ln(\sigma_{1})} \\ \vdots \\ \frac{\partial l}{\partial \ln(\sigma_{J})} \end{array}\right), \ {for} \ p\ge2 \ {and} \ J\ge1. $$

Note that the case where J=1 is equivalent to the case in “Formula for the gradient” section.

Taking the derivative of the log-likelihood with respect to βk,

$$ {}\begin{aligned} \frac{\partial l}{\partial \beta_{k}}=&\sum\limits_{j=1}^{J}\left\{\sum_{i=1}^{n_{j}}-\frac{x_{{ijk}} \phi\left(a_{{ij}}^{*}\right)}{\sigma_{j}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]} -\sum_{i\in S_{0j}}\frac{x_{{ijk}}\left[\phi\left(\nu_{{ij}}^{*}\right)-\phi\left(a_{{ij}}^{*}\right)\right]} {\sigma_{j}\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]}\right.\\ &\left.+\frac{1}{\sigma_{j}^{2}}\sum\limits_{i\in S_{1j}}x_{{ijk}}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right)\right\} \ \ {for} \ k = 1,\dots, p-1, \end{aligned} $$
(14)

and with respect to ln(σj),

$$ \begin{aligned} \frac{\partial l}{\partial \ln(\sigma_{j})}=&\sum\limits_{i=1}^{n_{j}}-\frac{a_{{ij}}^{*}\phi\left(a_{{ij}}^{*}\right)}{1-\Phi\left(a_{{ij}}^{*}\right)} -\sum\limits_{i\in S_{0j}}\frac{\left[\nu_{{ij}}^{*}\phi\left(\nu_{{ij}}^{*}\right)-a_{{ij}}^{*}\phi\left(a_{{ij}}^{*}\right)\right]}{\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)} -n_{1j}\\ &+\frac{1}{\sigma_{j}^{2}}\sum\limits_{i\in S_{1j}}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right)^{2} \ \ {for} \ j = 1,\dots, J. \end{aligned} $$
(15)

Therefore, the gradient vector is

$$ \nabla l(\boldsymbol{\theta})=\left(\begin{array}{cccc} \sum_{j=1}^{J}\sum_{i=1}^{n_{j}}-\frac{x_{ij1}\phi\left(a_{{ij}}^{*}\right)}{\sigma_{j}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]} -\sum_{i\in S_{0j}}\frac{x_{ij1}\left[\phi\left(\nu_{{ij}}^{*}\right)-\phi\left(a_{{ij}}^{*}\right)\right]}{\sigma_{j}\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]}\\ +\frac{1}{\sigma_{j}^{2}}\sum_{i\in S_{1j}}x_{ij1}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right) \\ \vdots \\ \sum_{j=1}^{J}\sum_{i=1}^{n_{j}}-\frac{x_{ij(p-1)}\phi\left(a_{{ij}}^{*}\right)}{\sigma_{j}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]} -\sum_{i\in S_{0j}}\frac{x_{ij(p-1)}\left[\phi\left(\nu_{{ij}}^{*}\right)-\phi\left(a_{{ij}}^{*}\right)\right]}{\sigma_{j}\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]}\\ +\frac{1}{\sigma_{j}^{2}}\sum_{i\in S_{1j}}x_{ij(p-1)}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right)\\ \sum_{i=1}^{n_{1}}-\frac{a_{i1}^{*}\phi\left(a_{i1}^{*}\right)}{1-\Phi\left(a_{i1}^{*}\right)} -\sum_{i\in S_{01}}\frac{\left[\nu_{i1}^{*}\phi\left(\nu_{i1}^{*}\right)-a_{i1}^{*}\phi\left(a_{i1}^{*}\right)\right]}{\Phi\left(\nu_{i1}^{*}\right)-\Phi\left(a_{i1}^{*}\right)} -n_{11}\\ +\frac{1}{\sigma_{1}^{2}}\sum_{i\in S_{11}}\left(y_{i1}-\mathbf{x}_{i1}^{T}\boldsymbol{\beta}\right)^{2}\\ \vdots \\ \sum_{i=1}^{n_{J}}-\frac{a_{{iJ}}^{*}\phi\left(a_{{iJ}}^{*}\right)}{1-\Phi\left(a_{{iJ}}^{*}\right)} -\sum_{i\in S_{0J}}\frac{\left[\nu_{{iJ}}^{*}\phi\left(\nu_{{iJ}}^{*}\right)-a_{{iJ}}^{*}\phi\left(a_{{iJ}}^{*}\right)\right]}{\Phi\left(\nu_{{iJ}}^{*}\right)-\Phi\left(a_{{iJ}}^{*}\right)} -n_{1J}\\ +\frac{1}{\sigma_{J}^{2}}\sum_{i\in S_{1J}}\left(y_{{iJ}}-\mathbf{x}_{{iJ}}^{T}\boldsymbol{\beta}\right)^{2} \end{array}\right). $$
(16)

Hessian for heteroskedastic model

The Hessian matrix for parameters θ in “Gradient for heteroskedastic model” section is derived by taking further partial derivatives. This matrix takes the form

$$\nabla^{2}l(\boldsymbol{\theta})=\left[\begin{array}{cccccc} \frac{\partial^{2} l}{\partial \beta_{1}^{2}} & \cdots & \frac{\partial^{2} l}{\partial \beta_{1} \partial\beta_{p-1}} & \frac{\partial^{2} l}{\partial \beta_{1} \partial\ln(\sigma_{1})} & \cdots & \frac{\partial^{2} l}{\partial \beta_{1} \partial\ln(\sigma_{J})}\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} l}{\partial\beta_{p-1} \partial \beta_{1}}& \cdots & \frac{\partial^{2} l}{\partial \beta_{p-1}^{2}} & \frac{\partial^{2} l}{\partial \beta_{p-1} \partial \ln(\sigma_{1})} & \cdots & \frac{\partial^{2} l}{\partial \beta_{p-1} \partial \ln(\sigma_{J})}\\ \frac{\partial^{2} l}{\partial\ln(\sigma_{1}) \partial \beta_{1}} & \cdots & \frac{\partial^{2} l}{\partial \ln(\sigma_{1}) \partial \beta_{p-1}} & \frac{\partial^{2} l}{\partial \ln^{2}(\sigma_{1})} & \cdots & \frac{\partial^{2} l}{\partial \ln(\sigma_{1})\partial \ln(\sigma_{J})} \\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} l}{\partial\ln(\sigma_{J}) \partial \beta_{1}} & \cdots & \frac{\partial^{2} l}{\partial \ln(\sigma_{J}) \partial \beta_{p-1}} & \frac{\partial^{2} l}{\partial \ln\left(\sigma_{J}\right)\partial \ln(\sigma_{1})} & \cdots &\frac{\partial^{2} l}{\partial \ln^{2}\left(\sigma_{j}\right)} \end{array}\right]. $$

Note that since the groups are assumed to be independent, \(\frac {\partial ^{2} l}{\partial \ln (\sigma _{j})\partial \ln (\sigma _{k})}=0\) for all jk, which reduces the Hessian matrix to

$$\nabla^{2}l(\boldsymbol{\theta})=\left[\begin{array}{cccccc} \frac{\partial^{2} l}{\partial \beta_{1}^{2}} & \cdots & \frac{\partial^{2} l}{\partial \beta_{1} \partial\beta_{p-1}} & \frac{\partial^{2} l}{\partial \beta_{1} \partial\ln(\sigma_{1})} & \cdots & \frac{\partial^{2} l}{\partial \beta_{1} \partial\ln(\sigma_{J})}\\ \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\ \frac{\partial^{2} l}{\partial\beta_{p-1} \partial \beta_{1}}& \cdots & \frac{\partial^{2} l}{\partial \beta_{p-1}^{2}} & \frac{\partial^{2} l}{\partial \beta_{p-1} \partial \ln(\sigma_{1})} & \cdots & \frac{\partial^{2} l}{\partial \beta_{p-1} \partial \ln(\sigma_{J})}\\ \frac{\partial^{2} l}{\partial\ln(\sigma_{1}) \partial \beta_{1}} & \cdots & \frac{\partial^{2} l}{\partial \ln(\sigma_{1}) \partial \beta_{p-1}} & \frac{\partial^{2} l}{\partial \ln^{2}(\sigma_{1})} & & 0 \\ \vdots & \ddots & \vdots & & \ddots & \\ \frac{\partial^{2} l}{\partial\ln\left(\sigma_{J}\right) \partial \beta_{1}} & \cdots & \frac{\partial^{2} l}{\partial \ln\left(\sigma_{J}\right) \partial \beta_{p-1}} & 0 & &\frac{\partial^{2} l}{\partial \ln^{2}\left(\sigma_{j}\right)} \end{array}\right]. $$

The individual components of this matrix are calculated as

$$ \begin{aligned} \frac{\partial^{2} l}{\partial \beta_{k}\partial\beta_{l}}=&\sum_{j=1}^{J}\left\{\sum\limits_{i=1}^{n_{j}}-\frac{x_{{ijk}}x_{{ijl}}\left\{a_{{ij}}^{*}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]\phi\left(a_{{ij}}^{*}\right)-\phi^{2}\left(a_{{ij}}^{*}\right)\right\}}{\sigma_{k}\sigma_{l}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}} \right.\\ &-\sum\limits_{i\in S_{0j}}\frac{x_{{ijk}}x_{{ijl}}\left\{ \begin{array}{c} \left[\nu_{{ij}}^{*}\phi\left(\nu_{{ij}}^{*}\right)-a_{{ij}}^{*}\phi\left(a_{{ij}}^{*}\right)\right]\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]\\ +\left[\phi\left(\nu_{{ij}}^{*}\right)-\phi\left(a_{{ij}}^{*}\right)\right]^{2} \end{array} \right\}}{\sigma_{k}\sigma_{l}\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}}\\ &\left. -\frac{1}{\sigma^{2}_{j}}\sum\limits_{i\in S_{1j}}x_{{ijk}}x_{{ijl}}\right\}\ {for} \ k= 1,\dots,p-1 \ {and} \ l = k,\dots,p-1, \end{aligned} $$
(17)
$$ \begin{aligned} &\frac{\partial^{2} l}{\partial \beta_{k} \partial \ln(\sigma_{j})}=\sum_{i=1}^{n_{j}}\frac{x_{{ijk}}\left\{\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]\phi\left(a_{{ij}}^{*}\right)\left[1-\left(a_{{ij}}^{*}\right)^{2}\right]+a_{{ij}}^{*}\phi^{2}\left(a_{{ij}}^{*}\right)\right\}}{\sigma_{j}\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}}\\ &-\sum_{i\in S_{0j}} \frac{x_{{ijk}}\left\{ \begin{array}{c} \left(\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right)\left(\phi\left(\nu_{{ij}}^{*}\right)\left[1-\left(\nu_{{ij}}^{*}\right)^{2}\right]- \phi\left(a_{{ij}}^{*}\right)\left[1-\left(a_{{ij}}^{*}\right)^{2}\right]\right)\\ -\left[\phi\left(\nu_{{ij}}^{*}\right)-\phi\left(a_{{ij}}^{*}\right)\right]\left[\phi\left(\nu_{{ij}}^{*}\right)\nu_{{ij}}^{*}-\phi\left(a_{{ij}}^{*}\right)a_{{ij}}^{*}\right] \end{array} \right\}}{\sigma_{j}\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}}\\ &-\frac{2}{\sigma_{j}^{2}}\sum_{i\in S_{1j}}x_{{ijk}}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right) \ {for} \ k = 1,\dots,p-1 \ {and} \ j = 1,\dots,J, \end{aligned} $$
(18)
$$ \begin{aligned} \frac{\partial^{2} l}{\partial \ln^{2}\left(\sigma_{j}\right)}=&\sum\limits_{i=1}^{n_{j}}\frac{a_{{ij}}^{*}\left\{\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]\phi\left(a_{{ij}}^{*}\right)\left[1-\left(a_{{ij}}^{*}\right)^{2}\right]+a_{{ij}}^{*}\phi^{2}\left(a_{{ij}}^{*}\right)\right\}}{\left[1-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}}\\ &-\sum\limits_{i\in S_{0j}}\frac{ \left\{ \begin{array}{c} \left(\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right)\left(\phi\left(\nu_{{ij}}^{*}\right)\left[\left(\nu_{{ij}}^{*}\right)^{3}-\nu_{{ij}}^{*}\right]-\phi\left(a_{{ij}}^{*}\right)\left[\left(a_{{ij}}^{*}\right)^{3}-a_{{ij}}^{*}\right]\right)\\ +\left[\phi\left(\nu_{{ij}}^{*}\right)\nu_{{ij}}^{*}-\phi\left(a_{{ij}}^{*}\right)a_{{ij}}^{*}\right]^{2} \end{array} \right\} } {\left[\Phi\left(\nu_{{ij}}^{*}\right)-\Phi\left(a_{{ij}}^{*}\right)\right]^{2}}\\ &-\frac{2}{\sigma_{j}^{2}}\sum_{i\in S_{1j}}\left(y_{{ij}}-\mathbf{x}_{{ij}}^{T}\boldsymbol{\beta}\right)^{2} \ {for} \ j = 1,\dots, J. \end{aligned} $$
(19)

Appendix B: Additional figures

Fig. 11
figure11

Average Log Mean Squared Error for δ from Six Different Estimation Methods in Two Population Model. The vertical dashed black line corresponds to the case when δ=0, i.e., μ1=μ2. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

Fig. 12
figure12

Performance Metrics for Common σ from Six Different Estimation Methods in Two Population Model. The vertical dashed black line corresponds to the case when δ=0, i.e., μ1=μ2. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

Fig. 13
figure13

Performance of Maximum Likelihood Estimate for σ as Function of Sample Size. The vertical dashed black line on the left figure corresponds to zero bias. GS = Gold Standard, i.e., uncensored observations with truncation adjustment; Uncens NT = Uncensored data with no truncation adjustment; DL = detection limit; Tobit = Tobit censored regression with no truncation adjustment; tcensReg = Censored regression with truncation adjustment

Abbreviations

CI:

confidence interval

CPD:

cycles per degree

DL:

Detection Limit

GS:

gold standard

IOL:

intraocular lens

ISO:

international standards organization

MSE:

mean squared error

tcensReg:

censored regression with truncation adjustment

TN:

truncated normal

Uncens NT:

uncensored with no truncation adjustment

References

  1. 1

    Hornung RW, Reed LD. Estimation of average concentration in the presence of nondetectable values. Appl Occup Environ Hyg. 1990; 5(1):46–51.

    CAS  Google Scholar 

  2. 2

    Lubin JH, Colt JS, Camann D, Davis S, Cerhan JR, Severson RK, Bernstein L, Hartge P. Epidemiologic evaluation of measurement data in the presence of detection limits. Environ Health Perspect. 2004; 112(17):1691–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  3. 3

    Schisterman EF, Vexler A, Whitcomb BW, Liu A. The limitations due to exposure detection limits for regression models. Am J Epidemiol. 2006; 163(4):374–83.

    PubMed  PubMed Central  Google Scholar 

  4. 4

    Helsel DR. Less than obvious-statistical treatment of data below the detection limit. Environ Sci Technol. 1990; 24(12):1766–74.

    CAS  Google Scholar 

  5. 5

    Analytical Methods Committee. Recommendations for the definition, estimation and use of the detection limit. Analyst. 1987; 112(2):199–204.

    Google Scholar 

  6. 6

    Zaugg SD, Sandstrom MW, Smith SG, Fehlberg KM. Methods of Analysis by the U.S. Geological Survey National Water Quality Laboratory-Determination of Pesticides in Water by C-18 Solid Phase Extraction and Capillary-Column Gas Chromatography/Mass Spectrometry with Selected-Ion Monitoring. Open-File Report 95-181, U.S. Geological Survey, Denver, Colorado. 1995. http://pubs.er.usgs.gov/publication/ofr95181.

  7. 7

    McDonald JF, Moffitt RA. The uses of tobit analysis. Rev Econ Stat. 1980; 62(2):318–21.

    Google Scholar 

  8. 8

    Greene WH. Econometric Analysis. In: Chap. 19. "Limited Dependent Variables, Truncation, Censoring and Sample Selection". 8th edn. United Kingdom: Pearson: 2018.

    Google Scholar 

  9. 9

    Hald A. Maximum likelihood estimation of the parameters of a normal distribution which is truncated at a known point. Skandinavisk Aktuarietidskrift. 1949; 32:119–34.

    Google Scholar 

  10. 10

    Gupta AK. Estimation of the mean and standard deviation of a normal population from a censored sample. Biometrika. 1952; 39(3):260–73.

    Google Scholar 

  11. 11

    Harter HL, Moore AH. Iterative maximum-likelihood estimation of the parameters of normal population from singly and doubly censored samples. Biometrika. 1966; 53(1-2):205–13.

    CAS  PubMed  Google Scholar 

  12. 12

    Tiku M. Estimating the mean and standard deviation from a censored normal sample. Biometrika. 1967; 54(1):155–65.

    CAS  PubMed  Google Scholar 

  13. 13

    Sarhan AE, Greenberg BG. Estimation of location and scale parameters by order statistics from singly and doubly censored samples. Ann Math Stat. 1956; 27(2):427–51.

    Google Scholar 

  14. 14

    Dixon WJ. Simplified estimation from censored normal samples. Ann Math Stat. 1960; 31(2):385–91.

    Google Scholar 

  15. 15

    Tobin J. Estimation of relationships for limited dependent variables. Econometrica. 1958; 26(1):24–36.

    Google Scholar 

  16. 16

    Buckley J, James I. Linear regression with censored data. Biometrika. 1979; 66(3):429–36.

    Google Scholar 

  17. 17

    Cohen Jr AC. Estimating the mean and variance of normal populations from singly truncated and doubly truncated samples. Ann Math Stat. 1950; 21(4):557–69.

    Google Scholar 

  18. 18

    Halperin M. Maximum likelihood estimation in truncated samples. Ann Math Stat. 1952; 23(2):226–38.

    Google Scholar 

  19. 19

    Amemiya T. Regression analysis when the dependent variable is truncated normal. Econometrica. 1973; 41(6):997–1016.

    Google Scholar 

  20. 20

    Owsley C, Sloane ME. Contrast sensitivity, acuity, and the perception of ‘real-world’ targets. Br J Ophthalmol. 1987; 71(10):791–6.

    CAS  PubMed  PubMed Central  Google Scholar 

  21. 21

    Pelli DG, Bex P. Measuring contrast sensitivity. Vis Res. 2013; 90:10–14.

    PubMed  Google Scholar 

  22. 22

    ISO 11979-7. Ophthalmic implants – Intraocular lenses – Part 7: Clinical investigations of intraocular lenses for the correction of aphakia. Vernier, Geneva, CH: Standard, International Organization for Standardization; 2018.

    Google Scholar 

  23. 23

    Lange K. Numerical Analysis for Statisticians. 2nd Edn.New York: Springer; 2010.

    Google Scholar 

  24. 24

    Ypma TJ. Historical development of the newton-raphson method. Soc Ind Appl Math. 1995; 37(4):531–51.

    Google Scholar 

  25. 25

    Broyden CG. The Convergence of a Class of Double-rank Minimization Algorithms 1. General Considerations. IMA J Appl Math. 1970; 6(1):76–90.

    Google Scholar 

  26. 26

    Fletcher R. A new approach to variable metric algorithms. Comput J. 1970; 13(3):317–22.

    Google Scholar 

  27. 27

    Goldfarb D. A family of variable-metric methods derived by variational means. Math Comput. 1970; 24(109):23–26.

    Google Scholar 

  28. 28

    Shanno DF. Conditioning of quasi-Newton methods for function minimization. Math Comput. 1970; 24(111):647–56.

    Google Scholar 

  29. 29

    Fletcher R, Reeves CM. Function minimization by conjugate gradients. Comput J. 1964; 7(2):149–54.

    Google Scholar 

  30. 30

    Henningsen A, Toomet O. maxLik: A package for maximum likelihood estimation in R. Comput Stat. 2011; 26(3):443–58.

    Google Scholar 

  31. 31

    Henningsen A. Estimating Censored Regression Models in R using the censReg Package. R Package vignettes. 2010; 5(2):12.

  32. 32

    Croissant Y, Zeileis A. Truncreg: Truncated Gaussian Regression Models. R package version 0.2-5. 2018. https://CRAN.R-project.org/package=truncreg.

  33. 33

    R Core Team. R: A Language and Environment for Statistical Computing. Vienna, Austria: R Foundation for Statistical Computing; 2019. https://www.R-project.org/, R Foundation for Statistical Computing.

    Google Scholar 

  34. 34

    Burkardt J. The Truncated Normal Distribution [PDF File]. https://people.sc.fsu.edu/~jburkardt/presentations/truncated_normal.pdf. Accessed 9 June 2020.

  35. 35

    Veall MR, Zimmermann KF. Pseudo-r2 measures for some common limited dependent variable models. J Econ Surv. 1996; 10(3):241–59.

    Google Scholar 

  36. 36

    McKelvey RD, Zavoina W. A statistical model for the analysis of ordinal level dependent variables. J Math Sociol. 1975; 4(1):103–20.

    Google Scholar 

  37. 37

    Lindstrom MJ, Bates DM. Newton-raphson and em algorithms for linear mixed-effects models for repeated-measures data. J Am Stat Assoc. 1988; 83(404):1014–22.

    Google Scholar 

  38. 38

    Zhang X, Wan AT, Zhou SZ. Focused information criteria, model selection, and model averaging in a Tobit model with a nonzero threshold. J Bus Econ Stat. 2012; 30(1):132–42.

    Google Scholar 

Download references

Acknowledgements

Not applicable.

Funding

Not applicable.

Author information

Affiliations

Authors

Contributions

JRW performed analyses, implemented software, and wrote first draft of the article. HWK contributed to analyses, oversaw initial methodology design, and reviewed manuscript. CMC provided additional methodological extensions, edited manuscript, and oversaw analyses. The author(s) read and approved the final manuscript.

Corresponding author

Correspondence to Justin R. Williams.

Ethics declarations

Ethics approval and consent to participate

Data analyzed for the motivating example were collected as part of clinical trials conducted by Alcon Laboratories, Inc. Details on these trials can be found on ClinicalTrials.gov using the identifiers NCT01510717 and NCT01424189. All participants or subjects’ legally authorized representatives were provided and signed informed consent to participate in these trials.

Consent for publication

Not applicable.

Competing interests

Hyung-Woo Kim is a former employee of Alcon Laboratories, Inc.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Williams, J.R., Kim, H. & Crespi, C.M. Modeling observations with a detection limit using a truncated normal distribution with censoring. BMC Med Res Methodol 20, 170 (2020). https://doi.org/10.1186/s12874-020-01032-9

Download citation

Keywords

  • Contrast sensitivity
  • limited dependent variables
  • Limited domain
  • Visual acuity