[R-br] dúvidas sobre utilização de disitribuição beta

Sexta Outubro 19 15:44:00 BRT 2012

Se for de ajuda, vou tentar explicar como são gerados os dados, quem sabe
ajude a visualizar melhor o problema.

Os valores vem de tipo um indice de reciprocidade de populações de uma
especie de planta.
Mais ou menos assim, salvo excessões, plantas tem estames(parte masculina)
e estigma(parte feminina).
No tentando devido a n fatores (como em nos humanos mesmo), cruzar com
parentes proximos é ruim, logo cruzar com vc mesmo é pessimo.
Dai algumas plantinhas tem uma "estrategia" para evitar isso, que é ter
estames em estigimas em alturas diferentes.
Ai uma planta so vai cruzar com outra. Tentando fazer uma ilustração
grosseira poderias pensar assim:
Note que no caso 1, uma plantinha vai mandar polen do estame baixinho da
planta 1 pro estigma baixinho da planta 2, e o mesmo pros altos, mandar via
polinizadores, o polen gruda na cabeça e no abdomen da abelhas e assim vai

#caso 1
plot(0,0,type="n")

segments(-0.3,-0.5,-0.3,0,col="pink",lwd=5)
segments(-0.5,-0.5,-0.5,1,col="blue",lwd=5)
text(-0.5, -0.6, "Plantinha 1")

segments(0.3,-0.5,0.3,0,col="blue",lwd=5)
segments(0.5,-0.5,0.5,1,col="pink",lwd=5)
text(0.5, -0.6, "Plantinha 2")

legend("topright",legend=c("Femea","Macho"),pch=19,col=c("pink","blue"))
title("Reciprocidade perfeita (valor indice perto de 1)")
#

no entando algumas vezes, por mutações as plantas ficam sem essa "ordem",
como no caso 2, o que gera ela polinizar ela mesma

#caso 2
plot(0,0,type="n")
segments(-0.3,-0.5,-0.3,0,col="pink",lwd=5)
segments(-0.5,-0.5,-0.5,0,col="blue",lwd=5)
text(-0.5, -0.6, "Plantinha 1")

segments(0.3,-0.5,0.3,1,col="blue",lwd=5)
segments(0.5,-0.5,0.5,1,col="pink",lwd=5)
text(0.5, -0.6, "Plantinha 2")

legend("topright",legend=c("Femea","Macho"),pch=19,col=c("pink","blue"))
title("Falha na reciprocidade (valor indice perto de 0)")
#

Agora um cara inventou um indice, que leva em conta todos os individuos de
uma população e da um valor, entre 0 e 1 de "nivel de reciprocidade"

dai olhando os dados, ele parecem com uma distribuição beta.
Eu não entendo perfeitamente a distribuição beta, mas ela tem dois
parametros.
e olhando alguns casos, variando os parametros da distribuição

#
par(mfrow=c(4,4))

for(i in 1:4) {
  for(j in 1:4) {
    hist(rbeta(1000,i,j),main=paste("shape1 =",i,"shape2 =",j),xlab="")
    função<-function(x) {dbeta(x,i,j)}
    curve(função(x)*100,0,1,add=T,lwd=2,col="red")
  }
}
#

Da pra ver alguns exemplos que se parecem exatamente com conjunto de dados:

#So que a hipotese a ser testada, é que dependendo do ambiente, as
populações tem uma tendencia a exibirem uma menor reciprocidade
#então pensando em populações que vivem em lagoas temporarios e populações
que vivem em lagoas permanentes (para uma planta aquatica)
#podemos fazer um exemplo que parece com os dados usando o rbeta():

set.seed(01)
dados<-data.frame(lagoa=rep(c("Temporario","Permanente"),each=50),indice=c(rbeta(50,4,1),rbeta(50,4,2)))
library(lattice)
histogram(~indice|lagoa,dados)

#usar uma regressão teria problemas devido a normalidade dos residuos
ajuste01<-lm(dados$indice~dados$lagoa)
summary(ajuste01)
par(mfrow=c(2,2))
plot(ajuste01)

#uma solução que vi em livros é a transformação  pelo arco sena da raiz
quadrada, mas dai é estimado parametros que são dificeis de entender
#dificil pensar no arcoseno da raiz quadrada de algo
histogram(~asin(sqrt(indice))|lagoa,dados)
ajuste02<-lm(asin(sqrt(dados$indice))~dados$lagoa)
summary(ajuste02)
par(mfrow=c(2,2))
plot(ajuste02)

#procurando alguma abordagem com GLM, pra descrever os dados usando
distribuição beta, vimos o betareg
library(betareg)
ajuste3<-betareg(indice~lagoa|lagoa,dados)
summary(ajuste3)
par(mfrow=c(2,2))
plot(ajuste03)

#o ajuste dele da 2 valores, uma media e um phi, que estava interpretando
como referentes ao shape 2 e shape 1 do comando rbeta()
#so que não tenho certeza se isso faz sentido

#mas o ajuste parece bem razoavel, mas não sei se essa abordagem usando
betareg é a correta pra esse tipo de dados:
round(tapply(dados$indice,dados$lagoa,mean),digits=3)
paste("estimado temporario
=",round(plogis(ajuste3$coefficients$mean[1]),digits=3))
paste("estimado permanente
=",round(plogis(ajuste3$coefficients$mean[1]+ajuste3$coefficients$mean[2]),digits=3))

#se alguem puder dar uma luz

Em 18 de outubro de 2012 10:54, Eliardo Costa <eliardocosta em gmail.com>escreveu:

> Olá Nicolay,
>
> Como seus dados apresentam 0's e 1's acredito que modelos de regressão
> beta inflacionados possam lhe ajudar, porém não sei se existe algo
> implementado no R. Pesquise sobre isso.
>
> Att, Eliardo.
>
> Em 18 de outubro de 2012 11:45, Nicolay Cunha <nicolaycunha em gmail.com>escreveu:
>
>> Olá a todos,
>>
>> Tenho uma variável dependente (Y) que varia de 0 a 1 e quero testar em
>> relação a duas variáveis categóricas (X1 com três níveis e X2 com 7
>> níveis). A categórica X2 transformo em números para ganhar graus de
>> liberdade.
>>
>> Os dados não atendem nenhum tipo de premissa de linearidade, portanto
>> tive que buscar outras alternativas para analisa-los.
>> Ao estudar sobre a melhor forma de ajustar esses dados, encontrei o
>> artigo (
>> http://cran.r-project.org/web/packages/betareg/vignettes/betareg.pdf)
>> sobre o tema, e os exemplos dados aparentam encaixar com o meu caso.
>>
>>
>> No entanto meu Y tem valores inteiros (0 e 1) e a distribuiçãoa aparenta
>> não aceitar isso, para tal transformei todos os 0s em 0.0001 e 1s em
>> 0.9999.  Ainda tenho dúvidas sobre a interpretação dos dados (inclusive o
>> papel do phi), utilização e escolha da melhor forma de gerar o modelo da
>> forma que estou propondo, logo gostaria de sugestões e orientações sobre
>> como proceder com essas análises.
>>
>>
>> Abaixo um CMR com um exemplo do que estou fazendo onde aparece uma
>> mensagem de aviso perdida ao estimar o modelo.
>>
>>
>> *"In betareg.fit(X, Y, Z, weights, offset, link, link.phi, type,
>> control) :
>>   no valid starting value for precision parameter found, using 1 instead"
>> *
>>
>>
>> #############################################################################
>>
>> X1<-factor(rep(c("A", "B", "C"),50)) ## categoria 1
>> X2<- rep(sample(1:7), len=150 ) ## categoria 2
>> Y<-rep(sample(0:100)/100, len=150) ## minha variável dependente em
>> proporção
>> Y[Y=='1']<-0.99999 # transformo 1 em 0.99999 pois a distribuição beta não
>> aceita 0 e 1, somente o intervalo entre eles
>> Y[Y=='0']<-0.00001 # transformo 0 em 0.00001
>>
>> library(betareg)
>> beta1<-betareg(Y~X1+X2, link="loglog")
>> coef(beta1)
>> summary(beta1)
>>
>> beta2<-betareg(Y~X1+X2|X2, link="loglog") ## adicionando X2 como um
>> regressor adicional
>> coef(beta2)
>> summary(beta2)
>>
>> AIC(beta1, beta2,  k = log(nrow(data.frame(X1,X2,Y))))
>>
>>
>> ##############################################################################
>>
>>
>>
>> Agradeço a ajuda!
>> Nicolay
>>
>>
>>
>>
>>
>>
>>
>> --
>> Nicolay Leme da Cunha
>>
>> Biólogo, Mestre, Doutorando em Ecologia e Conservação
>> Universidade Federal de Mato Grosso do Sul, 79070-900
>> Campo Grande, MS, Brasil
>> E-mail: nicolaycunha em gmail.com
>> lattes.cnpq.br/5916316648872099
>>
>>
>>
>>
>> --
>> Nicolay Leme da Cunha
>>
>> Biólogo, Mestre, Doutorando em Ecologia e Conservação
>> Universidade Federal de Mato Grosso do Sul, 79070-900
>> Campo Grande, MS, Brasil
>> E-mail: nicolaycunha em gmail.com
>> lattes.cnpq.br/5916316648872099
>>
>>
>> _______________________________________________
>> 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.
>>
>
>
> _______________________________________________
> 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.
>

-- 
Grato
Augusto C. A. Ribas

Site Pessoal: http://augustoribas.heliohost.org
Lattes: http://lattes.cnpq.br/7355685961127056
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121019/c68c5a39/attachment.html>