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