CODEWITHSUNDEEP
rcurve-fitting
user3036416
user3036416

Reputation: 1255

Error in nls singular gradient R

I have a data frame like this 

 time     WT   WTIC ATIC   RHIC esaIC          k  uIC
1 0.00 25.191 25.191 21.4 67.925 25.49 0.06152572 3.53
2 0.05 25.186 25.191 21.4 67.925 25.49 0.06152572 3.53
3 0.10 25.179 25.191 21.4 67.925 25.49 0.06152572 3.53
4 0.15 25.168 25.191 21.4 67.925 25.49 0.06152572 3.53
5 0.20 25.158 25.191 21.4 67.925 25.49 0.06152572 3.53
6 0.25 25.147 25.191 21.4 67.925 25.49 0.06152572 3.53

which I would like to fit using this non linear function

f <- function(x,a1,a2,a3,a4,a5,a6,par1,par2,par3) {  
      Tinf <- a2 - (par2*(1-a3/100)*a4)/(1+par2*a5*a4)
      kC <-par1*sqrt(a6)
      V <- kC + par3
      tau <- 1/(V*(1+par2*a5*a4))

      func <- a1 -(a1-Tinf)*(1-exp(-x/tau))
      return(func)
}

However, when using nls

nls(WT~f(time,WTIC,ATIC,RHIC,esaIC,k,uIC,par1,par2,par3), data=df, start=c(par1=1, par2=1,par3=1))

I get this error

Error in nlsModel(formula, mf, start, wts) : 
  singular gradient matrix at initial parameter estimates

I tried to change the starting values of the parameters but I still get the same error. Any help?

Upvotes: 1

Views: 1379

Answers (1)

mra68
mra68

Reputation: 2960

This "Singular Gradient" error is a drawback of the "nls" function. "optim" behaves better:

df <- read.table( text = 
 ' time     WT   WTIC ATIC   RHIC esaIC          k  uIC
  "0.00" "25.191" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"
  "0.05" "25.186" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"
  "0.10" "25.179" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"
  "0.15" "25.168" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"
  "0.20" "25.158" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"
  "0.25" "25.147" "25.191" "21.4" "67.925" "25.49" "0.06152572" "3.53"',
   header = TRUE )


f <- function(x,a1,a2,a3,a4,a5,a6,par1,par2,par3) {  
  Tinf <- a2 - (par2*(1-a3/100)*a4)/(1+par2*a5*a4)
  kC <-par1*sqrt(a6)
  V <- kC + par3
  tau <- 1/(V*(1+par2*a5*a4))

  func <- a1 -(a1-Tinf)*(1-exp(-x/tau))
  return(func)
}

#----------------------------------------------------------------------
# Essentially the same as f:

g <- function(v){f(v[1],v[2],v[3],v[4],v[5],v[6],v[7],v[8],v[9],v[10])}

#----------------------------------------------------------------------
# The function we want to minimize:

squaredError <- function(par)
{
  sum((df$"WT"-apply(cbind(df[,-2],par[1],par[2],par[3]),1,g))^2)
}

#----------------------------------------------------------------------
# Optimization of the parameters:

opt <- optim( par    = c(1,1,1),
              fn     = squaredError,
              method = "BFGS" )

#----------------------------------------------------------------
# Result:

opt

squaredError(opt$par + c( 1, 0, 0)*1e-3 )
squaredError(opt$par + c(-1, 0, 0)*1e-3 )
squaredError(opt$par + c( 0, 1, 0)*1e-3 )
squaredError(opt$par + c( 0,-1, 0)*1e-3 )
squaredError(opt$par + c( 0, 0, 1)*1e-3 )
squaredError(opt$par + c( 0, 0,-1)*1e-3 )

.

> opt
$par
[1] -0.04261273 -0.23600921  0.44504195

$value
[1] 4.572781e-05

$counts
function gradient 
     137      100 

$convergence
[1] 1

$message
NULL


> squaredError(opt$par + c( 1, 0, 0)*1e-3 )
[1] 4.581051e-05

> squaredError(opt$par + c(-1, 0, 0)*1e-3 )
[1] 4.583096e-05

> squaredError(opt$par + c( 0, 1, 0)*1e-3 )
[1] 4.900303e-05

> squaredError(opt$par + c( 0,-1, 0)*1e-3 )
[1] 4.939846e-05

> squaredError(opt$par + c( 0, 0, 1)*1e-3 )
[1] 4.57487e-05

> squaredError(opt$par + c( 0, 0,-1)*1e-3 )
[1] 4.575957e-05
>

Upvotes: 4

Related Questions