{'fraction_non_alphanumeric': 0.04149979975971165, 'fraction_numerical': 0.01566880256307569, 'mean_word_length': 4.328487596692451, 'pattern_counts': {'":': 0, '<': 32, '<?xml version=': 0, '>': 8, 'https://': 2, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 21, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Addressing Detection Limits with Semiparametric Cumulative Probability Models

Detection limits (DLs), where a variable is unable to be measured outside of a certain range, are common in research. Most approaches to handle DLs in the response variable implicitly make parametric assumptions on the distribution of data outside DLs. We propose a new approach to deal with DLs based on a widely used ordinal regression model, the cumulative probability model (CPM). The CPM is a type of semiparametric linear transformation model. CPMs are rank-based and can handle mixed distributions of continuous and discrete outcome variables. These features are key for analyzing data with DLs because while observations inside DLs are typically continuous, those outside DLs are censored and generally put into discrete categories. With a single lower DL, the CPM assigns values below the DL as having the lowest rank. When there are multiple DLs, the CPM likelihood can be modified to appropriately distribute probability mass. We demonstrate the use of CPMs with simulations and two HIV data examples. The first example models a biomarker in which 15% of observations are below a DL. The second uses multi-cohort data to model viral load, where approximately 55% of observations are outside DLs which vary across sites and over time.


Introduction
Detection limits (DLs) are not uncommon in biomedical research and other fields. For example, radiation doses may only be detected above a certain threshold (Wing et al., 1991), antibody concentrations may not be measured below certain levels (Wu et al., 2001), and X-rays may have lower limits of detection (Pan et al., 2017). In HIV research, viral load can only be detected above certain levels. To complicate matters, DLs often vary by assay and may change over time. For example, HIV viral load assays have had lower DLs at 400, 300, 200, 50, and 20 copies/mL depending on the commercial assay and year of application (Steegen et al., 2007).
Different types of analysis methods to handle DLs have been proposed for different purposes. In this manuscript, we will focus on studying the association between an outcome variable and covariates, where the outcome variable is subject to DLs. This is typically achieved with some sort of regression model. One common and simple method is to dichotomize the outcome as detectable or undetectable, and then to perform logistic regression (Jiamsakul et al., 2017). While it can be useful for some purposes, the dichotomization leads to information loss since the observed values inside the DLs are treated as if they are the same. Another common approach for handling DLs is substitution, where all nondetects are imputed with a single constant and a linear regression model is fit. The imputed constant may be, for example, the DL itself, DL/2, DL/ √ 2 (Hornung and Reed, 1990;Lubin et al., 2004;Helsel, 2011), or the expectation of the measurement conditional on being outside the DL under some assumed parametric model (Garland et al., 1993). For example, DL/2 corresponds to the expectation of a uniform distribution between 0 and the DL. Although simple, these substitution approaches typically result in biased estimation, underestimated variances, and thus sometimes wrong conclusions (Baccarelli et al., 2005;Fiévet and Della Vedova, 2010). In a third approach, one explicitly makes parametric assumptions on the distribution of the data, both within and outside the DLs. Parameters of interest can then be estimated by maximizing the censored data likelihood. Such a maximum likelihood approach is efficient and consistent when the distribution is correctly specified, but may perform poorly when distributional assumptions are incorrect.
To compound the problem, there is typically no way to examine model fit outside the DLs; goodness-of-fit of a parametric model inside DLs does not ensure goodness-of-fit outside the DLs (Baccarelli et al., 2005;Harel et al., 2014). A related, fourth approach for addressing DLs is to multiply impute values outside the DLs (Little and Rubin, 2019;Harel and Zhou, 2007). This approach may be computationally expensive, still requires parametric assumptions that can only be verified inside DLs, and may be particularly problematic with high rates of censoring or small sample sizes (Lubin et al., 2004;Zhang et al., 2009).
To avoid strong parametric assumptions, nonparametric methods such as Kaplan-Meier, score and rank-based methods have been proposed in two-sample comparisons (Helsel, 2011). Zhang et al. (2009) explored the use of the Wilcoxon rank sum test, other weighted rank tests, Gehan and Peto-Peto tests, and a novel nonparametric method for locationshift inference with DLs. Although attractive for two-sample tests, these nonparametric methods do not permit the inclusion of covariates.
In this manuscript, we propose a new approach for analyzing data subject to detection limits. Data with DLs effectively follow a mixture distribution, where those below a lower DL can be thought of as belonging to a discrete category, those above an upper DL belonging to another discrete category, while those inside the DLs are continuous. Whether discrete or continuous, the values are orderable. In earlier work, Liu et al. (2017) showed that continuous response variables can be modeled using a popular model for ordinal out-comes, namely the cumulative probability model (CPM), also known as the 'cumulative link model ' (Agresti, 2003). CPMs are a type of semiparametric linear transformation model, in which the continuous response variable after some unspecified monotonic transformation is assumed to follow a linear model, and the transformation is nonparametrically estimated (Zeng and Lin, 2007). These models are very flexible and can handle a wide variety of outcomes, including variables with DLs. Importantly, when fitting CPMs to data with DLs, minimal assumptions are made on the distribution of the response variable outside the DLs as these models are based on ranks, and values below/above DLs are simply the lowest/highest rank values. Because of their relationship to the Wilcoxon rank sum test (McCullagh, 1980), the CPM can be thought of as a semiparametric extension to permit covariates to the approaches that Zhang et al. (2009) found effective for handling DLs in two-sample comparisons. Finally, as will be shown, because CPMs model the conditional cumulative distribution function (CDF), it is easy to extract many different measures of conditional association from a single fitted model, including conditional quantiles, conditional probabilities, odds ratios, and probabilistic indexes, which permits flexible and compatible interpretation.
In Section 2, we review the CPM, illustrate its use for simple settings where there is only a single set of DLs, and then show how CPMs can be extended to address multiple DLs. We also propose a new method for estimating the conditional quantile from a CPM. In Section 3, we illustrate and demonstrate the advantages of the proposed approach using real data from two studies. The first study aims to measure the association between covariates and a biomarker whose values are below a DL in approximately 15% of observations. The second example is a large multi-cohort study of viral load (VL) after starting antiretroviral therapy among persons with HIV, where most observations are below DLs, but the DLs vary across sites and change over time. In Section 4, we demonstrate the performance of our method with simulations. The final section contains a discussion of the strengths and limitations of our method and future work.

Cumulative Probability Models
Transformation is often needed for regression of a continuous outcome variable Y to satisfy model assumptions, but specifying the correct transformation can be difficult. In a linear transformation model, the outcome is modeled as Y = H(β T X + ), where H(·) is an unknown monotonically increasing transformation, X is a vector of covariates, and follows a known distribution with CDF F . This linear transformation model can be equivalently expressed in terms of the conditional CDF, Let G = F −1 and α = H −1 ; α(·) is monotonically increasing but otherwise unknown. Then where G serves as a link function and the model becomes a cumulative probability model (CPM). The intercept function α(y) is the transformation of the response variable such that α(Y ) = β T X + . The β coefficients indicate the association between the response variable and covariates: fixing other covariates, a positive/negative β j means that an increase in X j is associated with a stochastic increase/decrease in the distribution of the response variable.
In the CPM (1), the intercept function α(y) can be nonparametrically estimated with a step function (Zeng and Lin, 2007;Liu et al., 2017). This allows great model flexibility.
Consider an iid dataset {(y i , x i ) : i = 1, . . . , n}. The nonparametric likelihood is where F (y − i |x i ) = lim t↑y i F (t|x i ). In nonparametric maximum likelihood estimation, the probability mass given any x will be distributed over the discrete set of observed outcome values. Thus we only need to consider functions for α(·) such that F (y|x i ) is a discrete distribution over the observed values. Let J be the number of distinct outcome values, denoted as a 1 < · · · < a J . Let S = {a 1 , . . . , a J }. These serve as the anchor points for the nonparametric likelihood. Let α j = α(a j ); then α 1 < · · · < α J . The nonparametric likelihood (2) can be written as Maximizing (4), we obtain the nonparametric maximum likelihood estimates (NPMLEs), (β,α), whereα = (α 1 , . . . ,α J−1 ). Note the multinomial form of the likelihood (4); because the probabilities in a multinomial likelihood add to one, α J is not estimated. Note also that the likelihood in (4) is identical to that of cumulative link models for ordinal data if the outcome Y is treated as ordinal with categories {a 1 , . . . , a J }. Liu et al. (2017) and Tian et al. (2020) have shown that CPMs can be fit to and work well for continuous and mixed types of responses. CPMs have also been shown to be consistent and asymptotically normal, with variance consistently estimated with the inverse of the information matrix under mild conditions including boundedness of the outcome variable (Li et al., 2022).
The NPMLEs and their estimated variances can be efficiently computed with the orm() function in the rms package in R (Harrell, 2020), which takes advantage of the tridiagonal nature of the Hessian matrix using Cholesky decomposition (Liu et al., 2017).
CPMs have several nice features. Some widely used regression methods model only one aspect of the conditional distributions (e.g., conditional mean for linear regression and conditional quantile for quantile regression). With the NPMLEs (β,α), we can estimate the conditional CDFs asF (y|x) = F (α j −β T x) where j is the index such that a j = max{a ∈ S : a ≤ y}; standard errors can be obtained by the delta method. Since conditional CDFs are directly modeled, other characteristics of the distribution, such as the conditional quantiles and conditional expectations, can be easily derived (Liu et al., 2017). Depending on the choice of link function, β may be interpretable; for example, with the logit link function, exp(β) is an odds ratio. Probabilistic indexes (De Neve et al., 2019), which are defined as Pr(Y 1 < Y 2 |X 1 , X 2 ), can also be easily derived; for example, with the logit

Single Detection Limits
In this subsection, we first present our method for the simple scenario that there is a single lower DL and/or a single upper DL. We will describe the general approach for multiple DLs in the next subsection.
Consider a dataset with a lower DL, l, and an upper DL, u. The outcome Y is observed if it is inside the DLs (i.e., l ≤ Y ≤ u) or censored if it is outside the DLs. The J distinct values of the observed outcomes are denoted as l ≤ a 1 < · · · < a J ≤ u. When there are no observations outside the DLs, these values are treated as ordered categories in CPMs and they are the anchor points in the nonparametric likelihood (4), and correspondingly there are J − 1 alpha parameters, α 1 < · · · < α J−1 . With observations outside the DLs, the likelihood (4) needs to be modified accordingly.
When there are observations below the lower DL, we do not know their values except that they are < l. As there is no way to distinguish them, we treat them as a single category, denoted as a 0 . Note that a 0 is not a value but a symbol for the additional category below a 1 . The nonparametric likelihood for a subject outcome censored at the lower DL l is where α 0 is the extra alpha parameter corresponding to category a 0 such that α 0 < α 1 .
Because a 1 , the previously lowest category, now has a category below it, the nonparametric likelihood for a subject with y i = a 1 becomes Similarly, when there are observations above the upper DL, they are also treated as a single category, denoted as a J+1 , which is a symbol for the additional category above a J .
The nonparametric likelihood for a subject censored at the upper DL u is Because a J is no longer the highest category, α J will need to be estimated, and the likelihood for a subject with y i = a J is now Put together, with observed data subject to a single lower DL and a single upper DL, the CPM likelihood is which is equivalent to (4) except with two new anchor points, a 0 and a J+1 . Therefore, (5) is maximized in an identical manner to (4), with outcomes below the lower DL and outcomes above the upper DL simply assigned to categories a 0 and a J+1 , respectively.
In summary, when there are data censored below the lower DL, we add a new anchor point a 0 < a 1 and a new parameter α 0 ; when there are data censored above the upper DL, we add a new anchor point a J+1 > a J and a new parameter α J . The alpha parameters to be estimated are (α 1 , . . . , α J−1 ) when there are no DLs or no data censored at DLs, (α 0 , α 1 , . . . , α J−1 , α J ) when both categories a 0 and a J+1 are added, (α 0 , α 1 , . . . , α J−1 ) when only a 0 is added, and (α 1 , . . . , α J−1 , α J ) when only a J+1 is added.
In practice, one can fit the NPMLE in these settings using the orm() function by setting outcomes below the lower DL to some arbitrary number < l and outcomes above the upper DL to some arbitrary number > u. Note that unlike single imputation approaches for dealing with DLs, the CPM estimation procedure is invariant to the choice of these numbers assigned to values outside the DLs. The CPM (1) assumes that after some unspecified transformation, the outcome follows a linear model both within and outside the DLs.
In contrast, parametric approaches to deal with DLs assume the full distribution of the outcome conditional on covariates is known, both within and outside DLs. Hence, CPMs make much weaker assumptions than fully parametric approaches.

Multiple Detection Limits
We now consider the general situation where data may be collected from multiple study Every subject has a vector X of covariates and three underlying random variables where Y is the true outcome and C L < C U are the lower and upper DLs.
When there is no upper DL, C U = ∞, and when there is no lower DL, C L = −∞. C L and C U are assumed to be independent of Y conditional on X; the vector X may contain variables for study sites or calendar time. This non-informative censoring assumption is typically plausible as DLs are determined by available equipment/assays.
We assume the CPM (1) holds for all subjects. Due to DLs, we may not always observe Y . Instead, we only observe (Z, ∆), where Z = max(min(Y, C U ), C L ) and ∆ is a variable indicating whether Y is observed or censored at a DL: ∆ = 1 and Z = Y if Y is observed, Given a dataset {(z i , δ i ; x i )} (i = 1, . . . , n), we first determine how many anchor points are needed to support the nonparametric likelihood of the CPM. Let J be the number of distinct values of z i among those with δ i = 1; they are denoted as a 1 < · · · < a J . For data without any DLs, these points are the anchor points for the nonparametric likelihood, and they are effectively treated as ordered categories in a CPM. Let S = {a 1 , · · · , a J } be the set of these values. When there are data with δ i = L, let l be the smallest z i with δ i = L.
Similarly, when there are data with δ i = U , let u be the largest we add a category into S below a 1 , denoted as a 0 ; note that it is not a value but a symbol for the additional category in S below a 1 . Similarly, if u ≥ a J , we add a J+1 into S, which is a symbol for the additional category above a J . Depending on the data, the number of ordered categories can be J, J + 1, or J + 2.
Consider the situation where both a 0 and a J+1 have been added to S (i.e., S = {a 0 , a 1 , . . . , a J , a J+1 }). When δ i = 1, the nonparametric likelihood for (z i , 1) is where j is the index such that a j = z i . When δ i = L, the nonparametric likelihood for where j is the index such that a j = max{a ∈ S : a < z i } when z i = l. When δ i = U , the where j is the index such that a j = min{a ∈ S : a > z i } when z i = u. The overall nonparametric likelihood is the product of these individual likelihoods over all subjects.
Note that if there are no uncensored observations between two lower (or upper) DLs, the two DLs are effectively treated as the same DL. A toy example to illustrate our definition is provided in Table S1 of the Supplementary Material. Slight modifications will be applied when no or only one additional category is added to S. When there is no need to add a 0 to S (i.e., when l > a 1 or there are no lower DLs), only the second row in the likelihood (7) for (z i , L) will be employed, and the likelihood for When there is no need to add a J+1 to S (i.e., when u < a J or there are no upper DLs), only the second row in the likelihood (8) for (z i , U ) will be employed, and the likelihood for Similar to the likelihood of CPM for data without DLs, the individual likelihoods presented above involve either one alpha parameter or two adjacent alpha parameters. As a result, the Hessian matrix continues to be tridiagonal, allowing us to use Cholesky decomposition to solve for the NPMLEs and efficiently estimate their variances. We have developed an R package, multipleDL available at https://github.com/YuqiTian35/multipleDL, which uses the optimizing() function in the rstan package to maximize the likelihood (Stan Development Team, 2020).

Interpretable Quantities and Conditional Quantiles
Interpretation of results after fitting CPMs to outcomes with DLs is similar to settings without DLs. Depending on the link function, β may be directly interpretable. The conditional CDF, probabilistic indexes, and conditional quantiles are also easily derived.
Note, however, that without additional assumptions on the distribution of the outcome outside DLs, conditional expectations cannot be estimated.
We now describe how to infer conditional quantiles from a CPM fitted on data with DLs.
The conditional CDF from a CPM for a given x can be computed asF where j is the index such that a j = max{a ∈ S : a ≤ y}; if there is no a ∈ S such that a ≤ y, thenF (y|x) = 0. For ease of presentation, we fix x and let P j =F (a j |x) (j = 0, 1, . . . , J, J + 1); for convenience, let P −1 = 0. Our goal is to define a quantile functionQ(p), where 0 < p < 1, for the conditional distribution given x.
The quantile function for a CDF F (·) is typically defined as Q(p) = inf{z : F (z) ≥ p}.
A plug-in estimator for an estimated CDF,F , isQ 0 (p) = inf{z :F (z) ≥ p}. When applied to our setting,Q 0 (p) = a j when P j−1 < p ≤ P j . This estimator may not be suitable for CPMs becauseF (·) is a step function and thereforeQ 0 (p) only takes values at the anchor points, which can be undesirable for continuous outcomes, especially when there is a large gap between adjacent anchor points.  Liu et al. (2017) proposed to estimate quantiles for CPMs with linear interpolation.
Specifically, given a fixed p, let j = j(p) be the index such that P j−1 < p ≤ P j . When p > P 0 , j ≥ 1 and defineQ 1 (p) = a j−1 +(p−P j−1 )/(P j −P j−1 )×(a j −a j−1 ), which is a linear interpolation between a j−1 and a j . When p ≤ P 0 ,Q 1 (p) is set to be a 0 . Recall that a 0 is not a value but a symbol for being below the lower DL, l; we thus relabel it as '<l', so when p ≤ P 0 ,Q 1 (p) = '<l'. For the linear interpolation between a 0 and a 1 , we set a 0 to be l.
Similarly, a J+1 is labeled '>u' and assigned the value u for the linear interpolation between a J and a J+1 .Q 1 (p) is illustrated as the dashed lines in Figure 1. An alternative definition is to interpolate between a j and a j+1 : Figure 1. For continuous data without DLs,Q 1 (p) andQ 2 (p) converge as the sample size increases. However, they are problematic for continuous data with DLs becausê Q 1 (p) < a J+1 for all p < 1 andQ 2 (p) > a 0 for all p > 0 even though there are non-zero estimated probabilities at the lower DL a 0 and upper DL a J+1 .
In situations where there is only a lower DL or an upper DL, our definition ofQ(p) is similar. Confidence intervals for the conditional quantiles can be estimated by applying a weighted linear interpolation to the confidence intervals of the conditional CDF similar to the above procedure (Liu et al., 2017).

Applications
In this section, we illustrate our method with two datasets, one from a biomarker study with a single lower DL and the other from a multi-center study with multiple DLs varying within and across centers.

Single Detection Limit
Our first example uses data from a study investigating the relationship between HIV, diabetes, obesity, and various biomarkers. Data were collected on 161 adults, some of whom were highly overweight (body mass index (BMI) ranged from 22 to 58 kg/m 2 ).
Several biomarkers were measured. Here, we focus on interleuken 4 (IL-4), a cytokine that is related to T-cell production and metabolism and has been seen to limit lipid accumulation in mice (Tsao et al., 2014). We examine the association between IL-4 and BMI, controlling for age, sex, HIV status, and diabetes status. Our measures of IL-4 had a single lower DL of 0.019 pg/ml, and 24 subjects (15%) had IL-4 values below the DL. The distribution of IL-4, which is right-skewed, is shown in Figure S1(A) in Supplemental Material S2.1.
We fit a CPM as described in Section 2.2, using the logit link; full results are in Table   S2 in the Supplementary Material S2.1. No transformation of IL-4 was needed. With the logit link function, the β parameters can be interpreted as log odds ratios. BMI was found to be negatively associated with IL-4 (p-value 0.023). Holding other covariates constant, a 5 kg/m 2 increase in BMI corresponded to a 22% decrease in the odds of having a higher IL-4 value (adjusted odds ratio 0.78, 95% confidence interval (CI) of (0.62, 0.97)).
The corresponding probabilistic index was 0.264, meaning that holding other variables of BMI could also be extracted from this analysis and is in Figure S3(B). The curve of conditional median as a function of BMI is similar to what was estimated with the CPM ( Figure S3(A)), but it is slightly lower and its confidence bands are tighter than those of the CPM. The tighter bands reflect the parametric assumption that the data are truly lognormally distributed. In contrast, the CPM does not require transformation of the data, and it non-parametrically estimates the best transformation. (iii) One could also directly estimate the conditional median as a function of BMI using quantile regression (Koenker and Hallock, 2001). This estimated curve is in Figure S3(C), which closely matches that estimated from the CPM. However, the confidence bands for median regression are wider than those of the CPM, and the 95% CI for the slope contains 0. One could argue that the CPM is assuming more than median regression (which only assumes a linear relationship between the median on the original outcome scale and the covariates); hence the narrower confidence bands. However, the CPM is able to yield several additional quantities (e.g., other quantiles, odds ratios, exceedance probabilities) from a single model that cannot be obtained from median regression. Also, the confidence bands obtained by the CPM do not go below the DL. (iv) Finally, one could dichotomize IL-4 into "undetectable" and "detectable" and fit a logistic regression model. However, logistic regression was not able to provide a stable estimation for this dichotomization. One could consider other dichotomiza-tions, but the choice is arbitrary. In fact, a beta coefficient in the CPM can be thought of as a weighted average of the log-odds ratios for logistic regression models that consider all possible orderable dichotomizations of the outcome.

Multiple Detection Limits
We illustrate our approach to handle multiple DLs with data from a multi-center HIV study. A traditional analysis in the HIV literature would dichotomize VL as detectable and undetectable and perform logistic regression (Jiamsakul et al., 2017). There are a few issues that make this analysis less than ideal. First, all VLs above the DL (nearly half of all observations) would be collapsed into a "detectable" category resulting in well-known loss of information due to dichotomizing continuous variables (Fedorov et al., 2009). Second, because the DL varies with time and by site, the analyst is forced to dichotomize at the largest DL (in this case 400 copies/mL) or else perform an analysis where values above the DL at one site are treated differently than they would be treated at another site. For example, a VL of 300 copies/mL measured in Mexico in 2005 would be measured as '<400' that same year in Peru; assigning this value as '<400' results in lost information but leaving it as "detectable" would make the outcome variable different across time and sites. A more parametric analysis might assume that the VL follows a specified distribution (e.g., lognormal distribution) and fit the censored data likelihood or multiply impute values below the DL from the assumed distribution to obtain estimated regression coefficients. However, distributional assumptions for values below the DL are strong and untestable, and given that over half of the response variables are below the DL, these assumptions would have a large impact on results.
In contrast, the CPM uses all available information (i.e., does not dichotomize the response variable) and makes much weaker assumptions than the fully parametric approaches. Similar to the parametric approaches, the CPM assumes non-informative censoring conditional on covariates (which is reasonable, given that DLs are determined by equipment/assays independent of true values) and that all observations follow a common distribution conditional on covariates, which permits borrowing information across sites  ART regimen, the time from ART initiation until the VL measurement, and calendar year are all associated with VL at 6 months. Holding other variables fixed, a 10-fold increase in VL at baseline is associated with a 44% increase in the odds of having a higher VL at 6 months (95% CI 34% to 54%).
Quantiles and cumulative probabilities are also easily extracted from the CPM. The first row of Figure 4 are the estimated conditional 50th and 90th percentiles of 6-month VL and the conditional probabilities for 6-month VL being greater than 1000 and 20 copies/mL as functions of age. The plots show that VL at 6 months is fairly similar across ages after fixing the other covariates. The smallest DL is 20 copies/mL, and all VL less than this DL belong to the smallest ordered category, which we label as '<20'. The second row of Figure  For comparisons, we also analyzed the data using competing approaches described earlier. First, we fit logistic regression to 6-month VL values dichotomized as <400 vs. ≥400 copies/mL, corresponding to the highest DL. Results are in Table S4 of the Supplementary Material S2.2. The CPM and the logistic regression model gave similar estimates of the beta coefficients (which are log odds ratios), although there were some differences in the estimates and the CIs from CPMs tend to be narrower, as expected. In logistic regression, the log odds ratios are based on the single undetectable vs. detectable dichotomization, while those in CPMs are based on dichotomizations at each response value. Second, we fit a full likelihood-based model assuming the outcome variable was normally distributed after log 10 (·) transformation. Note that even the log-transformed 6-month VL was still quite skewed (shown in Figure S4), and hence the assumptions of this fully parametric approach were questionable. The parameters in this approach and those from the CPM are not directly comparable because they are on different scales, however, the directions of associations were similar.

Simulations
Extensive simulations of CPMs with continuous data have been reported elsewhere (Liu et al., 2017;Tian et al., 2020). Here we present a limited set of simulations investigating the performance of CPMs with data subject to single and multiple DLs.

Single Detection Limits
Data were generated for sample sizes of n = 100 and n = 500 such that the outcome Y followed a normal linear model after log transformation in the following manner: ). In addition, we considered a scenario with a more complicated transformation: 6.
Note that the Y in scenario 6 is a monotonic transformation of that in scenario 2 with exactly the same censoring rate.
CPMs were fit to the observed data {X, Y } without any knowledge of the correct transformation or Y * . We simulated 10,000 replications under each scenario. Percent bias, root mean squared error (RMSE), and coverage of 95% CIs were estimated with respect to β, conditional medians for X = {0, 1}, and conditional CDFs at y = 1.5 for X = {0, 1}. Table 2 shows results under correctly specified models (i.e., probit link function and X correctly included). CPMs resulted in nearly unbiased estimation and good CI coverage. As the sample size increased, both the bias and RMSE decreased. Note that estimation of the condition medians was "perfect" in scenario 5 because the true conditional medians were below the lower DL due to the high censoring rate and the estimated conditional medians were always '<4', the lowest outcome category corresponding to below the DL. The estimate of β was more variable in scenario 5 because of the high censoring rate. The estimation of β in scenario 6, where data were generated from the complicated transformation, was exactly the same as that in scenario 2 because the same seed was used in all simulation scenarios and the order information above the DL was identical between scenarios 2 and 6. However, the conditional medians and CDFs depend on the scale of the outcome, and their estimates differed between scenarios 2 and 6. comparing CPMs with some widely used methods for handling DLs, specifically single imputation with l/2, single imputation with l/ √ 2, multiple imputation, and fully parametric maximum likelihood estimation (MLE). For all non-CPM approaches, we first correctly assumed that the outcome variable followed a log-normal distribution. With the imputation approaches, unobserved values were imputed, then a linear regression model was fit on the log-transformed outcome to obtain the β estimate, and median regression was used to estimate conditional medians. In multiple imputation, the correct tail distribution was used for imputing data and 10 iterations were performed for each data set. As expected, the MLE performed the best with the lowest bias and RMSE, and highest efficiency because the distributional assumptions matched the true distribution. The performance of multiple imputation was similar to that of the MLE, but with higher RMSE. As a semiparametric method, the CPM, also resulted in minimal bias and correct coverage but had slightly larger variance and RMSE. In contrast, the single imputation estimators were biased and tended to have poor coverage, especially for estimating β. We also evaluated the comparator methods under misspecification of the transformation. We simulated datasets with X ∼ N (5, 1), Y * = Xβ + , β = 1, ∼ N (0, 1), Y = Y * 2 , n = 1000, and l = 13.12 (approximately 17% censored). The non-CPM approaches assumed a normal linear model after an incorrectly specified log-transformation. As shown in the bottom half of Table   S5, only the CPM was able to properly estimate β and the conditional medians, because pre-transformation and strict distributional assumptions are not needed for fitting CPMs.
Finally, the Supplementary Material Table S6 shows the performance of CPMs for the data generated in scenario 2 under moderate and severe link function misspecification (i.e., fitting CPMs with logit and loglog link functions, respectively). Link function misspecification is equivalent to misspecification of the distribution of because F = G −1 . The performance of CPMs was reasonable with moderate link function misspecification with bias under 6% and coverage of 95% CI close to 0.95 with n = 100, although as low as 0.91 with n = 500. With severe link function misspecification, the performance of CPMs was noticeably worse, with bias as high as 12% and coverage as low as 0.60 for the conditional median at X = 1.

Multiple Detection Limits
To illustrate the use of CPMs with multiple detection limits, we simulated data from 3 study sites. The data were generated in a similar way as in Section 4.1, but different DLs were applied at different sites and the distribution of the covariate X was allowed to vary across sites in some scenarios. Specifically, we considered the following 5 scenarios: the quantiles at p = 0.03 (i.e., 3rd percentile) and CDFs at y = 0.05. Results from fitting the CPM based on 10,000 replications are shown in Supplementary Material S3. In summary, estimates had very low bias and confidence intervals had proper coverage in all simulation scenarios.

Discussion
In this paper, we have described an approach to address detection limits in response variables using CPMs. CPMs have several advantages over existing methods for addressing DLs. They make minimal distributional assumptions, they yield interpretable parameters, and they are invariant to the value assigned to measures outside DLs. Any values outside the lowest/highest DLs are simply assigned to the lowest/highest ordinal categories, and estimation proceeds naturally. CPMs are also easily extended to handle multiple DLs.
From simulation studies, we saw that CPMs performed well, even with high censoring rates and relatively small sample sizes. We also illustrated the use of CPMs with two quite different HIV datasets with censored response data. Similar datasets with limits of detection are quite common in biomedical research; the CPM is an effective analysis tool in these settings.
CPMs have some limitations. Although CPMs do not make distributional assumptions on the response variable, the link function must still be specified, which corresponds to making an assumption on the distribution of the response variable after an unspecified transformation. Performance can be poor with severe link function misspecification; however, CPMs appear to be fairly robust to moderate misspecification. In addition, because we do not make distributional assumptions outside DLs, we are not able to estimate conditional expectations after fitting a CPM; however, with DLs, conditional quantiles are probably more reasonable statistics to report anyway. The codes for applications and simulations are available at https://github.com/YuqiTian35/DetectionLimitCode.
Further research could consider extensions of CPMs to handle clustered or longitudinal data with DLs. It may be of interest to study the use of these models with right-censored failure time data (i.e., survival data), where each observation is potentially subject to a different censoring time; the current manuscript only considered situations with a relatively small number of potential censoring times (i.e., upper DLs).
Supplementary Material to "Addressing Detection Limits with Semiparametric Cumulative Probability Models" S1 Methods
Consider the logit link function, G(p) = log (p/(1 − p)). Then for j = 1, . . . , J −1, . Note that this equation also holds for j = 0 and j = J. Then Given data {(x i , y i )}, the log-likelihood is l(α, β) = i log p i,j(i) , where j(i) is the index such that a j(i) = y i . The score function with respect to β is Now consider the situation where there is a single binary X, coded as 0 and 1. Let n k be the number of observations with X = k (k = 0, 1), and n = n 0 + n 1 . Under the null of β = 0, the CPM estimate of α j isα j = G(P j ), whereP j is the CDF of the empirical distribution of {y i } at a j . As a result,γ i,j(i) =P j(i) . In this situation, the numerator of the score test statistic is Let k j be the number of observations with y i = a j , R j be the midrank for a j , and be the sum of the midranks for the observations with X = 1. Then Note that E(R 1 ) = n 1 (n + 1)/2 under the null of β = 0.

S1.2 Multiple Detection Limits
We show a toy example in Table S1 to illustrate the approach to handle multiple DLs described in Section 2.3.   The estimated transformation function. Figure S1 shows that the distribution of IL-4 is right-skewed and has a lower detection limit of 0.019. After log-transformation, the response variable is approximately normally distributed except for values below the DL. The results for fitting the CPM is shown in Table S2 and the estimated transformation function is in Figure S1. Based on the results of CPMs, we can derive conditional quantities. For example, Figure S2 shows the 90th percentile of IL-4 as a function of BMI, and the probabilities of IL-4 being greater than 0.019 (the DL) and greater than 0.05 as functions of BMI while setting other covariates to the median (for continuous covariates) or mode (for categorical covariates).
We fit median regression and the likelihood approach to the data and obtain the conditional medians. The comparison of the conditional medians is in Figure S3. Figure S4 shows the distribution of log-transformed 6-month VL and lower DLs.  The conditional median obtained by median regression.

S2.2 Multiple Detection Limits
Log10 transformed viral load    Figure S5 plots the estimated conditional quantiles and probabilities of the outcome being above 1000 and 20. The conditional median and its confidence intervals all fall into the smallest level "<20". For age, the general trend is the same as in the model without splines but for patients younger than 40 years old, age increase is correlated with a slight increase in 6-month viral load. The estimated conditional quantiles and probabilities as a function of prior AIDS events shown in the third row of Figure S5 are very similar to the results in the model without splines. In the second row of Figure S5, patients with large baseline viral load values tend to have much higher 6-month viral load measures. Results for analyzing the data using logistic regression and the full likelihood-based model are shown in Table S4.

S3 Simulation Results
We compared CPMs with some widely used methods for handling DLs under under scenario 2 with n = 1000 described in Section 4.1, specifically single imputation with l/2, single imputation with l/ √ 2, multiple imputation, and fully parametric maximum likelihood estimation (MLE). Results are in Table S5.
We show the performance of CPMs with moderate and severe link function misspecification in Table S6.
The simulation results for multiple detection limits are in Table S7.   The percent bias in scenario 4 is relatively high due to the small true values.