[R] Find position of asymptote

[Katrina Bennett]

Dear R Help,

Sorry I wasn't more clear before. Here is another crack at this.

What I am still trying to do is estimate the point on a line when the
slope changes or asymptotes. I have found some similar postings
talking about this but no answers.

https://stat.ethz.ch/pipermail/r-help/2003-January/028847.html for example
http://stackoverflow.com/questions/8245792/r-marking-slope-changes-in-loess-curve-using-ggplot2
(best one I see)

Based on the above, I've modified my code. I have provided below some
sample data so you can see the entire chain of work that summarizes
the issue I am trying to resolve.

dat <- c(75.44855206,NA,NA,NA,82.70745342,82.5335019,88.56617647,80.00128866,94.15418227,86.63987539,93.91052952,74.10612245,86.62289562,90.47961047,NA,NA,82.45320197,72.14371257,NA,71.44104803,72.59742896,68.36363636,NA,NA,61,NA,NA,71.26502909,NA,85.93333333,84.34248284,79.00522193,79.64223058,97.2074017,88.43700548,96.40413877,95.13511869,92.57379057,93.97498475,NA,97.55995131,89.53321146,97.21728545,93.21980198,77.54054054,95.85392575,86.25684723,97.55325624,80.03950617,NA,91.34023128,92.42906574,88.59433962,65.77272727,89.63772455,NA,NA,NA,NA,74.86344239,83.57594937,70.22516556,65.30543319,NA,NA,67.84852294,60.90909091,54.79303797,NA,52.18735363,33.47003155,NA,41.34693878,24.5047043,NA,NA,NA,NA,9.944444444,13.6875,NA,11.90267176,84.14285714,3.781456954,NA,1.432926829,4.26557377,1.823529412,0.444620253,4.711155378,NA,6.320284698,0.581632653,0.144578313,3.666666667,0,0,0,0,0,NA,0.032947462,0,0,10.54545455,0,NA,0.561007958,0.75,NA,0.048780488,0.74137931,NA,2.023339318,0,0,0,NA,NA,0.156950673,NA,0.283769634,32.81818182,NA,NA,0,NA,0,0,0,NA,0.212454212,3.120181406,NA,0.011811024,NA,0,0.120430108,5.928571429,1.75,0.679292929,0.97,NA,0,NA,NA,1,0.38547486,NA,1.460732984,0.007795889,0.05465288,0.004341534)

x.seq <- seq(1, 153,, 153)

dat.cl.2000 <- as.data.frame(cbind(x.seq, dat))

names(dat.cl.2000) <- c("x", "dat")Asym.2000 <-
mean(na.omit(dat.cl.2000$dat)[1:10])

na.dat.2000 <- which(!is.na(dat.cl.2000$dat))

loess.smooth.2000 <- loess.smooth(dat.cl.2000$x[na.dat.2000],
dat.cl.2000$dat[na.dat.2000], span=0.20)

diff.loess.smooth.2000 <- which(diff(loess.smooth.2000$y) ==
(min(diff(loess.smooth.2000$y))), arr.ind=T)

val.in.jday.2000 <- loess.smooth.2000$y[diff.loess.smooth.2000]

xmid.2000 <- val.in.jday.2000

scal.2000 <- (min(na.omit(dat.cl.2000$dat)) -
max(na.omit(dat.cl.2000$dat))) / (max(na.omit(dat.cl.2000$x)) -
min(na.omit(dat.cl.2000$x))) * 10

fit.dat.2000 <- nls(dat ~ SSlogis(x, Asym, xmid,scal), data =
dat.cl.2000, start=list(Asym=round(Asym.2000), xmid=round(xmid.2000),
scal=scal.2000), control = list(maxiter = 500, warnOnly = TRUE),
trace=TRUE)

year.ind <- 2000

fit.dat <- get(paste("fit.dat.", year.ind, sep=""))

a <- coef(fit.dat)[[1]]
m  <- coef(fit.dat)[[2]]
s  <- coef(fit.dat)[[3]]

sslogis <- expression(a/((1+exp((m-x)/s))))

sslogis.deriv1 <- D(sslogis, "x")

sslogis.deriv2 <- D(D(sslogis, "x"), "x")

fn.sslogis <- function (x, a, m, s) { a/((1+exp((m-x)/s))) }
fn.sslogis.deriv1 <- function (x, a, m, s) { a * (exp((m - x)/s) *
(1/s))/((1 + exp((m - x)/s)))^2 }
fn.sslogis.deriv2 <- function (x, a, m, s) {
-(a * (exp((m - x)/s) * (1/s) * (1/s))/((1 + exp((m - x)/s)))^2 -
  a * (exp((m - x)/s) * (1/s)) * (2 * (exp((m - x)/s) * (1/s) *
  (1 + exp((m - x)/s)))))/((((1 + exp((m - x)/s)))^2)^2) }

asym <- a

slope <- s

mid.point <- optimize(f=fn.sslogis.deriv1, interval=c(1:153), a=a,
m=m, s=s)$minimum #this will find the minimum point halfway through
melt

init <- optimize(f=fn.sslogis.deriv2, interval=c(1:mid.point), a=a,
m=m, s=s)$minimum #this is the start of the snowmelt season

term <- optimize(f=fn.sslogis.deriv2, interval=c(mid.point:153), a=a,
m=m, s=s, maximum=TRUE)$maximum #this is the end of the snowmelt
season

duration <- init-term


~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

Ok, this part of the code is where I really want to try and calculate
the points (init and term) where the curve asymptotes (using the first
derivative). On the below plot, this should be at ~40 and ~110 (these
represent days in my analysis).

plot.deriv1 <- fn.sslogis.deriv1(x.seq, a, m, s)
plot(plot.deriv1)

I found in the help pages above that I could do this based only on the
loess curve. I can take the diff of the loess prediction to calculate
a new function that has two local minimums.

loess.2000 <- loess(dat.cl.2000$x[na.dat.2000]~dat.cl.2000$dat[na.dat.2000])
loess.predict <- predict(loess.2000, newdata=x.seq)
loess.predict1 <- diff(loess.predict)

loess.predict2 <- c(NA,loess.predict1)
plot(x.seq, loess.predict2, main="diff(loess model)")

However, this still isn't providing me a means to get the points of
transition along this line.

Is there a better way to then pick off the change points or find the
asymptotes of a function in R?

Thank you.

Katrina