Cross-posted from Bioinformatics SE here. The program is used for regression of a bioassay standard curve, but the question I'm asking is statistical in nature.
I am attempting to write a tool to be used by my lab to automatically generate 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>
4
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)
793
794 if not in_bounds(x0, lb, ub):
--> 795 raise ValueError("`x0` is infeasible.")
796
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.leastsq 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!
(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)
print(len(x))
print(len(y_meas))
print(x)
print(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))
