how to get a geom_smooth for an nls assymptote function using SSasympOff

Question

I have the following dataframe:

df1<- structure(list(Site = structure(c(1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L), .Label = c("ALT01"), class = "factor"), Nets = 1:18, Cumulative.spp = c(12L,13L, 15L, 17L, 17L, 17L, 17L, 19L, 19L, 19L, 19L, 20L, 22L, 22L, 22L, 22L, 22L, 22L)), .Names = c("Site", "Nets", "Cumulative.spp"), row.names = c(NA, 18L), class = "data.frame")

and I am trying to get a ggplot2 plot with the geom_smooth response of this function:

Model1<-nls(Cumulative.spp ~ SSasympOff(Nets, A, lrc, c0), data = df1)

Typically if I had a model like this:

Model2 <- lm(Cumulative.spp ~ I(log(Nets), data = df1)

I tried two methods

Method 1

I would do this:

 library(ggplot2)

 ggplot(df1, aes(x=Nets, y = Cumulative.spp)) + geom_point() + geom_smooth(method="lm", formula=y~log(x), fill="blue", fullrange=T)

but when I try to do the same with the assymptote it does not work:

ggplot(df1, aes(x=Nets, y = Cumulative.spp)) + geom_point() + geom_smooth(method="nls", formula=y~SSasympOff(x, A, lrc, c0), color="blue", fullrange=T)

but I got this error and this plot:

Warning message:
Computation failed in `stat_smooth()`:
$ operator is invalid for atomic vectors

Method2

I tried predicting over the original dataframe to get a confidence interval and using geom_line over the predicted values and geom_ribbon on the interval, but when I do

predict(Model1, df1, interval = "confidence")

but I do not get the confidence interval, only the predicted values

any help would be appreciated

Why not just add the predicted values using `geom_line` and then the interval using `geom_ribbon`? — bouncyball, Dec 28 '16 at 14:51
hi @bouncyball I was trying to do that, but I don't get the interval, I tried `predict(Model1, df1, interval = "confidence")` as what I found in the `predict.nls` documentation, but it does not give me a confidence interval I will add this to the things I tried in the question. — Derek Corcoran, Dec 28 '16 at 14:58
you could perhaps use a bootstrap percentile method then...didn't realize getting confidence intervals would be so difficult with `nls`. My mistake — bouncyball, Dec 28 '16 at 15:00
@bouncyball hopefully you could discuss which approach you find to be better with me. — Derek Corcoran, Dec 28 '16 at 15:40

score 1 · Accepted Answer · answered Dec 28 '16 at 15:26

I thought since I suggested a bootstrap method I might demonstrate. In this case we'll be boostrapping the residuals (see Wikipedia for more information). I'm not too familiar with using nls, so someone may come along with a (valid) theoretical objection.

B <- 2500 # bootstrap iterations, big number
pred_mat <- matrix(0, nrow = 18, ncol = B) # initialize matrix
# extract residuals and predictions from original model
resids <- residuals(Model1)
preds <- predict(Model1)
df1$Pred <- preds
for(i in 1:B){
    # bootstrapped dependent variable
    new_y <- preds + sample(resids, replace = TRUE)
    df1$NewY <- new_y
    # fit model
    Model_Boot <- nls(NewY ~ SSasympOff(Nets, A, lrc, c0), data = df1)
    # extract predictions
    pred_mat[,i] <- predict(Model_Boot)
}

# add 2.5% and 97.5% percentile intervals to df1
df1 <- cbind(df1, t(apply(pred_mat, 1, FUN = function(x) quantile(x, c(.025, .975)))))
# rename appropriately
names(df1)[6:7] <- c('lower','upper')

# make plot
ggplot(df1, aes(x = Nets))+
    geom_point(aes(y = Cumulative.spp))+
    geom_line(aes(y = Pred))+
    geom_ribbon(aes(ymin = lower, ymax = upper),
                alpha = .2, fill = 'blue')

score 0 · Answer 2 · answered Dec 28 '16 at 15:38

0

I found this way to do it, which I think is a bit more simple than bouncyballs one, but I will let all of you judge the better answer, although I will upvote bouncyball's answer.

I found this old post about it, but I couldn't find the as.lm function in nls2, the I found this link with the function, and I decided to use the as.lm.nls function

ggplot(df1, aes(x=Nets, y = Cumulative.spp)) + geom_point() + geom_line(y = predict(as.lm.nls(Model1), interval = "confidence")[,1]) + geom_ribbon(ymax = predict(as.lm.nls(Model1), interval = "confidence")[,3], ymin = predict(as.lm.nls(Model1), interval = "confidence")[,2], fill = "blue", alpha = 0.5)

and I got this result faster than using the bootstrap approach

answered Dec 28 '16 at 15:38

Derek Corcoran

3,930
2
25
54

Unfortunately, the `as.lm.nls` function to linearize an nls which was in nls2 0.1.3 seemed to generate a lot of support requests because it does not work in the way that some people expect. This isn't a deficiency in how `as.lm.nls` is implemented but a misunderstanding of how this whole subject works on the part of a lot of people so I figured it would be best to remove the function. – G. Grothendieck Dec 29 '16 at 21:58
@G.Grothendieck I loved the way it worked, it is a lot easier than bootstrapping, I took it out from the link above. Can you recommend a citation on the mathematical process used to linearize the NLS, I would like to learn more about that. – Derek Corcoran Dec 29 '16 at 23:57
There is a book by Doug Bates that discusses this. – G. Grothendieck Dec 30 '16 at 02:13
Thanks @G.Grothendieck – Derek Corcoran Dec 30 '16 at 05:38

how to get a geom_smooth for an nls assymptote function using SSasympOff

2 Answers2