I have used nls to fit a sigmoidal curve to my data to calculate the value of x when y=0.5 - the x values I get look correct. I would also like the gradient at this point, so calculate the derivative as: grad<-params["scal"]*exp((x-params["xmid"])/params["scal"])/(1+exp((x-params["xmid"])/params["scal"]))^2
My y values are between 0-1 but I am getting very steep gradient values, up to ~3.5 so is my formula for the derivative incorrect? Visual inspection suggests gradients should be ~0.05.
Full code for fitting model to multiple Species within for loop:
results<-data.frame()
#iterate through each Species
Sps<-unique(data$Species)
for (sp in Sps) {
#subset data for current Species
subset_data<-filter(data, Species==sp)
#fit sigmoidal logistic curve for current Species
model<-nls(Y~SSlogis(X, Asym, xmid, scal), data=subset_data)
#extract parameters for current Species
params<-coef(model)
#calculate x when y=0.5 for current Species
x_when_y_0.5<-params["xmid"]+params["scal"]*log(0.5/(params["Asym"]-0.5))
#calculate gradient (derivative) when y=0.5 for current Species
grad_when_y_0.5<-params["scal"]*exp((bbday-params["xmid"])/params["scal"])/(1+exp((bbday-params["xmid"])/params["scal"]))^2
#store results in dataframe
results<-rbind(results, data.frame(Species=sp, X=x_when_y_0.5, Grad=grad_when_y_0.5))
}
SSlogis
(according to its documentation) uses the expression
## Logistic curve
Asym/(1+exp((xmid-input)/scal))
Wolfram Alpha gives us the derivative with respect to input
as
## Derivative of logistic curve:
Asym * exp((xmid - input) / scal) / (scal * (exp((xmid - input) / scal) + 1)^2)
The midpoint of this curve is when xmid = input
, so xmid - input = 0
simplifies this expression nicely:
## Derivative at midpoint:
Asym * exp(0) / (scal * (exp(0) + 1)^2)
## Derivative at midpoint:
Asym / (4 * scal)
So I would suggest
grad_at_midpoint = params["Asym"] / 4 / params["scal"]
Do note that, unless the fitted value of the asymptote is 1, the midpoint will not correspond to y = 0.5
. The midpoint will be when y = params["Asym"] / 2
. The x-value for the midpoint will be params["xmid"]
.