Interpretação saída do pacote contrast

Bom dia pessoal, Bom já sei que esse assunto já foi motivo de tópicos passados no nosso mail-list e no internacional ,tem material na internet etc, mais essa saída do pacote contraste não entrou na minha cabeça ainda não e venho pedir ajuda do grupo. Estou comparando minha variável resposta número de óvulos de fêmeas de um parasitóide com as variáveis explicativa tratamento (com 2 níveis: 1 casal e 10 casais) e geração (com 4 níveis: 3, 6, 9 e 12 gerações). Na saída da funcao contrast(), ele solta o trata- mento10 casais no intercepto(imagino) e não faz as comparacoes dentro de cada geração, ele só desdobra para o tratamento 1casal, mesmo assim sumiu com a geracao 12, onde pergunto onde que estou errando é na funcao ou na interpretação? Sendo os argumentos que usei:
##Contraste tratamentos x desenvolvimento oviposicao a mumia m.comple<-glm(des.ov.mu~tratamento*geracao,data=dados6, family="quasipoisson") anova.comple<-anova(m.comple,test="Chi") anova.comple Analysis of Deviance Table
Model: quasipoisson, link: log Response: des.ov.mu Terms added sequentially (first to last) Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 82 9.1366 tratamento 1 1.0969 81 8.0397 2.265e-07 *** geracao 3 2.2880 78 5.7518 4.445e-12 *** tratamento:geracao 3 2.7715 75 2.9802 1.328e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(m.comple)
10 casais
# ##Fazendo os contrastes require(contrast) contr<-contrast(m.comple, list(geracao=levels(geracao),tratamento="1casal"),list(geracao=levels(geracao),tratamento="10casal")) print(contr,X=TRUE) glm model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|) -0.298869262 0.03278233 -0.36417501 -0.23356352 -9.12 75 0.0000 0.006271251 0.03351287 -0.06048982 0.07303232 0.19 75 0.8521 -0.103435587 0.03532313 -0.17380286 -0.03306831 -2.93 75 0.0045 0.049914391 0.03237891 -0.01458771 0.11441650 1.54 75 0.1274 Contrast coefficients: (Intercept) tratamento1casal geracaog3 geracaog6 geracaog9 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 tratamento1casal:geracaog3 tratamento1casal:geracaog6 tratamento1casal:geracaog9 0 0 0 1 0 0 0 1 0 0 0 1
##
summary(m.comple)## Dando uma olhada no modelo
Call: glm(formula = des.ov.mu ~ tratamento * geracao, family = "quasipoisson", data = dados6) Deviance Residuals: Min 1Q Median 3Q Max -0.38017 -0.09964 0.00000 0.03447 0.87147 Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.09063 0.02250 92.934< 2e-16 *** tratamento1casal -0.29887 0.03278 -9.117 8.93e-14 *** geracaog3 -0.15385 0.03229 -4.765 9.04e-06 *** geracaog6 -0.15621 0.03313 -4.715 1.09e-05 *** geracaog9 -0.05985 0.03230 -1.853 0.067820 . tratamento1casal:geracaog3 0.30514 0.04688 6.509 7.62e-09 *** tratamento1casal:geracaog6 0.19543 0.04819 4.055 0.000121 *** tratamento1casal:geracaog9 0.34878 0.04608 7.570 7.83e-11 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 Null deviance: 9.1366 on 82 degrees of freedom Residual deviance: 2.9802 on 75 degrees of freedom AIC: NA Number of Fisher Scoring iterations: 4 Obrigado, Alexandre

Bom dia pessoal, Bom já sei que esse assunto já foi motivo de tópicos passados no nosso mail-list e no interna
cional,tem material na internet etc, mais essa saída do pacote contraste não entrou na minha cabeça ainda não e venho pedir ajuda do grupo. Estou comparando minha variável resposta número de óvulos de fêmeas de um parasitóide com as variáveis explicativa tratamento (com 2 níveis: 1 casal e 10 casais) e geração (com 4 níveis: 3, 6, 9 e 12 gerações). Na saída da funcao contrast(), ele solta o trata- mento10 casais no intercepto(imagino) e não faz as comparacoes dentro de cada geração, ele só des-
dobra para o tratamento 1casal, mesmo assim sumiu com a geracao 12, onde pergunto onde que estou errando é na funcao ou na interpretação? Sendo os argumentos que usei:
##Contraste tratamentos x desenvolvimento oviposicao a mumia m.comple<-glm(des.ov.mu~tratamento*geracao,data=dados6, family="quasipoisson") anova.comple<-anova(m.comple,test="Chi") anova.comple Analysis of Deviance Table
Model: quasipoisson, link: log
Response: des.ov.mu
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev P(>|Chi|) NULL 82 9.1366 tratamento 1 1.0969 81 8.0397 2.265e-07 *** geracao 3 2.2880 78 5.7518 4.445e-12 *** tratamento:geracao 3 2.7715 75 2.9802 1.328e-14 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#summary(m.comple) 10 casais
# ##Fazendo os contrastes require(contrast) contr<-contrast(m.comple, list(geracao=levels(geracao),tratamento="1casal"),list(geracao=levels(geracao),tratamento="10casal")) print(contr,X=TRUE) glm model parameter contrast
Contrast S.E. Lower Upper t df Pr(>|t|) -0.298869262 0.03278233 -0.36417501 -0.23356352 -9.12 75 0.0000 0.006271251 0.03351287 -0.06048982 0.07303232 0.19 75 0.8521 -0.103435587 0.03532313 -0.17380286 -0.03306831 -2.93 75 0.0045 0.049914391 0.03237891 -0.01458771 0.11441650 1.54 75 0.1274
Contrast coefficients: (Intercept) tratamento1casal geracaog3 geracaog6 geracaog9 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 0 1 0 0 0 tratamento1casal:geracaog3 tratamento1casal:geracaog6 tratamento1casal:geracaog9 0 0 0 1 0 0 0 1 0 0 0 1
## summary(m.comple)## Dando uma olhada no modelo
Call: glm(formula = des.ov.mu ~ tratamento * geracao, family = "quasipoisson", data = dados6)
Deviance Residuals: Min 1Q Median 3Q Max -0.38017 -0.09964 0.00000 0.03447 0.87147
Coefficients: Estimate Std. Error t value Pr(>|t|) (Intercept) 2.09063 0.02250 92.934< 2e-16 *** tratamento1casal -0.29887 0.03278 -9.117 8.93e-14 *** geracaog3 -0.15385 0.03229 -4.765 9.04e-06 *** geracaog6 -0.15621 0.03313 -4.715 1.09e-05 *** geracaog9 -0.05985 0.03230 -1.853 0.067820 . tratamento1casal:geracaog3 0.30514 0.04688 6.509 7.62e-09 *** tratamento1casal:geracaog6 0.19543 0.04819 4.055 0.000121 *** tratamento1casal:geracaog9 0.34878 0.04608 7.570 7.83e-11 *** --- Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Null deviance: 9.1366 on 82 degrees of freedom Residual deviance: 2.9802 on 75 degrees of freedom AIC: NA
Number of Fisher Scoring iterations: 4
Obrigado, Alexandre

Alexandre, Melhor forma de entender a saída e conhecer a partida. Assim, simule de um experimento onde você conhece o valor dos contrastes e veja o que a função retorna. Fiz isso com um exemplo normal, coloquei variância pequena só para não ter muita diferença entre parâmetros (theta) e estimativas. Veja se fica claro. da <- expand.grid(A=gl(2,1), B=gl(3,1), rept=1:4) X <- model.matrix(~A*B, da) colnames(X) theta <- c(10,1,0,-1,0,0) da$y <- X%*%theta+rnorm(nrow(X),0,0.001) m0 <- lm(y~A*B, da) require(contrast) c1 <- contrast(m0, list(A=levels(da$A), B="1"), list(A=levels(da$A), B="2")) c1 # fez B1 vs B2 dentro de A1 e dentro de A2 cbind(t(c1$X), theta) c1$X%*%theta c2 <- contrast(m0, list(A=levels(da$A), B="1"), list(A=levels(da$A), B="3")) # fez B1 vs B2 dentro de A1 e dentro de A2 c2 cbind(t(c2$X), theta) c2$X%*%theta c3 <- contrast(m0, list(A="1", B=levels(da$B)), list(A="2", B=levels(da$B))) # fez A1 vs A2 dentro de B1, B2 e B3 c3 cbind(t(c3$X), theta) c3$X%*%theta c4 <- contrast(m0, type="average", list(A="1", B=levels(da$B)), list(A="2", B=levels(da$B))) # fez A1 vs A2 na média dos níveis de B (efeito principal) c4 cbind(t(c4$X), theta) c4$X%*%theta À 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, agora ficou bem claro. Na verdade fica muito mais fácil de entender com dados simulados e no caso do pacote contrast os exemplos são baseados em dados de obesidade que achei de difícil compreensão. Abraco, Alexandre On 08/26/2011 02:03 PM, Walmes Zeviani wrote:
Alexandre,
Melhor forma de entender a saída e conhecer a partida. Assim, simule de um experimento onde você conhece o valor dos contrastes e veja o que a função retorna. Fiz isso com um exemplo normal, coloquei variância pequena só para não ter muita diferença entre parâmetros (theta) e estimativas. Veja se fica claro.
da <- expand.grid(A=gl(2,1), B=gl(3,1), rept=1:4) X <- model.matrix(~A*B, da) colnames(X) theta <- c(10,1,0,-1,0,0) da$y <- X%*%theta+rnorm(nrow(X),0,0.001)
m0 <- lm(y~A*B, da)
require(contrast)
c1 <- contrast(m0, list(A=levels(da$A), B="1"), list(A=levels(da$A), B="2")) c1 # fez B1 vs B2 dentro de A1 e dentro de A2 cbind(t(c1$X), theta) c1$X%*%theta
c2 <- contrast(m0, list(A=levels(da$A), B="1"), list(A=levels(da$A), B="3")) # fez B1 vs B2 dentro de A1 e dentro de A2 c2 cbind(t(c2$X), theta) c2$X%*%theta
c3 <- contrast(m0, list(A="1", B=levels(da$B)), list(A="2", B=levels(da$B))) # fez A1 vs A2 dentro de B1, B2 e B3 c3 cbind(t(c3$X), theta) c3$X%*%theta
c4 <- contrast(m0, type="average", list(A="1", B=levels(da$B)), list(A="2", B=levels(da$B))) # fez A1 vs A2 na média dos níveis de B (efeito principal) c4 cbind(t(c4$X), theta) c4$X%*%theta
À 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.
participantes (2)
-
ASANTOS
-
Walmes Zeviani