2
$\begingroup$

I'm writing a program to estimate the lower limit of detection for a nucleic acid assay. A typical analysis will have 5 concentrations with 10 independent replicates each, the dependent variable being whether or not the nucleic acid was detected. I calculate the probability of detection from the replicates (# positive / total, per concentration), then convert these to a probit score. Most of that is just background and not super important to the question.

I end up estimating the concentration at which the assay can detect nucleic acid 95% of the time by doing a linear regression on the results (x=concentration, y=probit score), and solving the linear equation for y=6.64... (the probit score corresponding to 95%).

How would I get the confidence interval for this extrapolated value? I already have the standard error for the slope and the intercept of the regression, and can compute the 95% CI for both, but I'm unsure how to report the CI for the estimated value/result.

Since it can be affected by both the slope and the intercept, do I just re-calculate the result using the steepest slope ($\hat{m} + \sigma_{\hat{m}}t$, where $\sigma_{\hat{m}}$ is the standard error for the slope and $t$ is the appropriate t-value for the desired confidence interval and degrees of freedom) and lowest intercept ($\hat{b}-\sigma_{\hat{b}}t$), set that as the upper bound of the confidence interval, then recalculate using the shallowest slope and highest intercept ($\hat{m}-\sigma_{\hat{m}}t$ and $\hat{b}+\sigma_{\hat{b}}t$)?

If it's relevant, the program is in Python (3.x) and I'm using scipy.stats for the regression (linregress()) and for finding the width of the confidence intervals (t.ins()). But if someone could just explain the math, I can write the code (I haven't found an obvious answer in the scipy.stats docs).

Cheers!

$\endgroup$

2 Answers 2

1
$\begingroup$

doing a linear regression on the results (x=concentration, y=probit score), and solving the linear equation for y=6.64... (the probit score corresponding to 95%)

probably isn't the best way to accomplish what you want.

A linear regression of probit score against concentration makes an implicit assumption that the variance in probit score values around their linearly predicted values is constant as a function of concentration. That's unlikely to be the case and might substantially bias the estimate of detection limit.

I think you can better accomplish what you want by doing a proper probit regression that takes the underlying binomial variance into account. That can be further enhanced by allowing for a nonlinear relationship with concentration. You can then plot predicted probabilities along with confidence intervals as a function of concentration. This page shows how to do probit regressions in several platforms including Python.

The Python example doesn't show plots, but plots (and the way to generate them) are illustrated in the R example. You use a "predict" function that returns means and confidence intervals based on the model in the probit scale over a range of concentrations, rescale to probability, and plot.

Some might prefer a logistic regression to the probit regression, but the differences are usually pretty minor.

Finally, the way that you are proceeding doesn't seem to be a standard way of setting detection limits. It's not clear from your description how you are taking into account the variance in blank values. You might want to see this article for guidance.

$\endgroup$
4
  • $\begingroup$ I'll have to give this a look. I know linear regression is the poor man's probit, but when I was asked to make an LLoD calculator using the probit method, the resources I found used this method. Seemed specious to me, but as I'm no statistician, I just went with it. Having plots shown in R is actually better; I like the ggplot package much better than matplotlib, and have a python wrapper for it. Re: logit vs. probit -- after researching both, I got the impression they were close enough that the difference would be smaller than the uncertainty in either of the models anyway $\endgroup$
    – W. MacTurk
    Sep 29, 2021 at 19:04
  • 1
    $\begingroup$ @W.MacTurk I've added a link to a more formal summary of how to define limit of detection, as it's not clear how you are taking variability in blanks into account--unless miraculously your blanks are always exactly 0. That's not been my experience. $\endgroup$
    – EdM
    Sep 29, 2021 at 19:14
  • $\begingroup$ I tried to edit my original response but went over the time limit. I'm sure the blank values are non-zero, but the report spit out by the machine marks anything under whatever its programmed LoB is as 0/not detected, so I have no real option to account for LoB variability. Incidentally, the site you linked first was one of the first I came across, but I found it difficult to follow and overly complex given I only have a single regressor. It may not be complex at all, but again, I'm no statistician, and a lot of the examples were over my head. $\endgroup$
    – W. MacTurk
    Sep 29, 2021 at 19:19
  • 1
    $\begingroup$ @W.MacTurk for linear regression that meets the necessary assumptions, see this Wikipedia entry. The sections on "Mean Response" and "Predicted Response" show how their variances change as quadratic functions of predictor value. The section on "Confidence Intervals" shows that you multiply the square root of variance by the appropriate t-value to get CI around the mean. Thus you can calculate 95% CI along the range of concentrations and back-calculate the concentrations that hit those CI at any given response value. $\endgroup$
    – EdM
    Sep 29, 2021 at 20:28
0
$\begingroup$

The scipy stats docs include a method for this:

https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html

Copied from there (but I would never use a lambda for this, personally):

# Two-sided inverse Students t-distribution
# p - probability, df - degrees of freedom
from scipy.stats import t
tinv = lambda p, df: abs(t.ppf(p/2, df))
ts = tinv(0.05, len(x)-2)
print(f"slope (95%): {res.slope:.6f} +/- {ts*res.stderr:.6f}")
slope (95%): 1.453392 +/- 0.743465
print(f"intercept (95%): {res.intercept:.6f}"
  f" +/- {ts*res.intercept_stderr:.6f}")
intercept (95%): 0.616950 +/- 0.544475

elsewhere, I find a different solution and created confidence intervals this way:

(probe_data is a dataframe and phenotype is a list of numbers (e.g. 0=control, 1=treatment group) to categorize the columns in the data. scipy returns a results dataframe with some stats, but not the confidence intervals.

alpha is typically 0.05, and results.rvalue is the regression coefficient and results.stderr is the standard error for the regression.

results = scipy.stats.linregress(phenotypes, probe_data)

# add in confidence intervals
r_z_value = np.arctanh(results.rvalue)
z_score = scipy.stats.norm.ppf(1-alpha/2)
# t_score = scipy.stats.t.ppf(1-alpha/2, n_samples)
lo_z = r_z_value - (z_score * results.stderr)
hi_z = r_z_value + (z_score * results.stderr)
ci_lower, ci_upper = np.tanh((lo_z, hi_z))

is that closer to what you are looking for?

$\endgroup$
2
  • 1
    $\begingroup$ i think OP wants to have a regression band around the estimated line, and not only confidence intervals around individual parameters $\endgroup$
    – rep_ho
    Sep 28, 2021 at 18:43
  • $\begingroup$ Yes, rep_ho is correct - I want the user to see the actual 95%CI band (and ideally have the width at a certain value listed directly) around the regression. I had already read the docs and implemented Marc's solution for want of anything better, but this IS an application to perform calculations -- i don't want the user to have to do more calculating afterwards;) $\endgroup$
    – W. MacTurk
    Sep 28, 2021 at 22:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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