
Caro Walmes, boa tarde. Li o comentário sobre curva exponencial no e-mail abaixo endereçado ao grupo. Pois bem, estou com um probleminha para resolver, a saber: Tenho estas duas séries: x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23) y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950 ,1000,1050,1100,1150,1200) Utilizei um script apresentado no Livro Conhecendo o R de Luiz A Pertenelli e Márcio P de Mello pág. 165 e não consegui resolver no R. Bom! Ajustei a curva y para exponencial Y no EXCEL conforme abaixo: Y<-c(187,206,227,250,275,303,334,368,405,446,491,541,596,657,723,797,877,966 ,1064,1172,1291,1422,1566) Por favor, como ajusto à exponencial a série y no R? Como chegar ao resultado Y no R? #Script utilizado que não funcionou: x<-c(1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23) y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950 ,1000,1050,1100,1150,1200) dados<-data.frame(x,y) plot(dados) funcao<-y~a*exp(b*x) exponencial<-nls(funcao, dados, start=c(a=1,b=1)) exponencial =====Erro no R ================== 10,11,12,13,14,15,16,17,18,19,20,21,22,23)
y<-c(100,150,200,250,300,350,400,450,500,550,600,650,700,750,800,850,900,950 ,1000,1050,1100,1150,1200)
dados<-data.frame(x,y)
plot(dados)
funcao<-y~a*exp(b*x)
exponencial<-nls(funcao, dados, start=c(a=1,b=1))
Erro em nls(funcao, dados, start = c(a = 1, b = 1)) : gradiente singular
exponencial
Erro: objeto 'exponencial' não encontrado
Onde estou errando? Abs. _____ <http://www.posicional.com/> assinatura_email_sites Cícero C Nunes Diretor Técnico Operacional Divisão Remuneração <mailto:c2n@posicional.com> c2n@posicional.com +55(16)3397.0226 +55(11)3280.0226 +55(11)9.9978.0026 SKYPE: consultormp Rua Padre Duarte, 1360 salas 6/7 Centro Araraquara, SP CEP.: 14.801-310 Antes de imprimir pense em sua responsabilidade com o meio ambiente Quem ama cuida! De: R-br [mailto:r-br-bounces@listas.c3sl.ufpr.br] Em nome de walmes . Enviada em: sexta-feira, 8 de novembro de 2013 18:33 Para: r-br@listas.c3sl.ufpr.br Assunto: Re: [R-br] Contrastes em modelos de regressão de poisson inflacionados por zeros Você não deve esperar ver uma curva exponencial para a média das contagens em função do tempo e uma sigmoide para a proporção de zeros não Poisson para os zeros se não simulou os dados com efeito de tempo e não declarou efeito de tempo na porção binomial. Para que isso aconteça a simulação tem que estar condizente, tem que ter os efeitos para que as curvas não sejam as "curvas nulas". #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(VGAM) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ # Dados artificiais. da <- expand.grid(tempo=rep(1:10), trat=gl(3,10)) ## simula porção Poisson X <- model.matrix(~trat+tempo, da) colnames(X) beta.pois <- c(-1,0.25,0.5,0.25) eta.pois <- X%*%beta.pois da$y.pois <- rpois(nrow(X), lambda=exp(eta.pois)) xyplot(y.pois~tempo|trat, da) ## simula porção binomial X <- model.matrix(~tempo, da) beta.bern <- c(0,0.3) eta.bern <- X%*%beta.bern da$y.bern <- rbinom(nrow(X), size=1, prob=1/(1+exp(-eta.bern))) xyplot(y.bern~tempo|trat, da) ## monta a contagem com inflação de zeros da$y <- with(da, ifelse(y.bern==0, 0, y.pois)) xyplot(y~tempo|trat, da) #------------------------------------------------------------------ # Modelo completo. compl.mod <- zeroinfl(y~trat+tempo|trat+tempo, data=da) coef(compl.mod) #------------------------------------------------------------------ # Predição do modelo considerando as duas porções. pred <- expand.grid(tempo=rep(1:10), trat=gl(3,1)) X <- model.matrix(~trat+tempo, pred) i <- grep("^count\\_", names(coef(compl.mod))) eta <- X%*%coef(compl.mod)[i] pred$m.pois <- exp(eta) i <- grep("^zero\\_", names(coef(compl.mod))) eta <- X%*%coef(compl.mod)[i] pred$p.zero <- exp(eta)/(1+exp(eta)) xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+ as.layer(xyplot(m.pois~tempo|trat, data=pred, type="l"))+ as.layer(xyplot(p.zero~tempo|trat, data=pred, type="l", lty=2, lwd=2))+ layer(panel.abline(h=1, lty=2)) xyplot(p.zero~tempo|trat, data=pred, type="l", lty=2, lwd=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 skype: walmeszeviani homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> linux user number: 531218 ==========================================================================