[R-br] multcomp/doBy comparação multipla com interação

Fernando Antonio de souza nandodesouza em gmail.com
Segunda Dezembro 3 20:07:55 BRST 2012 

Caro Walmes, muito obrigado pela atenção. Foi muito útil suas sugestões.
entretanto não estou entendendo algo e gostaria que vc me esclarecesse
Coloquei desta vez um CMR com um banco de dados em que a interação deu
significativa. Seguindo o comando que você colocou para no caso de
interação entretanto não apresenta nenhum efeito significante.. Não
entendo o que aconteceu. Você pode me explicar?


dput(abate[,1:8])structure(list(Animal = c(1L, 2L, 3L, 4L, 5L, 7L, 8L, 9L, 12L,
14L, 15L, 17L, 18L, 19L, 21L, 22L, 23L, 25L, 26L, 28L, 30L, 32L,
34L, 35L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 46L, 47L, 48L,
49L, 50L, 53L, NA, NA), Gest = structure(c(4L, 2L, 3L, 2L, 2L,
4L, 3L, 4L, 2L, 2L, 3L, 4L, 1L, 4L, 2L, 1L, 4L, 1L, 2L, 4L, 3L,
3L, 1L, 4L, 3L, 1L, 2L, 3L, 1L, 1L, 3L, 1L, 2L, 3L, 1L, 2L, 1L,
1L, NA, NA), .Label = c("0", "100", "130", "140"), class = "factor"),
    Manej = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
    2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L,
    1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, NA, NA
    ), .Label = c("1", "2"), class = "factor"), rep = c(2L, 3L,
    1L, 2L, 4L, 3L, 1L, 2L, 2L, 5L, 4L, 3L, 1L, 4L, 1L, 2L, 4L,
    4L, 3L, 1L, 3L, 5L, 6L, 1L, 2L, 5L, 1L, 2L, 4L, 3L, 3L, 3L,
    4L, 4L, 5L, 5L, 1L, 2L, NA, NA), PV = c(49000L, 40800L, 51800L,
    41900L, 35000L, 40500L, 35400L, 42300L, 36600L, 36600L, 38400L,
    46400L, 26900L, 50200L, 35500L, 35500L, 48400L, 40500L, 32500L,
    42800L, 36500L, 50000L, 31200L, 33500L, 46500L, 40900L, 36500L,
    31200L, 31500L, 35500L, 38700L, 30000L, 37900L, 42700L, 35500L,
    42000L, 34500L, 36700L, NA, NA), PCVZ = c(39841.2, 35457.1,
    42343.2, 32989.7, 26566.5, 30339.9, 25965.5, 28475.2, 30815.7,
    28162.6, 28103.5, 37905.8, 19174.5, 34319.4, 27963.2, 26609.4,
    40852.9, 34781.1, 23411.8, 43490.8, 28900.1, 38240.8, 22463.9,
    24631.9, 36483.9, 33688.2, 31838.7, 26048.1, 25646.6, 22815.6,
    31369.6, 22173.7, 30371.2, 36740.8, 25773, 38171.6, 26401.4,
    26711.2, NA, NA), CAR.PVZ = c(0.474383302, 0.518936969, 0.493585747,
    0.512281106, 0.508158771, 0.451550598, 0.427490324, 0.431954824,
    0.490009962, 0.511316427, 0.462575836, 0.485413842, 0.547602284,
    0.416673951, 0.507810265, 0.541162146, 0.474874489, 0.554899069,
    0.478391239, 0.494357427, 0.46020602, 0.489006506, 0.534190412,
    0.418156943, 0.482404567, 0.558058905, 0.515096408, 0.45684714,
    0.538083021, 0.547870755, 0.471794349, 0.518632434, 0.520229691,
    0.435483169, 0.52768401, 0.458456025, 0.537850266, 0.539099703,
    NA, NA), PCARC = c(18900L, 18400L, 20900L, 16900L, 13500L,
    13700L, 11100L, 12300L, 15100L, 14400L, 13000L, 18400L, 10500L,
    14300L, 14200L, 14400L, 19400L, 19300L, 11200L, 21500L, 13300L,
    18700L, 12000L, 10300L, 17600L, 18800L, 16400L, 11900L, 13800L,
    12500L, 14800L, 11500L, 15800L, 16000L, 13600L, 17500L, 14200L,
    14400L, NA, NA)), .Names = c("Animal", "Gest", "Manej", "rep",
"PV", "PCVZ", "CAR.PVZ", "PCARC"), class = "data.frame", row.names = c(NA,
-40L))

library(car)
library(doBy)
library(multcomp)


abate$Gest<-factor(abate$Gest)
abate$Manej<-factor(abate$Manej)
C<-lm(PCARC~Gest*Manej,abate)
Anova(C,type="III")
levels(abate$Manej)
X<-popMatrix(C, effect="Gest", at=list(Manej="2"))
cb <- combn(nrow(X), 2)
Xc <- X[cb[1,],]-X[cb[2,],]
summary(glht(C, linfct=Xc))
popMeans(C, effect="Gest", at=list(Manej="2"))


Em 3 de dezembro de 2012 16:39, Walmes Zeviani
<walmeszeviani em gmail.com>escreveu:

> Fernando,
>
> Parabéns, seu código é totalmente reproduzível (com excessão da função
> Anova() que tá ali sem carregar o pacote car, mas eu conhecia). CMR é meio
> raro na nessa lista. Para resolver seu problema você fazer assim.
>
> abate <-
> structure(list(Animal = c(1L, 2L, 3L, 4L, 5L, 7L, 8L, 9L, 12L,
> 14L, 15L, 17L, 18L, 19L, 21L, 22L, 23L, 25L, 26L, 28L, 30L, 32L,
> 34L, 35L, 37L, 38L, 39L, 40L, 41L, 42L, 43L, 44L, 46L, 47L, 48L,
> 49L, 50L, 53L, NA, NA), Gest = structure(c(4L, 2L, 3L, 2L, 2L,
> 4L, 3L, 4L, 2L, 2L, 3L, 4L, 1L, 4L, 2L, 1L, 4L, 1L, 2L, 4L, 3L,
> 3L, 1L, 4L, 3L, 1L, 2L, 3L, 1L, 1L, 3L, 1L, 2L, 3L, 1L, 2L, 1L,
> 1L, NA, NA), .Label = c("0", "100", "130", "140"), class = "factor"),
>     Manej = structure(c(1L, 1L, 1L, 2L, 2L, 2L, 2L, 2L, 1L, 1L,
>     2L, 1L, 2L, 2L, 2L, 2L, 1L, 1L, 2L, 1L, 2L, 1L, 2L, 2L, 1L,
>     1L, 1L, 2L, 2L, 1L, 1L, 2L, 1L, 1L, 2L, 2L, 1L, 1L, NA, NA
>     ), .Label = c("1", "2"), class = "factor"), rep = c(2L, 3L,
>     1L, 2L, 4L, 3L, 1L, 2L, 2L, 5L, 4L, 3L, 1L, 4L, 1L, 2L, 4L,
>     4L, 3L, 1L, 3L, 5L, 6L, 1L, 2L, 5L, 1L, 2L, 4L, 3L, 3L, 3L,
>     4L, 4L, 5L, 5L, 1L, 2L, NA, NA), PV = c(49000L, 40800L, 51800L,
>     41900L, 35000L, 40500L, 35400L, 42300L, 36600L, 36600L, 38400L,
>     46400L, 26900L, 50200L, 35500L, 35500L, 48400L, 40500L, 32500L,
>     42800L, 36500L, 50000L, 31200L, 33500L, 46500L, 40900L, 36500L,
>     31200L, 31500L, 35500L, 38700L, 30000L, 37900L, 42700L, 35500L,
>     42000L, 34500L, 36700L, NA, NA)), .Names = c("Animal", "Gest",
> "Manej", "rep", "PV"), class = "data.frame", row.names = c(NA,
> -40L))
>
> abate$Gest <- factor(abate$Gest)
> abate$Manej <- factor(abate$Manej)
> abate <- na.omit(abate)
>
> xtabs(~Gest+Manej, abate)
>
> C <- lm(PV~Gest*Manej, abate)
> anova(C)
> car::Anova(C,type="III")
>
> # interação não significativa à 5%!
> D <- lm(PV~Gest+Manej, abate)
>
> X <- doBy::popMatrix(D, effect="Gest")
>
> choose(nlevels(abate$Gest), 2)
> cb <- combn(nrow(X), 2)
> Xc <- X[cb[1,],]-X[cb[2,],]
> Xc
> rownames(Xc) <- apply(combn(levels(abate$Gest), 2), 2, paste, collapse="-")
>
> require(multcomp)
>
> summary(glht(D, linfct=Xc))
>
> # em caso de interação significativa fazer
> # isso aqui para cada nível que deseja fixar
>
> levels(abate$Manej)
> X <- doBy::popMatrix(C, effect="Gest", at=list(Manej="1"))
> X
>
> # o resto segue normal...
> # para automatizar os desdobramentos você pode usar a lapply() ou
> # similares.
>
> A opção de constrate nada mais é que uma restrição paramétrica no vetor de
> efeitos que te leva para uma particular solução ou conjunto de estimativas.
> Como você falou, o treatment anula (impõe/assume) que o efeito do primeiro
> nível de cada fator é zero, o soma zero impõe que a soma deles é zero. Em
> todas elas, o que você precisa estimar ao final são p-1 parâmetros/efeitos
> pois a restrição cuidou de um. Todos contra todos na saída do summary() não
> tem como especificar porque isso requer mais que p-1 hipóteses você só pode
> ir com p-1. Mas isso não é problema porque uma vez estimado os efeitos,
> qualquer função e quantidade de funções estimáveis a partir deles pode ser
> obtida.
>
> À 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 em ufpr.br
> skype: walmeszeviani
> twitter: @walmeszeviani
> homepage: http://www.leg.ufpr.br/~walmes
> linux user number: 531218
> ==========================================================================
>
> _______________________________________________
> R-br mailing list
> R-br em 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.
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20121203/eecd262e/attachment.html>

Mais detalhes sobre a lista de discussão R-br