[R-br] Problema com o pacote gamlss

Alexandre dos Santos alexandresantosbr em yahoo.com.br
Quarta Junho 15 12:06:31 BRT 2011


Felipe, 
   
Obrigado pelo help, mas quando simulo dados teóricos com a mesmas
características que os meus funciona, sendo:


library(emdbook) ##Para pegar a distribuição BB
logitinv=function(p) log(p/(1-p))
logit=function(x) exp(x)/(1+exp(x))

ndata=108##Número de parcelas

ntree=1+rpois(ndata,19)##Simulando o total de plantas

nhole=rnbinom(ndata,mu=5,size=0.3)##Simulando N insetos


##Simulando N de plantas mortas
ndeath=rbinom(ndata,1,prob=logit(-2+1*nhole))*
rbetabinom(ndata,ntree,prob=logit(-1+0.5*nhole),theta=100)
ndeath[is.na(ndeath)]=ntree[is.na(ndeath)]


dd=NULL
dd$y=cbind(ndeath,ntree)
dd$nhole=nhole

require(gamlss)
##Modelisando o número plantas mortas
mm=gamlss(formula=y~1+nhole,sigma.formula=~1,  
      nu.formula=~1+nhole,data=dd,
      family=ZIBB(mu.link="logit",sigma.link="log",nu.link="logit"))
#

Abraço,

Alexandre

-----Original Message-----
From: r-br-bounces em listas.c3sl.ufpr.br
[mailto:r-br-bounces em listas.c3sl.ufpr.br] On Behalf Of Felipe Emanoel
Barletta Mendes
Sent: quarta-feira, 15 de junho de 2011 16:39
To: r-br em listas.c3sl.ufpr.br; felipe em leg.ufpr.br
Subject: Re: [R-br] Problema com o pacote gamlss


Alexandre, na verdade você usa dd$y=cbind(ndeath,ntree-ndeath)
Eu testei ndeath, logo o que disse no e-mail anterior não procede.




Felipe Emanoel Barletta Mendes
> Bom dia Alexandre,
>
> Eu não tenho muita experiência com este modelo mas veja o teste que fiz
> com seus dados:
>
>
> ##### Análise Gráfica
> plot(density(dados$ndeath),main='Histograma e densidade
> estimada',ylim=c(0,.18)
> ,col='red')
> hist(dados$ndeath,freq=F,add=T,col="blue",border='white')
> abline(h=0,col='red')
> lines(density(dados$ndeath),col='red')
> summary(dados$ndeath)
>
>
>
> #### Estimacao numerica
> #
> f1 =
>
gamlss(dados$ndeath~1,family=ZIBB(mu.link="logit",nu.link="logit",sigma.link
="log"))
> Erro em gamlss(dados$ndeath ~ 1, family = ZIBB(mu.link = "logit", nu.link
> = "logit",  :
>   y values must be 0 <= y <= 1
>
>
> A mensagem sugere que a distribuição usada não é a ideal, você está
> modelando uma binomial e a variável resposta não segue distribuição
> binomial. Não seria este o motivo?
>
>
>
>
>
>
> Alexandre dos Santos
>>     Bom dia Pessoal,
>>
>>
>>
>>           Estou tentando modelizar a variável resposta ndeath (plantas
>> mortas)em função da variável nhole (numeros de insetos), com o uso de um
>> modelo hierárquico de Bernoulli e Beta binomial, para isso usei o pacote
>> gamlss, para estimar os parâmetros sendo:
>>
>>
>>
>> require(gamlss)
>>
>> mort_syn<-read.table(file="m_syn.txt", header=T)
>>
>> dd=NULL
>>
>> dd$y=cbind(ndeath,ntree-ndeath)
>>
>> dd$y
>>
>> dd$nhole=nhole
>>
>>
>>
>> mm=gamlss(formula=y~1+nhole,nu.formula=~1+nhole,sigma.formula=~1,
>>
>>       data=dd,
>>
>>       family=ZIBB(mu.link="logit",nu.link="logit",sigma.link="log"))
>>
>>
>>
>> Mais o objeto mm que deveria estimar o parâmetro de Nu esta dando NA,
>> sem
>> erros, sendo:
>>
>>> mm
>>
>>
>>
>> Family:  c("ZIBB", "Zero Inflated Beta Binomial")
>>
>> Fitting method: RS()
>>
>>
>>
>> Call:  gamlss(formula = y ~ 1 + nhole, sigma.formula = ~1, nu.formula =
>> ~1
>> +
>> nhole, family = ZIBB(mu.link = "logit", nu.link = "logit",
>>
>>     sigma.link = "log"), data = dd)
>>
>>
>>
>> Mu Coefficients:
>>
>> (Intercept)        nhole
>>
>>    -2.10530      0.02949
>>
>> Sigma Coefficients:
>>
>> (Intercept)
>>
>>      -3.422
>>
>> Nu Coefficients:
>>
>> (Intercept)        nhole
>>
>>      -2.367           NA
>>
>>
>>
>>  Degrees of Freedom for the fit: 4 Residual Deg. of Freedom   104
>>
>> Global Deviance:     420.808
>>
>>             AIC:     428.808
>>
>>             SBC:     439.537
>>
>>
>>
>>   Alguém que utiliza ou já utilizou o pacote gamlss, saberia me dizer o
>> que
>> esta acontecendo, os dados que utilizei estão anexo,
>>
>> Obrigado,
>>
>>
>>
>> Alexandre dos Santos
>>
>> Ingenieur forestier, Msc.
>>
>> INRA- Biostatistique et Processus Spatiaux (BioSP)
>>
>> Domaine Saint-Paul
>> Site Agroparc
>> 84914 -  Avignon - France
>> Tél. : +33 (0)6 87 95 16 29
>>
>>
>>
>> _______________________________________________
>> R-br mailing list
>> R-br em 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.
>
>
> Felipe E. Barletta Mendes
> (41)9189-5198
> (41)3025-2150
> (41)3328-7216
> http://www.leg.ufpr.br/doku.php/pessoais:felipe
>
> _______________________________________________
> R-br mailing list
> R-br em 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.
>


Felipe E. Barletta Mendes
(41)9189-5198
(41)3025-2150
(41)3328-7216
http://www.leg.ufpr.br/doku.php/pessoais:felipe

_______________________________________________
R-br mailing list
R-br em 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.



Mais detalhes sobre a lista de discussão R-br