Ajuda com modelo não linear
Boa noite a todosEstou fazendo doutorado na UEL no curso de Ciência animal. Ainda sou novato no uso do R e também não tenho base estatística aprofundada em modelos não lineares.E preciso de um help, estou adaptando uma fórmula de modelo não linear do proc do SAS para o R, mas li o help, busquei dicas no site http://leg.ufpr.br/~paulojus/embrapa/Rembrapa/Rembrapase30.html#x32-20100030 e fui seguindo o modelo explicado por esse autor, mas tem algumas coisas que não entendi e não estou conseguindo rodar o modelo.Primeiramente explicarei resumidamente o que se trata o modelo. Esse modelo é para poder determinar cinco outras variáveis, através da produção de gás pela fermentação de bacteriais do líquido ruminal in vitro, simulando assim o rumen do animal. As variáveis que tenho em mãos são o volume de gás produzido e o tempo da leitura, a partir deles preciso determinar:VcnF = volume de gás produzido pelo carboidratos não fibrosos (mL de gás)kdcnf = taxa de degradação dos carboidratos não fibrosos (%/hora)VcF = volume de gás produzido pelo carboidratos fibrosos (mL de gás)kdcf = taxa de degradação dos carboidratos fibrosos (%/hora)L= tempo de latência (h), este é o tempo que as bactérias gastam para começar a degradar o substrato.Mas ao rodar o modelo abaixo aparece o seguinte erro: fit30=nls(Y~ (VcnF/(1+exp(2-4*kdcnf*I(L-Tempo)))+ VcF/I(1+exp(2-4*kdcf*(L-Tempo)))), + data=cra30, + start=list(VcnF = 10,kdcnf = 0.01,L =0.1,VcF = 50,kdcf =0.01, Tempo=1), + alpha= (0.05))Error in qr.qty(QR, resid) : NA/NaN/Inf in foreign function call (arg 5)In addition: Warning message:In Ops.factor(lhs, rhs) : - not meaningful for factors Aproveitando o ensejo gostaria de perguntar como faço para limitar os parâmetros pois no proc do SAS está assim:VcnF = 10 to 200 by 50 kdcnf = 0.01 to 0.8 by 0.05 L = 0.1 to 5 by 0.2 VcF = 50 to 150 by 30 kdcf = 0.01 to 0.8 by 0.05; bounds L >=0; Caso alguém tenha materiais relativos a modelos não linear para que possa aprofundar melhor, também contribuirá bastante. Agradeço a todos, e espero que alguém possa me ajudar! Eduardo Lucas Terra PeixotoDiscente de doutorado - Ciência AnimalUniversidade Estadual de Londrina Cel. (43) 9653-5574
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
==========================================================================
participantes (2)
-
Eduardo Lucas Terra Peixoto -
Walmes Zeviani