20
$\begingroup$

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. Histogram of

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
}
$\endgroup$
6
  • 2
    $\begingroup$ Is the DL fixed or variable? How do the DLs compare to the quantified values of the variable? (If they are relatively small, the fixed-substitution method such as DL/2 introduces little error.) $\endgroup$
    – whuber
    Aug 8, 2020 at 14:51
  • 1
    $\begingroup$ @whuber Thanks for the questions! Detection limit is fixed, mercifully. But it is unfortunately not very small relative to the distribution of quantified values; I have added a plot to illustrate this. $\endgroup$
    – mkt
    Aug 8, 2020 at 18:48
  • 1
    $\begingroup$ Can you include your existing rstan code without the censoring? If so, we can probably tell you how to modify it. $\endgroup$
    – Ben
    Aug 12, 2020 at 7:15
  • 1
    $\begingroup$ @Ben Done, thanks for taking a look! $\endgroup$
    – mkt
    Aug 12, 2020 at 7:47
  • 1
    $\begingroup$ @YairDaon Well, throwing away data is inefficient and so I'd prefer not to do that. It's not 'bad data' or missing data; it's just censored. There is still useful information contained in there that I would like to make use of if possible. Throwing away these points is also truncating the axis which can lead to errors in inference. $\endgroup$
    – mkt
    Aug 12, 2020 at 8:00

4 Answers 4

6
+200
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ Thanks very much! I saw that chapter but was not sure if the logic extended to predictors as well. Would you be able to tell me how the regression model then gets specified in the code block? I've added my code in the question in case it helps. $\endgroup$
    – mkt
    Aug 12, 2020 at 7:52
  • $\begingroup$ You would first need to rewrite your code to make the interaction terms explicit. I don't see any interactions in the code you wrote. $\endgroup$
    – Tom Minka
    Aug 13, 2020 at 8:30
  • $\begingroup$ In the R code, I generated a column for the X1:X3 interaction by multiplying X1 & X3 to generate a new column X5; similarly for X2:X4 interaction. So the matrix of predictors X contains the interactions, and I just estimate the coefficients in the stan code. This is probably not great, but a consequence of the fact that stan syntax is somewhat unfamiliar to me, so I'm learning how to do that at the same time. $\endgroup$
    – mkt
    Aug 13, 2020 at 8:38
3
+400
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thank you, this is a great answer. I don't fully understand the for loop in the code yet or how to tweak it for a multiple regression (especially with interaction terms), but I will work on figuring that out. If you don't mind, I have one follow-up question: if one doesn't have a good idea of what the measurement error is, how would one estimate it for these purposes? That seems crucial. $\endgroup$
    – mkt
    Aug 13, 2020 at 15:44
  • $\begingroup$ I guess one way of doing this is to assign a prior to the measurement error, but I suspect it would be very poorly identified by the model (or not at all) - which may not be an issue, if the prior is meaningful (e.g. expert eliciation). With a lot of things I'm used to (say, blood pressure measurements in patients), we have at least a decent idea of how precise the devices used measure the things we are interested in. (Digression: A lot of modeling is generally ignoring measurement uncertainty in covariates, which I believe is an issue unless you have a randomized study...) $\endgroup$
    – Björn
    Aug 16, 2020 at 6:15
3
$\begingroup$

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.

$\endgroup$
1
  • 2
    $\begingroup$ Can you please edit to include the page of McElreath's Statistical Rethinking that discusses the hurdle model for chemical analyses? Or perhaps a relevant quotation from that resource (if it's brief enough for this format)? $\endgroup$
    – Sycorax
    Oct 17 at 14:50
0
$\begingroup$

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.

$\endgroup$
1
  • 1
    $\begingroup$ I wonder how well this approach would work: pretend the left-censored predictor is not censored but has a discontinuity in its relationship with Y at the censoring point, and have a second Y variance parameter for those observations to reflect incomplete conditioning on X (X < c instead of X = x). This would not force you to be right in selecting a distribution for X. $\endgroup$ Aug 16, 2020 at 11:21

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.