singular gradient matrix at initial parameter estimates

Question

I am fitting Lotka-Volterra model to some data and my simple code is :

    LCC<-c(11953.36,27413.50, 54823.75,  90904.73, 119250.34, 137679.14, 137750.79, 141107.43, 146356.56, 179280.00, 213037.07, 286138.03, 373380.14, 483474.16, 561950.42, 644919.07) 
AAR<-c(476068.4, 468621.0, 459684.1, 439587.9, 433429.0, 423889.7, 430239.4, 436614.3, 467925.9, 477742.2, 479049.2, 446933.5, 440711.2, 431425.5, 426563.1, 408660.5)
KAL<-c(888507, 887066, 875788, 872448, 855472, 817319, 787635, 793101, 800911, 787032, 763142, 733607, 706539, 673969, 656770, 672431)
PA<-c(457695.9,476068.4, 468621.0, 459684.1, 439587.9, 433429.0, 423889.7, 430239.4, 436614.3, 467925.9, 477742.2, 479049.2, 446933.5, 440711.2, 431425.5, 426563.1)
PK<-c(893515.0,888507, 887066, 875788, 872448, 855472, 817319, 787635, 793101, 800911, 787032, 763142, 733607, 706539, 673969, 656770)
PL<-c(11953.36,11953.36,27413.50, 54823.75,  90904.73, 119250.34, 137679.14, 137750.79, 141107.43, 146356.56, 179280.00, 213037.07, 286138.03, 373380.14, 483474.16, 561950.42)






LCCmodel<-nls(LCC~a1*PL/1+b1*PL+c1*PA+d1*PK, start=list(a1=-4.132e+05,b1=1.251e+00,c1=7.470e-01,d1=8.771e-02))

but it returns:

"Error in nls(LCC ~ a1 * PL/1 + b1 * PL + c1 * PA + d1 * PK, start = list(a1 = 122400,  :

singular gradient"

If i remove "PL" being multiplied with a1 from formula then it works. I can not skip it please help

Answer 1

There is a mistake in the formula. This can also be solved with DE:

xd=data.frame(PL,PA,PK)
library(NMOF)
yf= function(params,x){
  a1=params[1];b1=params[2];c1=params[3];d1=params[4]
  a1*x$PL/(1+b1*x$PL+c1*x$PA+d1*x$PK)
}

yf(rep(4,4),xd)

algo1 <- list(printBar = FALSE,
              nP  = 200L,
              nG  = 1000L,
              F   = 0.50,
              CR  = 0.99,
              min = rep(-5,4),
              max = rep(5,4))

OF1 <- function(Param, data) { 
  x <- data$x
  y <- data$y
  ye <- data$model(Param,x)
  aux <- y - ye; aux <- sum(aux^2)
  if (is.na(aux)) aux <- 1e10
  aux
}


data1 <- list(x = xd, y = LCC,  model = yf, ww = 1)
system.time(sol1 <- DEopt(OF = OF1, algo = algo1, data = data1))
sol1$xbest
OF1(sol1$xbest,data1)

plot(PL, LCC, pch=16, cex=0.8) 
lines(data1$x$PL,data1$model(sol1$xbest, data1$x),col=7,lwd=2)

#> sol1$xbest
#[1]  8.387138e-01  2.042065e-07 -2.004656e-06  7.055385e-07