17 Ago
2011
17 Ago
'11
11:53
Eduardo,
O grande gargalo com modelos de regressão não linear é o conhecimento que o
usuário tem do modelo. Isso porque o usuário deve passar valores iniciais
para o procedimento iterativo. Quando melhor são os valores maiores são as
chances de sucesso. Isso você pode extrair de uma pesquisa semelhante a sua
com dados de mesma natureza ou extrair essa informação do gráfico quando os
parâmetros têm interpretação cartesiana.
Aqueles limites que você declarou do SAS, penso eu, que são limites para um
procedimento brute-force que vai varrer dentro daqueles intervalos o melhor
valor inicial para os parâmetros, que depois serão passados para otimização.
Isso é comum é modelos com muitos parâmetros, primeiro faz se a busca num
conjunto discreto de valores e depois os melhores valores são passados para
o otimizador. No R você pode fazer isso com a função nls2::nls2(...,
algorithm="brute-force"), e passar os melhores valos para a nls().
Para o caso de encontrar chutes iniciais, você pode usar os recursos
gráficos interativos do R. O RStudio tem a função manipulate::manipulate(),
o pacote gWidgetsRGtk2 também oferece formas de implementar e o
playwith::playwith() também. O procedimento é plotar uma curva sobre os
dados e mover os deslizadores que alteram valores dos parâmetros de forma da
curva. Vai movendo até encontrar a melhor combinação e usá-la para passar
para a nls(). Veja estes exemplos
http://ridiculas.wordpress.com/2011/04/09/metodo-grafico-interativo-para-valores-iniciais-em-regressao-nao-linear/
http://ridiculas.wordpress.com/2011/07/01/metodo-grafico-para-valores-iniciais-em-regressao-nao-linear-com-rstudio/
Aqui eu fiz um ajuste com dados simulados pelo modelo que você passou, com
valores dos parâmetros dentro daqueles limites.
> # gerar dados para ilustrar exemplo
> VcnF <- 100
> kdcnf <- 0.1
> L <- 2
> VcF <- 75
> kdcf <- 0.9
> Tempo <- seq(0,4,l=30)
> Y <- (VcnF/(1+exp(2-4*kdcnf*(L-Tempo)))+ VcF/I(1+exp(2-4*kdcf*(L-Tempo))))
> plot(Y~Tempo)
>
> # soma erro aleatório
> set.seed(17082011)
> y <- Y+rnorm(length(Tempo),0,2)
> plot(y~Tempo)
>
> # ajusta modelo
> n0 <- nls(y~(VcnF/(1+exp(2-4*kdcnf*(L-Tempo)))+
+ VcF/I(1+exp(2-4*kdcf*(L-Tempo)))),
+ start=list(VcnF=100, kdcnf=0.4, L=2, VcF=75, kdcf=0.5))
> summary(n0)
Formula: y ~ (VcnF/(1 + exp(2 - 4 * kdcnf * (L - Tempo))) + VcF/I(1 +
exp(2 - 4 * kdcf * (L - Tempo))))
Parameters:
Estimate Std. Error t value Pr(>|t|)
VcnF 75.50288 20.46526 3.689 0.0011 **
kdcnf 0.10316 0.04780 2.158 0.0407 *
L 2.07184 0.05375 38.548 < 2e-16 ***
VcF 79.33454 9.54395 8.313 1.16e-08 ***
kdcf 0.82652 0.09224 8.960 2.80e-09 ***
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 1.771 on 25 degrees of freedom
Number of iterations to convergence: 13
Achieved convergence tolerance: 5.083e-06
>
> # gráfico da curva estimada sobre os valores observados
> c0 <- do.call(data.frame, as.list(coef(n0)))
> with(c0,
+ {
+ curve((VcnF/(1+exp(2-4*kdcnf*(L-x)))+
VcF/I(1+exp(2-4*kdcf*(L-x)))),
+ add=TRUE, col=2)
+ })
>
À disposição.
Walmes.
==========================================================================
Walmes Marques Zeviani
LEG (Laboratório de Estatística e Geoinformação, 25.450418 S, 49.231759 W)
Departamento de Estatística - Universidade Federal do Paraná
fone: (+55) 41 3361 3573
VoIP: (3361 3600) 1053 1173
e-mail: walmes@ufpr.br
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================