Desdobramento de interações em dados com distribuição de Poisson inflacionada por zeros
Bom dia Pessoal, Estou quebrando a cabeça para tentar desdobrar interações em dados com distribuição de Poisson inflacionada por zeros, na verdade é uma situação muito recorrente em entomologia, onde em algumas situações temos zero insetos e em outras muitos, considerando-se dados de contagens. Fiz uma pequena simulação abaixo para ilustrar o problema (Variáveis observadas Trat e Epoca e variável resposta y), se acho um modelo significativo que tem uma parte com distribuição de Poisson e outra com Binomial, como faço o desdobramento? Uma para cada distribuição ou tem alguma abordagem diferente? Segue CRM abaixo: ## Desdobramento de dados com distribuição de Poisson inflacionado por zeros --- require(pscl) y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2) trat <- as.factor(gl(3,600)) ##Criação dos tratamentos tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo dados<-as.data.frame(cbind(trat,tempo, y)) #------------------------------------------------------------------------------- # Análise de variância do dados inflacionados summary(m1 <- zeroinfl(y ~ trat | epoca, data = dados)) ## Modelo completo mnull <- update(m1, . ~ 1) ### Modelo nulo pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo ## Comparando com GLM Poisson sem inflação por zeros --------------------------- summary(p1 <- glm(y ~ trat + epoca, family = poisson, data = dados)) ##Teste de Vuong---------------------------------------------------------------- vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção # ## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial? Obrigado, -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Não estão sendo encontrados os objetos 'm1' e 'época' Bom dia Pessoal, Estou quebrando a cabeça para tentar desdobrar interações em dados com distribuição de Poisson inflacionada por zeros, na verdade é uma situação muito recorrente em entomologia, onde em algumas situações temos zero insetos e em outras muitos, considerando-se dados de contagens. Fiz uma pequena simulação abaixo para ilustrar o problema (Variáveis observadas Trat e Epoca e variável resposta y), se acho um modelo significativo que tem uma parte com distribuição de Poisson e outra com Binomial, como faço o desdobramento? Uma para cada distribuição ou tem alguma abordagem diferente? Segue CRM abaixo: ## Desdobramento de dados com distribuição de Poisson inflacionado por zeros --- require(pscl) y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2) trat <- as.factor(gl(3,600)) ##Criação dos tratamentos tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo dados<-as.data.frame(cbind(trat,tempo, y)) #------------------------------------------------------------------------------- # Análise de variância do dados inflacionados summary(m1 <- zeroinfl(y ~ trat | epoca, data = dados)) ## Modelo completo mnull <- update(m1, . ~ 1) ### Modelo nulo pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo ## Comparando com GLM Poisson sem inflação por zeros --------------------------- summary(p1 <- glm(y ~ trat + epoca, family = poisson, data = dados)) ##Teste de Vuong---------------------------------------------------------------- vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção # ## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial? Obrigado, -- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680
Olhando por cima, parece-me que você tem uma mistura de distribuições, não? Talvez o pacote flexmix possa te ajudar nesse aspecto. http://rss.acs.unt.edu/Rdoc/library/flexmix/html/flexmix.html abcs M 2013/9/22 Mauro Sznelwar <sznelwar@uol.com.br>
** Não estão sendo encontrados os objetos 'm1' e 'época'
Bom dia Pessoal,
Estou quebrando a cabeça para tentar desdobrar interações em dados com distribuição de Poisson inflacionada por zeros, na verdade é uma situação muito recorrente em entomologia, onde em algumas situações temos zero insetos e em outras muitos, considerando-se dados de contagens. Fiz uma pequena simulação abaixo para ilustrar o problema (Variáveis observadas Trat e Epoca e variável resposta y), se acho um modelo significativo que tem uma parte com distribuição de Poisson e outra com Binomial, como faço o desdobramento? Uma para cada distribuição ou tem alguma abordagem diferente? Segue CRM abaixo:
## Desdobramento de dados com distribuição de Poisson inflacionado por zeros ---
require(pscl)
y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2)
trat <- as.factor(gl(3,600)) ##Criação dos tratamentos
tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo
dados<-as.data.frame(cbind(trat,tempo, y))
#------------------------------------------------------------------------------- # Análise de variância do dados inflacionados
summary(m1 <- zeroinfl(y ~ trat | epoca, data = dados)) ## Modelo completo
mnull <- update(m1, . ~ 1) ### Modelo nulo
pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo
## Comparando com GLM Poisson sem inflação por zeros ---------------------------
summary(p1 <- glm(y ~ trat + epoca, family = poisson, data = dados))
##Teste de Vuong----------------------------------------------------------------
vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção
#
## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial?
Obrigado,
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680
_______________________________________________ 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.
-- Manoel Galdino https://sites.google.com/site/galdinomcz/
Bom dia Mauro, Na verdade troquei tempo por epoca e o m1 depende da instalação do pacote pscl, ficando o CRM corrigido: require(pscl) y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2) trat <- as.factor(gl(3,600)) ##Criação dos tratamentos tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo dados<-as.data.frame(cbind(trat,tempo, y)) #------------------------------------------------------------------------------- # Análise de variância do dados inflacionados summary(m1 <- zeroinfl(y ~ trat | tempo, data = dados)) ## Modelo completo mnull <- update(m1, . ~ 1) ### Modelo nulo pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo ## Comparando com GLM Poisson sem inflação por zeros --------------------------- summary(p1 <- glm(y ~ trat + epoca, family = poisson, data = dados)) ##Teste de Vuong---------------------------------------------------------------- vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção # Em 22/09/2013 19:34, Mauro Sznelwar escreveu:
Não estão sendo encontrados os objetos 'm1' e 'época'
Bom dia Pessoal,
Estou quebrando a cabeça para tentar desdobrar interações em dados com distribuição de Poisson inflacionada por zeros, na verdade é uma situação muito recorrente em entomologia, onde em algumas situações temos zero insetos e em outras muitos, considerando-se dados de contagens. Fiz uma pequena simulação abaixo para ilustrar o problema (Variáveis observadas Trat e Epoca e variável resposta y), se acho um modelo significativo que tem uma parte com distribuição de Poisson e outra com Binomial, como faço o desdobramento? Uma para cada distribuição ou tem alguma abordagem diferente? Segue CRM abaixo:
## Desdobramento de dados com distribuição de Poisson inflacionado por zeros ---
require(pscl)
y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2)
trat <- as.factor(gl(3,600)) ##Criação dos tratamentos
tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo
dados<-as.data.frame(cbind(trat,tempo, y))
#------------------------------------------------------------------------------- # Análise de variância do dados inflacionados
summary(m1 <- zeroinfl(y ~ trat | epoca, data = dados)) ## Modelo completo
mnull <- update(m1, . ~ 1) ### Modelo nulo
pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo
## Comparando com GLM Poisson sem inflação por zeros ---------------------------
summary(p1 <- glm(y ~ trat + epoca, family = poisson, data = dados))
##Teste de Vuong----------------------------------------------------------------
vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção
#
## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial?
Obrigado,
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br <mailto:alexandre.santos@cas.ifmt.edu.br> Lattes: http://lattes.cnpq.br/1360403201088680
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Antes de desdobrar a interação é preciso saber qual o objetivo da análise. É possível desdobrar de várias formas, fazendo contrastes específicos, comparações múltiplas. Eu gosto da abordagem de gráficos com os valores preditos (médias) e intervalos de confiança (ou bandas). Todavia existem pessoas da escola antiga que querer ver tabelas com letras ao lado de médias um em CV baixo, infelizmente. À 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 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
Walmes, Na verdade eu gostaria de saber quais tratamentos tiveram mesmo desempenho no tempo e representá-las através de curvas, mas não sei como lidar com a decomposição das somas de quadrados quando tenho duas distribuições (Poisson e Binomial), pois espero dois resultados: 1) Poisson: a quantidade de insetos representada para cada um dos tratamentos no tempo e; 2)Binomial: a ocorrência ou não dos insetos em cada um dos tratamentos no tempo. Teria alguma abordagem para sugerir? Obrigado, Alexandre CRM: ## Desdobramento de dados com distribuição de Poisson inflacionado por zeros --- require(pscl) y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2) trat <- as.factor(gl(3,600)) ##Criação dos tratamentos tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo dados<-as.data.frame(cbind(trat,tempo, y)) #------------------------------------------------------------------------------- # Análise de variância do dados inflacionados summary(m1 <- zeroinfl(y ~ trat | tempo, data = dados)) ## Modelo completo mnull <- update(m1, . ~ 1) ### Modelo nulo pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo ## Comparando com GLM Poisson sem inflação por zeros --------------------------- summary(p1 <- glm(y ~ trat + tempo, family = poisson, data = dados)) ##Teste de Vuong---------------------------------------------------------------- vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção # ## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial? Em 23/09/2013 09:26, walmes . escreveu:
Antes de desdobrar a interação é preciso saber qual o objetivo da análise. É possível desdobrar de várias formas, fazendo contrastes específicos, comparações múltiplas. Eu gosto da abordagem de gráficos com os valores preditos (médias) e intervalos de confiança (ou bandas). Todavia existem pessoas da escola antiga que querer ver tabelas com letras ao lado de médias um em CV baixo, infelizmente.
À 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 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br <mailto:walmes@ufpr.br> skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> linux user number: 531218 ==========================================================================
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Desdobrar as deviances tal como fazemos o desdobramento de somas de quadrados em modelos lineares gaussianos seria interessante mas eu acredito que os gráficos com os modelos ajutados (com bandas) já são suficientes para discutir os resultados. Isso é bem mais fácil de obter. Quando o leitor for ver os resultados ele vai gastar 5 segundos na tabela de desdobramento e 2 minutos vendo o seu gráfico. À 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 linux user number: 531218 ========================================================================== 2013/9/23 ASANTOS <alexandresantosbr@yahoo.com.br>
Walmes,
Na verdade eu gostaria de saber quais tratamentos tiveram mesmo desempenho no tempo e representá-las através de curvas, mas não sei como lidar com a decomposição das somas de quadrados quando tenho duas distribuições (Poisson e Binomial), pois espero dois resultados: 1) Poisson: a quantidade de insetos representada para cada um dos tratamentos no tempo e; 2)Binomial: a ocorrência ou não dos insetos em cada um dos tratamentos no tempo. Teria alguma abordagem para sugerir?
Obrigado,
Alexandre
CRM:
## Desdobramento de dados com distribuição de Poisson inflacionado por zeros ---
require(pscl)
y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2)
trat <- as.factor(gl(3,600)) ##Criação dos tratamentos
tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo
dados<-as.data.frame(cbind(trat,tempo, y))
#------------------------------------------------------------------------------- # Análise de variância do dados inflacionados
summary(m1 <- zeroinfl(y ~ trat | tempo, data = dados)) ## Modelo completo
mnull <- update(m1, . ~ 1) ### Modelo nulo
pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo
## Comparando com GLM Poisson sem inflação por zeros ---------------------------
summary(p1 <- glm(y ~ trat + tempo, family = poisson, data = dados))
##Teste de Vuong----------------------------------------------------------------
vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção
#
## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial?
Em 23/09/2013 09:26, walmes . escreveu:
Antes de desdobrar a interação é preciso saber qual o objetivo da análise. É possível desdobrar de várias formas, fazendo contrastes específicos, comparações múltiplas. Eu gosto da abordagem de gráficos com os valores preditos (médias) e intervalos de confiança (ou bandas). Todavia existem pessoas da escola antiga que querer ver tabelas com letras ao lado de médias um em CV baixo, infelizmente.
À 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 VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing listR-br@listas.c3sl.ufpr.brhttps://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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO)e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
_______________________________________________ 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.
Walmes, Atendi a sua sugestão e resolvi ajustar os modelos, no entanto, fiz as comparações múltiplas para saber a diferença entre os tratamentos e tudo bem, todos foram significativamente diferentes, porém na hora de fazer as curvas para cada tratamento no tempo não consigo obter os coeficientes do modelo para cada um dos tratamentos ou seja quando faço o summary() do modelo completo tenho apenas um coeficiente para trat e outro para tempo. Gostrai o que posso fazer para obter os coeficientes para cada tratamento? Segue CRM: ## Desdobramento de dados com distribuição de Poisson inflacionado por zeros --- require(pscl) require(multicomp) y1<- c(mapply(rpois, lambda=c(4,45,93), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2) trat <- as.factor(gl(3,600)) ##Criação dos tratamentos tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo dados<-as.data.frame(cbind(trat,tempo, y)) #------------------------------------------------------------------------------- # Análise de variância do dados inflacionados summary(m1 <- zeroinfl(y ~ trat*tempo | trat*tempo, data = dados)) ## Modelo completo mnull <- update(m1, . ~ 1) ### Modelo nulo pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo ## Comparando com GLM Poisson sem inflação por zeros --------------------------- summary(p1 <- glm(y ~ trat + tempo, family = poisson, data = dados)) ##Teste de Vuong---------------------------------------------------------------- vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção # ### Interação entre os tratamamentos no tempo----------------------------------- zi_m1<-zeroinfl(y ~ trat*tempo | trat*tempo, data = dados, link = "logit") summary(zinb_edanx2) ## Modelos completo ## Calculando os contrastes para o modelo inflacionado-------------------------- nr <- length(levels(trat)) contr <- matrix(0, nrow = nr, ncol = length(coef(zi_m1))) colnames(contr) <- names(coef(zi_m1)) rownames(contr) <- paste(levels(trat)[c(2, 3, 3)], levels(trat)[c(1, 1, 2)], sep = " - ") contr[,3:4] <- contrMat(numeric(nrow(contr)), type = "Tukey")[,-2] glht_zi <- glht(zi_m1, linfct = contr) ## Comparações multiplas------------------------------------------------------------------- summary(glht_zi) ## Realizando um ajuste para cada tratamento - Primeiro para Poisson sort(tapply(y,trat,mean)) levels(trat) zi_final<-zeroinfl(y ~ trat+tempo | trat+tempo, data = dados, link = "logit") summary(zi_final) ## Aqui não consigo descobrir os coeficientes para cada tratamento para fazer as curvas!!!!! Obrigado, Alexandre Em 26/09/2013 16:28, walmes . escreveu:
Desdobrar as deviances tal como fazemos o desdobramento de somas de quadrados em modelos lineares gaussianos seria interessante mas eu acredito que os gráficos com os modelos ajutados (com bandas) já são suficientes para discutir os resultados. Isso é bem mais fácil de obter. Quando o leitor for ver os resultados ele vai gastar 5 segundos na tabela de desdobramento e 2 minutos vendo o seu gráfico.
À 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 ==========================================================================
2013/9/23 ASANTOS <alexandresantosbr@yahoo.com.br <mailto:alexandresantosbr@yahoo.com.br>>
Walmes,
Na verdade eu gostaria de saber quais tratamentos tiveram mesmo desempenho no tempo e representá-las através de curvas, mas não sei como lidar com a decomposição das somas de quadrados quando tenho duas distribuições (Poisson e Binomial), pois espero dois resultados: 1) Poisson: a quantidade de insetos representada para cada um dos tratamentos no tempo e; 2)Binomial: a ocorrência ou não dos insetos em cada um dos tratamentos no tempo. Teria alguma abordagem para sugerir?
Obrigado,
Alexandre
CRM:
## Desdobramento de dados com distribuição de Poisson inflacionado por zeros ---
require(pscl)
y1<- c(mapply(rpois, lambda=c(5,20,45), MoreArgs=list(n=400)))##Criação da variável resposta Poisson y2<- c(mapply(rbinom, size=c(1,0,0), prob=c(0.5,1,1), MoreArgs=list(n=200)))##Criação da variável resposta Binomial y<-c(y1,y2)
trat <- as.factor(gl(3,600)) ##Criação dos tratamentos
tempo<- as.factor(rep(gl(6,100),3)) ### Criação da variável tempo
dados<-as.data.frame(cbind(trat,tempo, y))
#------------------------------------------------------------------------------- # Análise de variância do dados inflacionados
summary(m1 <- zeroinfl(y ~ trat | tempo, data = dados)) ## Modelo completo
mnull <- update(m1, . ~ 1) ### Modelo nulo
pchisq(2 * (logLik(m1) - logLik(mnull)), df = 2, lower.tail = FALSE) ## Teste de Chi o modelo completo foi significativo
## Comparando com GLM Poisson sem inflação por zeros ---------------------------
summary(p1 <- glm(y ~ trat + tempo, family = poisson, data = dados))
##Teste de Vuong----------------------------------------------------------------
vuong(p1, m1) ## O GLM Poisson estava tão mal ajustado que o Poisson iflacionado é a única opção
#
## E agora, como desdobrar isto? Faço um desdobramento para parte de poisson e outra para parte binomial?
Em 23/09/2013 09:26, walmes . escreveu:
Antes de desdobrar a interação é preciso saber qual o objetivo da análise. É possível desdobrar de várias formas, fazendo contrastes específicos, comparações múltiplas. Eu gosto da abordagem de gráficos com os valores preditos (médias) e intervalos de confiança (ou bandas). Todavia existem pessoas da escola antiga que querer ver tabelas com letras ao lado de médias um em CV baixo, infelizmente.
À 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 <tel:%28%2B55%29%2041%203361%203573> VoIP: (3361 3600) 1053 1173 e-mail: walmes@ufpr.br <mailto:walmes@ufpr.br> skype: walmeszeviani twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> linux user number: 531218 ==========================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br <mailto: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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone:(+55) 65 8132-8112 <tel:%28%2B55%29%2065%208132-8112> (TIM)(+55) 65 9686-6970 <tel:%28%2B55%29%2065%209686-6970> (VIVO) e-mails:alexandresantosbr@yahoo.com.br <mailto:e-mails:alexandresantosbr@yahoo.com.br> alexandre.santos@cas.ifmt.edu.br <mailto:alexandre.santos@cas.ifmt.edu.br> Lattes:http://lattes.cnpq.br/1360403201088680 ======================================================================
_______________________________________________ R-br mailing list R-br@listas.c3sl.ufpr.br <mailto: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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Alexandre, segue CMR. #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ ## Para ter curvas você precisa que tempo seja númerico com mais ## de 3 níveis, caso contrário não faz sentido. #------------------------------------------------------------------ # Dados artificiais. da <- expand.grid(trat=gl(2,1), tempo=1:10) X <- model.matrix(~trat+tempo, da); ncol(X) betas <- c(0.1,0.1,0.3) eta <- X%*%betas y1 <- rpois(da$trat, lambda=exp(eta)) y2 <- rbinom(y1, size=1, prob=0.7) da$y <- y1*y2 str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE) #------------------------------------------------------------------ # Ajuste do modelo. m0 <- zeroinfl(y~trat+tempo|trat, data=da) summary(m0) #------------------------------------------------------------------ # Predição do modelo considerando as duas porções. X <- model.matrix(~trat+tempo, da) i <- grep("^count\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.pois <- exp(eta) X <- model.matrix(~trat, da) i <- grep("^zero\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.zero <- exp(eta)/(1+exp(eta)) xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+ as.layer(xyplot(y.pois~tempo|trat, data=da, type="l"))+ as.layer(xyplot(y.zero~tempo|trat, data=da, type="l", lty=2, lwd=2))+ layer(panel.abline(h=1, lty=2)) # contínua: média da contagem ~ Poisson. # tracejada: probabilidade de um zero não Poisson. # abline: linha no 1, referência. #------------------------------------------------------------------ À 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 linux user number: 531218 ==========================================================================
Obrigado Walmes e uma última pergunta: A representação matemática dos meus modelos segundo a análise abaixo: Count model coefficients (poisson with log link): Estimate Std. Error z value Pr(>|z|) (Intercept) 0.870118 0.092920 9.364 < 2e-16 *** trat2 0.129159 0.029110 4.437 9.13e-06 *** tempo 0.701863 0.009766 71.871 < 2e-16 *** fica para tratamento 1: y = Intercept + e^(0.701863*tempo) e tratamento 2: y= 0.870118 + e^(0.701863*tempo) + 0.129159 ? Segue abaixo CRM para futuras consultas de membros da lista: #------------------------------------------------------------------ # Definições da sessão. rm(list=ls()) require(pscl) require(multcomp) require(lattice) require(latticeExtra) #------------------------------------------------------------------ ## Para ter curvas você precisa que tempo seja númerico com mais ## de 3 níveis, caso contrário não faz sentido. #------------------------------------------------------------------ # Dados artificiais. da <- expand.grid(trat=gl(4,1), tempo=1:10) X <- model.matrix(~trat+tempo, da); ncol(X) betas <- c(0.1,0.9,0.6,0.3,0.7) eta <- X%*%betas y1 <- rpois(da$trat, lambda=exp(eta)) y2 <- rbinom(y1, size=1, prob=0.7) da$y <- y1*y2 str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE) #------------------------------------------------------------------ # Ajuste do modelo. m0 <- zeroinfl(y~trat+tempo|trat, data=da) summary(m0) #------------------------------------------------------------------ # Predição do modelo considerando as duas porções. X <- model.matrix(~trat+tempo, da) i <- grep("^count\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.pois <- exp(eta) X <- model.matrix(~trat, da) i <- grep("^zero\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.zero <- exp(eta)/(1+exp(eta)) xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+ as.layer(xyplot(y.pois~tempo|trat, data=da, type="l"))+ as.layer(xyplot(y.zero~tempo|trat, data=da, type="l", lty=2, lwd=2))+ layer(panel.abline(h=1, lty=2)) # contínua: média da contagem ~ Poisson. # tracejada: probabilidade de um zero não Poisson. # abline: linha no 1, referência. #------------------------------------------------------------------ Em 30/09/2013 10:53, walmes . escreveu:
Alexandre, segue CMR.
#------------------------------------------------------------------ # Definições da sessão.
rm(list=ls()) require(pscl) require(multcomp) require(lattice) require(latticeExtra)
#------------------------------------------------------------------
## Para ter curvas você precisa que tempo seja númerico com mais ## de 3 níveis, caso contrário não faz sentido.
#------------------------------------------------------------------ # Dados artificiais.
da <- expand.grid(trat=gl(2,1), tempo=1:10) X <- model.matrix(~trat+tempo, da); ncol(X) betas <- c(0.1,0.1,0.3) eta <- X%*%betas y1 <- rpois(da$trat, lambda=exp(eta)) y2 <- rbinom(y1, size=1, prob=0.7) da$y <- y1*y2 str(da) xyplot(y~tempo|trat, data=da, jitter.x=TRUE)
#------------------------------------------------------------------ # Ajuste do modelo.
m0 <- zeroinfl(y~trat+tempo|trat, data=da) summary(m0)
#------------------------------------------------------------------ # Predição do modelo considerando as duas porções.
X <- model.matrix(~trat+tempo, da) i <- grep("^count\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.pois <- exp(eta)
X <- model.matrix(~trat, da) i <- grep("^zero\\_", names(coef(m0))) eta <- X%*%coef(m0)[i] da$y.zero <- exp(eta)/(1+exp(eta))
xyplot(y~tempo|trat, data=da, jitter.x=TRUE)+ as.layer(xyplot(y.pois~tempo|trat, data=da, type="l"))+ as.layer(xyplot(y.zero~tempo|trat, data=da, type="l", lty=2, lwd=2))+ layer(panel.abline(h=1, lty=2))
# contínua: média da contagem ~ Poisson. # tracejada: probabilidade de um zero não Poisson. # abline: linha no 1, referência.
#------------------------------------------------------------------
À 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 <tel:%28%2B55%29%2041%203361%203573> skype: walmeszeviani homepage: http://www.leg.ufpr.br/~walmes <http://www.leg.ufpr.br/%7Ewalmes> linux user number: 531218 ==========================================================================
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Para Poisson com função de ligação log o modelo fica exp(X%*%beta) que no caso fica * tratamento base: exp(b0+b1*trat) * tratamento outro: exp(b0+b0i+b1*trat); Para porção zero que considera Bernoulli com função de ligação logit fica exp(X%*%beta)/(1+exp(X%*%beta)), assim * tratamento base: exp(b0)/(1+exp(b0)) * tratamento outro: exp(b0+b0i)/(1+exp(b0+b0i)); As contas matriciais que fiz (veja onde tem o objeto eta) estão correspondendo à essas especificações que foram os modelos que eu declarei ao gerar os dados artificiais. À 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 linux user number: 531218 ==========================================================================
Walmes, Agradeço muito, mas ainda não entendi por que a variável tempo não é contemplada no modelo? Obrigado, Alexandre Em 30/09/2013 14:55, walmes . escreveu:
Para Poisson com função de ligação log o modelo fica exp(X%*%beta) que no caso fica * tratamento base: exp(b0+b1*trat) * tratamento outro: exp(b0+b0i+b1*trat); Para porção zero que considera Bernoulli com função de ligação logit fica exp(X%*%beta)/(1+exp(X%*%beta)), assim * tratamento base: exp(b0)/(1+exp(b0)) * tratamento outro: exp(b0+b0i)/(1+exp(b0+b0i)); As contas matriciais que fiz (veja onde tem o objeto eta) estão correspondendo à essas especificações que foram os modelos que eu declarei ao gerar os dados artificiais.
À 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 ==========================================================================
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
Nada impede de ser. Eu só quis que fosse para mostrar que os modelos não necessariamente precisam ter os mesmos termos. Declare o modelo que considerar mais correto para os seus dados e proceda com as adaptações correspondentes para gerar os resultados. À 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 linux user number: 531218 ==========================================================================
Muito obrigado pela ajuda Walmes, problema resolvido. Alexandre Em 30/09/2013 15:46, walmes . escreveu:
Nada impede de ser. Eu só quis que fosse para mostrar que os modelos não necessariamente precisam ter os mesmos termos. Declare o modelo que considerar mais correto para os seus dados e proceda com as adaptações correspondentes para gerar os resultados.
À 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 ==========================================================================
_______________________________________________ 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.
-- ====================================================================== Alexandre dos Santos Proteção Florestal IFMT - Instituto Federal de Educação, Ciência e Tecnologia de Mato Grosso Campus Cáceres Caixa Postal 244 Avenida dos Ramires, s/n Bairro: Distrito Industrial Cáceres - MT CEP: 78.200-000 Fone: (+55) 65 8132-8112 (TIM) (+55) 65 9686-6970 (VIVO) e-mails:alexandresantosbr@yahoo.com.br alexandre.santos@cas.ifmt.edu.br Lattes: http://lattes.cnpq.br/1360403201088680 ======================================================================
participantes (4)
-
ASANTOS -
Manoel Galdino -
Mauro Sznelwar -
walmes .