Единичная ошибка градиента во время начальной загрузки nls соответствует неверным данным

У меня есть набор данных, содержащий независимую переменную и набор зависимых переменных. Я хотел бы подогнать функцию к каждому набору независимых переменных, используя загрузочную нелинейную процедуру наименьших квадратов. В некоторых случаях независимые переменные имеют «хорошее качество», т. е. достаточно хорошо соответствуют функции. В других случаях они шумные.

Во всех случаях я могу использовать 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)

NLS соответствует данным mmFormula

Начальная загрузка отлично работает для данных хорошего качества:

# 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 для одновременного выполнения большого количества симуляций начальной загрузки?


r statistics-bootstrap nonlinear-functions
person Drew Steen    schedule 23.10.2012    source источник
comment
Я бы не стал подгонять модель Михаэлиса-Ментен только к трем различным значениям концентрации. Однако, возможно, вы могли бы улучшить предположение для начальных значений (в частности, KmGuess), сначала подобрав Lineweaver-Burk, используя lm.   -  person Roland    schedule 23.10.2012
comment
Да, я понимаю, что экспериментальная схема была далеко не оптимальной. Живи и учись. Использование Lineweaver-Burke в качестве начального предположения является хорошей идеей. Однако я не думаю, что начальное предположение является проблемой, потому что а.) nls подходит (без начальной загрузки) отлично работает с относительно плохими начальными предположениями, например. Км=100, Vмакс=0,5; б.) когда я изменяю функцию начальной загрузки на те же самые начальные предположения, я получаю ту же ошибку, и в.) я думаю, что плохие начальные предположения обычно вызывают ошибку сбоя сходимости, а не сингулярную ошибку градиента.   -  person Drew Steen    schedule 23.10.2012
comment
Ну, у вас есть некоторые данные, которые вообще не соответствуют модели. Иногда мне удавалось решать похожие проблемы (даже единичные ошибки градиента), используя разные начальные значения (с этим может помочь nls2). Также может помочь другой алгоритм оптимизации. Но если данные сильно нарушают модель, подгонка невозможна, что может произойти во время начальной загрузки.   -  person Roland    schedule 23.10.2012
comment
Но вот чего я не понимаю - все данные могут соответствовать модели. Это только повторно выбранные данные, которые не могут соответствовать модели.   -  person Drew Steen    schedule 23.10.2012
comment
Некоторые образцы начальной загрузки, вероятно, невозможно подогнать. Вы можете переписать свой бутстрап, чтобы исключить их. Если число исключенных выборок невелико, результат все еще может быть верным. Но мой совет не делать бутстрап на этих данных.   -  person Roland    schedule 23.10.2012
comment
Справедливо. Причина, по которой я хотел бы сделать их все, заключается в том, что я хочу впоследствии разработать критерий, который можно было бы отвергнуть. Также меня беспокоит, что я не понимаю ошибку.   -  person Drew Steen    schedule 23.10.2012
comment
давайте продолжим это обсуждение в чате   -  person Drew Steen    schedule 23.10.2012
comment
Возможно, вы могли бы загрузить остатки, чтобы лучше сохранить распределение x?   -  person Ben Bolker    schedule 24.10.2012

Ответы (1)


arrow_upward
3
arrow_downward

Это позволяет вам изучить, что происходит:

#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
person Roland    schedule 24.10.2012
comment
PS: избегайте использования subset внутри функций. Его справочная страница специально предостерегает от этого. Вместо этого используйте [. - person Roland; 24.10.2012