I am writing a Python script to automate a 5PL Standard Curve for ELISA data, using an XLSX file template. I am using pandas dataframes to hold the 96 wells as different rows, then in a separate dataframe, holding Standard Curve data (not all wells need a column for known concentration). I am also using scipy and numpy for analysis, then using matplotlib.pyplot for plotting the curve (to eventually be output to the XLSX), though I am open to suggestions.

So far, the program correctly takes the input, generates the aforementioned dataframes, and makes initial guesses for the 5PL's parameters. I was attempting to modify the 4PL code given here to work with my standard curve, using all actual standard replicates to optimize the curve, and obviously using 5PL instead of 4PL. However, when I use scipy's least_squares function for parameter optimization, I receive the following errors/traceback:

C:\Users\Jake\Anaconda3\lib\site-packages\ipykernel_launcher.py:6: RuntimeWarning: overflow encountered in power

ValueError                                Traceback (most recent call last)
<ipython-input-6-70b754f01d56> in <module>
      5 # Fit equation using least squares optimization
----> 6 plsq = least_squares(fun = residuals, jac = 'cs', x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))

~\Anaconda3\lib\site-packages\scipy\optimize\_lsq\least_squares.py in least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, loss, f_scale, diff_step, tr_solver, tr_options, jac_sparsity, max_nfev, verbose, args, kwargs)
    794     if not in_bounds(x0, lb, ub):
--> 795         raise ValueError("`x0` is infeasible.")
    797     x_scale = check_x_scale(x_scale, x0)

ValueError: `x0` is infeasible.

I was initially receiving these errors with the aforementioned article's suggested scipy.least_sq function, but I assumed it was minimizing exponents to 0 and encountering an overflow, so I swapped to scipy.least_squares instead due to its bounds argument. Code in question is below. Any help would be greatly appreciated! Thank you so much! And if this question is better suited to stats.SE, just let me know!

(I did find lmfit-py, though I was hoping to complete this project on my own if possible to get a better understanding of the regression. Still open to using it instead if that seems to be a better option.)

import matplotlib.pyplot as plt
from scipy.optimize import least_squares

def logistic5(x, A, B, C, D, E):
    """5PL logistic equation"""
    return D + ((A-D)/(np.power((1 + np.power((x / C), B)), E)))

def residuals(p, y, x):
    """Deviations of data from fitted 5PL Curve"""
    A, B, C, D, E = p
    err = y - logistic5(x, A, B, C, D, E)
    return err

def peval(x, p):
    """Evaluated value at x with current parameters"""
    A, B, C, D, E = p
    return logistic5(x, A, B, C, D, E)

x = []
y_meas = []
for ind, row in std_df.iterrows(): # Iterate through rows in std dataframe.
    for rep in std_df_rep_cols: # Iterate through replicate columns.
        if not pd.isnull(row[rep]): # If there is a value in the cell...
            x.append(row['conc']) # Add an iteration of the respective conc. to the x list.
            y_meas.append(row[rep]) # Add the cell's value to the y_meas list.

x = np.asarray(x)
y_meas = np.asarray(y_meas)


# Initial approximations of parameters
A = np.amax([np.amin(minBlk_ser)] + [0]) # Min asymptote
D = np.amax(y_meas) # Max asymptote
B = (D - A) / np.amax(x) # Steepness
C = np.amax(x) / 2 # Inflection Point, conc at which y = (D - A) / 2
E = 1 # Asymmetry factor
print(A, B, C, D, E)

p0 = [A, B, C, D, E] # List containing initial guesses as arg for least_squares
y_true = logistic5(x, A, B, C, D, E) # Assumed "true" curve based on initial params

# Fit equation using least squares optimization
plsq = least_squares(fun = residuals, x0 = p0, bounds = (0.0001, 100000), args=(y_meas, x))
  • $\begingroup$ Very interesting code. I can't see where you load import pandas or the data frame, but from memory iterrows and in particular pd is obviously a pandas command. The isnull is a bit hairy because most blanks will display an OD. Anyway with regards your bug .. there is a standard package used to input data for regression analysis, it is very famous (I would need to look it up) ... this would obviate using a list, which is what you are doing to define X0. $\endgroup$
    – M__
    Dec 5, 2019 at 12:13
  • $\begingroup$ Hi, @MichaelG.; thank you for the answer! I didn't include the pandas import/dataframe initialization in this question as they are earlier in the code. I can show them if it could help? The "isnull" is used not so much for blanks (measured and subtracted earlier in the code), but for missing replicates in the standard dataframe. It was a proactive choice in case, for instance, someone loaded wells with extra replicates for some (not all) standard concentrations and thereby left empty cells in the respective replicate column in std_df. Very nice catch! Is lmfit the package you're thinking of? $\endgroup$ Dec 5, 2019 at 22:39


Your Answer

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

Browse other questions tagged or ask your own question.