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 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@gmail.com lattes.cnpq.br/5916316648872099
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 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@gmail.com lattes.cnpq.br/5916316648872099
_______________________________________________ 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.
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@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@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 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@gmail.com lattes.cnpq.br/5916316648872099
_______________________________________________ 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.
-- Grato Augusto C. A. Ribas Site Pessoal: http://augustoribas.heliohost.org Lattes: http://lattes.cnpq.br/7355685961127056
Existe uma correção de borda para isto que funciona bem. O seu problema não parece ser com relação diretamente a isso, a mensagem diz respeito ao processo de maximização numerica, o betareg usa um esquema de busca de valores iniciais automatica quando esse procedimento não funciona ele solta a msg. Faz um ajuste preliminar e usa as estimativas deste ajuste novamente como chute inicial e verifique se os resultados estão sensíveis a isso. Vc pode disponibilizar seus dados ? -- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
Mas um duvida que fiquei é a seguinte: # 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") } } # Olhando o grafico acima, da pra comparar com como esta a distribuição dos dados e pensar num chute bom para o start, que da pra inserir pelo betareg.control(). Mas eu tenho dificuldade com matematica e não consegui fazer as estimativas de média e phi que ele estima no ajuste3 voltar a ser o valores que usei no shape1 e shape2 pra gerar os dados em rbeta(). O autor explica na parte 2 do artigo mas eu não entendi. Em 19 de outubro de 2012 16:06, Wagner Bonat <wbonat@gmail.com> escreveu:
Existe uma correção de borda para isto que funciona bem. O seu problema não parece ser com relação diretamente a isso, a mensagem diz respeito ao processo de maximização numerica, o betareg usa um esquema de busca de valores iniciais automatica quando esse procedimento não funciona ele solta a msg. Faz um ajuste preliminar e usa as estimativas deste ajuste novamente como chute inicial e verifique se os resultados estão sensíveis a isso. Vc pode disponibilizar seus dados ?
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
_______________________________________________ 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.
-- Grato Augusto C. A. Ribas Site Pessoal: http://augustoribas.heliohost.org Lattes: http://lattes.cnpq.br/7355685961127056
Agradeço a todos os comentários e ao Augusto por ajudar a explicar melhor todo o problema! Wagner, também sou bastante limitado em matemática por isso do tópico. Aqui estão os dados que estamos testando e o modelo. Ajuda sobre a melhor forma de análisa-los é muito bem vinda! Grato Nicolay _____________________________________________________________________________ #dados index<-c(0.749962502, 0.629968502, 0.990454495, 0.97559371, 0.959552189, 0.992212907, 0.991207709, 0.977094002, 0.960647357, 0.965756886, 0.710047142, 0.887653978, 0, 0.840074171, 0.991626994, 0.926527555, 0.855316235, 0.93723019, 0.779908921, 0.971128419, 0.84293306, 0.69501627, 0.927165576, 0.685086853, 0.97300735, 0, 0.985987221, 0, 0.856096633, 0, 0.749962502, 0.830983382, 0.868087428, 0, 0, 0.692213048, 0.746629335, 0.984632505, 0.969449446, 0.990320893, 0, 0, 0.793822786, 0.925777711, 0.832891421, 0.835937039, 0.083261143, 0.994466643, 0.990346881, 0.995345641, 0.995227688, 0.942101655, 0.987427629, 0.990576035, 0.926467349, 0.897074488, 0.534408541, 0.638856946, 0.996976761, 0.901154942, 0.464509646, 0.822662948, 0.963468058, 0.615255249, 0.884813406, 0.907454627, 0.950820473, 0.44401526, 0.952543568, 0.749516361, 0.691836864, 0.944174245, 0.983889943) periodicity<-structure(c(2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 1L, 2L, 1L, 1L, 1L, 2L, 2L, 1L, 1L, 1L, 1L, 1L, 2L, 2L, 2L, 2L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 3L, 2L, 2L, 1L, 1L, 2L), .Label = c("ephemerous", "permanent", "permanent_instable"), class = "factor") size<-c(1000, 200, 200, 2000, 100, 2000, 2000, 200, 200, 500, 100, 200, 25, 200, 50, 200, 100, 100, 100, 500, 100, 500, 100, 500, 500, 25, 200, 25, 100, 100, 25, 25, 200, 25, 25, 100, 200, 200, 500, 1000, 25, 50, 500, 500, 500, 100, 50, 500, 1000, 2000, 200, 2000, 2000, 500, 100, 50, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 2000, 500, 50, 50, 50, 200) #modelo library(betareg) index[index=='1']<-0.9999 index[findex=='0']<-0.0001 evenn_logit<-betareg(index~log(size)+periodicity) coef(evenn_logit) summary(evenn_logit) 2012/10/19 Augusto Ribas <ribas.aca@gmail.com>
Mas um duvida que fiquei é a seguinte:
# 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") } } #
Olhando o grafico acima, da pra comparar com como esta a distribuição dos dados e pensar num chute bom para o start, que da pra inserir pelo betareg.control(). Mas eu tenho dificuldade com matematica e não consegui fazer as estimativas de média e phi que ele estima no ajuste3 voltar a ser o valores que usei no shape1 e shape2 pra gerar os dados em rbeta(). O autor explica na parte 2 do artigo mas eu não entendi.
Em 19 de outubro de 2012 16:06, Wagner Bonat <wbonat@gmail.com> escreveu:
Existe uma correção de borda para isto que funciona bem. O seu problema não parece ser com relação diretamente a isso, a mensagem diz respeito ao processo de maximização numerica, o betareg usa um esquema de busca de valores iniciais automatica quando esse procedimento não funciona ele solta a msg. Faz um ajuste preliminar e usa as estimativas deste ajuste novamente como chute inicial e verifique se os resultados estão sensíveis a isso. Vc pode disponibilizar seus dados ?
-- Wagner Hugo Bonat LEG - Laboratório de Estatística e Geoinformação UFPR - Universidade Federal do Paraná
_______________________________________________ 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.
-- Grato Augusto C. A. Ribas
Site Pessoal: http://augustoribas.heliohost.org Lattes: http://lattes.cnpq.br/7355685961127056
_______________________________________________ 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.
-- 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 lattes.cnpq.br/5916316648872099
participantes (4)
-
Augusto Ribas -
Eliardo Costa -
Nicolay Cunha -
Wagner Bonat