nls() in R using entire matrix

Question

I have data which I want to fit to the following equation using R:

Z(u,w)=z0*F(w)*[1-exp((-b*u)/F(w))]

where z0 and b are constants and F(w), w=0,...,9 is a decreasing step function that depends on w with F(0)=1 and u=1,...,50.

Z(u,w) is an observed set of data in the form of a 50x10 matrix (u=50,...,1 down the side of the rows and w=0,...,9 along the columns). For example as I haven't explained that great, Z(42,3) will be the element in the 9th row down and the 4th column along.

Using F(0)=1 I was able to get estimates of b and z0 using just the first column (ie w=0) with the code:

n0=nls(zuw~z0*(1-exp(-b*u)),start=list(z0=283,b=0.03),options(digits=10))

I then found F(w) for w=1,...,9 by going through each columns and using the vlaues of b and z0 I found.

However, I was wanting to find a way to estimate all the 12 parameters at once (b, z0 and the 10 values of F(w)) as b and z0 should be fitted to all the data, not just the first column.

Does anyone know of any way of doing this? All help would be greatly appreciated!

Thanks James

Answer 1

This may be a case where the formula interface of the nls(...) function works against you. As an alternative, you can use nls.lm(...) in the minpack.lm package to perform non-linear regression with a programmatically defined function. To demonstrate this, first we create an artificial dataset which follows your functional form by design, with random error added (error ~ N[0,1]).

u  <- 1:50
w  <- 0:9
z0 <- 100
b  <- 0.02
F  <- 10/(10+w^2)
# matrix containing data, in OP's format: rows are u, cols are w
m  <- do.call(cbind,lapply(w,function(w)
                       z0*F[w+1]*(1-exp(-b*u/F[w+1]))+rnorm(length(u),0,1)))

So now we have a matrix m, which is equivalent to your dataset. This matrix is in the so-called "wide" format - the response for different values of w is in different columns. We need it in "long" format: all responses in a single column, with a separate columns identifying u and w. We do this using melt(...) in the reshape2 package.

# prepend values of u
df.wide <- data.frame(u=u, m)
library(reshape2)
# reshape to long format: col1 = u, col2=w, col3=z
df   <- melt(df.wide,id="u",variable.name="w", value.name="z")
df$w <- as.numeric(substr(df$w,2,4))-1

Now we have a data frame df with columns u, w, and z. The nls.lm(...) function takes (at least) 4 arguments: par is a vector of initial estimates of the parameters of the fit, fn is a function that calculates the residuals at each step, observed is the dependent variable (z), and xx is a vector or matrix containing the independent variables (u, v).

Next we define a function, f(par, xx), where par is an 11 element vector. The first two elements contain estimates of z0 and b. The next 9 contain estimates of F(w), w=1:9. This is because you state that F(0) is known to be 1. xx is a matrix with two columns: the values for u and w respectively. f(par,xx) then calculates estimate of the response z for all values of u and w, for the given parameter estimates.

library(minpack.lm)
# model function
f <- function(pars, xx) {
  z0 <- pars[1]
  b  <- pars[2]
  F  <- c(1,pars[3:11])
  u  <- xx[,1]
  w  <- xx[,2]
  z  <- z0*F[w+1]*(1-exp(-b*u/F[w+1]))
  return(z)
}
# residual function
resids <- function(p, observed, xx) {observed - f(p,xx)}

Next we perform the regression using nls.lm(...), which uses a highly robust fitting algorithm (Levenberg-Marquardt). Consequently, we can set the par argument (containing the initial estimates of z0, b, and F) to all 1's, which is fairly distant from the values used in creating the dataset (the "actual" values). nls.lm(...) returns a list with several components (see the documentation). The par component contains the final estimates of the fit parameters.

# initial parameter estimates; all 1's
par.start <- c(z0=1, b=1, rep(1,9))   
# fit using Levenberg-Marquardt algorithm
nls.out <- nls.lm(par=par.start, 
                  fn = resids, observed = df$z, xx = df[,c("u","w")], 
                  control=nls.lm.control(maxiter=10000, ftol=1e-6, maxfev=1e6))
par.final <- nls.out$par
results <- rbind(predicted=c(par.final[1:2],1,par.final[3:11]),actual=c(z0,b,F))
print(results,digits=5)
#               z0        b                                                                          
# predicted 102.71 0.019337 1 0.90456 0.70788 0.51893 0.37804 0.27789 0.21204 0.16199 0.13131 0.10657
# actual    100.00 0.020000 1 0.90909 0.71429 0.52632 0.38462 0.28571 0.21739 0.16949 0.13514 0.10989    

So the regression has done an excellent job at recovering the "actual" parameter values. Finally, we plot the results using ggplot just to make sure this is all correct. I can't overwmphasize how important it is to plot the final results.

df$pred <- f(par.final,df[,c("u","w")])
library(ggplot2)
ggplot(df,aes(x=u, color=factor(w)))+
  geom_point(aes(y=z))+ geom_line(aes(y=pred))