Teste t pareado como um caso especial de modelagem linear de efeitos mistos

20

Sabemos que um teste t emparelhado é apenas um caso especial de ANOVA unidirecional de medidas repetidas (ou dentro do sujeito), bem como modelo linear de efeitos mistos, que pode ser demonstrado com a função lme (), o pacote nlme em R como mostrado abaixo.

#response data from 10 subjects under two conditions
x1<-rnorm(10)
x2<-1+rnorm(10)

# Now create a dataframe for lme
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

Quando executo o seguinte teste t emparelhado:

t.test(x1, x2, paired = TRUE)

Eu obtive este resultado (você obterá um resultado diferente por causa do gerador aleatório):

t = -2.3056, df = 9, p-value = 0.04657

Com a abordagem ANOVA, podemos obter o mesmo resultado:

summary(aov(y ~ x + Error(subj/x), myDat))

# the F-value below is just the square of the t-value from paired t-test:
          Df  F value Pr(>F)
x          1  5.3158  0.04657

Agora posso obter o mesmo resultado no lme com o modelo a seguir, assumindo uma matriz de correlação simétrica definida positiva para as duas condições:

summary(fm1 <- lme(y ~ x, random=list(subj=pdSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.3142115  9 -0.7918878  0.4488
# xx2          1.3325786 0.5779727  9  2.3056084  0.0466

Ou outro modelo, assumindo uma simetria composta para a matriz de correlação das duas condições:

summary(fm2 <- lme(y ~ x, random=list(subj=pdCompSymm(form=~x-1)), data=myDat))

# the 2nd row in the following agrees with the paired t-test
# (Intercept) -0.2488202 0.4023431  9 -0.618428  0.5516
# xx2          1.3325786 0.5779727  9  2.305608  0.0466

Com o teste t emparelhado e a ANOVA de medidas repetidas unidirecionais, posso escrever o modelo tradicional de média de células como

Yij = μ + αi + βj + εij, i = 1, 2; j = 1, ..., 10

onde i indexa a condição, j indexa o sujeito, Y ij é a variável de resposta, μ é constante para o efeito fixo para a média geral, α i é o efeito fixo para a condição, β j é o efeito aleatório do sujeito após N (0, σ p 2 ) (σ p 2 é variação da população) e ε ij é residual após N (0, σ 2 ) (σ 2 é a variação dentro do sujeito).

Eu pensei que o modelo de média de células acima não seria apropriado para os modelos lme, mas o problema é que não consigo criar um modelo razoável para as duas abordagens lme () com a suposição da estrutura de correlação. O motivo é que o modelo lme parece ter mais parâmetros para os componentes aleatórios do que o modelo de média de células acima oferece. Pelo menos o modelo lme fornece exatamente o mesmo valor F, graus de liberdade e valor p, o que os gls não podem. Mais especificamente, gls fornece DFs incorretos devido ao fato de que ele não explica o fato de que cada sujeito tem duas observações, levando a DFs muito inflados. O modelo lme provavelmente é superparâmetro na especificação dos efeitos aleatórios, mas não sei qual é o modelo e quais são os parâmetros. Portanto, o problema ainda não foi resolvido para mim.

pólo azul
fonte
2
Não tenho certeza do que você está perguntando. O modelo que você anotou é precisamente o modelo para o modelo de efeitos aleatórios; a estrutura de correlação é induzida pelo efeito aleatório.
Aaron - Restabelece Monica
@ Aaron: o efeito aleatório βj no modelo de média das células deve seguir N (0, σp2). Minha confusão é: como esse termo (com apenas um parâmetro σp2) está associado à estrutura de correlação especificada pela simetria composta ou por uma simples matriz simétrica no modelo lme?
bluepole
Quando você calcula a correlação entre as duas observações sobre o mesmo assunto, a correlação é sigma_p ^ 2 / (sigma_p ^ 2 + sigma ^ 2) porque eles compartilham o mesmo beta_j. Veja Pinheiro / Bates, pág. Além disso, o modelo de efeito aleatório, como você o escreveu, é equivalente à simetria composta; outras estruturas de correlação são mais complexas.
Aaron - Restabelece Monica
@ Aaron: Obrigado! Eu já li o livro Pinheiro / Bates sobre isso e ainda não consegui descobrir os detalhes sobre os efeitos aleatórios. As páginas mais relevantes parecem ser o exemplo em P.160-161. Além disso, a saída de efeitos aleatórios de lme () com suposição de simetria composta parece não concordar com a correlação de σp2 / (σp2 + σ2) no modelo de média celular. Ainda perplexo com a estrutura do modelo.
22412 bluepole #
Bem, quase equivalente à simetria composta; no CS a correlação pode ser negativa, mas não com efeitos aleatórios. Talvez seja aí que a sua diferença surja. Consulte stats.stackexchange.com/a/14185/3601 para obter detalhes.
Aaron - Restabelece Monica

Respostas:

16

A equivalência dos modelos pode ser observada calculando a correlação entre duas observações do mesmo indivíduo, como segue:

YEuj=μ+αEu+βj+ϵEujβjN(0 0,σp2)ϵEujN(0 0,σ2)Cov(yEuk,yjk)=Cov(μ+αEu+βk+ϵEuk,μ+αj+βk+ϵjk)=Cov(βk,βk)=σp2Vumar(yEuk)=Vumar(yjk)=σp2+σ2σp2/(σp2+σ2)

Observe que os modelos, no entanto, não são exatamente equivalentes, pois o modelo de efeito aleatório força a correlação a ser positiva. O modelo CS e o modelo t-test / anova não.

EDIT: Existem outras duas diferenças também. Primeiro, os modelos CS e efeito aleatório assumem normalidade para o efeito aleatório, mas o modelo t-test / anova não. Em segundo lugar, os modelos CS e efeito aleatório são adequados usando a máxima verossimilhança, enquanto a anova é adequada usando quadrados médios; quando tudo está equilibrado, eles concordam, mas não necessariamente em situações mais complexas. Finalmente, eu seria cauteloso ao usar valores F / df / p dos vários ajustes como medidas de quanto os modelos concordam; veja a famosa mesa de Doug Bates em df's para mais detalhes. (END EDIT)

O problema com o seu Rcódigo é que você não está especificando a estrutura de correlação corretamente. Você precisa usar glscom a corCompSymmestrutura de correlação.

Gere dados para que haja um efeito de assunto:

set.seed(5)
x <- rnorm(10)
x1<-x+rnorm(10)
x2<-x+1 + rnorm(10)
myDat <- data.frame(c(x1,x2), c(rep("x1", 10), rep("x2", 10)), 
                    rep(paste("S", seq(1,10), sep=""), 2))
names(myDat) <- c("y", "x", "subj")

A seguir, veja como você ajustaria os efeitos aleatórios e os modelos de simetria composta.

library(nlme)
fm1 <- lme(y ~ x, random=~1|subj, data=myDat)
fm2 <- gls(y ~ x, correlation=corCompSymm(form=~1|subj), data=myDat)

Os erros padrão do modelo de efeitos aleatórios são:

m1.varp <- 0.5453527^2
m1.vare <- 1.084408^2

E a correlação e variação residual do modelo CS é:

m2.rho <- 0.2018595
m2.var <- 1.213816^2

E eles são iguais ao esperado:

> m1.varp/(m1.varp+m1.vare)
[1] 0.2018594
> sqrt(m1.varp + m1.vare)
[1] 1.213816

Outras estruturas de correlação geralmente não se ajustam a efeitos aleatórios, mas simplesmente especificando a estrutura desejada; Uma exceção comum é o modelo de efeito aleatório AR (1) +, que tem um efeito aleatório e correlação AR (1) entre observações no mesmo efeito aleatório.

EDIT2: Quando encaixo as três opções, obtenho exatamente os mesmos resultados, exceto que o gls não tenta adivinhar o df para o termo de interesse.

> summary(fm1)
...
Fixed effects: y ~ x 
                 Value Std.Error DF   t-value p-value
(Intercept) -0.5611156 0.3838423  9 -1.461839  0.1778
xx2          2.0772757 0.4849618  9  4.283380  0.0020

> summary(fm2)
...
                 Value Std.Error   t-value p-value
(Intercept) -0.5611156 0.3838423 -1.461839  0.1610
xx2          2.0772757 0.4849618  4.283380  0.0004

> m1 <- lm(y~ x + subj, data=myDat)
> summary(m1)
...
            Estimate Std. Error t value Pr(>|t|)   
(Intercept)  -0.3154     0.8042  -0.392  0.70403   
xx2           2.0773     0.4850   4.283  0.00204 **

(A interceptação é diferente aqui porque, com a codificação padrão, não é a média de todos os assuntos, mas a média do primeiro assunto.)

Também é interessante notar que o lme4pacote mais recente fornece os mesmos resultados, mas nem tenta calcular um valor-p.

> mm1 <- lmer(y ~ x + (1|subj), data=myDat)
> summary(mm1)
...
            Estimate Std. Error t value
(Intercept)  -0.5611     0.3838  -1.462
xx2           2.0773     0.4850   4.283
Aaron - Restabelecer Monica
fonte
Mais uma vez obrigado pela ajuda! Eu conheço essa parte da perspectiva do modelo de média celular. No entanto, com o seguinte resultado de lme () com simetria composta: Efeitos aleatórios: Fórmula: ~ x - 1 | subj Estrutura: Simetria composta StdDev xx1 1,1913363 xx2 1,1913363 Corr: -0,036 Residual 0,4466733. Ainda não consigo reconciliar esses números com o modelo de média de células. Talvez você possa me ajudar a resolver esses números?
bluepole
Além disso, algum pensamento sobre a formulação do modelo com outras estruturas de correlação, como matriz simétrica simples?
22612 bluepole #
Eu vejo! Eu deveria ter lido sua resposta no outro tópico com mais cuidado. Pensei em usar gls () antes, mas não consegui descobrir a especificação de correlação. É interessante que lme () com estrutura de simetria composta para o efeito aleatório ainda tenha o mesmo valor t, mas parece que as variações para os efeitos aleatórios não são diretamente interpretáveis. Eu realmente aprecio sua ajuda!
22412 bluepole #
Depois de pensar duas vezes, sinto que minha confusão original ainda não foi resolvida. Sim, os gls podem ser usados ​​para demonstrar a estrutura de correlação e as médias quadradas das colunas, mas o modelo abaixo dele não é exatamente o mesmo que o teste t pareado (ou ANOVA de medidas repetidas unidirecionais em geral), e essa avaliação é ainda suportado pelos DFs incorretos e pelo valor p de gls. Por outro lado, meu comando lme com simetria composta fornece os mesmos valores de F, DFs e p. A única coisa que me deixa intrigado é como o modelo lme é parametrizado conforme indicado no meu post original. Alguma ajuda por aí?
bluepole
Não sei como ajudá-lo. Você poderia escrever o que acha que são os dois modelos diferentes? Algo está errado na maneira como você pensa sobre um deles.
Aaron - Restabelece Monica
3

Você também pode considerar usar a função mixedno pacote afexpara retornar valores de p com a aproximação de Kenward-Roger df, que retorna valores de p idênticos como um teste t emparelhado:

library(afex)
mixed(y ~ x + (1|subj), type=3,method="KR",data=myDat) 

Ou

library(lmerTest)
options(contrasts=c('contr.sum', 'contr.poly'))
anova(lmer(y ~ x + (1|subj),data=myDat),ddf="Kenward-Roger")
Tom Wenseleers
fonte