r - 自举 nls 拟合不良数据期间出现奇异梯度错误

Question

我有一个包含一个自变量和一组因变量的数据集。我想使用自举非线性最小二乘法为每组自变量拟合一个函数。在某些情况下，自变量是“良好质量”的，即相当好地拟合函数。在其他情况下，它们很吵。

在所有情况下，我都可以使用 nls() 来估计参数。然而，当数据有噪声时， Bootstrap 会抛出错误 Error in nls(...) :奇异梯度。我可以理解为什么 nls 拟合噪声数据会失败，例如经过太多次迭代后未能收敛，但我不明白为什么它是一个奇异梯度误差，以及为什么我只能得到质量较差的重采样数据集。

代码:

require(ggplot2)
require(plyr)
require(boot)

# Data are in long form: columns are 'enzyme', 'x', and 'y'
enz <- read.table("http://dl.dropbox.com/s/ts3ruh91kpr47sj/SE.txt", header=TRUE)

# Nonlinear formula to fit to data
mmFormula <- formula(y ~ (x*Vmax) / (x + Km))

nls 完全能够拟合数据(即使在某些情况下，例如 a，我怀疑模型是否拟合数据。

# Use nls to fit mmFormula to the data - this works well enough
fitDf <- ddply(enz, .(enzyme), function(x) coefficients(nls(mmFormula, x, start=list(Km=100, Vmax=0.5))))

# Create points to plot for the simulated fits
xGrid <- 0:200
simFits <- dlply(fitDf, .(enzyme), function(x) data.frame(x=xGrid, y=(xGrid * x$Vmax)/(xGrid + x$Km)))
simFits <- ldply(simFits, identity) 

ggplot() + geom_point(data=enz, aes(x=x, y=y)) + geom_line(data=simFits, aes(x=x, y=y)) + 
  facet_wrap(~enzyme, scales="free_y") + aes(ymin=0)

引导对于高质量数据来说效果很好:

# Function to pass to bootstrap; returns coefficients of nls fit to formula
nlsCoef <- function(df, i) {
  KmGuess <- median(df$x)
  VmaxGuess <- max(df$y)
  dfSamp <- df[i,]
  nlsCoef <- coefficients(nls(mmFormula, dfSamp, start=list(Km=100, Vmax=0.5)))
}

eBoot <- boot(subset(enz, enzyme=="e"), nlsCoef, R=1000) #No error

但不适用于质量较差的数据

dBoot <- boot(subset(enz, enzyme=="d"), nlsCoef, R=10)
> Error in nls(mmFormula, dfSamp, start = list(Km = KmGuess, Vmax = VmaxGuess)) : 
   singular gradient

是什么导致了这个错误？鉴于我想使用 plyr 同时执行大量引导模拟，我该怎么办？

Answer 1

这可以让您检查发生了什么:

#modified function
#returns NAs if fit is not sucessfull
#does global assignment to store bootstrap permutations
nlsCoef <- function(df, i) {
  KmGuess <- median(df$x)
  VmaxGuess <- max(df$y)
  dfSamp <- df[i,]
  fit <- NULL
  try(fit <- nls(mmFormula, dfSamp, start=list(Km=100, Vmax=0.5)))
  if(!is.null(fit)){
    res <- coef(fit)
  } else{
    res <- c(Km=NA,Vmax=NA)
  }

  istore[k,] <<- i
  k <<- k+1
  res
}

n <- 100
istore <- matrix(nrow=n+1,ncol=9)
k <- 1

dat <- subset(enz, enzyme=="d")
dBoot <- boot(dat, nlsCoef, R=n) 

#permutations that create samples that cannot be fitted
nais <- istore[-1,][is.na(dBoot$t[,1]),]

#look at first bootstrap sample 
#that could not be fitted
samp <- dat[nais[1,],]
plot(y~x,data=samp)
fit <- nls(mmFormula, samp, start=list(Km=100, Vmax=0.5))
#error

您还可以使用自启动模型:

try(fit <- nls(y ~ SSmicmen(x, Vmax, Km), data = dfSamp))

这样，错误消息就会变得更加丰富。例如，一个错误是

too few distinct input values to fit a Michaelis-Menten model

这意味着，一些引导样本包含少于三个不同的浓度。但还有一些其他错误:

step factor 0.000488281 reduced below 'minFactor' of 0.000976562

您可以通过减小 minFactor 来避免这种情况。

以下内容令人讨厌。您可以尝试不同的拟合算法或起始值:

singular gradient matrix at initial parameter estimates

singular gradient