1
$\begingroup$

I have a bunch of curves that look like S-curves (stored in a db with Point format) and I am interested in finding this midpoint between both horizontal asymptotes for each of them. So, if X-axis represents time and Y-axis number of occurrences of a phenomenon, how can I project this midpoint over X-axis to know the date it happened?

I would like to do this operation in R software, please, can you help me with some guidelines about how to do it?

Thanks!

$\endgroup$

2 Answers 2

3
$\begingroup$

I'm pretty sure this is not a statistics question, so it will probably get moved. Also, it is quite vague, and I have no idea what 'stored in a db with Point format' means. I'll try to help you out though:

In theory you have a function $\pi(x)$ that is monotone increasing with range 0 to 1 and you want to find $x$ such that $\pi(x) = 0.5$. If you can compute the inverse of $\pi(x)$ analyticall then that $x$ would be $x= \pi^{-1}(0.5)$.

If you can't do it analytically you can always use a computer.

For example if $\pi(x) = \dfrac{e^{0.5+2x}}{1+e^{0.5+2x}}$ you can do this in R:

> x<- seq(-20,20,by=0.1)
> plot( exp(0.5+2*x)/(1+exp(0.5+2*x))~x)

enter image description here

> uniroot( function(x){exp(0.5+2*x)/(1+exp(0.5+2*x))}-0.5, c(-20,20))$root
[1] -0.2500002

enter image description here

$\endgroup$
3
  • 2
    $\begingroup$ I don't think there's any indication in the question that the S curves are between 0 and 1. In fact, since it says the y-values are "number of occurrences" I think it can't be the case $\endgroup$
    – Glen_b
    Commented Apr 11, 2014 at 15:25
  • 2
    $\begingroup$ Also, because the responses are counts, they likely have inherent variability. Thus more than an estimate of the midpoint (aka "LD50" in another setting) is needed: it would be good to derive a standard error of that estimate. This all suggests applying a poisson GLM, perhaps with logistic link, using time (or splines thereof) as explanatory variables, and then working out that middle time using methods of "inverse regression." Somehow it is starting to look like a statistics question :-). $\endgroup$
    – whuber
    Commented Apr 11, 2014 at 15:34
  • $\begingroup$ Well, in other words: my bunch of S curves are stored in a database. How I will retrieve and represent them is out of the question for simplicity. Yes, I need to find this 'midpoint' (maybe 'change from LOW to HIGH' is a better description) and Benjamin provided good ideas. As Glen_b said, values in my S curves are not bounded to [0, 1], but [0, infinity]. My concern is: if my S curves have the majority of points, say, on the left side, will 'uniroot' function still predict correctly the change LOW-to-HIGH? How this bias affect to the accuracy of the estimation? $\endgroup$
    – iamgin
    Commented Apr 16, 2014 at 9:04
2
$\begingroup$

This is what I did in the end, because it might be useful to others:

I used nls function to find a fit:

fit <- nls(Y ~ theta1/(1+exp(-(theta2 + theta3*X))), start=list(theta1=1, theta2=0.02, theta3=0.03), trace=TRUE)
p=predict(fit)

Theta1/2/3 correspond to the "upper limit" of the function, the "growth rate" and the "maximum growth rate". I guessed them taking a look to the data in the chart and trying. If you get an error such as "singular gradient" it is because your parameters are not correct for the function growing or maybe because your x/y values are too big for them and do not converge after several iterations. I think, not fully sure about what is going on.

To calculate the inverse of a function to be able to know the X value when Y is 0.5, I checked this post Solving for the inverse of a function in R. Like this I can do something like:

my_inverse = inverse(function(x) 0.74/(1+exp(4.14-0.47*x)))
res <- my_inverse(0.5)
res
$root
[1] 38.306

Where the coefficients are obtained from "fit" parameter. So, when Y=0.5, X=38.306. Maybe there is a easier method for the same, but I just could come up with this solution.

$\endgroup$

Your Answer

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

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