How get plot from nls in R?

Question

In R I use nls to do a nonlinear least-squares fit. How then do I plot the model function using the values of the coefficients that the fit provided?

(Yes, this is a very naive question from an R relative newbie.)

Answer 1

To give some context:

there are two main approaches to get the predicted values to plot. Choose either one depending on what you want to have control over:

• specifying a interpolation vector of x values (tt) and predicting the y values using the model with the predict function
• pulling the coefficients from the model and running the right side of our model expression on them

# Data:
x = c(0,20,40,60,80,100,120,140) # sample timepoints
C = c(147.8,78.3,44.7,29.5,15.2,7.8,3.2,3.9) # concentratios
data <- data.frame(C,x) # data we got 

# plot empirical data:
opar <- par(las = 1) # plot params to rotate ylabels 
plot(C~x, data = data, col = 2, pch=3) 

# modeling
mod <- formula(C ~ C0*exp(-k*x))
fit <- nls(mod, data, start=list(C0=100,k=0.1)) # fit model to data
# plot(profile(fit)) # profile iterations

## Plotting of the model:

# interpolation vector tt
t0 = min(data$x)
tmax = max(data$x)
niter = 301
tt <- seq(t0, tmax,length.out=niter)

#  using tt with predict to plot lines:
lines(tt, predict(fit, list(x = tt)),lty=4)

# using coefficents and curve (no tt vector nor niter needed) :
C0 <- coefficients(fit)["C0"] # extract coefficent initialC
k <-coefficients(fit)["k"] # extract coefficent k
curve(C0*exp(-k*x),t0,tmax,add = T,lty=4,col="red")

# or, use both tt and coefs with eval to get the lines:
# (this is ~more or less~ what happens under "the hood")
plot(tt,eval(C0*exp(-k*tt)),lty=4,col="blue",type = "l",add=T)

the last line is the most "bare-bones" R not using curve nor lines, and needs both the tt vector and coefficents C0 and k to be variables in the env.

Answer 2

I know what you want (I'm a Scientist). This isn't it, but at least shows how to use 'curve' to plot your fitting function over any range, and the curve will be smooth. Using the same data set as above:

nonlinFit <- nls(density ~ a - b*exp(-c*conc), data = DNase1, start = list(a=1, b=1, c=1) )

fitFnc <- function(x) predict(nonlinFit, list(conc=x))

curve(fitFnc, from=.5, to=10)

or,

curve(fitFnc, from=8.2, to=8.4)

or,

curve(fitFnc, from=.1, to=50) # well outside the data range

or whatever (without setting up a sequence of evaluation points first).

I'm a rudimentary R programmer, so I don't know how to implement (elegantly) something like ReplaceAll ( /. ) in Mathematica that one would use to replace occurrences of the symbolic parameters in the model, with the fitted parameters. This first step works although it looks horrible:

myModel <- "a - b*exp(-c*conc)"

nonlinFit <- nls(as.formula(paste("density ~", myModel)), data = DNase1, start = list(a=1, b=1, c=1) )

It leaves you with a separate 'model' (as a character string), that you might be able to make use of with the fitted parameters ... cleanly (NOT digging out a, b, c) would simply use nonlinFit ... not sure how though.

Answer 3

The function "curve" will plot functions for you.

Answer 4

coef(x) returns the coefficients for regression results x.

model<-nls(y~a+b*x^k,my.data,list(a=0.,b=1.,k=1))
plot(y~x,my.data)
a<-coef(model)[1]
b<-coef(model)[2]
k<-coef(model)[3]
lines(x<-c(1:10),a+b*x^k,col='red')

For example.

Answer 5

Using the first example from ?nls and following the example I pointed you to line by line achieves the following:

#This is just our data frame
DNase1 <- subset(DNase, Run == 1)
DNase1$lconc <- log(DNase1$conc)
#Fit the model
fm1DNase1 <- nls(density ~ SSlogis(lconc, Asym, xmid, scal), DNase1)

#Plot the original points
# first argument is the x values, second is the y values
plot(DNase1$lconc,DNase1$density)

#This adds to the already created plot a line
# once again, first argument is x values, second is y values
lines(DNase1$lconc,predict(fm1DNase1))

The predict method for a nls argument is automatically returning the fitted y values. Alternatively, you add a step and do

yFitted <- predict(fm1DNase1)

and pass yFitted in the second argument to lines instead. The result looks like this:

Or if you want a "smooth" curve, what you do is to simply repeat this but evaluate the function at more points:

r <- range(DNase1$lconc)
xNew <- seq(r[1],r[2],length.out = 200)
yNew <- predict(fm1DNase1,list(lconc = xNew))

plot(DNase1$lconc,DNase1$density)
lines(xNew,yNew)