8

独立変数と従属変数のセットを含むデータセットがあります。ブートストラップされた非線形最小二乗法を使用して、独立変数の各セットに関数を適合させたいと思います。場合によっては、独立変数は「高品質」です。つまり、関数に適度に適合します。他の場合、彼らは騒々しいです。

いずれの場合も、nls()パラメーターの推定値を取得するために使用できます。ただし、データにノイズが多い場合、ブートストラップはエラーをスローしますError in nls(...) : singular gradient。ノイズの多いデータへのフィッティングが失敗する理由は理解でき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)

mmFormulaのデータへのNLS適合

ブートストラップは、高品質のデータに対して正常に機能します。

# 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そして、多くのブートストラップシミュレーションを同時に実行するために使用したいので、どうすればよいですか？

4

1 に答える 1

3

これにより、何が起こるかを調べることができます。

#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))

そのエラーメッセージでもう少し有益になります。たとえば、1つのエラーは

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

これは、一部のブートストラップサンプルに含まれる濃度が3つ未満であることを意味します。しかし、他にもいくつかのエラーがあります。

step factor 0.000488281 reduced below 'minFactor' of 0.000976562

を減らすことで回避できる可能性がありますminFactor。

以下は厄介です。次の方法で、さまざまなフィッティングアルゴリズムまたは開始値を試すことができます。

singular gradient matrix at initial parameter estimates

singular gradient
于 2012-10-24T09:57:08.043 に答える