Comparações múltiplas usando MMC (mean–mean multiple comparisons) plots

*Bom dia pessoal, Estou tendo problemas em visualizar os resultados de uma comparação múltipla usando GLM, com o uso MMC plots do pacote HH. Bom, tenho quatro grupos funcionais de cupins que são humívoros (humi), ceifadores(ceifa), intermediários (inter) e xilófagos (xilo) que são meus tratamentos e como variável resposta o número de ocorrência de espécies de cada grupo (acorr). Para coleta da variável resposta, fizemos 6 transectos que são as repetições. Para ver se existe diferença entre os grupos funcionais fiz:* # mod2<-glm(ocorr~grupo, family="quasipoisson", data=trat1_90antes)## Um plot(mod2) prévio mostrou que quasipoisson foi a distribuição mais adequada aos dados que poisson summary(mod2) glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Deviance Residuals: Min 1Q Median 3Q Max -1.9148 -0.4946 0.1213 0.4178 1.3813 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 1.5041 0.1624 9.260 1.13e-08 *** grupohumi -1.0986 0.3249 -3.382 0.00296 ** grupointer -0.8979 0.3019 -2.974 0.00750 ** grupoxilo -0.6568 0.2780 -2.363 0.02838 * --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Dispersion parameter for quasipoisson family taken to be 0.7123377) Null deviance: 28.663 on 23 degrees of freedom Residual deviance: 16.889 on 20 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 5 anova(mod2, test="Chi") ## Anova do modelo Analysis of Deviance Table Model: quasipoisson, link: log Response: ocorr Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev Pr(>Chi) NULL 23 28.663 grupo 3 11.774 20 16.889 0.0008834 *** --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 plot(mod2) ##Teste de comparação múltipla summary(glht(mod2, linfct=mcp(grupo="Tukey"))) Simultaneous Tests for General Linear Hypotheses Multiple Comparisons of Means: Tukey Contrasts Fit: glm(formula = ocorr ~ grupo, family = "quasipoisson", data = trat1_90antes) Linear Hypotheses: Estimate Std. Error z value Pr(>|z|) humi - ceifa == 0 -1.0986 0.3249 -3.382 0.00405 ** inter - ceifa == 0 -0.8979 0.3019 -2.974 0.01491 * xilo - ceifa == 0 -0.6568 0.2780 -2.363 0.08255 . inter - humi == 0 0.2007 0.3794 0.529 0.95122 xilo - humi == 0 0.4418 0.3606 1.225 0.60632 xilo - inter == 0 0.2412 0.3401 0.709 0.89176 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 (Adjusted p values reported -- single-step method) > ##Resultados > tapply(trat1_90antes$ocorr, trat1_90antes$grupo,mean, na.rm = TRUE)##Media dos tratamentos ceifa humi inter xilo 4.500000 1.500000 1.833333 2.333333 > tapply(trat1_90antes$ocorr,trat1_90antes$grupo,sd,na.rm=TRUE)/sqrt(tapply(trat1_90antes$ocorr,trat1_90antes$grupo,length))###Erros padrão da media ceifa humi inter xilo 0.6708204 0.5000000 0.5426274 0.3333333 *Aqui tornou-se difícil saber a diferença entre os grupos funcionais em termos de atribui a, b, c, as médias, então resolvi fazer uma análise gráfica usando o pacote HH:* > ##Visualizando as comparações > contrast.1<- + if.R(r=glht(mod2, linfct=mcp(grupo="Tukey")), + s=multicomp(mod2, plot=FALSE)) > plot(contrast.1) *Que me deu um gráfico pouco informativo, porém seguindo o help do mmc {HH}, encontrei que o gráfico de comparações deveria ser feito usando as funções glht.mmc e multicomp.mmc, que requerem o uso de aov(), então fiz:* > > ## Contraste por mmc > contrast.aov <- aov(ocorr ~ grupo, data=trat1_90antes) > contrast.mmc <- + if.R(r=glht.mmc(contrast.aov, linfct=mcp(grupo="Tukey")), + s=multicomp.mmc(contrast.aov, plot=FALSE)) > plot(contrast.mmc) > # *Pois bem, agora tenho a comparação entre os grupos, mas não utilizando a distribuição de erros quasipoisson e não sei o que fazer, pois o pacote HH não me permite utilizar os objetos criados usando glm(), para se fazer as comparações múltiplas. Alguém do grupo teria alguma sugestão, Obrigado, Alexandre* -- Alexandre dos Santos Engenheiro Florestal, Dr. Universidade Federal de Lavras Departamento de Entomologia Laboratório de Entomologia Florestal Caixa Postal 3037 37200-000 - Lavras/MG Fone: +55 (35) 9223-0304

Alexandre, Você vai ter que abrir a função que faz o gráfico, ver o que ela recebe do objeto de classe aov() e criar uma nova função, com as devidas modificações, para que você passe os elementos um a um. Eu fiz isso uma vez mas perdi o script, mas eu lembro que não é difícil.Porém, antes de tentar isso pense que essas comparações são na escala do preditor linear e não na escala da variável, ou seja, apesar do gráfico ser elaborado e tal o usuário terá que reconhecer que o que está plotado são intervalos para diferenças na escala do preditor linear e fazer um exercício para converter esses valores na escala dos dados, no caso da Poisson, passando exp(). Possivelmente por essa razão que o desenvolvedor do pacote restringiu a função para classes lm() e aov(). À 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 twitter: @walmeszeviani homepage: http://www.leg.ufpr.br/~walmes linux user number: 531218 ==========================================================================

Walmes, Muito obrigado, vou seguir o seu conselho e quando conseguir posto o script. Alexandre Em 10/05/2012 16:01, Walmes Zeviani escreveu:
Alexandre,
Você vai ter que abrir a função que faz o gráfico, ver o que ela recebe do objeto de classe aov() e criar uma nova função, com as devidas modificações, para que você passe os elementos um a um. Eu fiz isso uma vez mas perdi o script, mas eu lembro que não é difícil.Porém, antes de tentar isso pense que essas comparações são na escala do preditor linear e não na escala da variável, ou seja, apesar do gráfico ser elaborado e tal o usuário terá que reconhecer que o que está plotado são intervalos para diferenças na escala do preditor linear e fazer um exercício para converter esses valores na escala dos dados, no caso da Poisson, passando exp(). Possivelmente por essa razão que o desenvolvedor do pacote restringiu a função para classes lm() e aov().
À 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> 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 Engenheiro Florestal, Dr. Universidade Federal de Lavras Departamento de Entomologia Laboratório de Entomologia Florestal Caixa Postal 3037 37200-000 - Lavras/MG Fone: +55 (35) 9223-0304
participantes (2)
-
ASANTOS
-
Walmes Zeviani