Obtendo os valores iniciais corretos para um modelo nls em R

12

Estou tentando ajustar um modelo simples de lei de energia a um conjunto de dados da seguinte maneira:

mydf:

rev     weeks
17906.4 1
5303.72 2
2700.58 3
1696.77 4
947.53  5
362.03  6

O objetivo é passar a linha de força e usá-la para prever revvlaues para as próximas semanas. Um monte de pesquisas me levou à nlsfunção, que eu implementei da seguinte maneira.

newMod <- nls(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))
predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))

Enquanto isso funciona para um lmmodelo, recebo um singular gradienterro, que eu entendo tem a ver com meus valores iniciais ae b. Tentei valores diferentes, chegando a traçar isso no Excel, passar sozinho, obter uma equação e usar os valores da equação, mas ainda assim recebo o erro. Eu olhei para um monte de respostas como essa e tentei a segunda resposta (não conseguia entender a primeira), mas sem resultado.

Eu realmente poderia usar alguma ajuda aqui sobre como encontrar os valores iniciais corretos. Ou, alternativamente, que outra função eu posso usar em vez de nls.

Caso você queira recriar mydfcom facilidade:

mydf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6)) 
NeonBlueHair
fonte
1
Embora declarado em termos de R (realmente precisa ser declarado em algum idioma), como encontrar valores iniciais apropriados para um ajuste de modelo não linear é suficientemente estatístico para ser abordado aqui, IMO. Não é realmente um Q de programação, por exemplo.
gung - Restabelece Monica

Respostas:

13

Esse é um problema comum nos modelos de mínimos quadrados não lineares; se seus valores iniciais estiverem muito longe do ideal, o algoritmo poderá não convergir, mesmo que ele possa se comportar bem próximo do ótimo.

Se você começar por tomar registros de ambos os lados e ajustar um modelo linear, você obter estimativas de e como a inclinação e interceptação (9,947 e -2,011) (edit: esse é o log natural)blog(a)b

Se você usar aqueles para orientar os valores iniciais para e tudo parece funcionar bem:bab

 newMod <- nls(rev ~ a*weeks^b, data=mydf, start = list(a=exp(9.947),b=-2.011))
 predict(newMod, newdata = data.frame(weeks=c(1,2,3,4,5,6,7,8,9,10)))
 [1] 17919.2138  5280.7001  2584.0109  1556.1951  1050.1230   761.4947   580.3091   458.6027
 [9]   372.6231   309.4658
Glen_b -Reinstate Monica
fonte
Isso é extremamente útil, muito obrigado! Eu tenho uma pergunta sobre como você conseguiu o seu valor "a" aqui. Tentei executar lm (log10 (rev) ~ log10 (semanas)) e depois usar a função "summary" e, enquanto eu recebo o mesmo valor "b", meu valor "a" sai para 4.3201. O que você fez de diferente para chegar a = 9.947?
NeonBlueHair
Observe que eu costumava expvoltar para valores não registrados, o que é uma pista indicando que eu usei a logfunção simples . Desde que você seja consistente com o log e o anti-log que você usa, obterá a mesma resposta para o valor inicial. Então você pode fazer a base 10 e eu posso fazer a base e tudo é igual. e
Glen_b -Reinstala Monica 27/11
Ah, você está totalmente certo. Erro amador da minha parte. Manteve o pensamento da notação matemática, esperando que "log" significasse base de log 10 e "ln" para log natural. Obrigado pelo esclarecimento.
NeonBlueHair
1
Para muitos matemáticos (e muitos estatísticos), um "registro" sem adornos é o registro natural, da mesma forma que um argumento sem adornos para uma função de pecado está em radianos. [Convenções conflitantes podem causar confusão, infelizmente, mas quando comecei a usar R, por exemplo, não pensei duas vezes no uso da função de log, já que R e compartilhamos a mesma convenção.]
Glen_b -Reinstate Monica
4

Experimentar

 newMod <- nls(rev ~ a*weeks^b, data=mydf, startlist(a=17919.2127344,b=-1.76270557120))

Me pediram para expandir um pouco essa resposta. Esse problema é tão simples que fico surpreso que o nls falhe. O problema real, porém, é com toda a abordagem R e a filosofia do ajuste não linear do modelo. No mundo real, seria escalado x para ficar entre -1 e 1 e y e y para ficar entre 0 e 1 (y = ax ^ b). Isso provavelmente seria suficiente para fazer com que os nls convergissem. É claro que, como Glen aponta, você pode ajustar o modelo log-linear correspondente. Isso se baseia no fato de existir uma transformação simples que lineariza o modelo. Isso geralmente não é o caso. O problema com rotinas R, como nls, é que elas não oferecem suporte para reparametrizar o modelo. Nesse caso, a reparameterização é simples, apenas redimensione / atualize x e y. No entanto, tendo se ajustado ao modelo, o usuário terá parâmetros a e b diferentes dos originais. Embora seja simples calcular os originais a partir deles, a outra dificuldade é que, em geral, não é tão simples obter os desvios padrão estimados para essas estimativas de parâmetros. Isso é feito pelo método delta, que envolve o hessiano da probabilidade logarítmica e de alguns derivados. O software não-linear de estimativa de parâmetros deve fornecer esses cálculos automaticamente, para que a reparameterização do modelo seja facilmente suportada. Outra coisa que o software deve suportar é a noção de fases. Você pode pensar em ajustar primeiro o modelo com a versão de Glen como fase 1. O modelo "real" se encaixa no estágio 2. a outra dificuldade é que, em geral, não é tão simples obter os desvios padrão estimados para essas estimativas de parâmetros. Isso é feito pelo método delta, que envolve o hessiano da probabilidade logarítmica e de alguns derivados. O software não-linear de estimativa de parâmetros deve fornecer esses cálculos automaticamente, para que a reparameterização do modelo seja facilmente suportada. Outra coisa que o software deve suportar é a noção de fases. Você pode pensar em ajustar primeiro o modelo com a versão de Glen como fase 1. O modelo "real" se encaixa no estágio 2. a outra dificuldade é que, em geral, não é tão simples obter os desvios padrão estimados para essas estimativas de parâmetros. Isso é feito pelo método delta, que envolve o hessiano da probabilidade logarítmica e de alguns derivados. O software não-linear de estimativa de parâmetros deve fornecer esses cálculos automaticamente, para que a reparameterização do modelo seja facilmente suportada. Outra coisa que o software deve suportar é a noção de fases. Você pode pensar em ajustar primeiro o modelo com a versão de Glen como fase 1. O modelo "real" se encaixa no estágio 2. O software não-linear de estimativa de parâmetros deve fornecer esses cálculos automaticamente, para que a reparameterização do modelo seja facilmente suportada. Outra coisa que o software deve suportar é a noção de fases. Você pode pensar em ajustar primeiro o modelo com a versão de Glen como fase 1. O modelo "real" se encaixa no estágio 2. O software não-linear de estimativa de parâmetros deve fornecer esses cálculos automaticamente, para que a reparameterização do modelo seja facilmente suportada. Outra coisa que o software deve suportar é a noção de fases. Você pode pensar em ajustar primeiro o modelo com a versão de Glen como fase 1. O modelo "real" se encaixa no estágio 2.

Ajustei o seu modelo com o AD Model Builder, que suporta fases de maneira natural. Na primeira fase, apenas a foi estimada. Isso coloca seu modelo no campo. Na segunda fase, a e b são estimados para obter a solução. O AD Model Builder calcula automaticamente os desvios padrão para qualquer função dos parâmetros do modelo por meio do método delta, de modo a incentivar a reparametrização estável do modelo.

Dave Fournier
fonte
2

O algoritmo de Levenberg-Marquardt pode ajudar:

modeldf <- data.frame(rev=c(17906.4, 5303.72, 2700.58 ,1696.77 ,947.53 ,362.03), weeks=c(1,2,3,4,5,6))

require(minpack.lm)
fit <- nlsLM(rev ~ a*weeks^b, data=modeldf, start = list(a=1,b=1))

require(broom)
fit_data <- augment(fit)

plot(.fitted~rev, data=fit_data)
Etienne Low-Décarie
fonte
1

Na minha experiência, uma boa maneira de encontrar valores iniciais para parâmetros de modelos NLR é usar um algoritmo evolutivo. Em uma população inicial (100) de estimativas aleatórias (pais) em um espaço de pesquisa, escolha as 20 melhores (descendentes) e use-as para ajudar a definir uma pesquisa em uma população subsequente. Repita até a convergência. Não há necessidade de gradientes ou hessianos, apenas avaliações de SSE. Se você não é muito ganancioso, isso geralmente funciona. Os problemas que as pessoas costumam ter é que estão usando uma pesquisa local (Newton-Raphson) para realizar o trabalho de uma pesquisa global. Como sempre, é uma questão de usar a ferramenta correta para o trabalho em questão. Faz mais sentido usar uma pesquisa global do EA para encontrar valores iniciais para a pesquisa local de Newton e, em seguida, deixar isso diminuir ao mínimo. Mas, como em todas as coisas, o diabo está nos detalhes.

Adrian O'Connor
fonte