How best to deal with a left-censored predictor (because of detection limits) in a multiple regression with interactions?

Question

Context: I'm new to Bayesian stats and am trying to fit a multiple regression with rstan. All variables are continuous and there is no hierarchical structure.

One of my predictors is left-censored because it falls below the detection limit for a chemical assay. What is the best way to deal with this in a multiple regression? So far, I can see a few possibilities:

1. A substitution rule, such as 'replace all values below the detection limit with a constant such as detection limit/2'. This is clearly not rigorous.
2. Multiple imputation, but (i) I don't know how to deal with the fact that values above the detection limit are likely to be generated by the imputation process, which I will know with high probability to be false, and (ii) I'm not sure how well multiple imputation plays with Bayesian approaches, since I can't think of a good way to aggregate the posterior distributions from fits to the different imputed datasets
3. Simulate values from a distribution that makes sense based on prior knowledge and the data, and randomly assign values below the detection limit to the relevant points. This suffers from similar problems to #2, since I would have to simulate many sets of values, model them separately, and then figure out how to integrate the posteriors.
4. Add a binary predictor indicating whether the censored predictor is above or below the detection limit. I'm not sure how best to integrate interactions with other variables in this case though. Would I need to build in interactions with both the binary and the censored predictor?

Am I missing better options? Are there useful Bayesian tricks that can help deal with this problem? I'm also open to non-Bayesian options.

The histogram below shows the distribution of values. The plot is on a log scale because that is most natural for this variable. For visual clarity, I have treated values below the detection limit (~25% of the data) as being 1/10 of the detection limit, and added a red line to separate them from the remaining points. Note that the red line is not the precise detection limit; the smallest quantified values to the right of the red line are at the putative limit. The fact that there are very few values exactly at the limit suggests that there may have been some variation in the detection limit between measurements, but I don't mind if that is ignored for the purposes of this question.

UPDATE:

Here's my rstan code, in case that is helpful. Betas 1 through 4 represent main effects, 5 & 6 are interaction terms (between 1 & 3 and 2 & 4). The censored predictor is therefore present in an interaction term as well, which is a complication I neglected to mention earlier.

data {
  int<lower=0> n;       // number of data items
  int<lower=0> k;       // number of predictors
  vector[n] Y;          // outcome vector
  matrix[n,k] X;        // predictor matrix
  int n2;               //the size of the new_X matrix
  matrix[n2,k] new_X;   //the matrix for the predicted values
}
parameters {
  real alpha; // intercept
  vector[k] beta; // coefficients for predictors
  real<lower=0> sigma; // error scale (cauchy truncated at zero)
}
model {
  beta[1] ~ normal(-0.75, 1);   //prior for beta
  beta[2] ~ normal(0, 3);   //prior for beta
  beta[3] ~ normal(0, 3);   //prior for beta
  beta[4] ~ normal(0, 3);   //prior for beta
  beta[5] ~ normal(0, 3);   //prior for beta
  beta[6] ~ normal(0, 3);   //prior for beta
  sigma ~ cauchy (0, 2.5);  //prior for sigma

  Y ~ normal(alpha + X * beta, sigma); // likelihood
}
generated quantities {
  vector[n2] y_pred;
  y_pred = new_X * beta; //the y values predicted by the model
}

Answer 1

rstan provides you with all the tools you need to solve this problem with Bayesian inference. In addition to the usual regression model of response $y$ in terms of predictors $x$, you should include a model of $x$ in the Stan code. This model should include the left-censoring. The Stan user manual chapter on censoring explains two different ways to do this in the Stan language. The first way is easier to incorporate into a regression model. The model for $x$ would look something like this (omitting the definition of N_obs and such):

data {
  real x_obs[N_obs];
}
parameters {
  real<upper=DL> x_cens[N_cens];
  real x[N];
}
model {
  x_obs ~ normal(mu, sigma);
  x_cens ~ normal(mu, sigma);
  x = append_array(x_obs, x_cens);
}

The key idea is that the censored data is represented by parameters whose upper limit is the detection limit. The censored data will be sampled alongside the other parameters in the model, so the posteriors you get will automatically integrate out the censored data.

Answer 2

Multiple imputation plays reasonably nicely with Bayesian inference. You just fit the Bayesian model on each imputation (making sure there's not too few, e.g. do at least 100 imputations or so) and then put the posterior samples together (=you use the mixture of the posteriors as the overall posterior). However, doing a good multiple imputation requires a multiple imputation tool that is aware of the left-censoring (if you ignore that, MI would more likely impute values like the non-censored observations). Technically, I think it would be valid to do multiple imputation and only select the imputation, for which values are below the limit of detection, but you very quickly get to where none of 1000s of imputations fulfill the criterion.

The substitution rule you mention apparently does not do too badly, if the censored quantity is the dependent variable in a model (see e.g. this paper for a list of references on the topic). How does it do for a covariate? No idea. I'd speculate it might be okay, if there's very few censored values. However, you have quite a few values that are censored.

The other obvious approach mentioned by Tom Minka is joint modeling of the covariate and the outcome of interest. I tried to really spell this out in Stan for an example like yours with a bit of made up data. I suspect that as usual my Stan program is not as efficiently written as it could be, but at least I hope it is reasonably clear.

library(rstan)

stancode = "
data {
  int<lower=0> N_obs; // Number of observation
  real y[N_obs]; // Observed y-values
  
  real x[N_obs]; // observed value or limit below which x is left-censored when x_censored=1
  int x_censored[N_obs]; // 1=left-censored, 0=not censored, 2=right-censored
  real measurement_error[N_obs]; // measurement error we know for the covariates
}
parameters {
  real mu; // intercept for the regression model for y
  real<lower=0> sigma; // residual SD for the regression model for y
  real beta; // regression coefficient for x in the regression model for y
  
  real x_randomeff[N_obs]; // A random effect we use to capture the underlying true value 
     // (obtained by multiplying by sigmax and adding mux - for more on the rationale for this parameterization look "non-centralized parameterization")
  real mux; // True population mean of the covariate values
  real<lower=0> sigmax; // True population SD of the covariate values
}
transformed parameters {
  real x_imputed[N_obs]; // Imputed values for x (or rather log(x))
  for (r in 1:N_obs){
    x_imputed[r] = mux + x_randomeff[r] * sigmax;
  }
}
model {
  // Specifying some wide weakly informative priors
  mu ~ normal(0, 100);
  sigma ~ normal(0, 100);
  beta ~ normal(0, 100);
  mux ~ normal(0, 10);
  sigmax ~ normal(0, 10);
  
  x_randomeff ~ normal(0,1);
  
  for (r in 1:N_obs){
    // Dealing with the covariate model
    if (x_censored[r]==1){
      target += normal_lcdf(x[r] | x_imputed[r], measurement_error[r]);
    } else if (x_censored[r]==2){
      target += normal_lccdf(x[r] | x_imputed[r], measurement_error[r]);
    } else {
      x[r] ~ normal(x_imputed[r], measurement_error[r]);
    }
    
    // fitting the regression model for y
    y[r] ~ normal(mu + x_imputed[r]*beta, sigma);
  }
  
}
"

sfit = stan(model_code = stancode,
         data=list(N_obs=12,
                   y=c(44, 40, 37, 33, 31, 27, 24, 19, 16, 13, 9, 6),
                   x=log( c(15,  7,  5,  3,  0.9, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5, 0.5) ),
                   x_censored = c(rep(0,5), rep(1, 7)),
                   measurement_error=rep(0.1, 12)),
         control=list(adapt_delta=0.95))

summary(sfit)$summary

As you can see the model even outputs what it imputed for the missing values. There's probably other ways of doing this, but this seemed reasonably intuitive to me. At the moment, I am using $log(x)\times \beta$ in the regression equation, but you could change that by exponentiating x_imputed[r].

Update: this paper just popped up in my Twitter feed.

Answer 3

In McElreath's Statistical Rethinking (2020) he gives an example almost exactly like what you are describing, where in chemical analyses there is a threshold below which something (e.g. the concentration of a specific compound) cannot be measured. In this case he discusses the use of a hurdle model. From what I'm reading on them, they could be applicable to your analysis, and they can also be fit relatively easily using Stan.

https://mc-stan.org/docs/2_20/stan-users-guide/zero-inflated-section.html

McElreath, R. (2020). Statistical rethinking: A Bayesian course with examples in R and Stan. CRC press.

Answer 4

Here's a somewhat related question: How small a quantity should be added to x to avoid taking the log of zero?

This looks like a very relevant paper that uses Bayesian regression with LOD censored predictors: https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6241297/

One simple and maybe less than ideal option is to add an indicator variable for whether an observation is below the LOD or not.