[R-br] Ajuda com os benditos contrastes.

[Walmes Zeviani]

Olympio,

Não é possível fazer testes para níveis de Trat1 dentro de Trat2 com a
formula assim ~Trat1/Trat2, isso por causa da codificação das matrizes
envolvidas. Você tem que usar Trat2/Trat1, ou Tra3*Trat2 (ou inverso), ou
ainda criar um novo fator combinando os níveis desses que pode ser
conseguido com a interaction(). No seu caso, você deve ajustar o modelo com
Trat2/Tempo/Trat1, o fator de níveis a serem comparados deve ser o último
na fórmula do modelo. Se for usar a contrast::contrast() o modelo tem que
ter estrutura cruzada e não aninhada. Veja,

Id: Unidade Amostral
#Estrutura dos dados:
set.seed(13032012)
d=expand.grid(Trat1=gl(3,18),Tempo=gl(2,1))
d$Trat2=gl(2,9,108)
d$Id=as.factor(rep(1:54,2))
d$Y=abs(rnorm(108,c(10,30,8),10))
str(d)

#Verificar estrutura para a unidade amostral i=3, por exemplo
d[d$Id==3,]: Tratamento1 =1, Tratamento2=1 Tempo=1, 2.

#Ajuste que fiz usando lme

library(nlme)

f=lme(Y~(Trat1/Trat2)*Tempo, random=~1|Id,d)
summary(f)

m0 <- lme(Y~Trat2/Tempo/Trat1, random=~1|Id,d)
fixef(m0)
anova(m0)

require(aod)

# 8 parâmetros = 8 g.l. da interação
grep("*.:Tempo*.:Trat1*.", names(fixef(m0)), value=TRUE)
idx <- grep("*.:Tempo*.:Trat1*.", names(fixef(m0)))

# última linha da anova(m0)
wald.test(Sigma=vcov(m0), b=fixef(m0), Terms=idx,
df=nrow(d)-length(fixef(m0)))

# assim é possível fazer testes de desdobramento de interação usando a
wald.test()
# é só saber que índices passar para Terms= para ter as hipóteses adequadas

require(contrast)

levels(d$Trat2)
levels(d$Trat1)
levels(d$Tempo)

c0 <- contrast(m0,
               list(Trat2="1", Tempo="1", Trat1="1"),
               list(Trat2="1", Tempo="1", Trat1="2"))

# modelo tem que estar na forma cruzada (*) e não aninhada (/)

m0 <- lme(Y~Trat2*Tempo*Trat1, random=~1|Id,d)

c0 <- contrast(m0,
               list(Trat2="1", Tempo="1", Trat1="1"),
               list(Trat2="1", Tempo="1", Trat1="2"))
c0
c0$X # vetor de contrastes

Depois é só você trocar os niveis na contrast() para demais comparações.
Você pode reunir os vetores de contraste e fazer esses contrastes com a
multcomp::glht() para ajustar o p-valor dos testes, visto que são muitas
hipóteses não ortogonais. Consulte a documentação dessas funções. Aconselho
você a sempre enviar perguntas para a R-br pois o sucesso de resposta é
maior, sua dúvida pode ajudar os demais além de ficar no arquivo da lista
para consulta futura.

À 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
twitter: @walmeszeviani
homepage: http://www.leg.ufpr.br/~walmes
linux user number: 531218
==========================================================================


2012/3/25 Olympio Neto <olympiotneto em gmail.com>

> Walmes, muito obrigado por você estar sempre disponível a ajudar... Olha,
> eu deveria te pagar, porque já me ajudou milhares de vezes com suas
> intervenções.
>
> Então, cara, eu tenho um problema. Vou colocar o que fiz aqui, usando
> aquele CMR que mandei.
>
> Id: Unidade Amostral
> #Estrutura dos dados:
> set.seed(13032012)
> d=expand.grid(Trat1=gl(3,18),Tempo=gl(2,1))
> d$Trat2=gl(2,9,108)
> d$Id=as.factor(rep(1:54,2))
> d$Y=abs(rnorm(108,c(10,30,8),10))
>
> #Verificar estrutura para a unidade amostral i=3, por exemplo
> d[d$Id==3,]: Tratamento1 =1, Tratamento2=1 Tempo=1, 2.
>
> #Ajuste que fiz usando lme
>
> library(nlme)
>  f=lme(Y~(Trat1/Trat2)*Tempo, random=~1|Id,d)
>
> Veja, e "f" é meu modelo final, o melhor e tal. Acontece que quero fazer
> comparações múltiplas considerando somente Trat1 e Tempo, desprezando Trat2
> dentro de Trat1. Estou me matando há 3 dias tentando isso e não consigo.
> Quando coloco Trat2, ela entra nos contrastes e não quero isso.
> Queria fazer isso, mas nem sei como te mostrar, desculpe, mas vai
> literalmente.
>
> Trat1:1 X Tempo:1 - Trat1:2 X Tempo:1
> Trat1:1 X Tempo: 1- Trat1:3 X Tempo:1
> Trat1:2 X Tempo:1 - Trat1:3 X Tempo:1
>
> Eu vi um exemplo parecido em
> http://www.ats.ucla.edu/stat/R/faq/testing_contrasts.htm, mas lá cria-se
> um novo modelo univariado com o comando interaction nos fatores. Poderia
> até fazer isso, mas isso iria alterar muito meu modelo, que já é
> complicado. Você vislumbra alguma forma de se resolver isso? Não queria
> perder a caracterização do aninhamento de Trat2 dentro de Trat1. Já tentei
> contrast, multcomp e tudo mais, mas sempre solicitam os níveis de saliva ou
> se coloco, eles aumentam as comparações.
>
> Desculpe te encher com esta dúvida idiota, mas os contrastes no R me
> arrebentam, e há anos. Já li um zilhão de livros do R, mas quando preciso
> ir às vias de fato, me mato.
>
> Abs e muito obrigado.
>
> Olympio
>
-------------- Próxima Parte ----------
Um anexo em HTML foi limpo...
URL: <http://listas.inf.ufpr.br/pipermail/r-br/attachments/20120326/6afc6f63/attachment.html>