nls() singular gradient matrix at initial parameter estimates error I am trying to run a nls() to interpolate height growth of a tree given a series of samples, where each sample has a value of age, and height. I am using the Richards (1959) method:
$$
\textrm{Height}=x\left(1-e^{-y\cdot\textrm{Age}}\right)^z
$$
my age and height values:
> age
[1]  0  3  6 17 29 47 57 63 72
> height
[1]   37   70  199  464  736 1008 1280 1462 1643
> 

the script I am using to calculate parameters x, y, and z, is:
e<-2.71828
growCurve<-nls(
  formula=height~x*(1-e^(-1*y*age))^z,
  start=c(x=0,y=0,z=0)
)
b <- environment(growCurve[["m"]][["incr"]])[["internalPars"]][1]
y <- environment(growCurve[["m"]][["incr"]])[["internalPars"]][2]
z <- environment(growCurve[["m"]][["incr"]])[["internalPars"]][3]

b
y
z

when I run this however, I get the error

singular gradient matrix at initial parameter estimates

I tried using another fix from this post, but I still get identical results so I am not really sure on how to proceed.

The problem seems to be in convergence of the z parameter so fix it at each of a sequence of values and optimize on the others for each and then take the result with the lowest residual sum of squares.  We use the plinear algorithm of nls which does not require starting values for the parameters  that enter linearly -- x enters linearly so its starting value can be omitted.  Since we are fixing z and since x enters linearly we only have to specify a starting value for y. To use plinear the right hand side should be a column that multiplies x but x itself is not specified (in a similar way to lm).
age <- c(0, 3, 6, 17, 29, 47, 57, 63, 72)
height <- c(37, 70, 199, 464, 736, 1008, 1280, 1462, 1643)

fo <- "height ~ cbind(x = (1 - exp(-y * age))^z)"
zseq <- seq(0, 2, 0.01)
sse <- sapply(zseq, function(z) {
  res <- try(nls(fo, start = list(y = 1), algorithm = "plinear"))
  if (inherits(res, "try-error")) NA else deviance(res)
})

z <- zseq[which.min(sse)]
z # z value which minimizes the SSE
## [1] 1

fm <- nls(fo, start = list(y = 1), algorithm = "plinear")
fm

giving
Nonlinear regression model
  model: height ~ cbind(x = (1 - exp(-y * age))^z)
   data: parent.frame()
        y    .lin.x 
3.384e-03 7.468e+03 
 residual sum-of-squares: 18004

Number of iterations to convergence: 9 
Achieved convergence tolerance: 1.103e-06

Plotting the solution
plot(height ~ age)
lines(fitted(fm) ~ age)

(continued after chart)

Note that a simple linear regression actually gives a better fit (lower residual sum of squares) and has fewer parameters (2 vs 3).
fm.lin <- lm(height ~ age)
deviance(fm.lin)
## [1] 13546.89

deviance(fm)
## [1] 18003.86

The following is just to show that the deviance function in R is equivalent to taking the residual sum of squares for both nls and lm.
sum(resid(fm.lin)^2)
## [1] 13546.89

sum(resid(fm)^2)
## [1] 18003.86

