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 ______________________________________________ R-help@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-help PLEASE do read the posting guide http://www.R-project.org/posting-guide.html and provide commented, minimal, self-contained, reproducible code.
