Plot a smooth and extrapolated curve using an nls model with several fitted parameters

Question

I feel that I am close to finding the answer for my problem, but somehow I just cannot manage to do it. I have used nls function to fit 3 parameters using a rather complicated function describing fertilization success of eggs (y-axis) in a range of sperm concentrations (x-axis) (Styan's model [1], [2]). Fitting the parameters works fine, but I cannot manage to plot a smoothed extrapolated curve using predict function (see at the end of this post). I guess it is because I have used a value that was not fitted on x-axis. My question is how to plot a smoothed and extrapolated curve based on a model fitted with nls function using non-fitted parameter on x-axis?

Here is an example:

library(ggplot2)

data.nls <- structure(list(S0 = c(0.23298, 2.32984, 23.2984, 232.98399, 2329.83993, 
23298.39926), fert = c(0.111111111111111, 0.386792452830189, 
0.158415841584158, 0.898648648648649, 0.616, 0.186440677966102
), speed = c(0.035161615379406, 0.035161615379406, 0.035161615379406, 
0.035161615379406, 0.035161615379406, 0.035161615379406), E0 = c(6.86219803476946, 
6.86219803476946, 6.86219803476946, 6.86219803476946, 6.86219803476946, 
7.05624476582978), tau = c(1800, 1800, 1800, 1800, 1800, 1800
), B0 = c(0.000102758645352932, 0.000102758645352932, 0.000102758645352932, 
0.000102758645352932, 0.000102758645352932, 0.000102758645352932
)), .Names = c("S0", "fert", "speed", "E0", "tau", "B0"), row.names = c(NA, 
6L), class = "data.frame")

## Model S

modelS <- function(Fe, tb, Be) with (data.nls,{
x <- Fe*(S0/E0)*(1-exp(-B0*E0*tau))
b <- Fe*(S0/E0)*(1-exp(-B0*E0*tb))
x*exp(-x)+Be*(1-exp(-x)-(x*exp(-x)))*exp(-b)})

## Define starting values

start <- list(Fe = 0.2, tb = 0.1, Be = 0.1)

## Fit the model using nls

modelS.fitted <- nls(formula = fert ~ modelS(Fe, tb, Be), data = data.nls, start = start, 
control=nls.control(warnOnly=TRUE,minFactor=1e-5),trace = T, lower = c(0,0,0), 
upper = c(1, Inf, 1), algorithm = "port")

## Combine model parameters

model.data <- cbind(data.nls, data.frame(pred = predict(modelS.fitted)))

## Plot

ggplot(model.data) +
geom_point(aes(x = S0, y = fert), size = 2) +
geom_line(aes(x = S0, y = pred), lwd = 1.3) +
scale_x_log10()

I have tried following joran's example here, but it has no effect, maybe because I did not fit S0:

r <- range(model.data$S0)
S0.ext <- seq(r[1],r[2],length.out = 200)
predict(modelS.fitted, newdata = list(S0 = S0.ext))
# [1] 0.002871585 0.028289057 0.244399948 0.806316161 0.705116868 0.147974213

Answer 1

You function should have the parameters (S0,E0,B0,tau,Fe,tb,Be). nls will look for the parameters in the data.frame passed to its data argument and only try to fit those it doesn't find there (provided that starting values are given). No need for this funny with business in your function. (with shouldn't be used inside functions anyway. It's meant for interactive usage.) In predict newdata must contain all variables, that is S0,E0,B0, and tau.

Try this:

modelS <- function(S0,E0,B0,tau,Fe, tb, Be) {
  x <- Fe*(S0/E0)*(1-exp(-B0*E0*tau))
  b <- Fe*(S0/E0)*(1-exp(-B0*E0*tb))
  x*exp(-x)+Be*(1-exp(-x)-(x*exp(-x)))*exp(-b)}

## Define starting values

start <- list(Fe = 0.2, tb = 0.1, Be = 0.1)

## Fit the model using nls

modelS.fitted <- nls(formula = fert ~ modelS(S0,E0,B0,tau,Fe, tb, Be), data = data.nls, start = start, 
                     control=nls.control(warnOnly=TRUE,minFactor=1e-5),trace = T, lower = c(0,0,0), 
                     upper = c(1, Inf, 1), algorithm = "port")

## Combine model parameters

model.data <- data.frame(
  S0=seq(min(data.nls$S0),max(data.nls$S0),length.out=1e5),
  E0=seq(min(data.nls$E0),max(data.nls$E0),length.out=1e5),
  B0=seq(min(data.nls$B0),max(data.nls$B0),length.out=1e5),
  tau=seq(min(data.nls$tau),max(data.nls$tau),length.out=1e5))
model.data$pred <- predict(modelS.fitted,newdata=model.data)


## Plot

ggplot(data.nls) +
  geom_point(aes(x = S0, y = fert), size = 2) +
  geom_line(data=model.data,aes(x = S0, y = pred), lwd = 1.3) +
  scale_x_log10()

Obviously, this might not be what you want, since the function has multiple variables and more than one vary in new.data. Normally one would only vary one and keep the others constant for such a plot.

So this might be more appropriate:

  S0 <- seq(min(data.nls$S0),max(data.nls$S0),length.out=1e4)
  E0 <- seq(1,20,length.out=20)
  B0 <- unique(data.nls$B0)
  tau <- unique(data.nls$tau)

model.data <- expand.grid(S0,E0,B0,tau)
names(model.data) <- c("S0","E0","B0","tau")

model.data$pred <- predict(modelS.fitted,newdata=model.data)



## Plot

ggplot(model.data) +
  geom_line(data=,aes(x = S0, y = pred, color=interaction(E0,B0,tau)), lwd = 1.3) +
  geom_point(data=data.nls,aes(x = S0, y = fert), size = 2) +
  scale_x_log10()