2

I have a dataset of 10 subjects, whose muscle response amplitude was measured across 11 different stimulus intensities (states) applied to their muscles. The output curve of the data looks like a sigmoidal function, with different stimulus intensities on the x-axis and resultant muscle amplitude outputs on the y-axis.

I would like to fit a Boltzmann sigmoid function to fit each subject's data. From the literature in my field, I would need to estimate my parameters of interest (i.e., slope, S50, and plateau from the function below) using a non-linear least squares fitting using R (typical once used is a Marquardt-Levenberg algorithm for least squares convergence).

From what I could find, I believe the Boltzmann function should look like:

y(x) = plateau / [1 + e^((S50-x)/slope)]

Plateau = maximum amplitude
S50 = stimulus intensity x required to evoke a response equal to half the maximum amplitude

Example data for one subject:

y = c(0.04, 0.053, 0.096, 0.219, 0.701, 1.333, 3.032, 3.556,4.33, 4.27, 4.167)
x = c(1,2,3,4,5,6,7,8,9,10,11)

Giving:

enter image description here

This post has someone outline how they do this with the "nls" package: Working with nls function produces a singular gradient

So I would have something like:

fit <- nls(y ~ plateau / (1 + exp((S50 - x)/slope), start = list( ??? ))

However, I'm still a little lost. How do I know what parameters to feed into the equation (which I think would go in the [???] section)? Do I need to specify plateau, S50, and slope? Isn't the whole point that I don't know these going in?

Is there a way to feed my data into some function to then get an output of the best non-linear equation that fits my data and gives me those S50, plateau, and slope parameters? Am I misunderstanding the process?

There is apparently a "gslnls" package available in R that uses the Marquardt-Levenberg algorithm to fit the data. I guess I am a little confused as to how to apply it, or if there is a simpler way to do this.

Any guidance would be greatly appreciated.

A.R.
  • 51
  • 5

1 Answers1

1

What goes inside the [???] are the initial guesses for the starting values for your three parameters. R and nls() does not know how your function behaves nor where are feasible solutions exist thus thus it needs some help on where to start looking. I recommend adding the trace=TRUE option to function to see where is it looking for the solution. Sometimes this takes some trial and error to get right!

Singular gradient error usually means one is entering an undefined region of the solution space, i.e. a value is going to infinity, square root of a negative number, the function is over defined with the number of adjustable parameters, etc.

Starting with the slope <1 seems to work.

y = c(0.04, 0.053, 0.096, 0.219, 0.701, 1.333, 3.032, 3.556,4.33, 4.27, 4.167)
x = c(1,2,3,4,5,6,7,8,9,10,11)
fit <- nls(y ~ plateau / (1 + exp((S50 - x)/slope)), start = list( plateau=1, S50=1, slope=.5), trace = TRUE)
Dave2e
  • 22,192
  • 18
  • 42
  • 50
  • Thanks! Just tried using ## fit2 <- nls(y ~ SSlogis(x, plateau, S50, slope)) ## and got the same outputs for plateau, slope, and S50 as with just the nls function with starting values. Is there a reason to use one over the other? SSlogis seems to take the guesswork out of it, but not sure if there might be scenarios where one works better than the other (at least, in the case of following a Boltzmann function). Also, do you happen to know if the nls function uses the Marquardt-Levenberg algorithm for least squares convergence? – A.R. Aug 04 '23 at 20:02
  • 1
    I assume since `nls()` is generic and can fit any model so it is just the go-to one people uses. I assume since there is this specific `SSlogis()` function to solve logistic equation, it is probably the better choice. I don't there is a concern of finding a local mininium so the starting point is not that critical. From the documentation, `nls()` uses the Gauss-Newton algorithm, though options to use others. – Dave2e Aug 04 '23 at 23:18