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
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@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@gmail.com
--
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@gmail.com
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.
_______________________________________________
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.