Ajuste de modelo não linear
Boa noite pessoal,o que eu poderia modificar neste ajuste para obter exito? Já variei os chutes iniciais e não consegui. Brasil<- c(13.63 ,13.93,14.09,15.12,15.98,16.08,16.45,17.96,18.83,19.28,19.18,20.89,21.87,21.85,23.43,25.38,25.64,26.99,28.68,31.92,32.92,33.32,34.37,34.99,36.08,37.72,40.41,42.81,44.18,44.91,48.97,51.59,53.44,56.49,62.2,66.49,67.01,68.41,74.14,77.58,80.39,85.71,89.58,96.3,100.9,102.8,107.69,111.51,111.42,116.59,118.85,121.42,127.53); anos <- c(1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013); fit<-nls(Brasil~a*(b^(anos)), start = list(a =18, b=1.4)) Error in nls(Brasil ~ a * (b^(anos)), start = list(a = 18, b = 1.4)) : número de iterações excedeu o máximo de 50 André Oliveira Souza. Graduação em Matemática, mestrado em estatística aplicada.Instituto Federal de Educação, Ciência e Tecnologia do Espirito Santo. IFES
O erro do R diz que tu excedeu o número máximo de iterações. Eu aumentei eles para 500 e, mesmo assim, não convergiu. Creio que o teu problema não tenha solução da maneira que foi proposto. Vai ser muito difícil encontrar uma função do tipo f(anos) = a*(b^anos) que funcione com um número tão grande para anos. Imagine uma função exponencial que cresce muito, muito, muito devagar até 1961 e, a partir daí, explode seu crescimento? Difícil, né? Para lidar com este problema, o que eu faria é o reduzir o valor dos anos, fazendo eles começarem em 1. Assim fica algo muito mais fácil de imaginar e interpretar: ##### Brasil <- c(13.63 ,13.93,14.09,15.12,15.98,16.08,16.45,17.96,18.83,19.28,19.18,20.89,21.87,21.85,23.43,25.38,25.64,26.99,28.68,31.92,32.92,33.32,34.37,34.99,36.08,37.72,40.41,42.81,44.18,44.91,48.97,51.59,53.44,56.49,62.2,66.49,67.01,68.41,74.14,77.58,80.39,85.71,89.58,96.3,100.9,102.8,107.69,111.51,111.42,116.59,118.85,121.42,127.53) anos <- c(1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013) anos2 <- anos-min(anos)+1 fit <- nls(Brasil ~ a*(b^(anos2)), start = list(a=12, b=1.4)) a <- coefficients(fit)[1] b <- coefficients(fit)[2] plot(Brasil ~ anos2) lines(a*(b^anos2), col="red", lty=2) ### Para encontrar os valores iniciais, eu tomei log(Brasil), ajustei uma reta a estes dados e tomei a exponencial dos valores encontrados. Aliás, linearizar a tua variável resposta é uma outra maneira de resolver este problema. Só tem que tomar cuidado com as inferências posteriores. 2016-10-15 22:56 GMT-03:00 Andre Oliveira via R-br <r-br@listas.c3sl.ufpr.br
:
Boa noite pessoal, o que eu poderia modificar neste ajuste para obter exito? Já variei os chutes iniciais e não consegui.
Brasil<- c(13.63 ,13.93,14.09,15.12,15.98,16. 08,16.45,17.96,18.83,19.28,19.18,20.89,21.87,21.85,23.43,25. 38,25.64,26.99,28.68,31.92,32.92,33.32,34.37,34.99,36.08,37. 72,40.41,42.81,44.18,44.91,48.97,51.59,53.44,56.49,62.2,66. 49,67.01,68.41,74.14,77.58,80.39,85.71,89.58,96.3,100.9,102. 8,107.69,111.51,111.42,116.59,118.85,121.42,127.53); anos <- c(1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971, 1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983, 1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995, 1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007, 2008,2009,2010,2011,2012,2013);
fit<-nls(Brasil~a*(b^(anos)), start = list(a =18, b=1.4))
*Error in nls(Brasil ~ a * (b^(anos)), start = list(a = 18, b = 1.4)) : número de iterações excedeu o máximo de 50*
André Oliveira Souza. Graduação em Matemática, mestrado em estatística aplicada.Instituto Federal de Educação, Ciência e Tecnologia do Espirito Santo. IFES
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br-guia) e forneça código mínimo reproduzível.
-- Marcus Nunes http://marcusnunes.me/
Obrigado Marcus! André Oliveira Souza. Graduação em Matemática, mestrado em estatística aplicada.Instituto Federal de Educação, Ciência e Tecnologia do Espirito Santo. IFES Em Domingo, 16 de Outubro de 2016 10:15, Marcus Nunes <marcus.nunes@gmail.com> escreveu: O erro do R diz que tu excedeu o número máximo de iterações. Eu aumentei eles para 500 e, mesmo assim, não convergiu. Creio que o teu problema não tenha solução da maneira que foi proposto. Vai ser muito difícil encontrar uma função do tipo f(anos) = a*(b^anos) que funcione com um número tão grande para anos. Imagine uma função exponencial que cresce muito, muito, muito devagar até 1961 e, a partir daí, explode seu crescimento? Difícil, né? Para lidar com este problema, o que eu faria é o reduzir o valor dos anos, fazendo eles começarem em 1. Assim fica algo muito mais fácil de imaginar e interpretar: #####Brasil <- c(13.63 ,13.93,14.09,15.12,15.98,16.08,16.45,17.96,18.83,19.28,19.18,20.89,21.87,21.85,23.43,25.38,25.64,26.99,28.68,31.92,32.92,33.32,34.37,34.99,36.08,37.72,40.41,42.81,44.18,44.91,48.97,51.59,53.44,56.49,62.2,66.49,67.01,68.41,74.14,77.58,80.39,85.71,89.58,96.3,100.9,102.8,107.69,111.51,111.42,116.59,118.85,121.42,127.53)anos <- c(1961,1962,1963,1964,1965,1966,1967,1968,1969,1970,1971,1972,1973,1974,1975,1976,1977,1978,1979,1980,1981,1982,1983,1984,1985,1986,1987,1988,1989,1990,1991,1992,1993,1994,1995,1996,1997,1998,1999,2000,2001,2002,2003,2004,2005,2006,2007,2008,2009,2010,2011,2012,2013) anos2 <- anos-min(anos)+1 fit <- nls(Brasil ~ a*(b^(anos2)), start = list(a=12, b=1.4)) a <- coefficients(fit)[1]b <- coefficients(fit)[2] plot(Brasil ~ anos2)lines(a*(b^anos2), col="red", lty=2)### Para encontrar os valores iniciais, eu tomei log(Brasil), ajustei uma reta a estes dados e tomei a exponencial dos valores encontrados. Aliás, linearizar a tua variável resposta é uma outra maneira de resolver este problema. Só tem que tomar cuidado com as inferências posteriores. 2016-10-15 22:56 GMT-03:00 Andre Oliveira via R-br <r-br@listas.c3sl.ufpr.br>: Boa noite pessoal,o que eu poderia modificar neste ajuste para obter exito? Já variei os chutes iniciais e não consegui. Brasil<- c(13.63 ,13.93,14.09,15.12,15.98,16. 08,16.45,17.96,18.83,19.28,19. 18,20.89,21.87,21.85,23.43,25. 38,25.64,26.99,28.68,31.92,32. 92,33.32,34.37,34.99,36.08,37. 72,40.41,42.81,44.18,44.91,48. 97,51.59,53.44,56.49,62.2,66. 49,67.01,68.41,74.14,77.58,80. 39,85.71,89.58,96.3,100.9,102. 8,107.69,111.51,111.42,116.59, 118.85,121.42,127.53); anos <- c(1961,1962,1963,1964,1965, 1966,1967,1968,1969,1970,1971, 1972,1973,1974,1975,1976,1977, 1978,1979,1980,1981,1982,1983, 1984,1985,1986,1987,1988,1989, 1990,1991,1992,1993,1994,1995, 1996,1997,1998,1999,2000,2001, 2002,2003,2004,2005,2006,2007, 2008,2009,2010,2011,2012,2013) ; fit<-nls(Brasil~a*(b^(anos)), start = list(a =18, b=1.4)) Error in nls(Brasil ~ a * (b^(anos)), start = list(a = 18, b = 1.4)) : número de iterações excedeu o máximo de 50 André Oliveira Souza. Graduação em Matemática, mestrado em estatística aplicada.Instituto Federal de Educação, Ciência e Tecnologia do Espirito Santo. IFES ______________________________ _________________ R-br mailing list R-br@listas.c3sl.ufpr.br https://listas.inf.ufpr.br/ cgi-bin/mailman/listinfo/r-br Leia o guia de postagem (http://www.leg.ufpr.br/r-br- guia) e forneça código mínimo reproduzível. -- Marcus Nunes http://marcusnunes.me/
André, Apenas para complementar a solução do Marcus, encaminho o código que faz uso de uma função com recursos do pacote rpanel para ajuste de modelo não linear de forma interativa. Essa função (na versão antiga) gerou o GIF que fica na capa da página da lista no Nabble (rs). O uso da interface não elimina o problema do chute inicial, só diminui porque ao invés de passar um valor você passa um intervalo e calibra visualmente. Além disso, tem a parte de análise de resíduos. Brasil <- c(13.63, 13.93, 14.09, 15.12, 15.98, 16.08, 16.45, 17.96, 18.83, 19.28, 19.18, 20.89, 21.87, 21.85, 23.43, 25.38, 25.64, 26.99, 28.68, 31.92, 32.92, 33.32, 34.37, 34.99, 36.08, 37.72, 40.41, 42.81, 44.18, 44.91, 48.97, 51.59, 53.44, 56.49, 62.2, 66.49, 67.01, 68.41, 74.14, 77.58, 80.39, 85.71, 89.58, 96.3, 100.9, 102.8, 107.69, 111.51, 111.42, 116.59, 118.85, 121.42, 127.53) anos <- 1961:2013 anos2 <- anos - min(anos) + 1 # Calibrando o chute inicial (por tentativa-erro faça aproximar). start <- list(a = 12, b = 1.049) plot(Brasil ~ anos2) with(start, curve(a * (b^(anos2)), xname = "anos2", add = TRUE, col = 2)) fit <- nls(Brasil ~ a * (b^(anos2)), start = start) # Verificando o ajuste. plot(Brasil ~ anos2) with(as.list(coef(fit)), curve(a * (b^(anos2)), xname = "anos2", add = TRUE, col = 4)) # rp.nls: função que usa recursos no pacote rpanel para ajuste de # modelos de regressão não linear. source(paste0("https://raw.githubusercontent.com/walmes/", "wzRfun/master/R/rp.nls.R")) library(rpanel) fit <- rp.nls(Brasil ~ a * (b^(anos2)), data = data.frame(Brasil = Brasil, anos2 = anos2), start = list(a = c(5, 20), b = c(0.95, 1.5))) summary(fit) # Função que era do pacote nls2() mas foi removida (guardei o fonte). source("~/Dropbox/CursoR/GeneticaEsalq/as.lm.R") library(nls2) mit <- as.lm(fit) # Diagnóstico do modelo. par(mfrow = c(2, 2)) plot(mit) layout(1) # ! Relação média variância ! # Linear! plot(log(Brasil) ~ anos2) À disposição. Walmes.
Olá Dr. Walmes, obrigado por compartilhar esse script do rpanel, muito legal! Mas ó, não consegui achar a função as.lm no caminho indicado. Daria para anexa-la aqui? Att. -----Mensagem Original----- De: "Walmes Zeviani via R-br" <r-br@listas.c3sl.ufpr.br> Enviada em: 17/10/2016 08:33 Para: "Andre Oliveira" <andreolsouza@yahoo.com.br>; "a lista Brasileira oficial de discussão do programa R." <r-br@listas.c3sl.ufpr.br> Assunto: Re: [R-br] Ajuste de modelo não linear André, Apenas para complementar a solução do Marcus, encaminho o código que faz uso de uma função com recursos do pacote rpanel para ajuste de modelo não linear de forma interativa. Essa função (na versão antiga) gerou o GIF que fica na capa da página da lista no Nabble (rs). O uso da interface não elimina o problema do chute inicial, só diminui porque ao invés de passar um valor você passa um intervalo e calibra visualmente. Além disso, tem a parte de análise de resíduos. Brasil <- c(13.63, 13.93, 14.09, 15.12, 15.98, 16.08, 16.45, 17.96, 18.83, 19.28, 19.18, 20.89, 21.87, 21.85, 23.43, 25.38, 25.64, 26.99, 28.68, 31.92, 32.92, 33.32, 34.37, 34.99, 36.08, 37.72, 40.41, 42.81, 44.18, 44.91, 48.97, 51.59, 53.44, 56.49, 62.2, 66.49, 67.01, 68.41, 74.14, 77.58, 80.39, 85.71, 89.58, 96.3, 100.9, 102.8, 107.69, 111.51, 111.42, 116.59, 118.85, 121.42, 127.53) anos <- 1961:2013 anos2 <- anos - min(anos) + 1 # Calibrando o chute inicial (por tentativa-erro faça aproximar). start <- list(a = 12, b = 1.049) plot(Brasil ~ anos2) with(start, curve(a * (b^(anos2)), xname = "anos2", add = TRUE, col = 2)) fit <- nls(Brasil ~ a * (b^(anos2)), start = start) # Verificando o ajuste. plot(Brasil ~ anos2) with(as.list(coef(fit)), curve(a * (b^(anos2)), xname = "anos2", add = TRUE, col = 4)) # rp.nls: função que usa recursos no pacote rpanel para ajuste de # modelos de regressão não linear. source(paste0("https://raw.githubusercontent.com/walmes/", "wzRfun/master/R/rp.nls.R")) library(rpanel) fit <- rp.nls(Brasil ~ a * (b^(anos2)), data = data.frame(Brasil = Brasil, anos2 = anos2), start = list(a = c(5, 20), b = c(0.95, 1.5))) summary(fit) # Função que era do pacote nls2() mas foi removida (guardei o fonte). source("~/Dropbox/CursoR/GeneticaEsalq/as.lm.R") library(nls2) mit <- as.lm(fit) # Diagnóstico do modelo. par(mfrow = c(2, 2)) plot(mit) layout(1) # ! Relação média variância ! # Linear! plot(log(Brasil) ~ anos2) À disposição. Walmes.
Desculpe minha distração. Eu deveria ter usado essa url: http://leg.ufpr.br/~walmes/cursoR/geneticaEsalq/as.lm.R. À disposição. Walmes.
Andre, segue minha contribuição com um modelo logístico com quatro parâmetros: phi1+(phi2-phi1)/1+exp((phi3-x)/phi4), em que: phi1 é a assíntota inferior; phi2 é a assíntota superior; phi3 é o ponto de inflexão; phi4 é um parâmetro de escala. Você pode tentar ajustar um modelo sem assíntota inferior e ver como o ajuste fica, note que assim o valor da função será zero somente em anos tendendo a menos infinito o que não faz sentido na prática. Pode tentar ajustar com a variável anos em sua escala original também. #============================================== Brasil <- c(13.63, 13.93, 14.09, 15.12, 15.98, 16.08, 16.45, 17.96, 18.83, 19.28, 19.18, 20.89, 21.87, 21.85, 23.43, 25.38, 25.64, 26.99, 28.68, 31.92, 32.92, 33.32, 34.37, 34.99, 36.08, 37.72, 40.41, 42.81, 44.18, 44.91, 48.97, 51.59, 53.44, 56.49, 62.2, 66.49, 67.01, 68.41, 74.14, 77.58, 80.39, 85.71, 89.58, 96.3, 100.9, 102.8, 107.69, 111.51, 111.42, 116.59, 118.85, 121.42, 127.53) anos <- 1961:2013 anos2 <- anos - min(anos) + 1 library(nlme) mod1 <- nls(Brasil~SSfpl(anos2,A,B,C,D)) parm <- coef(mod1) phi1 <- parm[1] phi2 <-parm[2] phi3 <- parm[3] phi4 <- parm[4] plot(Brasil~anos2) with(as.list(coef(mod1)), curve(phi1+(phi2-phi1)/(1+exp((phi3-anos2)/phi4)), xname = "anos2", add = TRUE, col = 4)) -- Atenciosamente, ================================== Felipe E. Barletta Mendes Estatístico(UFPR) - Conre3 9766-A Mestrando em Bioestatística(UEM) +55 (41)-92077191 +55 (41)-33287216 ===================================
participantes (5)
-
Andre Oliveira -
Felipe -
Graciliano -
Marcus Nunes -
Walmes Zeviani